The transition from evolutionary stability to branching: A catastrophic evolutionary shift

Evolutionary branching—resident-mutant coexistence under disruptive selection—is one of the main contributions of Adaptive Dynamics (AD), the mathematical framework introduced by S.A.H. Geritz, J.A.J. Metz, and coauthors to model the long-term evolution of coevolving multi-species communities. It has been shown to be the basic mechanism for sympatric and parapatric speciation, despite the essential asexual nature of AD. After 20 years from its introduction, we unfold the transition from evolutionary stability (ESS) to branching, along with gradual change in environmental, control, or exploitation parameters. The transition is a catastrophic evolutionary shift, the branching dynamics driving the system to a nonlocal evolutionary attractor that is viable before the transition, but unreachable from the ESS. Weak evolutionary stability hence qualifies as an early-warning signal for branching and a testable measure of the community’s resilience against biodiversity. We clarify a controversial theoretical question about the smoothness of the mutant invasion fitness at incipient branching. While a supposed nonsmoothness at third order long prevented the analysis of the ESS-branching transition, we argue that smoothness is generally expected and derive a local canonical model in terms of the geometry of the invasion fitness before branching. Any generic AD model undergoing the transition qualitatively behaves like our canonical model.

The consistency relations C1-C3 cannot determine the third and higher orders in the expansion (M2), because C3 gives only two constraints (C3a and C3b) among the k unknown functionss d k−d (w 1 , w 2 ), d = 1, . . . , k at order k. Again C2 is of no help in determining the unknown functions, it simply imposes the diagonal symmetry. To further constrain the coefficients of the expansion (M2) at order k ≥ 3, a specific class of ecological models must be considered to allow the direct computation of the directional derivatives (as done in Ref.20). Note that the monomorphic-dimorphic link is not fully exploited in C1, as it is valid also along the extinction boundaries (on which only one population is present). This is however of no help here, because the boundary in general is not straight (see Fig. 1e, f), so that, directionally, the link has only consequences at the singular point (x * , x * ).

Expansion of the resident-mutant coexistence region
The extinction boundary i of the resident-mutant coexistence region, along which only the resident x j is present (i = 1, 2, j = 2, 1) is defined by the invasion fitness of strategy x i being positive in the coexistence region and negative after crossing boundary i ( Fig. 1a-c).
The two boundaries are evidently related by the symmetry w.r.t. the diagonal x 1 = x 2 -one is obtained from the other by exchanging x 1 and x 2 in (M6)-so that below we focus on boundary 2. To approximate it locally to the singular point (x * , x * ), we rewrite it in polar coordinates (ε, θ ) as θ = θ 2 (ε), θ 2 (ε) being the function that gives the angle θ of the boundary point at distance ε from (x * , x * ). The function θ 2 (ε) is implicitly defined by Eq. (9) of the main text (the boundary definition in polar coordinates, reported below) which holds good for any (sufficiently small) ε ≥ 0. The approximation is in terms of an ε-expansion locally to ε = 0, i.e., θ 2 (ε) = θ 2 (0) + θ ′ 2 (0)ε + · · · + θ (k) 2 (0)ε k /k! + · · ·, to be used also for negative ε to describe the boundary across the diagonal x 1 = x 2 . The coefficients θ (k) 2 (0), k ≥ 0, of the expansion can be obtained by solving the ε-derivatives of Eq. (9) at ε = 0. The first derivative turns out to be the identity due to the fitness neutrality s x (x) = 0, whereas the second and third derivatives respectively solve for θ 2 (0) and θ ′ 2 (0). The result is reported in the Table (Extinction boundary 2), where only the monomorphic fitness derivatives with at least one order of derivation w.r.t. the mutant strategy are used-the pure x-derivatives ∂ x k s * are avoided by exploiting the fitness neutrality (see Table, Monomorphic fitness). In general, the k th -order coefficient θ (k) 2 (0) is determined by the monomorphic fitness derivatives up to order k + 2.
The first derivative θ ′ 2 (0) is nonzero close to the ESS-branching transition (under the coexistence and genericity conditions (2) and (7)) and determines the local curvature of the boundary-whether θ increases or decreases while moving away from (x * , x * ).

