Trans-heteroclinic bifurcation: a novel type of catastrophic shift

Global and local bifurcations are extremely important since they govern the transitions between different qualitative regimes in dynamical systems. These transitions or tipping points, which are ubiquitous in nature, can be smooth or catastrophic. Smooth transitions involve a continuous change in the steady state of the system until the bifurcation value is crossed, giving place to a second-order phase transition. Catastrophic transitions involve a discontinuity of the steady state at the bifurcation value, giving place to first-order phase transitions. Examples of catastrophic shifts can be found in ecosystems, climate, economic or social systems. Here we report a new type of global bifurcation responsible for a catastrophic shift. This bifurcation, identified in a family of quasi-species equations and named as trans-heteroclinic bifurcation, involves an exchange of stability between two distant and heteroclinically connected fixed points. Since the two fixed points interchange the stability without colliding, a catastrophic shift takes place. We provide an exhaustive description of this new bifurcation, also detailing the structure of the replication–mutation matrix of the quasi-species equation giving place to this bifurcation. A perturbation analysis is provided around the bifurcation value. At this value the heteroclinic connection is replaced by a line of fixed points in the quasi-species model. But it is shown that, if the replication–mutation matrix satisfies suitable conditions, then, under a small perturbation, the exchange of heteroclinic connections is preserved, except on a tiny range around the bifurcation value whose size is of the order of magnitude of the perturbation. The results presented here can help to understand better novel mechanisms behind catastrophic shifts and contribute to a finer identification of such transitions in theoretical models in evolutionary biology and other dynamical systems.

