1.. Introduction {#s1}
================

Bifurcations, both local and global, govern transitions between qualitatively different states in dynamical systems. Bifurcations have been experimentally reported in physics \[[@RSOS171304C1]--[@RSOS171304C4]\], chemistry \[[@RSOS171304C5],[@RSOS171304C6]\] and biology \[[@RSOS171304C7],[@RSOS171304C8]\]. Bifurcations in biological dynamical systems can govern extinction processes and population collapses \[[@RSOS171304C7],[@RSOS171304C9]--[@RSOS171304C14]\]. Such extinctions can be smooth or catastrophic when one or more bifurcation parameters are tuned.

Smooth transitions, found in some natural and social systems, cause soft changes between active and quiescent states, with a more easily reversed progressive deterioration \[[@RSOS171304C15]\]. These transitions can occur in chemical reactions, epidemic spreading, fixation of alleles in population genetics and computer virus propagation \[[@RSOS171304C16]--[@RSOS171304C20]\]. Typically, smooth transitions are governed by transcritical bifurcations, although other types of bifurcations such as pitchfork \[[@RSOS171304C21]\] 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. \[[@RSOS171304C22]\]) 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 \[[@RSOS171304C23]\].

Opposite to smooth shifts, other bifurcations can cause the so-called catastrophic transitions \[[@RSOS171304C9],[@RSOS171304C10]\]. 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 \[[@RSOS171304C24],[@RSOS171304C25]\], has been recently mapped in laboratory experiments with yeast \[[@RSOS171304C7]\]. Saddle--node bifurcations are also found in food-chain mathematical models \[[@RSOS171304C26]\] or in host--parasite theoretical models \[[@RSOS171304C27]\]. Also, global bifurcations such as saddle--nodes of periodic orbits \[[@RSOS171304C28]\], homoclinic bifurcations \[[@RSOS171304C29],[@RSOS171304C30]\], or the so-called chaotic crises \[[@RSOS171304C31]--[@RSOS171304C33]\] (which are mainly created by tangential heteroclinic events \[[@RSOS171304C34]\]) 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 \[[@RSOS171304C13],[@RSOS171304C14]\]. Sudden collapses can occur in socio-ecological systems \[[@RSOS171304C35]--[@RSOS171304C37]\] as well as in ecosystems. For instance, lakes, fisheries, coral reefs or savannahs can suffer abrupt collapses due to small changes in environmental conditions \[[@RSOS171304C30],[@RSOS171304C38]--[@RSOS171304C43]\]. Once such collapse occurs, recovery can be extremely problematic \[[@RSOS171304C15],[@RSOS171304C30]\].

In this paper, we carry out a general study of a new type of global bifurcation recently identified in the quasi-species equation \[[@RSOS171304C44]\]. 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 replication--mutation 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.

