Untangling Complex Systems, page 29
21
1
1/ r 12
b
b
f
1/ r
1
12
f
g
g
(iii)
(iv)
a
1
1/ r 21
a
1/ r 21
1
FIGURE 5.8 Phase diagrams for antagonistic relationships between species A and B. In (i), both r 12 and r 21 are less than 1; in (ii), both r 12 and r 21 are larger than 1; in (iii), r 12 > 1 and r 21 < 1; in (iv), r 12 < 1 and r 21 > 1. The black dots represent the possible steady states. In (ii), the continuous grey curve represents the stable eigenvector of
the saddle point and divides the positive quadrant of the a–b space into two regions: R1 and R2.
initial condition lies on the separatrix because random fluctuations will move the system out of the
separatrix (i.e., the eigenvector represented as a continuous grey curve in Figure 5.8[ii]) and one of the two species will disappear.
In situations (iii) and (iv), one species is stronger than the other. In (iii), r 12 > 1 and r 21 < 1, and species B finishes to dominate, whereas the other disappears. In (iv), r 12 < 1 and r 21 > 1, and it is species A to win and species B to die. In all cases, except (i), the competition finishes excluding one species.
5.6.2 muTualism
In a mutualistic relationship, the coefficients n and n appearing in [5.33], are both equal to 2. The A
B
differential equations describing how the number of individuals of species A and B change over time
are as follows:
dA
2
A
k
=
12
k 1 FA − k−1 A + k 12 AB = k 1 FA 1−
+
B
dt
KA k−1 K
A
[5.44]
dB =
2
B
k 21
k 2 FB − k 2 B + k 21 AB = k 2 FB 1−
−
A
dt
−
KB k−2 KB
To make easier the treatment of the system [5.44], we nondimensionalize it as we did in the case of the
antagonistic relationship. We introduce three dimensionless variables: a = AC A; b = BC B and τ = tCt, and we define the non-dimensionalizing constants A
1/
)
1/
)
C = ( K A , BC = ( K B and tC = k 1 F . We obtain:
da
k 12 K
=
a 1− a
B
+
b
dτ
k−1 K
A
[5.45]
db
k 2 F
k 21 K
=
b 1− b
A
+
a
dτ
k
1 F
k−2 KB
134
Untangling Complex Systems
Finally, fixing ρ = ( k
/ −
)
/ −
)
2 k
/ 1), r 12 = ( k 12 KB k 1 KA , and r 21 = ( k 21 KA k 2 KB , the two differential equations assume their final dimensionless form:
da = a(1− a+ r 12 b) = f
τ
d
[5.46]
db = ρ b(1− b+ r a
21 ) = g
dτ
The possible steady-state solutions are: ( a
) 0 0) (1 0) (0 )1 ((1 12) (1 12 21) (1 21)
ss , bss = ( ,
; , ; , ;
+ r
− r r , + r
1
( − r r ))
12 21
. As usual, their stability can be inferred by the linear analysis method. The Jacobian is
f
∂
f
∂
a
∂
b
∂
1− 2 ass + 12
r bss
12
r ass
J
ss
ss
=
=
[5.47]
g
∂
g
∂
ρ 21
r bss
ρ − 2ρ bss + ρ r a
21 ss
a
∂
b
∂
ss
ss
The steady state (0,0) is unstable, because both tr ( J ) and det ( J ) are positive.
The steady state (1,0) is a saddle point. In fact,
−1
r 12
J =
, tr ( J ) = ρ (1+ r 21) −1, det ( J ) = −ρ (1+ r 21) < 0 [5.48]
0
ρ (1+ ρ r 21)
The steady state (0,1) is also a saddle point. In fact,
1+ r 12
0
J =
, tr ( J ) = (1+ r 12 ) − ρ, det ( J ) = −ρ (1+ r 12 ) < 0 [5.49]
ρ r 21
−ρ
For the last steady state, we have
1+ r
12
r 12 ( r 12 + )
1
( r 12 r 21 − )
1
(1− r 12 r 21)
J =
ρ r
(
)
21 ( r 21 + )
1
ρ r 21 +1
(1− r
(
)
12 r 21 )
r 12 r 21 −1
1
( + r )
(
)
12 + ρ 1 + r
tr ( J
21
) =
(
[5.50]
r
)
12 r 21 −1
ρ 1
( + r )(
) ρ
(
)(
)
12
1+ r 21 − r 12 r 21 1+ r 12 1+ r
det ( J
21
) =
( r
)
12 r 21 −1 2
When r
( )
( )
12 r 21 < 1, det J
> 0 and tr J < 0. Therefore, the fourth steady state is stable. When r 12 r 21 > 1,
the fourth steady state does not lie in the positive quadrant of the a−b space and therefore, it does
not have physical meaning.
In synthesis, when r 12 r 21 > 1, the mutualism has large values for the kinetic constants k and in 12
k 21
the case of comparable carrying capacities. In such case, we have only three possible steady states
that are either unstable or saddle points. The populations of the two species grow unboundedly
The Emergence of Temporal Order in the Ecosystems
135
2
(i)
g
b 1
f
02 0
a 1
2
b
g
1
f
(ii)
0 0
a
1
2
FIGURE 5.9 Possible scenarios for a mutualistic relationship. Case (i) requires r 12 r 21 > 1; case (ii) holds when r 12 r 21 < 1 .
to infinity (see case (i) in Figure 5.9). On the other hand, when r 12 r 21 < 1, the mutualism is less strong. We have four steady states. The fourth one represents a stable solution, and all trajectories
of its phase space tend to ( a , b )
ss
ss = (>1, >1). In such steady state (see case (ii) in Figure 5.9), the
populations of the two species are larger than their maximum values achievable with a neutral
relationship.11
To see the dynamics of the populations of species A and B in the case of the other symbiotic
relationships, try exercises 5.7, 5.8, 5.9, and 5.10.
5.7 KEY QUESTIONS
• Describe the Lotka-Volterra model for the predator and prey relationship.
• Is the predator-prey relationship always oscillatory? Describe the dynamics in the case of
Type I, II and III functional responses.
• What does the Poincaré-Bendixson theorem state?
• Make examples of symbiotic relationships.
• Which are the possible outcomes of an antagonistic relationship?
• Why should human relationships be mutualistic?
5.8 KEY WORDS
Predator and prey; Functional and numerical responses; Symbiotic relationships.
5.9 HINTS FOR FURTHER READING
Besides Mathematical Biology by Murray (2002), additional books on mathematical modeling of
ecology are A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution by Otto and
Day (2007) and Elements of Mathematical Ecology by Kot (2001).
11 Ideally, any human relationship should be mutualistic. It is the only one that guarantees a reciprocal growth.
136
Untangling Complex Systems
5.10 EXERCISES
5.1. Apply the linear stability analysis to the improved Lotka-Volterra model where the growth
of the prey population follows the “logistic model” when there are not predators:
k
1 →
F + H ←
2 H
k−1
H + C
k
2→
C
2
C
k
3→
D
5.2. Apply the phase-space analysis to the dynamical system of exercise (5.1).
5.3. Solve the system of two differential equations [5.23] based on Type II functional response
and a numerical response proportional to the rate of prey consumption, using the follow-
ing values for the parameters (extracted from Harrison [1995] and regarding unicellular
Paramecium Aurelia and its unicellular predator Didinium nasutum):
K
k
1 F = 898; k 1 F = 1.85; k
12 4
. ; k
h = 25.5; K = 284.1; ′ h =
3 = 2.07.
Initial conditions: H
0 = 56 and C 0 = 23.
5.4. According to the results presented by Harrison (1995), when the prey Paramecium and the
predator Didinium are grown in a cerophyl medium thickened by methyl-cellulose, the
searching effectiveness of the Didinium dramatically reduces. Therefore, k becomes much
c
smaller and K much larger. The reduction of the nutrient supply for Paramecium may hinder
the growth rate k 1 F and diminish the carrying capacity K 1 F. For the following sets of parameters K
1 F = 400; k 1 F = 0.82; kh = 25.5; K = 350; k′ h = 12.4; k 3 = 2.07; H 0 = 87; C 0 = 10, determine the stationary state and its stability.
5.5. Consider the system of two differential equations [5.23], involving a Type II functional
response, and the values of the constants used in exercise 5.3:
K
898
1 85
25 5
284 1 k′
1 F =
; k 1 F = . ; kh = . ; K
=
. ; h = 12. ;
4 k
3 = .
2 0 .
7
Try to verify the Poincaré-Bendixson theorem by constructing a trapping region TR, i.e., a
closed set, such that the vector field points inward everywhere on its boundary.
5.6. Apply the linear stability analysis to confirm that for the system of two differential equa-
tions [5.32], there are three regions of stability as shown in Figure 5.7.
5.7. Predict the dynamics of a neutral relationship between two species, A and B, considering the general scheme [5.33].
5.8. Predict the dynamics of parasitism of species B towards species A. Considering the general scheme [5.33], fix n
A = 0 and nB = 2.
5.9. Predict the dynamics of amensalism between two species, A and B, assuming that
n
A = 0 and nB = 1 in the general scheme [5.33].
5.10. Predict the dynamics of commensalism involving two species, A and B, assuming that
n
A = 2 and nB = 1 in the general scheme [5.33].
The Emergence of Temporal Order in the Ecosystems
137
5.11 SOLUTIONS TO THE EXERCISES
5.1. The differential equations are
dH = k
2
1 FH − k−1 H
− k 2 HC
dt
dC = k 2 HC − k C
dt
3
The solutions of stationary states are:
H
(
2
2 1
−1 3 )
ss = k 3 k
/ 2 and Css = k k F − k k / k 2.
The Jacobian is:
k− k
− 1 3
−
k
k
3
J =
2
.
k− k
k
1 3
1 F −
0
k 2
tr ( J ) = − k− k
1 3 k 2 is always negative.
k
2
2
− k
1 3
k 1 Fk 3 k 2 k− k
det ( J ) = k
1 3
1 Fk 3 −
=
−
.
k 2
k 2
The eigenvalues are:
2
k
2
−
1
− k 3
k 1
− k 3
k k
1 3
2
−
±
4
− 1
k Fk −
tr ( J ) ± ( tr ( J )) − 4 det ( J )
k
k
3
k
λ
2
2
2
1,2 =
=
2
2
The determinant is positive if k 1 F > k− k
1 3 k 2 . If it is positive also the discriminant (Δ), the
eigenvalues are two real negative roots, and the stationary state is a stable node. When
Δ = 0, the stationary state is a stable star. Finally, if Δ < 0, the roots are two complex num-
bers having a negative real term. Therefore, the stationary state is a stable spiral, and the
predator and prey populations exhibit damped oscillations with the predator oscillations
lagging in phase behind the prey.
The determinant is null when k 1 F = k− k
1 3 k 2 . In this case, one eigenvalue is a negative real
number, and the other is 0: the system admits a line of stable fixed points. When the determi-
nant is negative, i.e., when k 1 F < k− k
1 3 k 2, the discriminant is positive, and the roots are real,
but one positive and the other negative. In this case, the stationary state is a saddle point.
5.2. In the system of two differential equations of exercise 5.1 there are only two variables, H, and C: dH = k
2
1 FH − k−1 H
− k 2 HC
dt
dC = k 2 HC − k C
dt
3
We can represent the evolution of the system in a phase space having C and H as coordi-
nates. The equations representing the null-clines are:
138
Untangling Complex Systems
C• = 0
( C is large and H is small)
H• = 0
H• < 0
( k 1 F/ k 2)
C• < 0
H• < 0
C
H• > 0
C• > 0
•
( H and C are large)
C < 0
H• > 0
( H and C are small)
C• > 0 ( H is large and C is small)
( k 3/ k 2)
H
( k 1 F/ k−1)
FIGURE 5.10 Phase plane for the “modified” Lotka-Volterra model when the det(J) is positive.
H = 0 when C = ( k 1 F − k−1 H ) k 2 representing a straight line with ( k
) as intercept
1 F/k 2
and –( k−
) as slope.
1 /k 2
C = 0 when H = ( k 3 k 2 ) representing a straight line parallel to the C axis.
When
k
( )
1 F > k− k
1 3 k 2 , i.e., det J
> 0, the phase space appears like that depicted in
Figure 5.10.
The intersection of the two nullclines is the stationary state that in the linear analysis of
stability we discovered it is stable (see exercise 5.1). This result is confirmed qualitatively
by the vector field shown in Figure 5.10. The signs of
H and
C in the different regions of
1
1/ r 12
b
b
f
1/ r
1
12
f
g
g
(iii)
(iv)
a
1
1/ r 21
a
1/ r 21
1
FIGURE 5.8 Phase diagrams for antagonistic relationships between species A and B. In (i), both r 12 and r 21 are less than 1; in (ii), both r 12 and r 21 are larger than 1; in (iii), r 12 > 1 and r 21 < 1; in (iv), r 12 < 1 and r 21 > 1. The black dots represent the possible steady states. In (ii), the continuous grey curve represents the stable eigenvector of
the saddle point and divides the positive quadrant of the a–b space into two regions: R1 and R2.
initial condition lies on the separatrix because random fluctuations will move the system out of the
separatrix (i.e., the eigenvector represented as a continuous grey curve in Figure 5.8[ii]) and one of the two species will disappear.
In situations (iii) and (iv), one species is stronger than the other. In (iii), r 12 > 1 and r 21 < 1, and species B finishes to dominate, whereas the other disappears. In (iv), r 12 < 1 and r 21 > 1, and it is species A to win and species B to die. In all cases, except (i), the competition finishes excluding one species.
5.6.2 muTualism
In a mutualistic relationship, the coefficients n and n appearing in [5.33], are both equal to 2. The A
B
differential equations describing how the number of individuals of species A and B change over time
are as follows:
dA
2
A
k
=
12
k 1 FA − k−1 A + k 12 AB = k 1 FA 1−
+
B
dt
KA k−1 K
A
[5.44]
dB =
2
B
k 21
k 2 FB − k 2 B + k 21 AB = k 2 FB 1−
−
A
dt
−
KB k−2 KB
To make easier the treatment of the system [5.44], we nondimensionalize it as we did in the case of the
antagonistic relationship. We introduce three dimensionless variables: a = AC A; b = BC B and τ = tCt, and we define the non-dimensionalizing constants A
1/
)
1/
)
C = ( K A , BC = ( K B and tC = k 1 F . We obtain:
da
k 12 K
=
a 1− a
B
+
b
dτ
k−1 K
A
[5.45]
db
k 2 F
k 21 K
=
b 1− b
A
+
a
dτ
k
1 F
k−2 KB
134
Untangling Complex Systems
Finally, fixing ρ = ( k
/ −
)
/ −
)
2 k
/ 1), r 12 = ( k 12 KB k 1 KA , and r 21 = ( k 21 KA k 2 KB , the two differential equations assume their final dimensionless form:
da = a(1− a+ r 12 b) = f
τ
d
[5.46]
db = ρ b(1− b+ r a
21 ) = g
dτ
The possible steady-state solutions are: ( a
) 0 0) (1 0) (0 )1 ((1 12) (1 12 21) (1 21)
ss , bss = ( ,
; , ; , ;
+ r
− r r , + r
1
( − r r ))
12 21
. As usual, their stability can be inferred by the linear analysis method. The Jacobian is
f
∂
f
∂
a
∂
b
∂
1− 2 ass + 12
r bss
12
r ass
J
ss
ss
=
=
[5.47]
g
∂
g
∂
ρ 21
r bss
ρ − 2ρ bss + ρ r a
21 ss
a
∂
b
∂
ss
ss
The steady state (0,0) is unstable, because both tr ( J ) and det ( J ) are positive.
The steady state (1,0) is a saddle point. In fact,
−1
r 12
J =
, tr ( J ) = ρ (1+ r 21) −1, det ( J ) = −ρ (1+ r 21) < 0 [5.48]
0
ρ (1+ ρ r 21)
The steady state (0,1) is also a saddle point. In fact,
1+ r 12
0
J =
, tr ( J ) = (1+ r 12 ) − ρ, det ( J ) = −ρ (1+ r 12 ) < 0 [5.49]
ρ r 21
−ρ
For the last steady state, we have
1+ r
12
r 12 ( r 12 + )
1
( r 12 r 21 − )
1
(1− r 12 r 21)
J =
ρ r
(
)
21 ( r 21 + )
1
ρ r 21 +1
(1− r
(
)
12 r 21 )
r 12 r 21 −1
1
( + r )
(
)
12 + ρ 1 + r
tr ( J
21
) =
(
[5.50]
r
)
12 r 21 −1
ρ 1
( + r )(
) ρ
(
)(
)
12
1+ r 21 − r 12 r 21 1+ r 12 1+ r
det ( J
21
) =
( r
)
12 r 21 −1 2
When r
( )
( )
12 r 21 < 1, det J
> 0 and tr J < 0. Therefore, the fourth steady state is stable. When r 12 r 21 > 1,
the fourth steady state does not lie in the positive quadrant of the a−b space and therefore, it does
not have physical meaning.
In synthesis, when r 12 r 21 > 1, the mutualism has large values for the kinetic constants k and in 12
k 21
the case of comparable carrying capacities. In such case, we have only three possible steady states
that are either unstable or saddle points. The populations of the two species grow unboundedly
The Emergence of Temporal Order in the Ecosystems
135
2
(i)
g
b 1
f
02 0
a 1
2
b
g
1
f
(ii)
0 0
a
1
2
FIGURE 5.9 Possible scenarios for a mutualistic relationship. Case (i) requires r 12 r 21 > 1; case (ii) holds when r 12 r 21 < 1 .
to infinity (see case (i) in Figure 5.9). On the other hand, when r 12 r 21 < 1, the mutualism is less strong. We have four steady states. The fourth one represents a stable solution, and all trajectories
of its phase space tend to ( a , b )
ss
ss = (>1, >1). In such steady state (see case (ii) in Figure 5.9), the
populations of the two species are larger than their maximum values achievable with a neutral
relationship.11
To see the dynamics of the populations of species A and B in the case of the other symbiotic
relationships, try exercises 5.7, 5.8, 5.9, and 5.10.
5.7 KEY QUESTIONS
• Describe the Lotka-Volterra model for the predator and prey relationship.
• Is the predator-prey relationship always oscillatory? Describe the dynamics in the case of
Type I, II and III functional responses.
• What does the Poincaré-Bendixson theorem state?
• Make examples of symbiotic relationships.
• Which are the possible outcomes of an antagonistic relationship?
• Why should human relationships be mutualistic?
5.8 KEY WORDS
Predator and prey; Functional and numerical responses; Symbiotic relationships.
5.9 HINTS FOR FURTHER READING
Besides Mathematical Biology by Murray (2002), additional books on mathematical modeling of
ecology are A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution by Otto and
Day (2007) and Elements of Mathematical Ecology by Kot (2001).
11 Ideally, any human relationship should be mutualistic. It is the only one that guarantees a reciprocal growth.
136
Untangling Complex Systems
5.10 EXERCISES
5.1. Apply the linear stability analysis to the improved Lotka-Volterra model where the growth
of the prey population follows the “logistic model” when there are not predators:
k
1 →
F + H ←
2 H
k−1
H + C
k
2→
C
2
C
k
3→
D
5.2. Apply the phase-space analysis to the dynamical system of exercise (5.1).
5.3. Solve the system of two differential equations [5.23] based on Type II functional response
and a numerical response proportional to the rate of prey consumption, using the follow-
ing values for the parameters (extracted from Harrison [1995] and regarding unicellular
Paramecium Aurelia and its unicellular predator Didinium nasutum):
K
k
1 F = 898; k 1 F = 1.85; k
12 4
. ; k
h = 25.5; K = 284.1; ′ h =
3 = 2.07.
Initial conditions: H
0 = 56 and C 0 = 23.
5.4. According to the results presented by Harrison (1995), when the prey Paramecium and the
predator Didinium are grown in a cerophyl medium thickened by methyl-cellulose, the
searching effectiveness of the Didinium dramatically reduces. Therefore, k becomes much
c
smaller and K much larger. The reduction of the nutrient supply for Paramecium may hinder
the growth rate k 1 F and diminish the carrying capacity K 1 F. For the following sets of parameters K
1 F = 400; k 1 F = 0.82; kh = 25.5; K = 350; k′ h = 12.4; k 3 = 2.07; H 0 = 87; C 0 = 10, determine the stationary state and its stability.
5.5. Consider the system of two differential equations [5.23], involving a Type II functional
response, and the values of the constants used in exercise 5.3:
K
898
1 85
25 5
284 1 k′
1 F =
; k 1 F = . ; kh = . ; K
=
. ; h = 12. ;
4 k
3 = .
2 0 .
7
Try to verify the Poincaré-Bendixson theorem by constructing a trapping region TR, i.e., a
closed set, such that the vector field points inward everywhere on its boundary.
5.6. Apply the linear stability analysis to confirm that for the system of two differential equa-
tions [5.32], there are three regions of stability as shown in Figure 5.7.
5.7. Predict the dynamics of a neutral relationship between two species, A and B, considering the general scheme [5.33].
5.8. Predict the dynamics of parasitism of species B towards species A. Considering the general scheme [5.33], fix n
A = 0 and nB = 2.
5.9. Predict the dynamics of amensalism between two species, A and B, assuming that
n
A = 0 and nB = 1 in the general scheme [5.33].
5.10. Predict the dynamics of commensalism involving two species, A and B, assuming that
n
A = 2 and nB = 1 in the general scheme [5.33].
The Emergence of Temporal Order in the Ecosystems
137
5.11 SOLUTIONS TO THE EXERCISES
5.1. The differential equations are
dH = k
2
1 FH − k−1 H
− k 2 HC
dt
dC = k 2 HC − k C
dt
3
The solutions of stationary states are:
H
(
2
2 1
−1 3 )
ss = k 3 k
/ 2 and Css = k k F − k k / k 2.
The Jacobian is:
k− k
− 1 3
−
k
k
3
J =
2
.
k− k
k
1 3
1 F −
0
k 2
tr ( J ) = − k− k
1 3 k 2 is always negative.
k
2
2
− k
1 3
k 1 Fk 3 k 2 k− k
det ( J ) = k
1 3
1 Fk 3 −
=
−
.
k 2
k 2
The eigenvalues are:
2
k
2
−
1
− k 3
k 1
− k 3
k k
1 3
2
−
±
4
− 1
k Fk −
tr ( J ) ± ( tr ( J )) − 4 det ( J )
k
k
3
k
λ
2
2
2
1,2 =
=
2
2
The determinant is positive if k 1 F > k− k
1 3 k 2 . If it is positive also the discriminant (Δ), the
eigenvalues are two real negative roots, and the stationary state is a stable node. When
Δ = 0, the stationary state is a stable star. Finally, if Δ < 0, the roots are two complex num-
bers having a negative real term. Therefore, the stationary state is a stable spiral, and the
predator and prey populations exhibit damped oscillations with the predator oscillations
lagging in phase behind the prey.
The determinant is null when k 1 F = k− k
1 3 k 2 . In this case, one eigenvalue is a negative real
number, and the other is 0: the system admits a line of stable fixed points. When the determi-
nant is negative, i.e., when k 1 F < k− k
1 3 k 2, the discriminant is positive, and the roots are real,
but one positive and the other negative. In this case, the stationary state is a saddle point.
5.2. In the system of two differential equations of exercise 5.1 there are only two variables, H, and C: dH = k
2
1 FH − k−1 H
− k 2 HC
dt
dC = k 2 HC − k C
dt
3
We can represent the evolution of the system in a phase space having C and H as coordi-
nates. The equations representing the null-clines are:
138
Untangling Complex Systems
C• = 0
( C is large and H is small)
H• = 0
H• < 0
( k 1 F/ k 2)
C• < 0
H• < 0
C
H• > 0
C• > 0
•
( H and C are large)
C < 0
H• > 0
( H and C are small)
C• > 0 ( H is large and C is small)
( k 3/ k 2)
H
( k 1 F/ k−1)
FIGURE 5.10 Phase plane for the “modified” Lotka-Volterra model when the det(J) is positive.
H = 0 when C = ( k 1 F − k−1 H ) k 2 representing a straight line with ( k
) as intercept
1 F/k 2
and –( k−
) as slope.
1 /k 2
C = 0 when H = ( k 3 k 2 ) representing a straight line parallel to the C axis.
When
k
( )
1 F > k− k
1 3 k 2 , i.e., det J
> 0, the phase space appears like that depicted in
Figure 5.10.
The intersection of the two nullclines is the stationary state that in the linear analysis of
stability we discovered it is stable (see exercise 5.1). This result is confirmed qualitatively
by the vector field shown in Figure 5.10. The signs of
H and
C in the different regions of