Global and local bifurcations are extremely important since they govern the transitions between different qualitative regimes in dynamical systems. These transitions or tipping points, which are ubiquitous in nature, can be smooth or catastrophic. Smooth transitions involve a continuous change in the steady state of the system until the bifurcation value is crossed, giving place to a second-order phase transition. Catastrophic transitions involve a discontinuity of the steady state at the bifurcation value, giving place to first-order phase transitions. Examples of catastrophic shifts can be found in ecosystems, climate, economic or social systems. Here we report a new type of global bifurcation responsible for a catastrophic shift. This bifurcation, identified in a family of quasi-species equations and named as transheteroclinic bifurcation, involves an exchange of stability between two distant and heteroclinically connected fixed points. Since the two fixed points interchange the stability without colliding, a catastrophic shift takes place. We provide an exhaustive description of this new bifurcation, also detailing the structure of the replication-mutation matrix of the quasispecies equation giving place to this bifurcation. A perturbation analysis is provided around the bifurcation value. At this value the heteroclinic connection is replaced by a line of fixed points in the quasi-species model. But it is shown that, if the replication-mutation matrix satisfies suitable conditions, then, under a small perturbation, the exchange of heteroclinic connections is preserved, except on a tiny range around the bifurcation value whose size is of the order of magnitude of the perturbation. The results presented here can help to
Smooth transitions, found in some natural and social systems, cause soft changes between active and quiescent states, with a more easily reversed progressive deterioration [15]. These transitions can occur in chemical reactions, epidemic spreading, fixation of alleles in population genetics and computer virus propagation [16][17][18][19][20]. Typically, smooth transitions are governed by transcritical bifurcations, although other types of bifurcations such as pitchfork [21] or saddle-node bifurcations (whenever the node and the saddle collide at a point located at the boundary of the admissible domain, see e.g. [22]) can give place to continuous population declines until extinction is achieved. The transcritical bifurcation typically involves the collision of two fixed points and an exchange of stability between them [23].
Opposite to smooth shifts, other bifurcations can cause the so-called catastrophic transitions [9,10]. For instance saddle-node bifurcations (whenever the saddle and the node collide at the interior of the admissible domain and the orbits that where attracted by the node jump to a different, distant, attractor). This bifurcation, typical of catalytic systems [24,25], has been recently mapped in laboratory experiments with yeast [7]. Saddle-node bifurcations are also found in food-chain mathematical models [26] or in host-parasite theoretical models [27]. Also, global bifurcations such as saddle-nodes of periodic orbits [28], homoclinic bifurcations [29,30], or the so-called chaotic crises [31][32][33] (which are mainly created by tangential heteroclinic events [34]) can give place to discontinuous shifts. Catastrophic shifts are of extreme importance in nature since they involve a radical change in the steady state of a system once the bifurcation threshold is crossed [13,14]. Sudden collapses can occur in socio-ecological systems [35][36][37] as well as in ecosystems. For instance, lakes, fisheries, coral reefs or savannahs can suffer abrupt collapses due to small changes in environmental conditions [30,[38][39][40][41][42][43]. Once such collapse occurs, recovery can be extremely problematic [15,30].
In this paper, we carry out a general study of a new type of global bifurcation recently identified in the quasi-species equation [44]. This bifurcation, named as trans-heteroclinic, involves a catastrophic shift. Identifying the mechanisms responsible for catastrophic regime shifts and distinguishing them from their smoother counterparts constitutes a timely subject with enormous important applications. We start by introducing the mathematical model where this bifurcation was found to provide some background for what is followed. Next, we provide a description of this new type of transition involving a discontinuous shift. For values of the parameter below the critical one there exists a heteroclinic connection from a point, say Q 1 , to a distant point, say Q 2 . After crossing the critical value the heteroclinic connection is recovered, but the stability roles of Q 1 and Q 2 are exchanged. We provide the conditions (within the replicationmutation matrix) needed for this bifurcation to take place in the general quasi-species replicator model. Then, we develop a perturbation analysis around the bifurcation value.

Background
In the following lines we introduce the mathematical model that displays a trans-heteroclinic bifurcation. This model was obtained from the well-known Eigen's quasi-species equation [45] (see electronic supplementary material, §1), and it considers a heterogeneous pool of tumour cell phenotypes competing with healthy cells [44]. Assuming a constant population, Eigen's model can be represented with the general equation:ẋ  dynamic equations describing self-replicating species in a chemostat. A is an n × n matrix with all entries non-negative and corresponds, for Eigen's model, to the replication-mutation matrix, which includes all the terms tied to the replication and mutation rates of the sequences x i=1...n . The entries of matrix A, which determine the shape of the fitness landscape and how the quasi-species move across this landscape, can also be represented as a ij = f i Q ji , being f i the fitness of sequence i (i.e. replication rates) and Q ji the mutation matrix [46]. Φ is the so-called outflow term, given by a linear function of x that introduces competition and keeps the population constant so that u T x = 1 is preserved for all t. That is We remark that an eigenvector v of A with non-negative components and normalized such that u T v = 1, determines a fixed point of the quasi-species equation. The corresponding eigenvalue coincides with Φ(v). It is well known that explicit solutions can be obtained for Eigen's quasi-species model [47,48] by introducing a vector y ∈ R n defined as where the integration is done along a solution x(τ ). Then the equation for y becomesẏ = Ay, easy to solve, and x(t) is recovered as y(t)/u T y(t).
Equation (2.1) provides a general model for investigating replication-mutation dynamics. This equation can be adapted to investigate the dynamics of different systems such as the one displayed in figure 1, which is here introduced as a model example of the trans-heteroclinic bifurcation [44]. Specifically, this model has been recently used to explore the population dynamics of healthy cells (labelled x 0 ) competing with tumour cells phenotypes (x 1...7 ). These cells can carry mutations or present anomalies in genes involved in cell growth (compartment R in figure 1b, i.e. proto-oncogenes and tumour suppressor genes), in genes keeping genomic stability (compartment S), and in housekeeping genes (compartment H, see [44] for details). The model parameters are given by the replication rate of cells r > 0; by the increase in proliferation of tumour cells due to anomalies in compartment R, given by δ r > 0; by the tumour cells mutation rates 0 < μ < 1; and by the increase in mutation rates (δ μ ) of tumour cells with anomalies in compartment S, with 0 < δ μ < 1 − μ. As was done in [44], we introducer = r + δ r , μ = μ + δ μ , α 1 =μr/2 and α 2 = μr/2. According to the previous terms, the dynamical system introduced in [44] and built from equation (2.1) can be written in a compact form aṡ The model above displays a threshold behaviour at a critical mutation rate, μ c = δ r /(r + δ r ), resulting in a discontinuous or jump transition. This catastrophic shift, caused by a trans-heteroclinic bifurcation, is displayed in figure 1c, where the equilibrium populations of cells are plotted against μ. When μ < μ c the whole population is dominated by tumour cells. However, at μ = μ c a discontinuous transition occurs, resulting in an abrupt and catastrophic extinction of the tumour cells and a dominance of healthy cells (that persists for μ > μ c ). The topological changes of equation (2.2) are represented in figure 2, below (a), at (b), and above (c) the bifurcation. When μ < μ c , the fixed point (1, 0, . . . , 0) is unstable and the flow achieves a globally asymptotically stable fixed point that involves the dominance of tumour cell phenotypes (red circle in figure 2a). Just at the bifurcation value μ = μ c the heteroclinic connection between the two alternative states (dominance of tumour cells versus dominance of healthy cells) is a line of fixed points, and the asymptotic states largely depend on the initial conditions (figure 2b). Finally, once μ > μ c , the fixed point (1, 0, . . . , 0) becomes a global attractor (blue circle in figure 2c) due to the stability exchange between the two equilibria. Figure 2 also displays time series associated to each one of the different scenarios found below, at and above the bifurcation.
We note that the same catastrophic shift is found in a quasi-species mathematical model describing the dynamics of the so-called survival-of-the-flattest effect (see [49]). The trans-heteroclinic bifurcation for this dynamical system is discussed in electronic supplementary material, §2.
Next, we provide a description of the trans-heteroclinic bifurcation and the structure of the replication-mutation matrix giving place to this type of catastrophic shift.

Structure of the replication-mutation matrix to find the trans-heteroclinic bifurcation
We consider the general quasi-species model given by equation (2.1). According to the Perron-Frobenius theorem the matrix A has a positive dominant eigenvalue, λ 1 , such that all the other eigenvalues, λ j , j > 1, satisfy matrices, respectively, with 0 < k < n, then λ 1 has multiplicity 1 and there exists a non-trivial eigenvector, say v 1 , with eigenvalue λ 1 , which has all the components positive. Hence, if A is irreducible and we consider an initial vector y(0) = x(0) with non-zero component in v 1 , the solution tends to the fixed point associated to this eigenvector, and it is an attracting fixed point of the system.
We would like to have a line of fixed points. So, we should assume that A is reducible, with a block triangular structure as mentioned before. The matrix, after a suitable permutation, can have a 2 × 2 block structure or a higher number of blocks. Then it is possible to have a dominant eigenvalue λ 1 with algebraic multiplicity equal to 2. But if in the associated subspace the matrix has a non-diagonal Jordan structure (that is, the geometrical multiplicity is 1) there appear terms of the form t exp(λ 1 t) in the solution, which are the dominant ones (assuming that the coefficient is non-zero). In this case, one has again a single attracting fixed point.
The matrix A can depend continuously on several parameters such as replication and mutation rates. For concreteness we assume all of them fixed except one of them, say μ, and we shall consider A(μ). Then the eigenvalues depend on μ in a continuous way.
The next proposition gives sufficient conditions so that the quasi-species equation has a transheteroclinic bifurcation at some critical value μ = μ c .  x 4 x 0 x 4 x 6 x 7 x 5 x 0 x 0 x 4 x 4 x 6 x 7 x 5 x 6 x 7 x 5 x 0 x 4 x 6 x 7 (a) Orbits attracted by the fixed point P * 3 (red solid circle, which is stable with μ c > μ = 3 10 ) with x 0 = 0 (see [44] for notation) before the trans-heteroclinic bifurcation. (b) Dynamics at the bifurcation value μ = μ c = 1 3 , for which the final state depends on the initial conditions (the small red dots indicate the asymptotic population state). (c) Phase portrait above μ c , with μ c < μ = 11 30 . Here the orbits achieve the fixed point P * 2 with x 0 = 1 (see [44]), which is globally asymptotically stable (blue solid circle). Notice that in all panels, the trajectories are attracted by the heteroclinic connection. Stable (unstable) fixed points are indicated with solid (open) circles. On each panel one has taken nine initial conditions on top and below the line joining P * 2 and P * 3 . The time series in the second (resp. third) row correspond to initial data taken below (resp. on top) of the line joining P * 2 and P * 3 in the related phase portrait. Here the dynamics of healthy (x 0 variable, thick blue line) and tumour cells populations (thin lines in red, except for lethal phenotypes, indicated with blue thin lines) are represented. All of the variables are plotted although some of them appear overlapped or go rapidly to extinction. connection, from Q 2 (μ) to Q 1 (μ) if μ > μ c and from Q 1 (μ) to Q 2 (μ) if μ < μ c . In this case we say that a transheteroclinic bifurcation appears for the critical value μ c .
Proof. From now on all the eigenvectors v of A are assumed to be normalized, that is, u T v = 1.
Consider first the case μ = μ c . For any c real, cv 1 (μ c ) + (1 − c)v 2 (μ c ) is an eigenvector of A for the eigenvalue λ 1 := λ 1 (μ c ) = λ 2 (μ c ), giving rise to a fixed point. Then we obtain a line of fixed points. Moreover, for any initial condition x(0) with non-zero component in v 1 or v 2 , the solution can be written as where Y(t) is a linear combination of functions of the form e λ j t , for j > 2 (eventually multiplied by a polynomial in t if λ j is a multiple eigenvalue) and so, e −λ 1 t Y(t) tends to 0 as t goes to ∞. Therefore, if c 2 1 + c 2 2 > 0, the solution tends to a fixed point in the line before. Assume μ near μ c . The invariant line contains two fixed points Q 1 (μ) and Q 2 (μ) associated to v 1 (μ) and v 2 (μ) respectively. Given an initial condition on the line x = cv 1 , this solution tends to Q 1 (μ) as t goes to ∞ and to Q 2 (μ) when t goes to −∞. So, there is a heteroclinic connection from Q 2 (μ) to Q 1 (μ). In a similar way, if λ 1 (μ) < λ 2 (μ), a heteroclinic connection from Q 1 (μ) to Q 2 (μ) exists.

Remark 2.2.
Assume that for μ = μ c the matrix A has the structure B 0 C D , with all of the entries nonnegative and the dominant eigenvalues of B and D coinciding and equal to λ 1 . Then, for λ 1 there exists u 1 , eigenvector of B, and u 2 , eigenvector of D, with non-negative components. It is clear that v 1 = 0 u 2 is an eigenvector of A with non-negative components. To get a second independent eigenvector of A with eigenvalue λ 1 , one should look for v 2 = u 1 w 2 , where w 2 must satisfy Cu 1 + Dw 2 = λ 1 w 2 . In general, if such a vector w 2 exists, it could have some negative components. This is why we ask explicitly for the existence of v 2 with all components non-negative in condition (1).

Remark 2.3.
Beyond the condition on the structure of A, we require the condition that the dominant eigenvalues of B and D coincide. Furthermore, the existence of w 2 requires (λ 1 I − D)w 2 = Cw 1 . As λ 1 is a simple eigenvalue of D, the kernel of (λ 1 I − D) has dimension one. Hence, the condition reduces to ask that Cw 1 has no component in that kernel. Summarizing, the request reduces to two scalar conditions (beyond the positivity or zero values of the components of w 2 ). Therefore, the condition for the transheteroclinic bifurcation has codimension two in the equations of the family (2.1), if we assume that the matrix A has the structure given in remark 2.2, as it happens for the system (2) in the present study.

Remark 2.4.
A trans-heteroclinic bifurcation can also appear under more degenerate conditions. The matrix A could have a dominant eigenvalue with algebraic and geometric multiplicity equal to 3 or higher, or to have, for instance, a dominant eigenvalue with algebraic multiplicity 4 and geometric multiplicity 2, with two 2 × 2 Jordan blocs, such that there are two independent solutions of the form t exp(λ 1 t). This would also lead to a line of fixed points.
For simplicity in the exposition let us denote as L the line of fixed points which appears when μ = μ c , and L(μ) the invariant line for μ near μ c . Note that the linear part of the equation for variables y(t) do not depend at all on any variable of block z, while variables z(t) depend on both block variables y and z. The dynamical system given by equations (2.4) and (2.5) can be considered as a model of two different ecological communities which define each block. Considering that B, C and D need to have non-negative entries one can expect different possibilities for trans-heteroclinic bifurcations to occur in ecological communities in a chemostat. Keeping the competition term indicated by φ, and naming the entries of block C as c ij , some general cases where this bifurcation could be found can be listed:

Ecological interpretation of matrix A and extension to other biological systems
-Case c ij = 0 ∀i, j. Purely competitive communities. For this system, the dynamics corresponds to two communities with Malthusian replication competing for resources. Hence, according to our results, the trans-heteroclinic bifurcation could underlie the competitive exclusion principle ('Gause's principle') among communities in chemostats (see [50]). -Case c ij > 0 ∀i, j. Competitive communities with asymmetric linear facilitation. The case where all of the elements of matrix C are positive may correspond to a facilitation process (e.g. production of nutrients) carried out by community y on z (and not the other way around, i.e. commensalism). This facilitation can allow community z to better reproduce.
The previous two cases consider all of the elements c ij to be zero or positive to illustrate different types of ecological interactions among the species of each block (or community). All other cases given by all possible combinations of these values of the elements c ij might define different types of interactions: competition (c ij = 0) and/or facilitation (c ij > 0) between the species of two communities in chemostats' ecosystems.

The effect of a perturbation
The next question to analyse will be what happens to the solutions of the initial quasi-species model when a sufficiently small, arbitrary and smooth perturbation is added to the system. Let us consider the perturbed systemẋ where, as before, A depends on a parameter μ and G is assumed to be normalized (e.g. G = 1) and with bounded differential DG < M for some fixed M > 0. The bounds are assumed to hold in the full simplex n i=1 x i = 1, x j ≥ 0, j = 1, . . . , n. Anyway, for a local study of the trans-heteroclinic bifurcation, it would be enough to assume that these bounds hold in a vicinity of L.
Note that if we want to keep the condition u T x = 1, the vector field G(x) must be replaced by G(x) − u T G(x)x and to be normalized after this modification.
First we analyse the stability properties of the fixed points for the unperturbed vector field F(x) = Ax − Φ(x)x, under the hypothesis of proposition 2.1. We shall denote the eigenvalues of A by λ i , skipping the explicit dependence on μ for simplicity.
To prove the claim we perform a change of variables to put A in Jordan normal form. Let z = P −1 x where P is an n × n matrix such that P whereΦ(z) = Φ(Pz) = u T APz = u T PJz. Using that the first two columns of P are the normalized eigenvectors v 1 , v 2 we get that the first two components of u T P are equal to 1. Thereforê for some constants β j , j = 3, . . . , n. The differential of the vector field in z is M(z) := J −Φ(z)I n − zDΦ(z) where I n denotes the identity in dimension n.
In the new variable z the coordinates of Q 1 are z 1 = 1, z j = 0, j = 2, . . . , n. It is easy to check that at Q 1 the matrix M becomes an upper triangular matrix ⎛ and then the eigenvalues are −λ 1 , The eigenvalues at Q 2 are computed in the same way taking into account that Q 2 has coordinates z 1 = 0, z 2 = 1, z j = 0, j = 3 . . . , n.
When we return to the value μ = μ c with λ 2 = λ 1 , instead of a single fixed point one has the line of fixed points L. A fixed point Q in that line has z coordinates as z 1 = α, Therefore, the eigenvalues of the differential are 0, −λ 1 , λ j − λ 1 , j = 3, . . . , n. Beyond the 0 due to the existence of the line of fixed points, L, all the other eigenvalues have negative real part if the assumptions of proposition 2.1 hold. The previously stated facts say that L is a normally hyperbolic invariant manifold (hereafter NHIM). That is, at the points of L the space R n splits as a one-dimensional space associated to the eigenvalue 0 (i.e. in the direction of L) and an (n − 1)-dimensional space associated to all the eigenvalues with negative real part.
It is well known (see, e.g. [51] or [52]) that NHIM persist under sufficiently small perturbation. This shows that for the system F + εG there exists an invariant curveL ε , close to L and also normally hyperbolic. And the same happens if one takes into account variations of μ close to μ c . We summarize these results as follows: Proposition 2.5. Let us consider the system (2.6) with A(μ) satisfying the hypothesis of proposition 2.1, the normalization of G and a bound on DG , as mentioned. Then, if ε is sufficiently small, for μ in a vicinity of μ c , there exists an attracting invariant curveL ε (μ), close to L(μ).

Robustness of the trans-heteroclinic bifurcation under perturbations
The final question to consider is what happens to the dynamics in the new invariant attracting curvê L ε (μ) for μ near μ c . That is, we want to see the robustness of the trans-heteroclinic bifurcation when a perturbation is added to the quasi-species model. Proposition 2.6. Let us consider the system (2.6) with A(μ) satisfying the hypothesis of proposition 2.1. Let μ 1 = μ c − and μ 2 = μ c + for some > 0 small enough. Assume that K < |λ 1 (μ j ) − λ 2 (μ j )|, j = 1, 2, for some given positive K.
Remark 2.7. For a small fixed value of ε and small enough, that is for μ close enough to μ c , the dynamics on the invariant curveL ε (μ) can have other fixed points or, eventually, no fixed points at all.
Proof. For simplicity we skip the explicit dependence on μ. If ε = 0 and μ = μ 2 , the fixed points Q 1 and Q 2 are the only ones which remain from the line of fixed points L. As has been mentioned, the eigenvalues of DF(Q 1 ) are λ i − λ 1 for i = 2, 3, . . . , n and −λ 1 , showing that Q 1 is an attractor. The eigenvalues at DF(Q 2 ) are λ i − λ 2 for i = 1, 3, . . . , n and −λ 2 , showing that Q 2 has a one-dimensional unstable manifold. It gives the connection. The same happens for μ = μ 1 , exchanging the role of the points. . Blue (red) color corresponds to unstable (stable) fixed points. Hence, the dynamics on the invariant line, projected on the x 1 variable and for every fixed value of μ, is going away from the points in blue and tending to the points in red. Several bifurcations such as saddle-node (s-n) or transcritical are indicated. Now let us consider the case ε = 0. By assumption all the eigenvalues of DF(Q 1 ) and DF(Q 2 ) are bounded away from zero by K , for > 0 small enough. Here we use also the fact that the eigenvalues λ 1 and λ 2 are at a finite distance of λ j , J ≥ 3, according to assumption (2) in proposition 2.1. Hence, the differential matrices are regular, and the implicit function theorem applied to F + εG ensures that there exist single fixed points Q 1 (ε) and Q 2 (ε) for F + εG if ε is sufficiently small, in some neighbourhoods around Q 1 and Q 2 . In the part ofL ε (μ) between Q 1 (ε) and Q 2 (ε) which excludes these neighbourhoods the field F is bounded away from zero and, hence, no new fixed points can appear for F + εG on it when ε = 0 is sufficiently small because of the bound on G . Furthermore, because of the bound on DG the dominant eigenvalues of the differential at Q 1 (ε) and Q 2 (ε) change only slightly and the dimensions of the manifolds at these points, are preserved if ε is sufficiently small.
Finally, let us consider the case close to zero (including the case μ = μ c , i.e. = 0). By the properties of NHIM under perturbation, the line L will be replaced by a nearby (at distance O(ε)) invariant curveL ε . The dynamics onL ε depends on the perturbation G and can be of any of the types which appear in onedimensional flows: all the points move (no fixed points), either in one direction or the other; several fixed points with alternative stability properties show up; some points can be at a saddle-node bifurcation, and the like.
As a conclusion, except for a narrow interval of μ around μ c whose size is of the order of ε, the change from heteroclinic connection from Q 2 (ε) to Q 1 (ε) to a heteroclinic connection from Q 1 (ε) to Q 2 (ε), persists. As an example we look at the model considered here given by equations (2.2) and perturb with a vector field εG such that all the components of G are zero except the first one, which is taken as g(x 1 ) = cos(15x 1 ). As done in some examples in [44] we take the parameter values r = 0.1 and δ r = 0.05, which give μ c = 1 3 and δ μ = 0.05. For coherence with the present notation we rename the variables x 0 , . . . , x 7 used in [44] as x 1 , . . . , x 8 , we also assume that the condition u T x = 1 is satisfied for all time if it is satisfied at t = 0 and we use the notation P * 2 (ε), P * 3 (ε) for the fixed points. We recall that for ε = 0 the first components of P * 2 and P * 3 are, respectively, x 1 = 1 and x 1 = 0. The perturbation G can be modified, as explained before, to satisfy u T x = 1 for all time or, alternatively, we redefine the function Φ in the equations as Φ = r(x 1 + x 3 ) +r(x 5 + x 7 ) + εg(x 1 ), wherer = r + δ r . In particular this implies that P * 2 is not changing when μ and ε are changed. Figure 3 shows the fixed points in this case as a function of μ, setting ε = 10 −4 . Specifically, this figure shows that for μ − μ c < −0.003366 there are only two fixed points P * 2 and P * 3 (ε) inL ε and, so, there is only a heteroclinic connection from P * 2 to P * 3 (ε). This connection subsists until μ − μ c ≈ −0.00161, together with other new fixed points. One can also see several saddle-node (s-n) bifurcations and a transcritical bifurcation at x 1 = 1. Several solutions have no meaning in the context of the biological model because some of the variables can be negative or greater than 1. It is easy to check (both theoretically and numerically) that the size of the domain in μ for which the heteroclinic connection between P * 2 and P * 3 (ε) do not exists is O(ε). Figure 4 displays for μ = μ c the fixed points on the lineL ε which in the present example and for μ = μ c coincides with L. The plot shows the variables (x 1 , x 8 ). As before the blue (red) points correspond to unstable (stable) fixed points. Only six of them play a role if we are interested in x 1 ∈ [0, 1]. From every unstable point there are heteroclinic connections to the two nearby stable points.

Discussion
In a recent article [44], we identified a novel bifurcation involving a catastrophic shift. This bifurcation, which is global and was named as a trans-heteroclinic bifurcation, was found in the quasi-species equation [45], used to investigate the dynamics of competition between healthy and tumour heterogeneous cells given by different phenotypes [44]. Starting from this context, and to provide a complete description of this new bifurcation, we have here developed a thorough analysis of the mechanisms responsible for this bifurcation, paying special attention to the structure of the replication-mutation matrix of Eigen's quasi-species model where this type of bifurcation is found.
The trans-heteroclinic bifurcation causes a catastrophic shift since the exchange of stability between two fixed points is produced without their collision, in contrast with the transcritical and the saddlenode bifurcations [23]. Hence, systems with trans-heteroclinic transitions are not bistable because each one of the two equilibria, when stable, is globally asymptotically stable. We have also analysed the effect of small perturbations at the bifurcation value, where a line of fixed points lives within the simplex.
For a general ordinary differential equation, the fact that it has a line of fixed points is a codimension infinity bifurcation, i.e. is extremely degenerate. However, in the present case, inside the quasi-species models and assuming that the matrix A is irreducible, is just a codimension two bifurcation (see remark 2.4).
We have noticed that this type of transition is also found in a quasi-species mathematical model describing the so-called survival-of-the-flattest effect (see [49]). The survival-of-the-flattest effect, described experimentally in viroids [53] and in competing computer programs [54], involves the survival of slow-reproducing replicons that are robust against deleterious mutations, and the vanishing of fastreproducing ones that are not robust to deleterious mutations. Specifically, the model studied in [49] presented a discontinuous transition when increasing the mutation rate that involved the outcompetition of the fast-replicating, fit quasi-species by the slow-replicating, flat one. The replication-mutation matrix for this model fulfils the conditions for the presence of the trans-heteroclinic bifurcation discussed in this article (see electronic supplementary material, §2).
It is known that transients typically suffer long delays at the vicinity of bifurcations. For instance, near a saddle-node bifurcation occurring at a critical value μ c and assuming that these fixed points exist for μ < μ c , the delay, named delayed transition, is τ = O((μ c − μ) −1/2 ). This scaling law, identified in different nonlinear systems such as charge density waves and hypercycles [23,25,55], has been experimentally characterized in electronic circuits [2]. For transcritical bifurcations one has τ = O(|μ − μ c | −1 ), and the delay is similar for a pitchfork bifurcation, when an attracting point, until μ = μ c , becomes unstable for μ > μ c and two nearby stable points appear. In this case τ = O((μ − μ c ) −1 ). Pitchfork bifurcations have been recently identified in experiments with elastic structures [4]. The transheteroclinic bifurcation has the same type of behaviour: τ = O(|μ − μ c | −1 ), both for μ < μ c and for μ > μ c . The derivation follows immediately from [44] (see also electronic supplementary material, §3). Summarizing, the bifurcation described in this article, identified in a family of quasi-species equations [45], involves a new type of catastrophic (abrupt) transition arising within the framework of evolutionary biology. The quasi-species model provides a very powerful mathematical framework to investigate the evolutionary dynamics of biological information. Initially conceived to explore the dynamics of prebiotic replicons [45], quasi-species theory has been widely used to investigate the dynamics of RNA viruses [56][57][58][59][60], cancer evolution [44,61,62] and compositional assemblies [63]. The results of this article provide the conditions for which catastrophic shifts due to trans-heteroclinic bifurcations should be expected in quasi-species systems, also introducing a novel type of bifurcations in nonlinear dynamical systems. Indeed, our results can be extended to any system with (exponentially) self-replicating species in a chemostat due to the generality of equation (