3/9 3 Expansion of the dimorphic invasion fitness
We now assume that the dimorphic fitness s x 1 ,x 2 (y) is three-times differentiable at (x 1 , x 2 , y) = (x * , x * , x * ), in the sense specified in Sect.1. As stated in the main text, we recall once more that this is so far shown to be the case only for the class of unstructured ecological models under stationary coexistence (done in Ref. 20 by direct computation of the directional expansion (M2)), though we expect the assumption to hold good for any class of ecological models yielding a smooth monomorphic fitness for a one-dimensional strategy. We hence assume the existence of the 3 rd -order local expansion for (∆x 1 , ∆x 2 ) within the coexistence region, and we determine the expansion's coefficients by applying the consistency relations C1-C3. Applying C2 implies the symmetry constraints s * We therefore eliminate the unknown coefficients with d 1 < d 2 and write the constraints implied by C1 and C3 in the 12 unknowns with d 1 ≥ d 2 at orders k = 1, 2, 3, plus the zero-order coefficient s * (see Table, Dimorphic fitness). So doing, we also eliminate the relation C3b, as it is implied by C2 and C3a.
Here is where we really need the differentiability of the dimorphic fitness. To exploit C1 ′ and constrain the coefficients of the expansion (M9), we need to impose the (ε, ∆y)-derivatives of C1 ′ at (ε, ∆y) = (0, 0). Such derivatives formally involve the partial derivatives of the dimorphic fitness at (x 1 , x 2 , y) = (x * , x * , x * ), that are not defined. Equivalently, we can substitute the truncated expansion (M9) in the left-hand side of C1 ′ and then differentiate. So doing, we obtain the following linear 4/9 constraints, organized by the order of the considered derivative: ε 2 : 2 2(∂ xy s * ) 2 + 2∂ xy s * ∂ y 2 s * + (∂ y 2 s * ) 2 s * 200 − 2(2∂ xy s * + ∂ y 2 s * )∂ y 2 s * s * 110 + (2∂ xy s * + ∂ y 2 s * )(∂ y 2 s * ) 2 = 0, (M12d) ε∆y : ε 2 ∆y : 6 2(∂ xy s * ) 2 + 2∂ xy s * ∂ y 2 s * + (∂ y 2 s * ) 2 s * 201 − 6(2∂ xy s * + ∂ y 2 s * )∂ y 2 s * s * 111 ε∆y 2 : 2∂ xy s * s * 102 + ∂ y 2 s * ∂ xy 2 s * = 0, (M12i) Of course the constraints in (M12) include those in (M10), obtained by C1 at order 0 and with the pure ∆y-derivatives. For the εand mixed-derivatives, the results at 1 st order (M12b,c) are exploited to obtain the constraints at 2 nd order and the results at 1 st and 2 nd orders (M12b-f) are exploited to obtain the constraints at 3 rd order. This allows to eliminate θ ′ 2 (0) at 2 nd order and θ ′′ 2 (0) at 3 rd order, so that only the expressions for θ 2 (0) and θ ′ 2 (0) (Table, Extinction boundary 2) are involved and substituted where needed. Also the coexistence condition (2) is taken into account to remove the denominators coming from θ 2 (0) and θ ′ 2 (0). The constraints implied by C1 ′ up to 2 nd order (M12a-f) are obviously redundant w.r.t. those implied by C1-C3 in (M10) and (M11). By contrast, each of the new 3 rd -order constraints (M12g-i) equivalently solves, together with (M11f-i), for the five 3 rd -order coefficients in (M9). The results are reported in the Table (Dimorphic fitness, 3 rd order). Note that, under our smoothness assumption for the dimorphic fitness, the 3 rd -order coefficients determine (according to the linear combinations in Table, Directional derivatives) the directional functionss 30 (w 1 , w 2 ),s 21 (w 1 , w 2 ),s 12 (w 1 , w 2 ) appearing at order 3 in the directional expansion (M2). The results indeed coincide with those found for the class of unstructured ecological models under stationary coexistence. 20 Unfortunately, the constraints implied by properties C1 ′ -C3 at 4 th order are not enough to solve for the 15 4 th -order coefficients of the expansion (M9). Out of the 16 constraints, only 14 are independent. In general, we have (k + 1)(k + 2)/2 coefficients at order k and the number of constraints, counting redundancies, is k + 1 for C1 ′ and C3a and (k/2)(k/2 + 1) (k even) or (k + 1) 2 /4 (k odd) for C2. Thus, even counting redundancies, the number of unknowns exceeds the number of constraints starting from order 6. This does not necessarily mean that the geometry of the dimorphic fitness, locally to the singularity (x 1 , x 2 , y) = (x * , x * , x * ), is not fully determined by the local geometry of the monomorphic fitness. The two fitness functions are linked to each other by the underlying ecological model, that is a much stronger link than C1 ′ . Only by exploiting the full ecological link we can then answer the above question, starting from order 4. Whether the answer is yes or no remains an open theoretical issue of AD.

