Theory of Stable Dynamical Systems
Gassmann, Fritz
Group Dynamics of Innovative Systems, Paul Scherrer Institute
CH-5232 Villigen PSI, Switzerland
Tel. +41 310 26 47; Fax. +41 310 44 16
E-mail gassmann@psi.ch
Abstract
During development of a dynamical model to simulate the central T-cell subsystem of the human
immune system, its extraordinary stability lead to the assumption that it might be deduced solely from
stability considerations. We demonstrate that linear stability conditions together with additional
general requirements indeed define a low number (of the order of 10) of dynamical systems from about
10" possibilities for systems with up to four components. At least two of the most simple linearly stable
systems play central roles in biology, indicating that stable dynamical subsystems resulting from
evolutionary processes might be understood on a mathematical basis. We expect that our dynamical
subsystems are generic and will be found as basic building blocks in biological, social or technical
systems.
Key words: linearly stable dynamical systems, theoretical analysis, generic
structures, immune system, evolution, structural biology.
Introduction
The immune system! fulfills the difficult task to attack and eliminate a large
variety of unknown invaders (antigens such as viruses, bacteria) without harming the
own tissue. Based on extensive research over the last decades, a detailed but still
growing model of the dynamics of the immune system could be developed. At the
heart of this very complicated structure are B- and T-cells, maturing in the bone
marrow and in the Thymus, respectively. They interact mediated by chemical
substances (cytokines) such as interleukins, interferons, etc. Growth of a specific
subgroup of T-cells can be triggered by antigen presenting cells displaying short
strings of foreign peptides (resulting from cutting an antigenic protein into small
pieces) that match receptors located on the surface of the T-cells. With this
mechanism, only T-cells directed against the particular antigens, are activated and
multiply to initiate a process finally eliminating the antigens specifically. The vast
number of T-cell subsystems (estimates give 10’..10'' different T-cell receptors each
specific for a single antigen) is permanently present and waits to be activated by an
antigen attack. If one specific subsystem gets triggered, cytokine-levels change due to
the immune response. As cytokines act on all subsystems in parallel, they are all
influenced by the excursion of one single subsystem. It is of vital importance that
these subsystems remain in a stable state under any conditions. As soon as antigens
are eliminated, the subsystems involved in the elimination process, from any state of
excitation, have to relax to their quiescent and stable ground state. These
considerations illustrate that stability is one of the most important issues for the
immune system.
A central part of the immune system, being in the focus of ongoing research on
the development of allergies’, is the interaction of two types of helper T-cells called
Thl and Th?’ that interact by secretion of cytokines as shown in Fig. 1.
. +
+ interferon-y
antigen decay
level decay
interleukin
+
’ decay
>
decay
Fig. 1: Interaction of helper T-cells of type 1 and 2 mediated by interferon-y (IFN-y) and interleukins
(IL-4, 5, 13).
Examining this symmetric interaction scheme, the fundamental questions arises, if
other structures would be possible. From the modeling point of view, the main
question will be how to translate the interaction scheme into differential equations,
i.e., into a dynamical model. In the present paper, we answer both questions on the
basis of generic stability considerations and other generic assumptions. Our result will
be that there is one unique structure that can be interpreted as two cell-types
interacting with two substances (the one shown in Fig. 1) and fulfilling all imposed
generic conditions. Our investigation explains why evolutionary processes, starting
from the most simple systems, necessarily had to find this structure.
With the same set of generic assumptions, we will analyze systems with 1, 2, 3
and 4 components. However, the method we develop in the present paper can be
applied also to systems with more than four components.
Systems with 1 component
A dynamical model for a 1-component system can be written as follows:
dx x
art A(T) + 9(T)x-= ()
f(T) and g(T) are nonlinear functions of external parameters T. In many cases, T is,
e.g., temperature influencing a system, but T can be considered as a vector including
more than one external influence (or force). For our investigation, we linearize f and
g and write f(T)=a,u and g(T)=a,v with the external forces uy, v being positive
(including zero). Due to our linearization, the parameters a; vary according to the
operating point (i.e., the linearized domain around T). Further, parameter values in
most biological or social systems are subject to many influences and therefore can
vary in a considerable range. We assume that the parameters a, are either positive or
negative, but we exclude their sign to change. This assumption is important, because
it defines whether we call two systems with the same structure but different parameter
values to be different or the same. In many cases, stability of a system will change if
the sign of certain parameters change. For similar reasons, we restrict all variables
to positive values. For many systems, this restriction arises quite naturally, e.g., for
variables describing concentrations or populations. The term -x/t describes a decay
with a positive time constant t we always include into our equations. For biological,
biochemical or social systems, such a decay term describes lifetimes of cells,
biomolecules or individuals. In general, such systems are called dissipative. Finally,
we impose our stability condition that there must be a linearly stable equilibrium
state for all , v >0. The equilibrium condition reads:
= ore (2.1)
—-a,v
T
Our condition for x to be positive, irrespective of the norms of all a, implies aj20, a,20
and a,<0. To ensure linear stability, we claim the derivative of the r.h.s. of (1) towards
x to be negative:
va,-1<0 (2.2)
T
giving once more the same condition a,<0. According to (2.1), ay defines the value of
x for uw=v=0. By a linear transformation x—>x-x,), the equilibrium for vanishing
external forces can always be shifted to any point. In this report, we assume zero
to be the equilibrium point for vanishing external forces and therefore, the
parameter a) will be omitted. Further, a term a,vx with negative a, is an inhibiting
external influence and drives the system into a "dead" state for large v. To avoid such
"dead" states, we do not allow the variables to vanish with large external forces.
We therefore skip the term a,;vx with negative a,. The resulting absolutely stable
system is therefore:
dx x
a ep a,20 GB)
The graphical representation of the relaxation dynamics (3) is given in Fig. 2.
Fig. 2: Relaxation dynamics (3) for a
population or substance x with an external
+ influence (effect of 4). Decay refers to the
external negative term -x/t with positive time constant
influence ™
decay
Though we imposed with (2.2) only linear stability, i-e., relaxation of the system
to its equilibrium after small perturbations, the analytical solution of the dynamical
equation (3) for constant 4 shows that the system returns to the same equilibrium
point from any initial state and therefore, its equilibrium is not only linearly but also
absolutely stable. In some nonlinear dissipative systems we will analyze below, linear
stability will not exclude the existence of limit cycles and chaotic dynamics, and
therefore, our linearly stable dissipative systems we will identify later might not
always show stable behavior after large perturbations.
Systems with 2 components
A system with components x and u that is controlled by an external influence (or
force) u can be written as follows, if only bilinear or linear terms in all variables
are taken into account:
o 0
ox a "4a,x|1]4a,,u]°)-*
dt sa ra 22") 4 r,
u
u (4.1)
o 0
du 1 Oo] u
aa ah +a,,ul |-—
i 1|
u
u
u and 7; are positive, aj, a,; can be positive or negative. The brackets [] give different
possibilities for the respective terms. In general, for n components (variables), there
are
(° + 2) +(n 4) «= {+H combinations with repetitions (4.2)
ry !
For n=2, we can combine the different terms and the equations in 300 ways. To
specify the system, a first set of generic requirements is applied:
We require a coupling of each equation to at least one other equation for p=0
to avoid trivial systems, i.e., we do not choose n+1 zeros in any equation.
We avoid combinations with a,x-x/t, because such terms can be combined
into one term x/T'.
Neither equation should read dx/dt=x{...}, because this would imply an
equilibrium state x=0 independent of the other variables and of the external
influences. We reject such equations due to only weak interactions with other
components of the system. In fact, only the relaxation time constant would be
influenced by other components or by external forces and the component
governed by such an equation would vanish and decouple from the remaining
components: A n-component system would collapse to a simpler (n-1)-
component system.
The combination a,wuta,u = a,(u+a,/a,)u = a,u'u does not introduce new
dynamics and so can be considered as equivalent to a,p'u (the linear
transformation u'= u+a,/a, eliminates the term a,u).
¢ We exclude equations of the form u(a, ,x+a,,) to simplify our analysis. It has
to be investigated later, if interesting systems are excluded by this
prescription.
¢ Mathematical isomorphisms (systems becoming identical after exchanging
Takii
names of variables, e.g., x and u) are excluded.
ing into account the above stated requirements, (4.1) reduces to the following
equations for the variables x and u:
1) OX 4 sea
af 7 SHRM
2) ox = au+aux-~
dt ity 6.1)
dx x
3) oie apx+ au-—
1
dx x
4) — =aput+aux-—
) gp Panu ey :
d
1) Ye autax—
dt a
2) Ye aytaux-
dt tT, (6.2)
3) a0 utax-2
dt | aa
2
du u
4) — =ayx+aux-—
yap 7 ext a .
2 1
investigated. We introduce the nomenclature i,k with k2i (i, k=1...4) to refer to the
system consisting of equation i for x and k for u. To demonstrate our stability
analysis, we begin with the most simple linear system 1,1. The matrix of partial
derivatives (we call it stability matrix) becomes:
Combinations of (5.1) and (5.2) give [2 )-$4-10 systems to be
‘ 1 (6)
It follows from basic mathematics‘ that a stationary solution is linearly stable (i.e.
stable against small perturbations) if the real parts of all eigenvalues of the stability
matrix are negative. Calculation of the eigenvalues > leads to the following
polynomial P(A):
-[»+ +)p+t)-aa, (7)
T, T,
The two eigenvalues are defined by P(A)=0. t, are positive decay time constants, but
all a, are not restricted at the moment to positive or negative values. If the term d=a,a,
would vanish, the zero points of P(A) would be obvious: they would be located in the
left complex half plane. However, d>0 shifts the graph of the polynomial downwards,
and moves its right zero point towards positive real numbers as shown in Fig. 3. This
means that critical combinations of the parameters a,, and 1, exist for which one
eigenvalue gets zero, i.e., beyond these critical combinations, one real eigenvalue
becomes positive. As we look for systems that are linearly stable for any combination
of parameter values (however, without changing their signs), such runaway solutions
for certain combinations are forbidden and we impose, e.g., a,<0 and a,>0 giving d<0.
P(
F >)
Fig. 3: Two polynomials P(A) with
d=a,a, = 0 and d>0. At a critical
value d., A,=0, and for d>d., the
stationary state of the system
becomes unstable. For d<0, the
function is shifted upwards and the
real X approach each other until
they merge. From this point on, its
real parts are constant and their
imaginary parts separate vertically
and symmetrically to the A-axis on
the point-dashed-line (we assume
the imaginary axis drawn as
vertical axis on the same plot).
Because also imaginary zero points would have negative real parts, system 1,1 is
acceptable from the stability point of view. However, the stationary solution for x
reads
x
= UT
a,+a,a,T,
1
1-a,a,1,T,
a,>0, a,<0,a,>0,a,>0 (8)
and would give negative values for -a,a;t,>a,. To avoid this possibility, we set a,=0
and find the first system 1,1 meeting all our requirements:
system 1,1:
M -ayrau-™ a,>0,a,<0 (9)
4
du u
— =a,x- a,>0
dt * ti, am
external
influence
decay decay
Fig. 4.1: The linearly stable 2-component system 1,1 according to equation (9).
Analogous analysis of system 1,2 shows, that both conditions (no positive real
eigenvalues, positive variables) can be met for all parameter values if a,>0, a,>0,
a,>0, a,<0. As a last check, we have to search for complex eigenvalues crossing the
imaginary axis towards the right complex half plane. For this aim, we eliminate a,u in
P(A) using the first equation of (5.1) to get
P(r) -[ + r+ t-ax|-a(%-anl (10.1)
qT, T, qT,
To find eigenvalues on the imaginary axis, we evaluate P(A) for A=in giving the
following equations for the real and the imaginary part:
y= 2a -a,(-a,}u
- T,t. Tt.
F ' 4 ' (10.2)
if 2 +—+(-a, » =0
tT, Tv,
The imaginary part of equation (10.2) implies n=0 and the real part equation becomes
Que hag ca)|% 4 as 10.3)
TT, T,
From the dynamical equation for x follows (for the stationary solution)
x
—-ay=au>0 (10.4)
Tt,
and this implies that (10.3) can never be satisfied. Therefore, no complex eigenvalues
can cross the imaginary axis to the right complex half plane with any combinations of
the parameter values, and system 1,2 is linearly stable and meets all our conditions:
system 1,2:
dx x
a eee a,>0,a,>0 (11.1)
1
ju u
=~ =ai+aux-— a,>0,a,<0
dt 9
external +
influence <«—____._ &xternal
influence
decay decay
Fig. 4.2: The linearly stable 2-component system 1,2 according to equation (11.1) can be interpreted as
a cell population u and a substance x produced by the cells that has a inhibitory effect on the
cells.
System 1,2 shown in Fig. 4.2 is often found in a biological context: cells u
stimulated by an external influence produce a cytokine x inhibiting their growth. In
this context, the external influence on the cytokine x might come from other cells
producing the same cytokine. The external influence on the cells u could be a stimulus
for cell multiplication.
Investigation of all 10 combinations of equations (5.1) and (5.2) give the results
summarized in Tab. 1. The two additional systems 2,2 and 2,4 that were found are
given below as equations (11.2), (11.3) and in graphical form in Fig. 5. The two
systems 1,3 and 2,3 with a,=0 are special cases of systems 1,1 and 2,1 (2,1 is identical
with 1,2 after exchanging x and u) and therefore, they give no new dynamics.
system 2,2:
May raxu-* a,>0,a,<0 (11.2)
4
OY aut aux- a,>0,a,<0
dt 7s
system 2,4:
dx x
ap aN aU a,>0,a,<0 (11.3)
4
oo = aux aux a,>0,a,<0
2
+
< external
influence
external
influence
decay decay
Fig. 5: The linearly stable 2-component systems 2,2 (without broken line) and 2,4 (including broken
line) according to (11.2) and (11.3).
The analysis for 2-component systems shows that mainly combinations of
equations | and 2 lead to systems meeting our conditions. Though for systems with
3 or more components, valid systems containing equations of type 3 and 4 can
surely be found, we restrict our further analysis to systems combined with
equations of type 1 and 2 only. With this restriction, we exclude systems with
nonlinear external influences (such as, e.g., system 2,4) described by ax, a,wy, au
and av. Our motivation to postpone this class of systems to a later analysis is to
simplify analytical calculations as well as the structure of the present report.
system
positive
u, xX
stability
Li
yes
yes
12
P(A) = (- + +
atteax
yes
yes
13
1,4
-[+4 ne t-ax
}-a(ansaa)
no
2,2
P(r)=[ae tau
Tv
ye
A+ —-a,x
yes
yes
2,3
2,4
P(2) = - + sau
Ty
A+ A - ax] - ax(ayu + au)
t
yes
yes
33
P(A) -[r+2-an
%
neta
no
34
P(a) = [-+t-a0
tT,
xed nax
T,
J-alansau)
no
44
) -(au+axj(au+au}
no
Table 1: Polynomials P(A) for all 10 combinations of equations (5.1) and (5.2). The signs of the
parameters a, were chosen to avoid negative variables for the equilibrium state (y
u, x"). If this first requirement could be met, the P(A) was checked for the absence of eigenvalues with
in column "positive
positive real parts (yes in column "stability"). Systems 1,1 / 1,2 / 1,3 / 2,2 / 2,3 / 2,4 are retained as
valid systems. As we exclude "dead" equations with a;<0, a; is set to zero in systems 1,3 and 2,3. With
this restriction, these two systems are special cases of systems 1,1 and 1,2.
11
Systems with 3 components
We reverse now the procedure followed in the previous section and begin with
the stability requirements. We will demonstrate first that a stability matrix leading to a
polynomial P(A) of the form (7) must have the generic form:
*- z0
Oz (12)
zoe
In every row and every column, there is only one matrix element different from zero
(z) in addition to the diagonal elements (¢). All elements z or * can be constants or
functions of those variables that are not restricted by zeros in the same row. With the
most simple settings z= a; and * =-1/1 , the eigenvalue problem is the following:
1 a 0
— 1
% 4 § 5
0 -— a, |*/nl=An (13)
T,
4 c
3 "3
To calculate the eigenvalues A, we solve the equations
n ~2eae4)
1 1
te daft) (14)
Ib 2
stent]
3 3
giving the polynomial P(A)
- T(-+2)-TT as)
with the correct form (7). Introducing a constant b instead of zero in the first row of
the stability matrix (13) would give
P(
(16)
For b=0 (16) reduces to (15). However, for b different from zero, an additional
relation between the parameters has to be true for preventing real parts of eigenvalues
to become positive:
12
ab < (7)
TT,
As we look for systems that are stable irrespective of the values of the involved
parameters, inequalities such as (17) have to be avoided.
As the choice of the matrix (12) is only one of different possibilities, we
investigate variants arising from permutations of the variables: Exchanging the names
of u and v in a system with variables x, u, v transforms (12) into the following matrix:
- 02
ze? 0 (18)
Ooze
In the first row are the partial derivations of the x-equation in the order x, u, v. As u
and v are exchanged, the second and third positions have to be exchanged. The
second and third row refer to the u- and y-equations that change place according to
the transformation u'=v, v'=u. The transformed system has to be arranged in the order
x, u', v'. Further permutations of (18) do not result in new structures.
We explained above why we accept only matrices with a maximum of one non-
zero position in addition to the diagonal element in every row. We explain now, why
we require the same condition for the columns. The following equivalent stability
matrices fulfill the row conditions, but not the column conditions:
*0z *z0 *zo0 *zo0 *0z *0z
z* O}=|z ¢ O|/=!z ¢ O}=]0 © zj/>/0 © zl=|0 « 2{(19)
zoe zoe Oze Ooze O0ze zoe
As all 6 matrices are equivalent, it is enough to analyze the first one. The two zeros in
the second column mean that the x and v-equations are not influenced by the u-
equation, i.e., the x-v-system can be solved independently, and the system can be
reduced to a 2-component system influencing a 1-component system. It is therefore
enough to investigate systems with a stability matrix of the form (18), i.e., we have to
investigate the following differential equations (excluding nonlinear external
influences and terms of the form u(a,,x+a,)):
OK 4 L+a,v 1%
dt Nes x] T,
du 1] u
at =ajut aa] - L (20)
Wea veagul ||
dt v] T
Tab. 2 gives an overview of the results we obtained for the remaining
on 4-3-2
—— =4 combinations.
3 3-2-1
system positive stability signs of a; valid
variables system
1,1,1 yes no a,>0, a,>0, a,<0 no
a,=a,=0, a,>0
1,1,2 yes yes a,>0, a,>0, a,<0 yes
a,=a,=0, a>0
1,2,2 yes yes a,>0, a,<0, a,<0 yes
a,=0, a;>0, a,>0
2,2,2 yes yes a,<0, a,<0, a,<0 yes
a,>0, a;>0, a,>0
Tab. 2: Overview of our results from the analysis of systems with 3 components combined from
equations of type 1 and 2 according to (20).
We identified the three systems 1,1,2 / 1,2,2 / 2,2,2 as valid systems according to our
requirements. The special conditions a,=a,=0 we imposed on system 1,1,2 and a,=0
for system 1,2,2 helped to simplify the analysis, but they could turn out to be
generalizable to a,>0, a;>0. These three systems are given by equations (21), (22),
(23) and in graphical form in Figs. 6 and 7.
dx
x
3
—=av-— system 1,1,2:
dt 7° ot,
WU ax a, >0,a, >0,a, <0
dt Ty
Wm aut auv-~ a,>0O
3
XX awe. system 1,2,2:
dt A
= ayes a.xu- 2 a, >0,a,<0,a, <0
T,
Wm aur auv- a, >0,a,>0
3
* = a+ a,VX— ue system 2,2,2:
4
WU ea ueaxu- 2 a, <0,a,<0,a, <0
dt T
WV auaie. a,>0,a, >0,a, >0
dt t
Q)
(22)
(23)
14
external
influence
external
influence ~
external
influence
external
influence
external
influence
external
influence
Fig. 6: System 1,1,2
according to (21).
Fig. 7.1: System 1,22
according to (22).
Fig. 7.2: System 2,2,2
according to (23).
15
We would like to add some comments concerning two other systems:
System 1,1,1 referenced in Tab. 2 is excluded due to stability reasons. Its polynomial
P(A)
P(A) = [: + +h + +f. + 4) -a,a,a, (24.1)
tT, T, T3
gives for A=in the following conditions resulting from P(in)=0:
WY +=-a.a,+[]* a, <0,a, > 0,a, >0
w= So
iT,
by
(24.2)
We can therefore separate the t-terms from the a-terms giving
(Bee}Ls}Te-se4>2 es
The negative product term is cancelled by the same positive term from the product
{}{} and so, the I.h.s. of (25) is a sum of positive terms. For any choice of the 1; it is
possible to choose a set of a; that meet the condition (25), i.e., complex eigenvalues
are found on the imaginary axis and therefore, system 1,1,1 gets oscillatory unstable
for many parameter combinations.
System 1,2,3 with a,>0, a,<0, a,<0, a,=0, a,>0, a;>0 is an example of a valid system
with nonlinear external influence we excluded from further investigations to simplify
the present analysis.
16
Systems with 4 components
We follow our procedure developed for 3 component systems and begin with the
stability matrix and its 9 possible permutations. There are two symmetry groups
shown in (26): an asymmetric group with six elements and a symmetric group with
three elements.
asymmetric group
* 02z0 - z0
ee“ NOOO
Oo ON
.
COON Oo
* NOOO
ooo.
oo oO N
*NOO
.
.
ooo i=}
“NOOO
(26)
ON O
Oo
.
0
zoe
0
symmetric group
* 02z0 - 0 Zz
0*-0z O° 0
=>
zoe 0 Ooze 0
0zo0-s. z0o0-*s.
Within the two symmetry groups, all elements are connected by permutations of the
equations as explained above (see (12), (18) and (19)).
The symmetric group consists of systems composed of two non-interacting 2-
component systems. This property can be understood with the third element, where
the two non-connected subsystems are indicated with broken lines. The symmetric
group can therefore be reduced to 2-component systems we discussed above.
The asymmetric group consists of six elements that show up as symmetric pairs
with respect to the main diagonal. The symmetric pairs (located vertically one on top
of the other in (26)) are connected by the two permutations 1-3 and 2-4. A 2-3
permutation and a 1-2 permutation connects the pairs in (26) horizontally. As all six
elements are connected by permutations, they represent the same system, and we can
restrict our investigation to the first element:
-02z0
07 02z
Q7)
Ooze 0
zoos.
17
Due to our above explained restrictions, there remain only | 4 =e
wa 5:4:3:2 _¢
combinations of the following equations to be investigated:
dx x dx x
1) aoe 2) qe ant aK
1 4
d d
1) Gen att av- = 2) Gen att avy- 2
2 ? 28
du u du u ee)
1) at OF a 2) at tt ay
3 3
dv v dv v
1) qe ae 2) Ge 7 at ax "|,
Systems combined of equations of type 1 and 2 have the following stability
matrix:
Aran) 0 2 [| 0
T, 1 x
0 14 a,Vv o 0 a,* |
. 1 0 ” -
0 a,: mae: ay: 0
u T; 1
leading to the eigenvalue-polynomial:
P(0) = [a-t-au file oe fn il ae SD
coffe
The system 1,1,1,1 (all upper values in brackets [] in (30)) gives the same structure
for P(A) as the above investigated system 1,1,1 and has to be excluded due to the
possibility of complex eigenvalues with positive real parts.
System 2,1,1,1 (equivalent to system 1,1,1,2) leads to
P(in) = [n + Lo aun + +Yin + +fin + 4) -a,a,a,a,x = 0 (31)
t. tT. Tt. T,
1 2 3 ‘4
18
The real and imaginary parts give the conditions
a) a (32)
2 _ bk TTT, apt GT,
1
S=-au
int T
For the special case a,=a,=a,=0, we can find the following relation from the four
equations defining system 2,1,1,1
u
a,a,a,X = , a
TyTgVy
(33)
to eliminate x in (32). With a,<0 and all other parameters as well as u,x positive, the
second equation of (32) always gives two real yn (symmetric to zero) that can be
inserted into the first equation giving:
Rek-wzs} -2e{St-of
-{St-affy -auy Alps -au$4|
i opt T
G4)
A combination of parameters such that (34) holds would mean that two symmetric
imaginary eigenvalues would be about to enter the right complex half plane and the
system would show oscillatory instability. To exclude such instability, we have to
show that (34) never holds. On the I.h.s. of (34), we have 49+25=74 positive terms,
on the r.h.s. we have 5-7-9=315 positive terms. Using basic mathematics, it can be
shown that all 74 terms show up among the 315 terms also on the r.h.s. (a practical
way to prove this is to arrange all terms as powers of (-a,u) and to prove the
cancellation separately for power 0, 1, 2, and 3). The result is that (34) becomes the
form
0 =sum of positive terms (35)
what obviously never can be true. System 2,1,1,1 has therefore no oscillatory
instability and the last point we have to show is that its variables in the stationary
states are always positive. We set the time derivatives of the dynamical equations to
zero and solve the equations analytically giving:
X= ay {te vis 4aa.cu} (36)
OL
O =-8,a,a,a,T,T,1,T, > 0
From (36) follows x>0 for all possible parameter values and from the remaining linear
equations we get also y, u and v positive. All our requirements for a system to be
accepted are met by system 2,1,1,1 according to equations (37). Fig. 8 shows the
system in graphical form.
19
system 2,1,1,1: a, <0,a,>0,a,>0,a,>0
dx x
—=aut+aux-— a,>0,a,=a,=a,=0 (37)
dt T,
dy y du u_ dv v
ay apo oN, yp RO
dt » at Tt, dt Te,
Fig. 8: Linearly stable 4-
component system
2A,11.
*
a ~
S/o \ .
/
external +
influence p> decay
a) /
\ S
Se 4
+
The difference to our stable 2-component system 1,2 (Fig. 4.2) is the series of three
substances acting on one another without being influenced directly by external forces
instead of only one substance that can be externally influenced. Further investigations
have to show if the condition a;=a,=a,=0 (i.e. no external influence on y, u, v) is
necessary or could be omitted.
The symmetric system 2,2,1,1 (equivalent to system 1,1,2,2) has the following P(A):
P(r) = [: + A. aula + A. av] + +f. + 4) -aa, wl (38)
T, T, T; T, T3T,
In (38) the product a,agxy has been replaced by uv/(t,t,) according to the dynamical
equations for the special case a,=a,=0. From (38) follows that a,, a, be negative, and
the stationary states obey the relations
a ant
X= y= a,,a,,a,,a, > 0
Tay Tau 1 gs Ags By
71 t
aa,t a,a,T, (39)
us Acts v= StH aa, <0
A 4 TRALACT IML 47 au
%; 7 au
20
The u-equation can be solved giving
ae ibs votes
“2,
‘4 T, (40)
b 1 T,
= -4,0,a, t+ aaa ta, C= aau—
TT, 2
Only the solution with the positive sign in front of the square root leading to u>0 has
been retained, i.e., there is always exactly one solution for any arbitrary choice of the
parameters (meeting the conditions a,, a,<0). The polynomial (38) leads to the
following conditions for P(in)=0:
nt-vif(4-ay] tenet t}s(2-aal tots 1 jt paecet. tho
T T Ts Uy Te TT Ty 1t,| TT, | T TT,
1 1 -av-au {1 1\(-au. -av
n'fav-eus St} -¥ +— 4 (2-41 a+ |
iJ Bek TTT, TT q 3 UG
(41)
u and v have to be inserted into (41) from (39) and (40), leading to an extremely
complicated expression. However, analogous to the procedure explained before, we
can prove that no solution exists by reducing (41) to an equation of type (35).
Therefore, system 2,2,1,1 is a second system with 4 components that meets all our
conditions. The system is given in graphical form in Fig. 9.
system 2,2,1,1:
dx x dy y
—=au+avx-— —=au+auy-— a,>0,a,>0,a,=-a (42)
dt t, dt ty
du u dv v
— =a,x-— —=ay-— a, <0,a,<0,a,>0,a,>0
dt, dt a q 7 \s 6 a
+
+
- u
external decay
influence decay
. Vv
+
+ decay
-——+|
decay
Fig. 9: Linearly stable 4-component system 2,2,1,1 according to equations (42).
21
The structure of this system is equivalent to the T1-T2-cell system, the basic building
block of the immune system, shown in Fig. 1. Here, x, y are T1 and T2 helper cells, u,
v are interferon-y and an indicator for Interleukins 4, 5, 13, respectively. Fig. 10
shows the same system in a form more familiar for scientists working in the field of
system dynamics.
ate decay ine Tt
a -
ae
_changerste "| Tt helper |_ seca 7
cels
i *
change rate
deday rate
change rate
Peal interteukin
interleukin
decavtine
—~ “lesteukin
x
YT2 helper
change rate T2 calls decay rate T2
28 WW Revine rz
Fig. 10: Linearly stable 4-component system 2,2,1,1 in system dynamics view.
Summary of imposed conditions and results
For didactic reasons, we introduced our conditions step by step at the place where
we needed them. To avoid the impression of arbitrariness, we sum up here the entire
list of conditions to show once again their generic character.
* we impose the stability condition that for all external influences 1>0 and for
all parameter values (not changing signs), a stationary state exists and is
stable.
¢ the stationary state with all variables being zero must result from u=0.
¢ equilibrium values for all variables must be positive for all parameter
combinations and all u>0.
* each component decays, described by a positive time constant.
* quadratic terms like x’ are excluded, only bilinear terms are allowed, such as
ux, vx, etc. (linearization around equilibrium).
¢ neither equation should read dx/dt=x{...}, because this would imply an
equilibrium state with x=0 only or equilibrium states for all x (with {}=0). For
22
both cases, the terms within the brackets {} would have to fulfill conditions
that cannot be true for arbitrary choices of the parameters.
¢ We require a coupling of each equation to at least one other equation for u=0
to avoid trivial systems.
¢ the combination a,wuta,u = (a,u+a,)ju = au does not introduce new
dynamics and so can be considered as equivalent to a,u'u.
¢ mathematical isomorphisms (systems getting identical after exchanging its
variables) are considered as one and the same system.
All of the above stated conditions arise quite naturally in a biological or sociological
context. However, we omitted systems with nonlinear external influences of type aux
and equations of the form u(a,,x+a,) to simplify our analysis. It is an open question,
how many interesting systems are lost by these restrictions.
The main result of our mathematical analysis is a reduction of the number of
valid systems from about 10'° to 10. From these 10 most simple systems, at least 2 are
of central importance in a biological context. The 10 linearly stable systems are
summarized in "stenographic" form in Fig. 11.
The number 10'° for possible system configurations is derived from the 4-
component systems according to relation (4.2):
(4+2)1+(4— i} _ ("| _ 720°
4 4
=1.12-10" (43)
1-component system
>
2-component systems
system 1,2 system 1,1 system 2,2 system 2,4
3-component systems
KO
system 1,1,2 system 1,2,2 system 2,2,2
4-component systems Legend:
input +
' + output
input—4
component
— positive, linear input
positive, nonlinear input
A)
—e negative, linear input
— negative, nonlinear input
system 2,1,1,1 system 2.2.11 ------ external input not yet analyzed
Fig. 11: Overview of the ten most simple linearly stable systems found among about 10"° different
possibilities. Components with inhibitory nonlinear inputs (open circles) may be interpreted as
cells, the other components as chemical substances. Broken lines refer to external inputs that
were neglected in this report to simplify calculations. All components are assumed to have a
decay-term driving the respective variables to zero when all external influences (forces) are
absent. Systems with nonlinear external influences of type a,ux and equations of the form
u(a);X+a,) were omitted in 3- and 4-component systems. Also systems 2,2,2,1 and 2,2,2,2
were not yet analyzed.
24
Discussion
We have analyzed interactions in systems with up to four components based on a
set of conditions we imposed on them. Most of our conditions are based on general
requirements we think are mandatory for a lot of systems where stability in a very
large parameter space is an important issue. However, our set of restrictions might be
too strong for certain classes of systems (only systems with linear external influences,
no equations of the form u(a,,x+a,,)). Therefore, our resulting ten systems have to be
considered as a subset of all conceivable linearly stable systems with up to four
components. A first extension of this minimal subset might follow from the inclusion
of additional external influences (indicated by broken lines in Fig. 11). Then, more
systems would be found if we considered also nonlinear external influences. In
addition, weaker conditions might enlarge the class of linearly stable systems. We
would like to point out that our stability analysis only ensures linear stability against
small amplitude oscillations around stationary states. Nonlinear systems can display
large amplitude nonlinear oscillations (called limit cycles) not detected with our linear
analysis, i.e., some of our systems might show this type of instability and will have to
be cancelled from the list of absolutely stable systems.
Though the number of possibilities grows quicker than exponentially with
increasing number of components n (see (4.2) and (43)), the number of most simple
stable system configurations found for n=1,2,3,4 was not increasing: we found one to
four systems for every n. The next step to analyze the about 10'° 5-component
systems seems feasible and might indicate how a search for systems with more
components might go on. It is an interesting open question, if solutions can be found
for every n or if there exists an upper limit for the number of components. A different
but also important question concerns the stability of systems resulting from the
combination of many systems with four and less components. In the immune system,
the helper cell subsystems are built upon specific T1 and T2 cells resulting in a large
number (in the order of 10’ to 10'') of parallel 2,2,1,1-subsystems that are
interconnected by a few substances (such as interferons and interleukins). What are
the conditions that an interconnection of stable systems remains stable?
Considerable activity in the field of system dynamics was devoted to define
generic structures over the last decade® resulting in the definition of a series of useful
tools (e.g. generic system archetypes, generic infrastructures) supporting model
conceptualization, systemic thinking, and communication to share dynamic insights.
In contrast, our classification is based on mathematical concepts and has to be
considered as a theory aimed towards understanding of system structure rather than a
set of practical rules helping model building and communication. We think that the
combination of both, theory and practical rules, are the ingredients having the
potential to boost the development of the field of system dynamics. Our approach
towards a theory of stable systems seems to have a considerable potential for further
development.
25
Conclusions
We have developed a theory of stable dynamical systems based on imposed
generic conditions and some simplifications. We found the surprising result that our
requirements were strong enough to reduce the variety of possible systems from the
order of 10'° to 10. Several additional candidates are open to be analyzed further. It is
interesting that two of the ten most simple systems play basic roles in biology, where
stability in extremely varying conditions is mandatory. We think that the here
presented theory of stable systems reveals important generic dynamical structures that
can be used as building blocks for larger systems. It is open for further investigations,
if the theory can be extended to systems with five or more components and to
interactions of different systems. A general theory of stable dynamical systems would
boost development of the interdisciplinary field of system dynamics. Further
applications for such a theory could be found in the upcoming new research field
called structural biology dealing with models to simulate the very complex
biochemical reactions within a living cell.
Acknowledgement
I thank the immunology group at Children's Hospital in Zurich (Switzerland),
lead by Roger P. Lauener, for many interesting discussions that triggered my
wondering why extremely large systems can remain stable under constantly varying
conditions. I am indebted to Silvia Ulli-Beer, who introduced me in the views and
conceptions of the system dynamics community and to David Andersen for his
motivating interest in my research.
References:
"Abbas, A. K. & Lichtman, A. H. 2006: Basic Immunology — Functions and Disorders of the
Immune System, Saunders, Philadelphia (PA), 322 p.
Lauener, R. P. et al. 2002: Expression of CD14 and Toll-like receptor 2 in farmers! and non-
farmers’ children, Lancet 360, 465f.
* Jiang, H. & Chess, L. 2006: Regulation of immune responses by T cells, N. Engl.J.Med. 354, 1166-
1176.
Hsu, Sze-Bi 2006: Ordinary differential equations with applications, World Scientific, Series on
applied mathematics 16, 244 p.
Wolstenholme, E. F. 2003: Towards the definition and use of a core set of archetypal structures in
system dynamics, Sys. Dyn. Rev. 19, 7-26.