2.. Results {#s2}
===========

2.1.. Background {#s2a}
----------------

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 \[[@RSOS171304C45]\] (see electronic supplementary material, §1), and it considers a heterogeneous pool of tumour cell phenotypes competing with healthy cells \[[@RSOS171304C44]\]. Assuming a constant population, Eigen's model can be represented with the general equation: $$\overset{˙}{\textbf{x}} = A\textbf{x} - \mathit{\Phi}(\textbf{x})\textbf{x},$$where $\textbf{x} = (x_{1},\ldots,x_{n}) \in \mathbb{R}^{n}$ is the population vector with all components non-negative and normalized initially by **u**^*T*^**x**(0)=1 where $\textbf{u}^{T} = (1,1,\ldots,1) \in \mathbb{R}^{n}$. Equation ([2.1](#RSOS171304M2x1){ref-type="disp-formula"}) can be further generalized to the 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 \[[@RSOS171304C46]\]. *Φ* 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 $\mathit{\Phi}(\textbf{x}) = \textbf{u}^{T}A\textbf{x} = \sum\limits_{j = 1}^{n}(\sum\limits_{i = 1}^{n}a_{i,j})x_{j}$.

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 \[[@RSOS171304C47],[@RSOS171304C48]\] by introducing a vector $\textbf{y} \in \mathbb{R}^{n}$ defined as $$\textbf{y} = \exp(f(t))\textbf{x},\quad f(t):=\int_{0}^{t}\mathit{\Phi}(\textbf{x}(\tau))\, d\tau,$$where the integration is done along a solution **x**(*τ*). Then the equation for **y** becomes $\overset{˙}{\textbf{y}} = A\textbf{y}$, easy to solve, and **x**(*t*) is recovered as **y**(*t*)/**u**^*T*^**y**(*t*).

Equation ([2.1](#RSOS171304M2x1){ref-type="disp-formula"}) 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](#RSOS171304F1){ref-type="fig"}, which is here introduced as a model example of the trans-heteroclinic bifurcation \[[@RSOS171304C44]\]. 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 1](#RSOS171304F1){ref-type="fig"}*b*, i.e. proto-oncogenes and tumour suppressor genes), in genes keeping genomic stability (compartment *S*), and in housekeeping genes (compartment *H*, see \[[@RSOS171304C44]\] 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 \[[@RSOS171304C44]\], we introduce $\hat{r} = r + \delta_{r},$ $\hat{\mu} = \mu + \delta_{\mu}$, $\alpha_{1} = \hat{\mu}r/2$ and $\alpha_{2} = \mu\hat{r}/2$. According to the previous terms, the dynamical system introduced in \[[@RSOS171304C44]\] and built from equation ([2.1](#RSOS171304M2x1){ref-type="disp-formula"}) can be written in a compact form as $$\left. \begin{matrix}
{\overset{˙}{x}}_{0} & {= x_{0}(r - \mathit{\Phi}),} \\
{\overset{˙}{x}}_{2} & {= x_{2}((1 - \hat{\mu})r - \mathit{\Phi}),} \\
{\overset{˙}{x}}_{4} & {= x_{4}((1 - \mu)\hat{r} - \mathit{\Phi}),} \\
{\overset{˙}{x}}_{6} & {= x_{6}((1 - \hat{\mu})\hat{r} - \mathit{\Phi}) + \alpha_{1}x_{2} + \alpha_{2}x_{4},} \\
{\overset{˙}{x}}_{1} & {= - x_{1}\mathit{\Phi},} \\
{\overset{˙}{x}}_{3} & {= \alpha_{1}x_{2} - x_{3}\mathit{\Phi},} \\
{\overset{˙}{x}}_{5} & {= \alpha_{2}x_{4} - x_{5}\mathit{\Phi},} \\
{\qquad\text{and}\qquad\quad{\overset{˙}{x}}_{7}} & {= \hat{\mu}\hat{r}x_{6} - x_{7}\mathit{\Phi},} \\
\end{matrix} \right\}$$where $\mathit{\Phi} = r(x_{0} + x_{2}) + \hat{r}(x_{4} + x_{6})$. 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 1](#RSOS171304F1){ref-type="fig"}*c*, 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~). Figure 1.The trans-heteroclinic bifurcation in a cancer quasi-species model. (*a*) Schematic diagram of the model, in which different tumour phenotypes compete with healthy cells (see \[[@RSOS171304C44]\]). (*b*) Structure of the sequence space for the cancer phenotypic model. Here tumour phenotypes include cells with anomalies in replication-related genes (proto-oncogenes and tumour suppressor genes, compartment *R*); anomalies in genes preserving genomic stability (compartment *S*); and anomalies in housekeeping genes (compartment *H*). Alterations in these compartments are assumed to define different cancerous phenotypes. (*c*) Catastrophic transition governed by the trans-heteroclinic bifurcation when increasing the mutation rate of tumour cells above a critical threshold *μ*~c~ \[[@RSOS171304C44]\]. This bifurcation separates two different scenarios: (i) persistence of tumour cells and vanishing of the healthy ones when *μ*\<*μ*~c~ and (ii) extinction of tumour cells and dominance of healthy cells when *μ*\>*μ*~c~.

The topological changes of equation ([2.2](#RSOS171304M2x2){ref-type="disp-formula"}) are represented in [figure 2](#RSOS171304F2){ref-type="fig"}, 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 2](#RSOS171304F2){ref-type="fig"}*a*). 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 2](#RSOS171304F2){ref-type="fig"}*b*). Finally, once *μ*\>*μ*~c~, the fixed point (1,0,...,0) becomes a global attractor (blue circle in [figure 2](#RSOS171304F2){ref-type="fig"}*c*) due to the stability exchange between the two equilibria. [Figure 2](#RSOS171304F2){ref-type="fig"} also displays time series associated to each one of the different scenarios found below, at and above the bifurcation. Figure 2.Topological changes associated to the trans-heteroclinic bifurcation. Phase portraits for system equations ([2.2](#RSOS171304M2x2){ref-type="disp-formula"}) with *r*=0.1, *δ*~*r*~=0.05 and *δ*~*μ*~=0.05 displayed in the projection (*x*~0~,*x*~4~). (*a*) Orbits attracted by the fixed point $P_{3}^{\ast}$ (red solid circle, which is stable with $\mu_{c} > \mu = \frac{3}{10}$) with *x*~0~=0 (see \[[@RSOS171304C44]\] for notation) before the trans-heteroclinic bifurcation. (*b*) Dynamics at the bifurcation value $\mu = \mu_{c} = \frac{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 $\mu_{c} < \mu = \frac{11}{30}$. Here the orbits achieve the fixed point $P_{2}^{\ast}$ with *x*~0~=1 (see \[[@RSOS171304C44]\]), 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}^{\ast}$ and $P_{3}^{\ast}$. The time series in the second (resp. third) row correspond to initial data taken below (resp. on top) of the line joining $P_{2}^{\ast}$ and $P_{3}^{\ast}$ 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.

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 \[[@RSOS171304C49]\]). 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.

2.2.. Structure of the replication--mutation matrix to find the trans-heteroclinic bifurcation {#s2b}
----------------------------------------------------------------------------------------------

We consider the general quasi-species model given by equation ([2.1](#RSOS171304M2x1){ref-type="disp-formula"}). 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 λ~1~≥*Re*(λ~*j*~). Furthermore, if the matrix *A* is irreducible, that is, there does not exist a permutation matrix *P* so that $PAP^{T} = \begin{pmatrix}
B & 0 \\
C & D \\
\end{pmatrix}$, where *B* and *D* are square *k*×*k* and (*n*−*k*)×(*n*−*k*) 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(\lambda_{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 trans-heteroclinic bifurcation at some critical value *μ*=*μ*~c~.

Proposition 2.1*Assume that the matrix* *A*(*μ*), *with all entries non-negative, depends on a parameter* *μ*. *Assume that there exists a critical value* *μ*~c~ *such that the following conditions hold*: (1) *For* *μ*=*μ*~c~ *the matrix* *A*(*μ*) *has a dominant eigenvalue with algebraic and geometric multiplicity equal to 2* $$\lambda_{1}(\mu_{c}) = \lambda_{2}(\mu_{c}) > \text{Re}(\lambda_{j}(\mu_{c})),\quad j = 3,\ldots,n$$*and two independent eigenvectors* ***v***~1~(*μ*~c~),***v***~2~(*μ*~c~) *for the eigenvalue* λ~1~(*μ*~c~)*, with non-negative components.*(2) *For* *μ* *close to* *μ*~c~ $$\lambda_{1}(\mu) > \lambda_{2}(\mu) > \text{Re}(\lambda_{j}(\mu)),\quad j = 3,\ldots,n,\quad\text{if}\,\mu > \mu_{c}$$*and* $$\lambda_{2}(\mu) > \lambda_{1}(\mu) > \text{Re}(\lambda_{j}(\mu)),\quad j = 3,\ldots,n,\quad\text{if}\,\mu < \mu_{c}.$$*Then, for any* *μ* *in a small neighbourhood of* *μ*~c~, *there is an attracting invariant line* *L*(*μ*) *such that for* *μ*=*μ*~c~ *all points in the line are fixed, and for* *μ*≠*μ*~c~*, there are in* *L*(*μ*) *two fixed points* *Q*~1~(*μ*),(*μ*)*, and a heteroclinic 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 trans-heteroclinic 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.Let **v**~1~(*μ*),**v**~2~(*μ*) be eigenvectors for λ~1~(*μ*),λ~2~(*μ*), respectively. The line defined by **x**=*c***v**~1~(*μ*)+(1−*c*)**v**~2~(*μ*), $c \in \mathbb{R}$ is invariant. Notice that if we take initial conditions in that line, that is, **x**(0)=*c***v**~1~(*μ*)+(1−*c*)**v**~2~(*μ*) then the solution is given by $$\textbf{x}(t) = \frac{c\, e^{\lambda_{1}t}\textbf{v}_{1}(\mu) + (1 - c)\, e^{\lambda_{2}t}\textbf{v}_{2}(\mu)}{c\, e^{\lambda_{1}t} + (1 - c)\, e^{\lambda_{2}t}}.$$Clearly it belongs to the line for any real *t*. Furthermore, if **v**~1~(*μ*),**v**~2~(*μ*) have non-negative components and *c*∈\[0,1\], then **x**(*t*) has non-negative components for all *t*.Consider first the case *μ*=*μ*~c~. For any *c* real, *c***v**~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 $$\textbf{x}(t) = \frac{c_{1}\textbf{v}_{1}(\mu_{c}) + c_{2}\textbf{v}_{2}(\mu_{c}) + e^{- \lambda_{1}t}\textbf{Y}(t)}{c_{1} + c_{2} + e^{- \lambda_{1}t}\textbf{u}^{T}\textbf{Y}(t)},$$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 $\infty$. Therefore, if $c_{1}^{2} + 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**=*c***v**~1~(*μ*)+(1−*c*)**v**~2~(*μ*), 0\<*c*\<1, the solution is ([2.3](#RSOS171304M2x3){ref-type="disp-formula"}). If λ~1~(*μ*)\>λ~2~(*μ*), this solution tends to *Q*~1~(*μ*) as *t* goes to $\infty$ and to *Q*~2~(*μ*) when *t* goes to $- \infty$. 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.2Assume that for *μ*=*μ*~c~ the matrix *A* has the structure $\begin{pmatrix}
B & 0 \\
C & D \\
\end{pmatrix}$, with all of the entries non-negative 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 $\textbf{v}_{1} = \begin{pmatrix}
0 \\
\textbf{u}_{2} \\
\end{pmatrix}$ is an eigenvector of *A* with non-negative components. To get a second independent eigenvector of *A* with eigenvalue λ~1~, one should look for $\textbf{v}_{2} = \begin{pmatrix}
\textbf{u}_{1} \\
\textbf{w}_{2} \\
\end{pmatrix}$, where **w**~2~ must satisfy *C***u**~1~+*D***w**~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.3Beyond 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~=*C***w**~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 *C***w**~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 trans-heteroclinic bifurcation has codimension two in the equations of the family ([2.1](#RSOS171304M2x1){ref-type="disp-formula"}), if we assume that the matrix *A* has the structure given in remark [2.2](#RSOS171304ST2){ref-type="statement"}, as it happens for the system (2) in the present study.

Remark 2.4A 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(\lambda_{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~.

2.3.. Ecological interpretation of matrix *A* and extension to other biological systems {#s2c}
---------------------------------------------------------------------------------------

Let us consider the variable $x = (\begin{matrix}
y \\
z \\
\end{matrix})$ which evolves according to $$\frac{d\mathbf{y}}{dt} = B\mathbf{y} - \phi\mathbf{y},$$and $$\frac{d\mathbf{z}}{dt} = C\mathbf{y} + D\mathbf{z} - \phi\mathbf{z}.$$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](#RSOS171304M2x4){ref-type="disp-formula"}) and ([2.5](#RSOS171304M2x5){ref-type="disp-formula"}) 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: --- 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 \[[@RSOS171304C50]\]).--- 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.

2.4.. The effect of a perturbation {#s2d}
----------------------------------

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 $$\overset{˙}{\textbf{x}} = A\textbf{x} - \mathit{\Phi}(\textbf{x})\textbf{x} + \varepsilon G(\textbf{x}),$$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 $\sum\limits_{i = 1}^{n}x_{i} = 1,x_{j} \geq 0,j = 1,\ldots,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**)=*A***x**−*Φ*(**x**)**x**, under the hypothesis of proposition [2.1](#RSOS171304ST1){ref-type="statement"}. We shall denote the eigenvalues of *A* by λ~*i*~, skipping the explicit dependence on *μ* for simplicity.

We recall that for *μ*≠*μ*~c~ there are two fixed points *Q*~1~, *Q*~2~ associated to the eigenvectors **v**~1~,**v**~2~.

We claim that the eigenvalues of the linearization of the vector field *F* at *Q*~1~, *DF*(*Q*~1~), are −λ~1~,λ~2~−λ~1~,λ~*j*~−λ~1~, *j*=3,...,*n*. The eigenvalues of *DF*(*Q*~2~) are −λ~2~,λ~1~−λ~2~,λ~*j*~−λ~2~ *j*=3,...,*n*. Then, if λ~1~\>λ~2~ the point *Q*~1~ is attracting and *Q*~2~ has a one-dimensional unstable manifold and a stable one of dimension *n*−1. If λ~1~\<λ~2~, they exchange their roles.

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*^−1^*AP*=*J* $$J = \begin{pmatrix}
\lambda_{1} & 0 & 0 \\
0 & \lambda_{2} & 0 \\
0 & 0 & J_{3} \\
\end{pmatrix},$$where *J*~3~ contains the remaining Jordan blocks, each one of them being an upper triangular matrix. The system for the new variable **z** becomes $$\overset{˙}{\textbf{z}} = J\textbf{z} - \hat{\mathit{\Phi}}(\textbf{z})\textbf{z},$$where $\hat{\mathit{\Phi}}(\textbf{z}) = \mathit{\Phi}(P\textbf{z}) = \textbf{u}^{T}AP\textbf{z} = \textbf{u}^{T}PJ\textbf{z}$. 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. Therefore $$\hat{\mathit{\Phi}}(\textbf{z}) = \lambda_{1}z_{1} + \lambda_{2}z_{2} + \sum\limits_{j = 3}^{n}\beta_{j}z_{j},$$for some constants *β*~*j*~, *j*=3,...,*n*. The differential of the vector field in **z** is $\mathcal{M}(\textbf{z}):=J - \hat{\mathit{\Phi}}(\textbf{z})I_{n} - \textbf{z}D\hat{\mathit{\Phi}}(\textbf{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 $$\begin{pmatrix}
{- \lambda_{1}} & \lambda_{2} & \ldots \\
0 & {\lambda_{2} - \lambda_{1}} & \ldots \\
0 & 0 & {J_{3} - \lambda_{1}{Id}_{n - 2}} \\
\end{pmatrix},$$and then the eigenvalues are −λ~1~,λ~2~−λ~1~,λ~*j*~−λ~1~, *j*=3,...,*n*.

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~=*α*, *z*~2~=1−*α*, *z*~*j*~=0, *j*=3,...,*n*, for some real *α*. Then M becomes $$\begin{pmatrix}
{- \alpha\lambda_{1}} & {- \lambda_{2}} & \ldots \\
{- (1 - \alpha)\lambda_{1}} & {- (1 - \alpha)\lambda_{1}} & \ldots \\
0 & 0 & {J_{3} - \lambda_{1}{Id}_{n - 2}} \\
\end{pmatrix}.$$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](#RSOS171304ST1){ref-type="statement"} 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 $\mathbb{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. \[[@RSOS171304C51]\] or \[[@RSOS171304C52]\]) that NHIM persist under sufficiently small perturbation. This shows that for the system *F*+*εG* there exists an invariant curve ${\hat{L}}_{\varepsilon}$, 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](#RSOS171304ST1){ref-type="statement"}*, 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 curve* ${\hat{L}}_{\varepsilon}(\mu)$*, close to* *L*(*μ*).

2.5.. Robustness of the trans-heteroclinic bifurcation under perturbations {#s2e}
--------------------------------------------------------------------------

The final question to consider is what happens to the dynamics in the new invariant attracting curve ${\hat{L}}_{\varepsilon}(\mu)$ 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](#RSOS171304ST1){ref-type="statement"}. *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*.*Then for* *ε* *sufficiently small, in* ${\hat{L}}_{\varepsilon}(\mu)$ *there exists a heteroclinic connection from* *Q*~1~(*ε*) *to* *Q*~2~(*ε*) *if* *μ*=*μ*~1~ *and, a heteroclinic connection from* *Q*~2~(*ε*) *to* *Q*~1~(*ε*) *if* *μ*=*μ*~2~.

Remark 2.7For a small fixed value of *ε* and Δ small enough, that is for *μ* close enough to *μ*~c~, the dynamics on the invariant curve ${\hat{L}}_{\varepsilon}(\mu)$ 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.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](#RSOS171304ST1){ref-type="statement"}. 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 of ${\hat{L}}_{\varepsilon}(\mu)$ 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 $\mathcal{O}(\varepsilon)$) invariant curve ${\hat{L}}_{\varepsilon}$. The dynamics on ${\hat{L}}_{\varepsilon}$ depends on the perturbation *G* and can be of any of the types which appear in one-dimensional 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](#RSOS171304M2x2){ref-type="disp-formula"}) 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 \[[@RSOS171304C44]\] we take the parameter values *r*=0.1 and *δ*~*r*~=0.05, which give $\mu_{c} = \frac{1}{3}$ and *δ*~*μ*~=0.05. For coherence with the present notation we rename the variables *x*~0~,...,*x*~7~ used in \[[@RSOS171304C44]\] 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}^{\ast}(\varepsilon),P_{3}^{\ast}(\varepsilon)$ for the fixed points. We recall that for *ε*=0 the first components of $P_{2}^{\ast}$ and $P_{3}^{\ast}$ 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 $$\mathit{\Phi} = r(x_{1} + x_{3}) + \hat{r}(x_{5} + x_{7}) + \varepsilon g(x_{1}),$$where $\hat{r} = r + \delta_{r}$. In particular this implies that $P_{2}^{\ast}$ is not changing when *μ* and *ε* are changed.

[Figure 3](#RSOS171304F3){ref-type="fig"} 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}^{\ast}$ and $P_{3}^{\ast}(\varepsilon)$ in ${\hat{L}}_{\varepsilon}$ and, so, there is only a heteroclinic connection from $P_{2}^{\ast}$ to $P_{3}^{\ast}(\varepsilon)$. 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}^{\ast}$ and $P_{3}^{\ast}(\varepsilon)$ do not exists is $\mathcal{O}(\varepsilon)$. Figure 3.Fixed points and bifurcations identified for the perturbed model. In the horizontal axis we plot $\mu - \frac{1}{3}$ and in the vertical axis the *x*~1~ coordinate of the fixed points $(x_{1}^{\ast})$. 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.

[Figure 4](#RSOS171304F4){ref-type="fig"} displays for *μ*=*μ*~c~ the fixed points on the line ${\hat{L}}_{\varepsilon}$ 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. Figure 4.Invariant line at bifurcation value for the perturbed model. The invariant curve for $\mu = \frac{1}{3} = \mu_{c}$ seen in the (*x*~1~,*x*~8~) variables. Blue (red) points correspond to unstable (stable) equilibria. The area framed with the orange lines is the biologically-meaningful phase space given by the simplex.

3.. Discussion {#s3}
==============

In a recent article \[[@RSOS171304C44]\], 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 \[[@RSOS171304C45]\], used to investigate the dynamics of competition between healthy and tumour heterogeneous cells given by different phenotypes \[[@RSOS171304C44]\]. 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 saddle--node bifurcations \[[@RSOS171304C23]\]. 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](#RSOS171304ST4){ref-type="statement"}).

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 \[[@RSOS171304C49]\]). The survival-of-the-flattest effect, described experimentally in viroids \[[@RSOS171304C53]\] and in competing computer programs \[[@RSOS171304C54]\], involves the survival of slow-reproducing replicons that are robust against deleterious mutations, and the vanishing of fast-reproducing ones that are not robust to deleterious mutations. Specifically, the model studied in \[[@RSOS171304C49]\] 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 $\tau = \mathcal{O}({(\mu_{c} - \mu)}^{- 1/2})$. This scaling law, identified in different nonlinear systems such as charge density waves and hypercycles \[[@RSOS171304C23],[@RSOS171304C25],[@RSOS171304C55]\], has been experimentally characterized in electronic circuits \[[@RSOS171304C2]\]. For transcritical bifurcations one has $\tau = \mathcal{O}(|\mu - \mu_{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 $\tau = \mathcal{O}({(\mu - \mu_{c})}^{- 1})$. Pitchfork bifurcations have been recently identified in experiments with elastic structures \[[@RSOS171304C4]\]. The trans-heteroclinic bifurcation has the same type of behaviour: $\tau = \mathcal{O}(|\mu - \mu_{c}|^{- 1})$, both for *μ*\<*μ*~c~ and for *μ*\>*μ*~c~. The derivation follows immediately from \[[@RSOS171304C44]\] (see also electronic supplementary material, §3).

Summarizing, the bifurcation described in this article, identified in a family of quasi-species equations \[[@RSOS171304C45]\], 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 \[[@RSOS171304C45]\], quasi-species theory has been widely used to investigate the dynamics of RNA viruses \[[@RSOS171304C56]--[@RSOS171304C60]\], cancer evolution \[[@RSOS171304C44],[@RSOS171304C61],[@RSOS171304C62]\] and compositional assemblies \[[@RSOS171304C63]\]. 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 ([2.1](#RSOS171304M2x1){ref-type="disp-formula"}).

Supplementary Material
======================

###### Supporting Information

We thank the members of the Complex Systems Lab for useful comments and the people from the UB-UPC dynamical systems groups for discussions.

Data accessibility {#s4}
==================

All data are including in the article and in the electronic supplementary material.

Authors\' contributions {#s5}
=======================

J.S. built the mathematical model; C.S. and R.M. carried out the analysis of the mathematical model and the related theoretical calculations; J.S. and C.S. produced and analysed the numerical data; J.S., C.S. and R.M. wrote the article. All authors gave final approval for publication.

Competing interests {#s6}
===================

We have no competing interests.

Funding {#s7}
=======

The research leading to these results has received funding from 'la Caixa' Foundation and from a MINECO grant awarded to the Barcelona Graduate School of Mathematics (BGSMath) under the 'María de Maeztu' Program (grant MDM-2014-0445). R.M. and C.S. have been supported by the Spanish grant MTM2013-41168-P. C.S. has also been partially supported by the grants 2014 SGR 1145 (Catalonia) and MTM2016-80116-P (Spain). J.S. has been partially funded by the CERCA Programme of the Generalitat de Catalunya.

[^1]: Electronic supplementary material is available online at <https://doi.org/10.6084/m9.figshare.c.3983904>.