The canonical model of the ESS-branching transition
Using the results derived in Sects. 2 and 3, we now derive the canonical model (8). We consider the continuous evolutionary dynamics ruled by the so-called AD canonical equation 11, 12 in the limit of rare and infinitesimally small mutational steps. We note however that the assumption of rare mutation can be relaxed 14 and that the dynamics of model (8) (the direction of evolution in Eq. (8a), the ecological scaling rates in Eqs. (8b,c), and the fitness gradients in Eqs. (8d,e)) inform as well about the adaptive dynamics driven by sufficiently small but finite strategy mutations. 23,24 In the simple setting of unstructured ecological models under stationary coexistence, the AD canonical equation readṡ before branching (monomorphic phase) anḋ after branching (dimorphic phase), where µ(x) and σ (x) 2 respectively scale with the frequency and variance of mutations in strategy x (half of which are at disadvantage and go extinct),n(x) andn i (x 1 , x 2 ) are the resident equilibrium densities (see 5/9 Fig. 1c), andẋ is the time-derivative of x on an evolutionary time scale that is fully separated from the ecological one. In more complex ecological settings-including structured population models and non-stationary attractors of coexistence-the fitness gradients in Eqs. (M13) and (M14) are scaled by the average birth outputs shown by the populations at the ecological regime, 11, 12, 18 in lieu of the resident equilibrium densities. This ecological scaling has been shown equivalent to four times the effective population size, 35 as defined in population genetics.
As an approximation of the stochastic dynamics driven by finite mutations, Eqs. (M13) and (M14) give better results the more mutations are small. The convergence of the stochastic dynamics to the deterministic solution, as mutational steps become infinitesimal, is however non-uniform, 12 i.e., closer to a singular strategy the approximation requires smaller mutations. As a consequence, branching must be discussed by resorting on finite mutations close to the singular strategy 36 (see also the discussion on finite mutations following Fig. 2 in the main text).
To derive the canonical model (8), we first get rid of the evolutionary scaling 1 2 µ(x i )σ (x i ) 2 in Eq. (M14). Locally to the singular point (x * , x * ) this can be done in two steps. A near-identity coordinate transformation, z i = z i (x 1 , x 2 ), i = 1, 2 (near-identity meaning ∂ z i /∂ x j = 1 if i = j, 0 otherwise), whose expansion can be set to eliminate all the derivatives of µ and σ in the expansion of the scaling factor around x 1 = x 2 = x * ; a time-scaling τ = 1 2 µ(x * )σ (x * ) 2 t, τ being the new time. For simplicity, we keep on using variables x i (actually ∆x i ) and t for the new variables and time.
Second, we approximate the ecological scaling factor in Eq. (M14). To avoid specific assumptions at the ecological level, the idea is to replace the resident equilibrium densityn i (x 1 , x 2 ) with a quantity that is determined by the geometry of the monomorphic fitness and that qualitatively behaves as the birth output of population i, i = 1, 2, locally to the singular point (x * , x * ). So doing, we lose the quantitative mapping between our canonical model and the AD canonical equation (M14), that can however be saved only working with a specific class of ecological models.
The birth output of population i is positive in the resident-mutant coexistence region and vanishes while approaching the extinction boundary i. It is therefore sign-equivalent, within the coexistence region (boundaries included), to the invasion fitness of strategy x i that defines boundary i in (M6). However, while the monomorphic fitness is smooth and quadratically vanishes at the singular point (x * , x * ), the birth output i is discontinuous at (x * , x * )-population i being absent on boundary i and present along boundary j, i = j.
The unfolding parameter of the ESS-branching transition-that we move across zero-is ∂ y 2 s * . Four other parameters are left in model (8): ∂ xy s * , ∂ x 2 y s * , ∂ xy 2 s * , ∂ y 3 s * , i.e., all second and third independent derivatives of the monomorphic fitness-the pure x-derivatives ∂ x 2 s * and ∂ x 3 s * being related to the others by the fitness neutrality (see Table, Monomorphic fitness). The first and last parameters are constrained by the coexistence and genericity conditions (2) and (7), whereas ∂ x 2 y s * and ∂ xy 2 s * play no role in the transition. In fact, though they both appear in η i (∆x 1 , ∆x 2 ) (see Eqs. (8b,c)) and only ∂ xy 2 s * appears in s i (∆x 1 , ∆x 2 ) (see Eqs. (8d,e)), at the transition (∂ y 2 s * = 0) the curvature of the extinction boundaries and the flow nullclines are unaffected (see the expression of θ ′ 2 (0) in Table, Extinction boundary 2, where the effect of ∂ x 2 y s * and ∂ xy 2 s * is modulated by ∂ y 2 s * , and note that ∂ xy 2 s * is multiplied by ∂ y 2 s * in s i (∆x 1 , ∆x 2 )), whereas the flow of model (8) is perturbed only in the cubic (∆x 1 , ∆x 2 )-terms. With ∂ xy s * constrained in sign, ∂ y 3 s * is the only relevant coefficient of the canonical model.
Finally, we note that model (8) serves as canonical model for the ESS-branching transition also in the presence of other evolving (one-dimensional) strategies (e.g. the transition to a higher polymorphism or the coevolution in a multi-species model). Whether the other strategies are at singularity or slowly evolving, 37 our analysis in Ref. 20, Appx. C, shows that they remain selectively neutral during the initial phase of branching in strategy x.

