Appendix A
Stability analysis of the system without evolutionary adaptation
The local stability of the system is described by equation (1) when Gi = 0 by linearizing the dynamics near the nontrivial equilibrium. Local stability was judged by whether the characteristic equation of their Jacobian matrix satisfied the Routh–Hurwitz criteria. The nontrivial equilibrium was calculated as,
(A-1a)
(A-1b)
where asterisks relate to the equilibrium. The feasible condition when the equilibrium of the two species was positive was calculated as follows:
Case I
Case II
The two cases defined above correspond to facultative and obligate mutualism, respectively.
Under the equilibrium condition, the Jacobian matrix was calculated as follows:
.(A-2)
The characteristic equation for determining the eigenvalues is, where
(A-3a)
(A-3b)
The equilibrium point is locally stable if w1, w2 > 0. The equilibrium of Case I was always stable and that of Case II was always unstable.
Appendix B
Stability analysis of the system with evolutionary adaptation
The local stability of the system was described by equations (1) and (2) by linearizing the dynamics near the nontrivial equilibrium. Local stability was judged by whether the characteristic equation of their Jacobian matrix satisfied the Routh–Hurwitz criteria.
The Jacobian matrix (not shown because of its complexity) under the equilibrium condition was calculated as follows:
(B-1a)
(B-1b)
(B-1c)
(B-1d)
where primes denote the first derivatives. For the specific functions used in this study,and
(B-2)
For the convenience of the notations,.
By using (B-1c), (B-1d) and (B-2b), the equilibrium condition can be re-written as
(B-3)
Consider a special case where all parameters are symmetrical. In such a case, if u* = v*, it is clear that the equilibrium condition (B-3) is met.
The equations (B-1c) and (B-1d) must be positive because of the shapes of the functions (B-2b) which is met when the equilibrium population sizes are positive. In the symmetrical case,the feasible condition when the equilibrium of the two species is positive is as follows;
Case I
Case II
whereand . These two cases indicate that coevolution goes towards facultative mutualism when the benefits of mutualism are smaller, whereas it goes towards obligate mutualism when the benefits of mutualism are larger.
The characteristic equation for determining the eigenvalues isThe equilibrium point is locally stable if .and according to the Routh–Hurwitz criteria. The full stability condition is too complex to show except for a limited case where rX(u*) = rY(v*) = r, and In this case,
w1 =,
w2 =,
w3 =,
w4 =.
Linear cost function
When the cost function is linear (), all stability conditions except for are satisfied. The stability condition is then
(B-4)
Since , condition (B-4) is always met if mutualism at equilibrium is facultative. For obligate mutualism, the stability condition is reduced to
(B-5)
In the limit of GX GY, the left-hand side of condition (B-5) tends towards infinity, hence the condition can never be satisfied. Thus, in the symmetric case (GX=GY), the system can never be stable when mutualism is obligate.
Nonlinear cost function
When the cost function is nonlinear, some stability conditions can be demonstrated as w1 > 0, w3 > 0 and w4 > 0,
,(B-8a)
,(B-8b)
(B-8c)
These equations suggest that the shape of the cost function affects stability. When the cost function is concave, > 0, the system is likely to be unstable. If the concave function in two species in the model simulation was used, the system often collapsed, causing an overflow. For example, even when< 2 (i.e. r > 0), (B-8a) is not held if Gi and/or is large. When the cost function is convex, < 0, the system is likely to be stable. In fact, the first condition (B-8a) is always held when < 0. However, the other conditions also depend on whether is >2 or not. When < 2 (i.e. r > 0), the conditions (B-8b) and (B-8c) are always held. In contrast, when > 2 (i.e. r << 0), larger or Gi is also required for stability (see also Fig. 2).
Although the last stability condition was too complex to evaluate the stability condition, a necessary condition for stability can be shown as w2> 0.
Case1: < 2 (i.e. r > 0)
.(B-9a)
Case2: > 2 (i.e. r < 0)
.(B-9b)
These equations suggest that < 0 is required for stability and when < 0, the first condition is always held. The second condition may be held when the speed of evolution is faster and/or is larger.
Although I could not derive all stability conditions, the analytical results can be complemented by a numerical stability analysis (Fig. 2, Fig. D3). In addition, in Fig. D7 and D8, the numerical results for stability if the symmetry assumption for parameters is relaxed are demonstrated.
Appendix C
Stability analysis of the system with a slow evolutionary adaptation
When the evolutionary dynamics are much slower than the population dynamics (<1), the coevolutionary dynamics are approximated by averaged dynamics when the population dynamics without evolution is stable:
(C-1a)
,(C-1b)
whereWi(WX= rX() – X*+ and WY= rY() – Y*+ ) is the fitness of the mutant traits and , respectively.
The Jacobian matrix (not shown because of its complexity) of the system (C-1) under the equilibrium condition can be calculated as
(C-2a)
(C-2b)
The characteristic equation for determining the eigenvalues isThe equilibrium point is locally stable if w1, w2 > 0 according to the Routh–Hurwitz criteria. For simplicity, a symmetrical case was examined where rX(u*) = rY(v*) = r, and
The stability conditions are
(C-3a) and
(C-3b)
Because an assumed stable equilibrium is required in population dynamics without evolution in this analysis, the stability conditions are reduced to
(C-4a)
(C-4b)
This suggests that stability does not depend on the speed of evolution and that a strongly convex cost function and/or a large value of the shape parameter of the mutualism function tends to stabilize the system.
Appendix D
Supplemental Figures
Fig. D1. Function of mutualistic benefits. The figure is an example of species X. Different colours indicate the difference in parameter values of and
Fig. D2. Bifurcation diagrams of population and trait dynamics of species X (a) and Y (b) in relation to their ratios of adaptation speed. The maximum, mean and minimum oscillations values are plotted. GX is fixed at 0.05. Parameter values are the same as those in Fig. 1.
Fig. D3. The relationship between adaptation speed and local stability of the equilibrium. Within the shaded region, the non-trivial equilibrium is locally unstable. The stability is evaluated by the signs of real part dominant eigenvalues. Parameter values are ri = 1, = 1.99 and = 10.
Fig. D4. The relationship between adaptation speed and the amplitude of population oscillation of X varying with . Parameter values are the same as those in Fig. 1.
Fig. D5. The relationship between adaptation speed and the amplitude of population oscillation of X varying with Parameter values are the same as those in Fig. 1.
Fig. D6. Dynamics of population sizes and trait values depending on different initial settings. Dotted lines indicate the dynamics in a case where the initial values are close to equilibrium. Solid lines indicate the dynamics in a case where the initial values are not close to equilibrium. GX = 0.05 and GY = 0.1. Other parameter values are same as those in Fig 1.
Fig. D7. The relationship between adaptation speed and the amplitude of population oscillations of X varying with Dashed lines indicate the boundary which separates the local stability of the equilibrium. Note that the oscillation can occur even when the equilibrium is locally stable (left panel). Parameter values are same as those in Fig. 1.
Fig. D8. The relationship between adaptation speed and the amplitude of population oscillations of X varying with Parameter values are the same as those of Fig. 1 and other information is the same as in Fig. D7.
Fig. D9. The relationship between the asymmetry of between species and the amplitude of population oscillations of X (a) (b)Gi = 0.01. Other parameter values are the same as those in Fig. 1.
Fig. D10. The parameter dependence on the equilibrium of population sizes and trait values. Parameter values are ri = 1, = 1, = 30 and = 2 in (a); ri = 1, = 1, = 30 and = 2 in (b); ri = 1, = 1 and = 2 in (c). Black points in (c) indicates that the values are the same between species.