Unfolding of the ESS-branching transition
We analyze in this last section the dynamics of the canonical model (8), under the coexistence and genericity conditions (2) and (7), while varying the unfolding parameter ∂ y 2 s * across zero. Model (8) is defined only within the resident-mutant coexistence region defined byñ i (∆x 1 , ∆x 2 ) ≥ 0, i = 1, 2, see Eqs. (M16a,b) and (8b,c). However, for the purpose of the mathematical analysis, we extend the validity of the model's equations in a full neighborhood of (∆x 1 , ∆x 2 ) = (0, 0).
Note that E1 and E2 are symmetric boundary equilibria, respectively lying on the extinction boundaries 1 and 2, while E3 lies on the diagonal (and is therefore unfeasible for the evolutionary dynamics). The four equilibria are all involved in the ESS-branching transition occurring at ∂ y 2 s * = 0, as they collide at (0, 0) at the transition. Under conditions (2) and (7), equilibria E0-E3 intersect transversely as ∂ y 2 s * moves across zero. The transition classifies as a non-simple branch point bifurcation [38][39][40] (not to be confused with the branching point of AD!), i.e., the transversal intersection of more than two ∂ y 2 s * -parameterized equilibrium branches. This bifurcation generically requires the continuation problem 41 that defines the intersecting equilibrium branches to have a nullspace with dimension larger than two at the bifurcation. Specifically, the continuation problem is defined in the space (∆x 1 , ∆x 2 , ∂ y 2 s * ) by C(∆x 1 , ∆x 2 , ∂ y 2 s * ) := η 1 (∆x 1 , ∆x 2 , ∂ y 2 s * ) s 1 (∆x 1 , ∆x 2 , ∂ y 2 s * ) η 2 (∆x 1 , ∆x 2 , ∂ y 2 s * ) s 2 (∆x 1 , ∆x 2 , ∂ y 2 s * ) = 0, where ∂ y 2 s * is explicitly mentioned as an argument of η i and s i . The Jacobian of the (vector-valued) function C w.r.t.
(∆x 1 , ∆x 2 , ∂ y 2 s * ) is indeed a (2 × 3) null matrix at the bifurcation (easy to check), i.e., the nullspace is three-dimensional. Due to the diagonal symmetries of the canonical model (8), this bifurcation occurs with codimension-one, i.e., moving a single unfolding parameter (see Ref. 38, Sect. 8.2, for further detail). Two cases can be distinguished, namely ∂ y 3 s * > 0 and ∂ y 3 s * < 0, whose unfoldings are pictured in Fig. M1 (top and bottom panels, respectively). The movements and stability of the four equilibria, as ∂ y 2 s * goes from negative to positive, are evident