Dynamical effects of loss of cooperation in discrete-time hypercycles

Journal Pre-proof Dynamical effects of loss of cooperation in discrete-time hypercycles Júlia Perona, Ernest Fontich, Josep Sardanyés

PII: DOI: Reference:

S0167-2789(19)30757-2 https://doi.org/10.1016/j.physd.2020.132425 PHYSD 132425

To appear in:

Physica D

Received date : 2 December 2019 Revised date : 13 February 2020 Accepted date : 14 February 2020 Please cite this article as: J. Perona, E. Fontich and J. Sardanyés, Dynamical effects of loss of cooperation in discrete-time hypercycles, Physica D (2020), doi: https://doi.org/10.1016/j.physd.2020.132425. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier B.V.

Journal Pre-proof

1

Dynamical effects of loss of cooperation in discrete-time hypercycles

2

J´ ulia Perona,1 Ernest Fontich,1, 2, ∗ and Josep Sardany´es3, 2, † 1

3 4 5

Departament de Matem` atiques i Inform` atica, Gran Via 585, 08007 Barcelona. Universitat de Barcelona (UB), Spain 2 Barcelona Graduate School of Mathematics (BGSMath), Campus de Bellaterra, Edifici C, 08193 Bellaterra, Spain 3 Centre de Recerca Matem` atica (CRM), Campus de Bellaterra, Edifici C, 08193 Bellaterra, Spain Hypercycles’ dynamics have been widely investigated in the context of origins of life, especially using time-continuous dynamical models. Different hypercycle architectures jeopardising their stability and persistence have been discussed and investigated, namely the catalytic parasites and the short-circuits. Here we address a different scenario considering RNA-based hypercycles in which cooperation is lost and catalysis shifts to density-dependent degradation processes due to the acquisition of cleaving activity by one hypercycle species. That is, we study the dynamical changes introduced by a functional shift. To do so we use a discrete-time model that can be approached to the time continuous limit by means of a temporal discretisation parameter, labelled C. We explore dynamical changes tied to the loss of cooperation in two-, three-, and four-member hypercycles in this discrete-time setting. With cooperation, the all-species coexistence in two- and three-member hypercycles is governed by an internal stable fixed point. When one species shifts to directed degradation, a transcritical bifurcation takes place and the other hypercycle members go to extinction. The asymptotic dynamics of the four-member system is governed by an invariant curve in its cooperative regime. For this system, we have identified a simultaneous degenerate transcritical–Neimark-Sacker bifurcation as cooperation switches to directed degradation. After these bifurcations, as we found for the other systems, all the cooperative species except the one performing degradation become extinct. Finally, we also found that the observed bifurcations and asymptotic dynamical behaviours are independent of C. Our results can help in understanding the impact of changes in ecological interactions (i.e., functional shifts) in multi-species systems and to determine the nature of the transitions tied to co-extinctions and out-competition processes in both ecosystems and RNA-based systems.

lP repro of

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

Keywords: Bifurcations, Discrete-time hypercycles, Functional shifts, Origins of Life, Nonlinear dynamics, Theoretical ecology

27

Jou

rna

28

∗

†

corresponding author. E-mail: [email protected] corresponding author. E-mail: [email protected]

Journal Pre-proof 2 I.

29

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85

lP repro of

32

Hypercycles [1] are nonlinear dynamical systems formed by n polynucleotides with catalytic activity. Hypercycles have been mainly studied within the framework of prebiotic evolution and origins of life, providing a potential solution to the so-called prebiotic information crisis [1–3, 5]. The generality of hypercycle (replicator) equations has also allowed to employ this model in neural networks [6, 7], virus replication [8–10, 12], immune system [13], or ecosystem dynamics [5, 18], among others. Also, parallelisms about the error threshold and hypercycles have been discussed within the framework of the emergence of language [19]. Interestingly, hypercycles have been experimentally built using coiled-coil peptides [20], yeast cell populations [21], and cooperative engineered bacteria growing with catalytic parasites [22]. It has been argued that hypercycle species may need two minimal conditions in order to be evolutionary stable [3, 4], namely: (i) catalytic replication and (ii) capability of information storage. These two properties are found simultaneously in RNA-based replicons such as ribozymes (ribonucleic acid enzymes). RNAs with loop and stem structures, similar to those of modern tRNAs [23], are known to be stable against hydrolysis [24] also having replicability potential [25, 26]. Indeed, smaller functional RNAs have been found in viroids [27] and other RNAs [28]. Ribozymes are short RNA molecules able to catalyse specific biochemical reactions, similar to the action of protein enzymes [26, 29]. Hence, ribozymes have been considered as potential candidates forming the first autonomous, self-replicating molecular systems involved in the origins of life [3, 4, 25, 30–33]. Some of these hypothetic prebiotic RNAs were supposed to participate in ribosome-free translation of an appropriate messenger [31, 34]. Different activities have been described for natural and in vitro (e.g., peptide-bond formation [35]) evolved ribozymes. Certain introns can catalyse their own excision (self-cleavage) from single-stranded RNA (ssRNA) [26] and ligase reactions by RNA catalysts can occur even with short RNA sequences [36]. Moreover, the same RNA sequences can catalyse trans-esterification reactions for elongation of one monomer [26], ligation of two independent ssRNAs [37, 38], and cleavage of RNA into smaller sequences [26–28, 39] (see [31, 40] for reviews). Despite the functional properties of ribozymes, RNA-catalysed self-replication from RNA templates seems to be quite limited. However, recent experiments evolving catalysts at sub-zero temperatures have revealed that the combination of RNAs with cold-adaptative mutations with a previously described 50 extension operating at ambient temperatures enabled catalysing the synthesis of an RNA sequence longer than itself (adding up to 206 nucleotides) [41]. Moreover, recent experiments have shown the spontaneous formation of catalytic cycles and networks from mixtures of RNA fragments able to self-assemble into self-replicating ribozymes [42], providing evidences for selective advantage of cooperative systems composed by ribozymes. From the modelling point of view, hypercycles have been mainly investigated with continuous time approaches, for both well-mixed i.e., ordinary differential equations (ODEs) [1, 43–46, 48] and spatially-resolved [49, 51, 52] systems. ODEs typically reveal that the asymptotic coexistence for hypercycles with n = 2, 3, 4 species is typically governed by an interior stable equilibrium [1, 46, 53]. More specifically, the case n = 2 has a stable node [46], while cases n = 3 and n = 4 are governed by stable foci with fast and hardly damped oscillations [1, 45, 53], respectively. Moreover, a multitude of analyses (both numerical and analytical) have revealed that for n > 4, populations undergo selfsustained oscillations in its cooperative regime [1, 43, 44, 54]. To date, very few works have investigated hypercycles considering discrete time (e.g., using difference equations or maps [55, 56]), being mainly analysed with cellular automata models [47, 48, 50]. Specifically, the system studied by Hofbauer and others [55, 56] revealed that discretetime hypercycles with n = 2, 3 have an interior stable fixed point governing coexistence dynamics, while the case n = 4 involves oscillating coexistence governed by an invariant curve. In this article we consider the discrete hypercycle model developed in [55] to investigate the impact on the dynamics and the bifurcations when one of the species shifts from cooperative to antagonistic interactions. By the cyclic character of the system we can assume the species that shifts is the first one. To date, different architectural changes having a negative impact on hypercycles have been thoroughly investigated. These include the so-called catalytic parasites [47–49, 52] and short-circuits [50, 51], suggested to impair hypercycle’s stability thus constraining the increase of information. The case we investigate in this article is different since the cyclic structure of the hypercycle is maintained but a given replicator instructs the degradation of the next species of the system, instead of providing catalytic aid. This new system is inspired in the existence of ribozymes with trans-cleaving functions. For example, minimal trans-cleaving RNA hammerheads were generated several decades ago [57, 58]. Also, both in vitro and in vivo hammerhead ribozymes with trans-cleaving activity against viroids have been described more recently [39]. As mentioned, we are interested in the dynamics when a given species shifts from cooperative to antagonistic interactions i.e., density-dependent degradation, focusing on small hypercycles with n = 2, 3, 4 species. Although we are not modelling this functional shift explicitly by considering mutations in the catalytic motifs and their change to cleaving motifs, we investigate this shift by taking a replication constant both either positive (catalysis) or negative (cleavage). The paper is organised as follows. In Section II we introduce the studied model [55], showing its relation with the ODEs model as the discretisation time parameter C → ∞. Then, we compute the fixed points and the

rna

31

Jou

30

INTRODUCTION

Journal Pre-proof 3

87 88 89 90 91 92 93 94 95

eigenvalues for the general model in Section III A. In Section III B we analytically prove that, for any number of species, when the first species shifts to directed degradation, the asymptotic behaviour is a fixed point in the corner of the phase space, involving the out-competition of all other species providing catalytic aid. In Section III C we investigate the dynamics for the studied hypercycles with n = 2, 3, 4 species. In particular, we analytically determine the basin of attraction in the domain of the system. In Section III D we analytically obtain the rates of convergence for those cases where the ω−limit is a fixed point, showing the relevant parameters in the asymptotic expression. Numerical computations confirm the analytic findings, in particular, we illustrate the linear dependence of the number of iterations to the fixed points with the parameter C. Finally, in Section III E we provide a numerical study of the invariant curves found for the case n = 4 and ki > 0. The bifurcations tied to the functional shifts are also discussed in Sections III C and III E. Finally, Section IV is devoted to final conclusions. II.

96

lP repro of

86

MATHEMATICAL MODEL

In this section we introduce the discrete-time hypercycle model proposed by Hofbauer [55], that will be employed in this work to determine the impact of functional shifts in hypercycles. Let xi denote the concentration of the i-th species, Si , and ki the kinetic constants that quantify the strength of catalysis that the i − 1 species provides to the i-th species. For notational convenience the subindices i are modulo n, i.e., x0 = xn and also xn+1 = x1 . The system is determined by an n-dimensional function F : Rn → Rn , Fi (x) being the concentration xi in the next generation, i.e., Fi (x) represents the concentration after one unit of time. This function considers replication rate of Si to be proportional to the amount of Si−1 , according to the product xi xi−1 (catalytically-assisted replication), taking into account that the (i − 1)-th species contributes to the replication of i-th one. We write Fi (x) ∼ xi (C + ki xi−1 ),

C > 0.

Next, we P determine the proportionality factor A(x) imposing the total population to be a constant. So if n we want i=1 Fi (x) = 1. This means ! n n X X A(x)xi (C + ki xi−1 ) = A(x) C + ki xi xi−1 = 1. i=1

Pn

i=1

xi = 1

i=1

We introduce

φ(x) =

n X

ki xi xi−1

i=1

and then A(x) has to be equal to (C + φ(x))−1 . Therefore, we have the following discrete-time system:

rna

97

Fi (x) =

C + ki xi−1 xi , C + φ(x)

1 ≤ i ≤ n.

(1)

The dynamics of Map (1) spans the following (n − 1)-simplex:

Jou

n n o X S n−1 = x = (x1 , ..., xn ) ∈ Rn | xi = 1 and xi ≥ 0 for i = 1, · · · , n . i=1

To compare this map with an analogous continuous time model we rewrite the i-th component of F as follows: Fi (x) − xi =

so that

C + ki xi−1 ki xi−1 − φ(x) xi − xi = xi C + φ(x) C + φ(x)

Fi (x) − xi C = xi (ki xi−1 − φ(x)) . C −1 C + φ(x)

98 99

Interpreting now C −1 as the time interval between two generations, the Map (1) can be seen as the Euler C −1 step of the differential equation x˙ i = xi (ki xi−1 − φ(x)),

1 ≤ i ≤ n,

(2)

Journal Pre-proof 4 where we have used that lim

C→∞

101 102 103 104 105 106 107 108 109

because φ(x) is bounded. The term φ(x) is equivalent to the dilution outflow used in time-continuous models, which introduces competition between replicators also ensuring a constant population. We note that for large values of C, the discrete system introduced above will have similar properties to system Eqs. (2). As mentioned, the main goal of our article is to investigate how dynamics change considering the range e.g., k1 ≥ −1 and ki > 0 for 2 ≤ i ≤ n. Due to the cyclic structure of hypercycles, setting a negative k1 is the same that fixing any other single value of ki6=1 ≥ −1 and all others to kj6=i > 0 (i.e., in this study we will focus on the change of sign of one parameter). A positive value of k1 means that species x1 receives catalytic aid from species xn . For k1 = 0 no interaction happens between x1 and xn , and the hypercycle becomes a catalytic chain (see [4]). For k1 < 0, species xn degrades species x1 (i.e., by trans-cleaving ribozymes activity). Since we admit k1 ≥ −1, in order to have C + φ(x) > 0 in S n−1 when k1 < 0 we should take C > −k1 /4. This leads us to assume C > 1/4 in all cases. III.

110

111 112 113 114 115 116

lP repro of

100

C =1 C + φ(x)

RESULTS AND DISCUSSION

In the next Sections we will characterise the dynamics of Map (1). In Section A we will study the fixed points and their local stability. Section B discusses the behaviour of the system setting a negative k1 value. In Section C we analyse the particular cases of two- three- and four-species systems, focusing on the dynamics and the bifurcations identified in the studied hypercycles. Section D provides analytical and numerical results of the rates of convergence to the point attractors. Finally, Section E provides a study on the invariant curves and the bifurcations for the case n = 4 when k1 = 0. A.

117

Fixed poins and eigenvalues

We begin studying the fixed points ofP Map (1). In this work we assume C > 1/4, k1 ≥ −1 and ki > 0 for 2 ≤ i ≤ n. n Let ∆n−1 be the hyperplane {x ∈ Rn | j=1 xj = 1}. Note that S n−1 ⊂ ∆n−1 . To understand the bifurcation that occurs when k1 = 0 we consider F defined in ∆n−1 \ {x ∈ Rn | C + φ(x) = 0}. For k1 6= 0 we introduce the quantity M=

119

Proposition 1. We have

Pn (a) If k1 6= 0, −1/( j=2

1 kj ),

˜ n−1 = {x ∈ ∆n−1 | xi 6= 0, ∀i}. then F has a unique fixed point pn in ∆

rna

118

pn = (p1 , ..., pn )

121

122 123 124

with

pi =

1 , ki+1 M

1 ≤ i ≤ n.

(3)

The point pn ∈ S n−1 \ ∂S n−1 if and only if k1 > 0. Moreover, when k1 → 0, pn converges to (0, ..., 0, 1). If Pn ˜ n−1 . k1 = 0 or k1 = −1/( j=2 k1j ) then F has no fixed points in ∆

˜ n−1 . Then, x is a fixed point if and only if ki xi xi−1 = 0 for all i. If k1 > 0 the previous (b) Let x ∈ ∆n−1 \ ∆ conditions are also equivalent to φ(x) = 0. The points q n,m = (q1m , . . . , qnm ) such that qim = δm,i , 0 ≤ m ≤ n, are always fixed points (here δ is the Kronecker delta).

Jou

120

n X 1 . k j=1 j

Proof. (a). We assume xi 6= 0 for all i. From the condition Fi (x) = xi we get C + ki xi−1 =1 C + φ(x)

and hence ki xi−1 = φ(x) for all i. Then k2 x1 = k3 x2 = ... = kn xn−1 = k1 xn . If k1 = 0 there are no fixed points in ˜ n−1 . When k1 6= 0 we can write xi = (k1 /ki+1 )xn and determine the value of xn imposing the condition that the ∆ point is in ∆n−1 : k1 xn

n−1 X j=1

1 + xn = 1. kj+1

Journal Pre-proof 5 Directed degradation

Heterocatalysis

a

1

1

S1

0.8

S1

0.6

0.8

1 0.8 0.6

0.8

0.4

0.4

0.2

0.6

0.2

Predation, Theoretical ecology, x2

0

0.4

0.6

0 100

50

150

Predation, Theoretical ecology, x2 0.4

0.8 0.6

100

50

0.6

0.4

0.2

0.2

S2

0

0.8

0.4

lP repro of

0.2

0

1

0.2

0 0

0 0

50

0.2

100

0.4

0.6

0.8

= (x1 , ..., x

b

0

S2

1

0

0

0

50

0.2

100

0.4

0.6

0.8

= (x1 , ..., x

1

0.8

S1

S1 0.6

0.8

x3

x3

0.6

0.4

0.4

0.2

S3

S2

S3

0.2 0

0.2

S2

0

0

0.4

0.05

= (x1 , ..., x

c

S1

S1

0.2

0.15

1

S3

0.5

S3

0.1

= (x1 , ..., x

0.8

0.6

x4

0.4

x4

0.4

S2

0.3

S2

0.2

0.2 0

0.1

0.15

0.2

0

0.3

0.25

S4

= (x1 , ..., x

S4

0.05

0.1

0.15

0.2

= (x1 , ..., x

125

If k1 = −1/

Pn−1

1 j=1 kj+1 , there n−1 ˜ n−1

\∆

˜ n−1 . Otherwise k1 xn = 1/M and we get (3). is not a fixed point in ∆

. There exists l such that xl = 0 and xl+1 6= 0. If x is a fixed point we have

Jou

(b). Let x ∈ ∆

rna

FIG. 1. Schematic diagram of the studied hypercycles formed by (ribozyme) species Si (with i = 2, ..., 4) and their dynamical outcomes displayed in phase portraits considering cooperation (heterocatalysis represented with solid black arrows, setting k1 = 0.5 and ki6=1 = 1), and emergence of directed degradation (trans-cleaving activity, indicated by dashed red arrows, using k1 = −0.5 and ki6=1 = 1). (a) Two-member hypercycle: the insets display time series for x1 (black) and x2 (red) using the same initial conditions of the orbits of the phase portrait. The insets for directed degradation show time series also for x1 (black) and x2 (red) using different initial conditions that achieve the stable fixed point q 2,2 = (0, 1) (small orange dot). (b) Three-member hypercycle with a stable focus as coexistence attractor (fixed point p3 ). The three-species system with directed degradation displays a stable fixed point at q 3,3 = (0, 0, 1). (c) Four-member hypercycle with oscillatory coexistence governed by an attracting invariant curve (shown in black). Two different initial conditions are shown: one spiralling towards (green iterations) the periodic attractor and another spiralling outwards (blue iterations) displayed in a two-dimensional projection. Directed degradation for this case has a single point attractor at q 4,4 = (0, 0, 0, 1). In all panels we have set C = 10.

xl+1 =

C + kl+1 xl xl+1 C + φ(x)

and hence φ(x) = kl+1 xl = 0. Therefore

xi =

126 127

C + ki xi−1 1 xi = xi + ki xi xi−1 , C C

for all i and thus ki xi xi−1 = 0. Conversely, ki xi xi−1 = 0 for all i implies φ(x) = 0 and one immediately gets that x is a fixed point. To obtain the eigenvalues of DF (pn ) and study the stability of the fixed point pn obtained in (a) of Proposition 1 it is convenient to use the following baricentric variables. ki+1 xi yi = Pn , j=1 kj+1 xj

1 ≤ i ≤ n.

Journal Pre-proof 6 When k1 > 0, this change of variables sendsPS n−1 to S n−1 bijectively and, more generally, for k1 ∈ R, the points Psends n n of S n−1 except the ones on the hyperplane j=1 kj+1 xj = 0 to S n−1 . Whenever defined, i.e., when j=1 kj+1 xj 6= 0, the differential of the change has rank n − 1 and, actually it is a (local) diffeomorphism from S n−1 to S n−1 . This means that we can compute the eigenvalues of DF at pn in baricentric coordinates. In such coordinates F reads C+

j=1 yj−1 yj

Ψ(y)

n

yi ,

where

Ψ(y) =

n X j=1

1 yj , kj+1

lP repro of

Fi (y) =

yi−1 C + Ψ(y) Pn 1

and the fixed point p located at (1/n, ..., 1/n). It is not difficult to compute the partial derivatives and obtain ∂Fi n 2 (p ) = 1 − , ∂yi n(CM + 1) 128

∂Fi n 1 2 (p ) = − , ∂yi−1 CM + 1 n(CM + 1)

Then, the differential DF (pn ) is a circulant matrix  c0 c1 cn−1 c0  cn−2 cn−1  c1

with

c0 = 1 −

c2

∂Fi n −2 (p ) = , ∂yl n(CM + 1)

l 6= i, i − 1.

 . . . cn−1 . . . cn−2   . . . cn−3   ... . . . c0

(4)

2 1 2 −2 , cn−1 = − , and ci = for 1 ≤ i ≤ n − 2. n(CM + 1) CM + 1 n(CM + 1) n(CM + 1)

It is known [60] that the eigenvalues of a circulant matrix as (4) are λm =

n−1 X

cj e2πijm/n ,

j=0

where i denotes the imaginary unit

√

0 ≤ m ≤ n − 1,

−1, with corresponding eigenvectors (1, e−2πim/n , ..., e−2πi(n−1)m/n ).

In our case

λm = 1 +

130 131 132 133

0 ≤ m ≤ n − 1.

The eigenvalue λ0 corresponds to the eigenvector (1, 1, ..., 1) which is transversal to S n−1 . The other eigenvalues Pn−1 correspond to eigenvectors tangent to S n−1 . Indeed, when m 6= 0, l=0 e−2πilm/n = 0. To compute the eigenvalues of DF (q n,n ) we first look for the linearisation of F (in the original coordinates) at q n,n = (0, 0, ..., 1). To do so we translate it to the origin by means of the change of coordinates xn = ξn + 1, xi = ξi , 1 ≤ i ≤ n − 1. In these variables the map is expressed as:

rna

129

1 e2πim/n , CM + 1

Jou

C + k1 (ξn + 1) F˜1 (ξ) = ξ1 , ˜ C + φ(ξ) C + ki ξi−1 F˜i (ξ) = ξi , 2 ≤ i ≤ n − 1, ˜ C + φ(ξ)

C + kn ξn−1 (ξn + 1) − 1, F˜n (ξ) = ˜ C + φ(ξ)

134 135 136 137

Pn ˜ where φ(ξ) = j=1 kj ξj ξj−1 + k1 ξ1 + kn ξn−1 . From these expressions we readily obtain   1 + kC1 0 . . . 0  0 1 . . . 0 . DF (q n,n ) =    ... −k1 0 . . . 1 C

The eigenvalues are 1+k1 /C and 1. The eigenvalue 1+k1 /C corresponds to the eigenvector (1, 0, ..., 0, −1). The eigenvalue 1 corresponds to the (linearly independent) eigenvectors (0, 1, −1, 0, ..., 0), (0, 1, 0, −1, ..., 0), ..., (0, 1, 0, ..., 0, −1) and (0, ..., 0, 1). All these vectors are tangent to S n−1 except the last one. Proceeding in an analogous way we can check that the eigenvalues of DF (q n,i ) are 1 + ki+1 /C and 1.

Journal Pre-proof 7 B.

138

140 141 142 143

144 145

In this section we will prove that for k1 ≤ 0 the dynamics achieves the fixed point q n,n . This involves that the species that performs directed degradation will outcompete all of the others. Let us go back to Map (1). As mentioned, by the cyclic structure of the map we only deal with the case k1 ≤ 0. But, by the symmetry, in the same way we have that if kj ≤ 0 and ki6=j > 0 the dynamics achieves q j−1,j−1 . We now assume that −1 ≤ k1 ≤ 0, and ki > 0 for 2 ≤ i ≤ n. These conditions ensure that for x ∈ S n−1 both C + ki xi−1 and C + φ(x) are positive.

lP repro of

139

When k1 ≤ 0 the basin of attraction of q n,n contains S n−1 \ ∂S n−1

Proposition 2. Assume C > 1/4, −1 ≤ k1 ≤ 0, and ki > 0 with 2 ≤ i ≤ n. If x0 ∈ S n−1 \ ∂S n−1 then {F m (x0 )} converges to q n,n = (0, 0, ..., 1). m m 0 0 Proof. We write xm = (xm / ∂S n−1 , 0 < x0i < 1 for all i. Moreover, by the form of F , 1 , ..., xn ) = F (x ). Since x ∈ m m 0 < xi < 1 for all m and i. First, we check that {x1 } is strictly decreasing and converges to 0. Indeed, since k1 ≤ 0 m m m m m and xm 1 < 1, k1 xn x1 ≥ k1 xn and since ki > 0 for 2 ≤ i ≤ n, φ(x ) > k1 xn . Then

0<

C + k1 xm n <1 C + φ(xm )

xm+1 = 1

and

C + k1 xm n xm < xm 1 , C + φ(xm ) 1

m ≥ 1.

By compactness of S n−1 there is a subsequence {xmk } of {xm } which converges to some x ˜ = (˜ x1 , ..., x ˜n ) ∈ S n−1 . m Note that, by monotonicity, {x1 } converges to x ˜1 . We assume that x ˜1 > 0 to get a contradiction. Taking limit in k +1 xm = 1

we get

k C + k1 xm n xmk C + φ(xmk ) 1

C + k1 x ˜n =1 C + φ(˜ x)

146 147 148

which implies k1 x ˜n = φ(˜ x), or equivalently k1 x ˜n (1 − x ˜1 ) = k2 x ˜2 x ˜ 1 + k3 x ˜3 x ˜ 2 + · · · + kn x ˜n x ˜n−1 . The left hand side is less or equal than zero while the right hand one is bigger or equal than zero. Therefore ki x ˜i x ˜i−1 = 0, 2 ≤ i ≤ n, and, in particular, x ˜2 x ˜1 = 0 which gives x ˜2 = 0. From

and

rna

m C + k1 xm xm+1 n x1 1 m+1 = C + k xm xm x2 2 1 2

0<

we have that

n

xm 1 xm 2

o

C + k1 xm n <1 C + k2 xm 1

is strictly decreasing, in particular is bounded from above. Then

149

Jou

k xm = 1

k xm mk 1 m k x2 x2

converges to 0 which provides the desired n contradiction. o xm Now, we claim that, for 1 ≤ i ≤ n − 1, xmi is strictly monotone for m ≥ M for some M (depending on i) and i+1 xi → 0. Indeed, by the previous arguments the statement is true for i = 1. We assume it is true for 1 ≤ i ≤ n − 2. Let ki+1 xm i , m→∞ ki+2 xm i+1

γi = lim

Note that 0 ≤ γi ≤ ∞. If γi > 1, or γi = 1 and

n

xm i xm i+1

o

1 ≤ i ≤ n − 2.

is decreasing, then

ki+1 xm i > 1, ki+2 xm i+1

Journal Pre-proof 8 150

n m o i+1 ˜. ˜ for some M ˜ and then xm is strictly increasing for m ≥ M for m ≥ M n m o xi+2 x If γi < 1, or γi = 1 and xmi is increasing, then i+1

152

ˆ for some M ˆ and then for m ≥ M

n

xm i+1 xm i+2

o

ˆ. is strictly decreasing for m ≥ M

lP repro of

151

ki+1 xm i < 1, ki+2 xm i+1

Now, to prove that {xm i+1 } converges to zero we distinguish two cases: γi > 0 and γi = 0. When γi > 0 there exists m0i such that ki+1 xm γi i > , m ki+2 xi+1 2

and then, from

xm i+1 <

153

we get xm i+1 → 0. When γi = 0, there exists m ˜ 0i such that

1 ki+1 xm i < ki+2 xm 2 i+1 Obviously, for m ≥ m ˜ 0i ,

m ≥ m0i ,

ki+1 2 m x , ki+2 γi i

for

m≥m ˜ 0i .

C + (1/2)ki+2 xm C + ki+1 xm i+1 i < < 1. m C + ki+2 xi+1 C + ki+2 xm i+1

If we assume that {xm i+1 } does not converge to 0, then there exists ε > 0 and infinitely many indices m such that ki+2 xm > ε and therefore infinitely many factors i+1 C + ki+1 xm C + (1/2)ε i . < C + ki+2 xm C +ε i+1

This means that, given q,

156

157

rna

155

with qm

Finally, since xm ∈ S n−1 , xm n → 1. C.

i+2

i+2

Case studies: Hypercycles with n=2, n=3, and n=4 members 1.

Jou

154

 q q xm C + (1/2)ε m xi+1 i+1 < , m > q, xm C +ε xqi+2 i+2 n m o xi+1 xm i+1 m → ∞ as m → ∞. Clearly, xm → 0 and xm i+1 < xm gives that {xi+1 } converges to 0.

Case n=2

In this case the model is essentially one dimensional. When ki > 0 it has a unique inner fixed point   k1 k2 p2 = , , k1 + k2 k1 + k2

and the fixed points q 2,1 = (1, 0) and q 2,2 = (0, 1). The eigenvalue at p is 1+

158 159 160 161 162

1 CM C(k1 + k2 ) e2πi/2 = = < 1. CM + 1 CM + 1 C(k1 + k2 ) + k1 k2

The eigenvalues at q 2,1 and q 2,2 are 1 + k2 /C and 1 + k1 /C, respectively. Actually, p2 attracts all points of S 1 \ ∂S 1 . When k1 → 0 with k2 fixed, the fixed point p2 tends to q 2,2 and they undergo a transcritical bifurcation. When k1 ≤ 0 all points of S 1 \ ∂S 1 tend to q 2,2 . The bifurcation diagram obtained by iteration of Map (1) and tuning −1 ≤ k1 ≤ 1 is displayed in Fig. 2(a). Here, for 0 < k1 ≤ 1 the coexistence equilibrium is given by the fixed point p2 . At k1 = 0 the points p2 and q 2,2 collide in a transcritical bifurcation. Then, for negative values of k1 the point q 2,2 is stable.

Journal Pre-proof 9 Functional shift

a

directed degradation

b

cooperation

6

6×10

S1

0.5

S2

S2 0 -1

0

-0.5

Iterations to attractor

S1

k1

6

4×10

6

2×10

0

lP repro of

Stable equilibrium

1

6

4×10

6

3×10

6

2×10

0 0

1

0.5

4×10

4

4

8×10

C

FIG. 2. (a) Bifurcation diagram obtained by iteration of Map (1) when n = 2 using k1 as a control parameter with k2 = 1 and C = 10. Black and red lines denote the equilibrium population of species x1 and x2 respectively. For 0 < k1 ≤ 1 the dynamics is attracted by the fixed point p2 , while for −1 ≤ k1 < 0 the stable fixed point is q 2,2 , involving the persistence of the second replicator and the extinction of S1 . At k1 = 0 the fixed points p2 and q 2,2 collide in a transcritical bifurcation. (b) Linear dependence of parameter C on the number of iterations needed to achieve the attractors fixing k2 = 1 and: (upper panel, for attractor p2 ) k1 = 1 (black), k1 = 0.75 (red), k1 = 0.5 (blue), k1 = 0.25 (green); (lower panel, for attractor q 2,2 ) we have used the same values of k1 than in the upper panel but with negative sign. We consider δ = 10−6 . In all panels we have used x1 (0) = 0.75, x2 (0) = 0.25 as initial conditions. 2.

163

When ki > 0 the inner fixed point is given by 3

p =

and the corresponding eigenvalues are

λ1,2 = 1 +

|λ1,2 |2 = 1 + 166 167 168 169 170 171 172 173 174 175

176

1 eiθ1,2 , CM + 1

2 cos θ1,2 + CM + 1



1 CM + 1

θ1 =

2



,

2π , 3

=1−

θ2 =

1 CM + 1

4π . 3



1−

1 CM + 1



< 1.

The other fixed points, according to Proposition 1, satisfy φ(x) = 0. The only possibilities are q 3,1 = (1, 0, 0), q 3,2 = (0, 1, 0) and q 3,3 = (0, 0, 1). They have an eigenvalue of modulus greater than 1. The point p3 is an attractor. In Ref. [55] it is proved, by using a strict Lyapunov function, that S 2 \ ∂S 2 is the basin of attraction of p3 . When k1 → 0 with k2 and k3 fixed, p3 tends to q 3,3 and they undergo a (degenerate) transcritical bifurcation. At the bifurcation, the two eigenvalues are 1. A special feature is that at the bifurcation there is a segment of fixed points {x2 = 0, x1 + x3 = 1} with q 3,3 in an extreme of it. After the bifurcation, i.e., when k1 < 0, p3 is outside S 2 , it is unstable. Moreover, when k1 ≤ 0, q 3,3 attracts all points of S 2 \ ∂S 2 . The dynamics for n = 3 is displayed in Fig. 3(a) by means of a bifurcation diagram built iterating Map (1). Here, similarly to the case n = 2, the hypercycle persists for 0 < k1 ≤ 1 because the point p3 is stable. At k1 = 0, there is a degenerate transcritical bifurcation between the points p3 and q 3,3 , and for negative values of k1 the third member outcompetes all other species i.e., the fixed point q 3,3 attracts all points of S 2 \ ∂S 2 .

Jou

165 164

1 1 1 , , k2 M k3 M k1 M

rna

We have



Case n=3

3.

Case n=4

For ki > 0 the dynamics is governed by an invariant curve [55, 56] that allows the coexistence of all of the species by means of an oscillatory regime (see Figs. 1, 4(a), and 6). When ki > 0, the inner fixed point is given by   1 1 1 1 4 p = , , , , k2 M k3 M k4 M k1 M

Journal Pre-proof 10

a

b S1

Iterations to attractor

6

0.8 S3

S2

0.6

S1

0.4 0.2

S3

S2

0 -1

0

-0.5

5

6×10

0 9 1×10

8

8×10

8

4×10

0 0

1

0.5

k1

1×10

lP repro of

Stable equilibrium

1

3

2×10

3

4×10

C

FIG. 3. (a) Bifurcation diagram obtained by iteration of Map (1) when n = 3 using k1 as a control parameter with k2 = 1, k3 = 0.5, and C = 10. Here we show equilibria for variables x1 (black), x2 (red), and x3 (green). For positive k1 the dynamics achieve the fixed point p3 . At k1 = 0 there is a degenerate transcritical bifurcation between the fixed points p3 and q 3,3 . For negative k1 the fixed point q 3,3 is an attractor. (b) Linear dependence of parameter C on the number of iterations needed to achieve the attractor q 3,3 fixing k2 = 1 and: (upper panel, for attractor p3 ) k1 = 1 (black), k1 = 0.75 (red), k1 = 0.5 (blue), k1 = 0.25 (green); (lower panel, for attractor q 3,3 ) here we have used k1 = −1. Due to the extremely long transients obtained for k1 < 0 we here consider δ = 10−5 and a shorter range for C. Here the four values of k1 < 0 give place to very similar transient times, which are displayed overlapped and also have a linear dependence on C. In all panels we have used x1 (0) = 0.5, x2 (0) = 0.35, and x3 (0) = 0.15 as initial conditions.

and its eigenvalues are

1 eiθj CM + 1

λj = 1 + We have 2

2

|λ1 | = |λ3 | = 1 +

180 181 182 183 184

185

186 187 188 189 190 191 192 193 194 195 196

2

>1

θj = ei2πj/4 ,

and

2

|λ2 | =

1 ≤ j ≤ 3.



1 1− CM + 1

2

< 1.

Moreover, on ∂S 3 we have the fixed points q 4,i , with qi4,i = δij (δ being the Kronecker delta) and the segments of fixed points {(α, 0, 1 − α, 0)| α ∈ [0, 1]}, {(0, α, 0, 1 − α)| α ∈ [0, 1]}. When k1 = 0 we also have the segment of fixed points {(α, 0, 0, 1 − α)| α ∈ [0, 1]}. When k1 → 0, p tends to q 4,4 and at the bifurcation value k1 = 0 all eigenvalues are equal to 1. At the bifurcation and after it, i.e. when k1 ≤ 0, q 4,4 attracts all points of S 3 \ ∂S 3 . Figure 4 displays how local maxima and minima obtained from time series for the dynamics on the invariant curve change at decreasing k1 from 1 to 0. Notice that the invariant curve shrinks (see also Fig. 6(a)), finally collapsing at k1 = 0 (the stability of the invariant curve as well as the bifurcations occurring at crossing k1 = 0 are discussed in Section III E below).

rna

178 179

1 CM + 1

Jou

177



with

D.

Rates of convergence to the point attractors

In this section we study the rates of convergence of the point attractors of the system. For that, given an initial condition x0 ∈ S n−1 \ ∂S n−1 , we compute the number of iterations m to arrive to a ball of radius δ centered at the attractor. We have several cases depending on m and on whether the attractor is the inner point or a vertex. Also, the computations are different if the fixed point is hyperbolic or not. For n = 2 and n = 3, if k1 > 0, the attractors are p2 and p3 , respectively, which are hyperbolic. If k1 ≤ 0, the attractors are q 2,2 and q 3,3 . The point q 2,2 is hyperbolic if k1 < 0 while both q 2,2 and q 3,3 have eigenvalues equal to one in the other cases. When n = 4, if k1 > 0, there is not an attracting fixed point. If k1 ≤ 0, the attractor is q 4,4 , which has eigenvalues equal to 1. Here attractor is understood as a fixed point, which attracts all points of the interior of the simplex. Notice that, in some cases, they have eigenvalues equal to 1. Together with the analytical derivations developed along this section, we also provide numerical results computing the number of iterations to achieve the attractors, showing their linear dependence with the discretisation parameter C (and with replication constants, see below). Specifically, Fig. 2(b) displays this linear relation between

Journal Pre-proof

S2

S1

S3

S4

0.2

S2

Local minima x_1(n) on the invariant S4 curve

P* = (0, 0, 0, 1)

0 x_3(n)

0.6 0.4 0.2 0 -1

P* = (0, 0, 0, 1)

-0.5

0

≥ k1

0.5

Equilibrium coordinate, x4

0.4

k2 x ˜2 Equilibrium coordinate,

S3

0.6

x_2(n)

0.4

lP repro of

k3 x ˜3 x Equilibrium coordinate,

Local maximaSx_1(n) on the invariant curve 1

0.6

1

Local maxima and minima Local maxima and minima

Local maxima and minima Local maxima and minima

Equilibrium coordinate, x ˜2 x ˜1

11

0.2

P* = (0, 0, 0, 1)

0

1

P* = (0, 0, 0, 1)

0.8

0.6 0.4

x_4(n)

0.2

1

0 -1

-0.5

0

≥ k1

0.5

1

FIG. 4. Bifurcation diagram for the four-member hypercycle obtained by iteration of Map (1) using −1 ≤ k1 ≤ 1 as control parameter, setting k2,3,4 = 1 and using the initial condition x0 (0) = x1 (0) = x2 (0) = 0.025 and x4 (0) = 0.925. The black and red dots display, respectively, the local maxima and minima of each variable obtained from time series once the dynamics has settled on the invariant curve for k1 > 0 (right y-axis). For 1 ≤ k1 < 0, the equilibrium of each coordinate is also displayed (left y-axis). Here the only species that persists is x4 . In all panels we set C = 10.

S1

S3

rna

Iterations to attractor

9

2×10

9

1×10

S2

S4

8

5×10

Jou

0 0

200

400

C

600

800

1000

FIG. 5. Number of iterations to reach the attractor q 4,4 for n = 4 and their relation with C, setting 1 ≤ C ≤ 103 and k2 = 1, k3 = 0.5, k4 = 1, with; k1 = −1 (black line), k1 = −0.8 (red) dots, and k1 = −0.6 (black squares). Here we have used δ = 10−6 and the initial conditions x0 (0) = x1 (0) = x2 (0) = 0.025 and x4 (0) = 0.925. See Section III D for details.

197 198 199 200 201 203 202 204

C and the iterations to the coexistence attractor p2 (upper panel) and to the out-competition attractor q 2,2 (lower panel). Also, we have found the linear relation between C and the iterations to the coexistence attractor p3 (upper panel in Fig. 3(b)) and the out-competition one q 3,3 (lower panel in Fig. 3(b)). Finally, Fig. 5 also displays the linear relation between constant C and the out-competition attractor q 4,4 . Specifically, we have obtained that for the Pn 1 2 3 points p , p the times are proportional to C( i=1 ki ), n = 2, 3, and for the points q n,n , n = 2, 3, 4, the times are proportional to C/kn . Next, we describe in detail the computation of the number of iterations in the more involved cases i.e., for p3 when

Journal Pre-proof 12 205 206

n = 3 and q 4,4 when n = 4. The other cases are studied using the same ideas in a much simpler way. For the latter we will just make some comments on the variations on the arguments and give the results. Convergence to p3

1.

207

lP repro of

When n = 3 the system is essentially two dimensional. We use the variables x1 , x2 to describe S 2 . The eigenvalues at p3 have already been computed and are √ 3 1 ¯1, +i , λ2 = λ λ1 = 1 − 2(CM + 1) 2(CM + 1) with

|λ1 | = |λ2 | = 1 − 208 209 210

1 1 < 1. + CM + 1 (CM + 1)2

Since p3 is a hyperbolic attractor (without resonances) we can apply Poincar´e’s theorem [61] and get that the system is locally conjugated to its linear part L by an analytical conjugation h defined in a neighbourhood of 0 sending 0 to p3 and satisfying Dh(0) = Id . Specifically, we have F ◦h=h◦L

in a neighbourhood of 0.

From basic algebra we know that there exists a non-singular matrix B such that B   cos ϕ1 − sin ϕ1 ˜ L = |λ1 | , ϕ1 = arg λ1 . sin ϕ1 cos ϕ1 Clearly,

212 213 214 215 216

˜ = h ◦ B and we have We take h

  m m cos(mϕ1 ) − sin(mϕ1 ) ˜ L = |λ1 | . sin(mϕ1 ) cos(mϕ1 )

˜ =F ◦h◦B =h◦L◦B =h◦B◦L ˜ ◦ L. ˜=h ˜ F ◦h

(6)

˜ is defined in a ball of radius r, Br (0). Using (6) we can extend the domain of h ˜ to a neighbourhood We assume that h ˜ ˜ defined on U of 0 such that h(U) is contained in the image by F of its domain of invertibility. Indeed, we start with h −1 ˜ ˜ ˜ ˜ Br (0) and we inductively use h = F ◦ h ◦ L to extend, at step j, the domain of h from B|λ1 |−j+1 r (0) to B|λ1 |−j r (0). This can be done while F −1 exists. Then, eventually we have to stop at some step j0 . Notice that if the parameter ˜ can be extended to R2 . C is big enough, F is globally invertible in the simplex and in such case the domain of h
˜ B|λ|−j0 r (0) . Let x0 ∈ S 2 \ ∂S 2 . Since p3 is a global attractor (Theorem 3 of [55]) there exists We denote U = h m0 ≥ 1 such that F m0 (x0 ) ∈ U. We can write

rna

211

(5)

˜ with LB = L

−1

˜ 0 )), F m (x0 ) = F m−m0 (F m0 (x0 )) = F m−m0 (h(y

0

for some y ∈ B|λ|−j0 r (0). Then

m ≥ m0

˜ 0 )) − h(0)k ˜ ˜◦L ˜ ˜ m−m0 (y 0 ) − h(0)k. kF m (x0 ) − p3 k = kF m−m0 (h(y = kh

Jou

˜ ˜ m−m0 y 0 is very close to 0 and we have Dh(0) = B, Since we look for m such that L

˜◦L ˜ ˜ m−m0 y 0 k . kh ˜ m−m0 (y 0 ) − h(0)k ˜ m−m0 (y 0 )k. kB −1 k−1 kL . kBk kL

˜ m−m0 y 0 k = δ is equivalent to Moreover, kL

If C is big,

log |λ1 | =

217 218 219

m=

log δ − log ky 0 k + m0 . log |λ1 |

−1 1 3 1 1 + − +O 2 CM 2 (CM ) 3 (CM )3



1 (CM )4



,

˜ −1 (F m0 (x0 ))k−1 ) + m0 . and then m is of the order CM (log δ −1 − log kh −1 m ˜ (F 0 (x0 )) should be interpreted as a constant depending on the initial Here, and in the following cases, log h condition.

Journal Pre-proof 13

221 222

Convergence to p2 (when k1 > 0)

2.

220

In this case p2 is a hyperbolic fixed point (k1 > 0) and the corresponding eigenvalue is λ = (CM )/(1 + CM ). Using the same strategy as before, we obtain

3.

223

lP repro of

˜ −1 (F m0 (x0 )) log δ − log h + m0 = log λ   i h 1 ˜ −1 (F m0 (x0 )))−1 1 + O + m0 . = CM log δ −1 − log(h CM

m≈

Convergence to q 2,2 (when k1 < 0)

The eigenvalue corresponding to q 2,2 is 1 + k1 /C < 1. Similarly as before we now have   i C h ˜ −1 (F m0 (x0 )))−1 1 + O k1 log δ −1 − log(h m≈ + m0 . (−k1 ) C 4.

224

225 226

Convergence to q 4,4 (when k1 ≤ 0)

The point q 4,4 is not hyperbolic and this fact forces to introduce several technicalities. We start with a lemma which provides control on the convergence of some sequences. Lemma 3. Let γ > 0, m0 ≥ 0 and {z m } be a sequence of positive numbers. If z m+1 ≥

If z m+1 ≤

227

zm 1+γz m

zm ≥

z m0 , 1 + (m − m0 )γz m0

m ≥ m0 .

zm ≤

z m0 , 1 + (m − m0 )γz m0

m ≥ m0 .

for m ≥ m0 then

zm 1+γz m

for m ≥ m0 then

The same statement is true with strict inequalities with the conclusions for m > m0 .

rna

Proof. Let {ξ m } be the auxiliary sequence defined by ξ m0 = z m0 and ξ m+1 =

ξm , 1 + γξ m

m ≥ m0 .

m0

Jou

ξ m We easily check by induction that ξ m = 1+(m−m ≥ ξ m for all m ≥ m0 . Indeed, when m . We claim that z 0 )γξ 0 t m = m0 this is obviously true. Assuming it is true for m ≥ m0 , and using that ϕ(t) = 1+γt is strictly increasing in (0, ∞) we have

z m+1 ≥

228

229 230 231 232 233 234 235 236 237

zm ξm ≥ = ξ m+1 . 1 + γz m 1 + γξ m

Then the result is obtained. The second part follows in the same way. m Let (x01 , x02 , x03 , x04 ) ∈ S 3 \ ∂S 3 . We already know from the proof of Proposition 2 that the sequences {xm 1 /x2 }, m m m m m m m }, {x } converge to 0 and {x /x } and {x /x } are strictly monotone from some index on, that {x }, {x {xm 4 } 3 4 1 2 3 2 3 converges to 1. In the next claims we will use a small constant ε > 0 and an integer m0 sufficiently big. They will be the ones needed for certain conditions on sequences to be met, and may be different at different places. We will require a finite m m m (small) number of such conditions. Given ε ∈ (0, 1) there exists m0 such that xm 1 < ε, x2 < ε, x3 < ε and x4 > 1 − ε m m m for m ≥ m0 . Since {xm /x } is strictly decreasing, x /x < β for some β > 0. 1 1 1 2 m≥0 1 2 C+k3 xm m m m m 2 Moreover, since xm < 1 and hence k3 xm 3 /x4 → 0, {x3 /x4 } is strictly decreasing for m ≥ m0 , then C+k4 xm 2 < k4 x3 3 for m ≥ m0 .

Journal Pre-proof 14 238

m Claim 4. {xm 1 /x3 }m≥0 converges to 0.

Proof. First we consider the case k1 = 0. We have that C xm xm+1 1 1 m+1 = C + k xm xm x3 3 2 3 xm

243

Assume ε is small enough so that (k3 /(β1 C))ε < 1. Using that log

240 241

244

245

log

C C+(k3 /β1 )xj1

lP repro of

242

m 1 and hence { xm }m≥0 is strictly decreasing. We know that xm 1 ≤ β1 x 2 . 3 m m m m To get a contradiction we assume that limm→∞ x1 /x3 = β2 > 0. Then xm 1 > β2 x3 for m ≥ 0 and φ(x ) = m+1 C 1 m m m m m m m m m k2 x2 x1 + k3 x3 x2 + k4 x4 x3 < (k2 + k3 /β2 + k4 /β2 )x1 . Then x1 = C+φ(xm ) x1 > 1+γ1 xm x1 , where γ1 = 1 0 0 −1 (k2 + k3 /β2 + k4 /β2 )/C for all m ≥ 0. By Lemma 3, xm > x /(1 + mγ x ) . Then 1 1 1 1 ! ! m−1 m−1 0 X Y x x01 xm C C 1 1 = = exp . (7) log m 0 j j x3 x3 C + k3 x2 x03 j=0 C + k3 x2 j=0

239

1 1+t

x01

249 250 251 252 253

1

xm

Since

C+k1 (1−ε) C

< 1, we also have

xm 1 xm 3

→ 0.

m Claim 5. {xm 2 /x3 }m≥0 converges to 0.

2

1+A 1 < , 1+B 1 + B − 2A

Jou

C + k2 xj1 C+

257

258 259

(8)

On the other hand, using that if A > 0 and −1 + 2A < B < 1 + 2A

we have

256

0 (k3 /k4 )xm 2 . 0 1 + (m − m0 )γ3 xm 2

rna

3

xm 3 ≥

255

m ≥ m0 .

m m m Proof. We assume that lim xm 2 /x3 = β3 > 0. By the condition k3 x2 < k4 x3 for m ≥ m0 we have β3 ≤ k4 /k3 . m m m m m Then x2 > (β3 −ε)x3 for m > m0 . Moreover, since x1 /x3 < ε for m ≥ m0 we also have that xm 1 < (ε/(β3 −ε))x2 . C+k2 xm m+1 1 m m m m m m m m m m 1 Then φ(x ) ≤ k2 x2 x1 + k3 x3 x2 + k4 x4 x3 < (k2 + k3 + k4 /(β3 − ε))x2 . Then x2 = C+φ(xm ) x2 > 1+γ3 xm x2 , 2 where γ3 = (k2 + k3 + k4 /(β3 − ε))/C for all m ≥ m0 . m x2 0 m m By Lemma 3, xm 2 ≥ 1+(m−m )γ xm0 for m ≥ m0 . Moreover, using again that k3 x2 < k4 x3 , 0

254

<

to −∞ when m → ∞ and hence x1m → 0 which is a contradiction. 3 When k1 < 0, we use that given ε > 0 there exists m0 such that if m ≥ m0 then xm 4 > 1 − ε. Then

247

248

C C+k3 xj2

< −(log 2)(k3 /(β1 C))xj1 < −(log 2)(k3 /(β1 C)) 1+jγ1 x0 for m ≥ 0 and therefore the sum in (7) diverges

xm+1 C + k1 (1 − ε) xm 1 1 , m+1 ≤ C xm x3 3 246

< −(log 2)t for t ∈ (0, 1) we have log

k3 xj2

≤

C + εk2 xj3

C + k3 (β3 −

ε)xj3

≤

1 1 + γ3 xj3

(9)

,

j ≥ m0 ,

with γ3 = (k3 (β3 − ε) − 2εk2 )/C and ε so small that γ3 > 0. Then ! ! m−1 m−1 Y C + k2 xj xj0 X xm C + k3 xj2 xj20 2 1 2 = = exp log . j xm xj30 C + k3 xj2 xj30 3 j=j0 C + k3 x2 j=j0 Assume j0 is big enough so that γ3 xj30 < 1. Using that log log 1+γ1 xj 3 3 xm 2 hence xm 3

<

−(log 2)γ3 xj3 .

1 1+t

< −(log 2)t for t ∈ (0, 1) we have log

(10) C+k2 xj1 C+k2 xj2

≤

Taking into account (8) we get that the sum in (10) diverges to −∞ when m → ∞ and

→ 0 which is a contradiction.

Journal Pre-proof 15 To estimate the distance from F m (x0 ) to q 4,4 we write m m m 2 m 2 m 2 m 2 m m m 2 k(xm 1 , x2 , x3 , x4 ) − (0, 0, 0, 1)k = (x1 ) + (x2 ) + (x3 ) + (x1 + x2 + x3 ) " !2 !2 # m xm xm xm xm xm 1 2 1 x2 1 2 2 = 2(xm + 3 ) 1+ m + m + m m + m x3 x3 x3 x3 x3 xm 3

so that the asymptotic behaviour depends on how {xm 3 } tends to 0.

261

Claim 6. Given ε > 0 there exists m0 ≥ 1 such that

lP repro of

260

0 0 xm xm m 3 3 0 ≤ x3 ≤ 0 , m 1 + (m − m0 )((k4 + ε ν2 )/C)x3 1 + (m − m0 )((k4 + ε ν1 )/C)xm 3 262

m ≥ m0 ,

(11)

(12)

where ν1 = k1 − 2k3 − k4 and ν2 = k3 + εk2 .

m m m m m m Proof. By the previous claims we have that xm 1 < εx3 and x2 < εx3 for m ≥ m0 . Also x1 , x2 , x3 < ε for m ≥ m0 . First we establish the bounds m φ(xm ) ≥ εk1 xm 3 + k4 (1 − ε)x3 ,

and

m m φ(xm ) ≤ ε2 k2 xm 3 + ε k3 x3 + k4 x3 ,

Then, using (9),

m ≥ m0 ,

m ≥ m0 .

C + εk3 xm 1 C + k3 xm 3 2 ≤ ≤ C + φ(xm ) C + (k4 + ε(k1 − k4 ))xm 1 + ((k + εν1 )/C)xm 4 3 3 which gives

xm+1 ≤ 3

1 xm 3 , 1 + ((k4 + εν1 )/C)xm 3

and by Lemma 3, we obtain the right hand side inequality of the claim. On the other hand

263

rna

C 1 C + k3 xm 2 ≥ = , C + φ(xm ) C + (k4 + ε(k3 + εk2 ))xm 1 + ((k + εν2 )/C)xm 4 3 3 and, using Lemma 3 again, we obtain the other inequality.

0 n−1 With the information on the rate of convergence of {xm \ ∂S n−1 , the 3 } we can now estimate, given x ∈ S m 0 4,4 number of iterations m for F (x ) to arrive to a distance δ from q . The condition for m is obtained putting xm 3 =δ in (12). From this we get

264 265 267 266

268

Jou

C C 0 0 (1/δ − 1/xm (1/δ − 1/xm 3 ) + m0 < m < 3 ) + m0 . k4 + εν2 k4 + εν1

That is, apart from a transitory, the number of iterations for xm 3 to get δ is essentially proportional to C/k4 , and by (11), the number√ of iterations for xm to arrive to a neighbourhood of q 4,4 of radius δ is given by the previous formula changing δ by δ/ 2. 5.

Convergence to q 3,3 (when k1 ≤ 0).

Following the same scheme as before, to estimate the distance from F m (x0 ) to q 3,3 we write # "  xm 2 xm 1 1 m m m 2 m 2 m 2 m m 2 m 2 k(x1 , x2 , x3 ) − (0, 0, 1)k = (x1 ) + (x2 ) + (x1 + x2 ) = 2(x2 ) 1 + m + x2 xm 2

Journal Pre-proof 16 a

b

0.8

= 10. x4

0.7 0.6 0.5

= 10. x4

x1,2,3

lP repro of

0.25 0.2

0.15

0.1

0.05

0

= (x1

5000

10000

15000

20000

iterations

FIG. 6. (a) Evolution of the invariant curve in a projection of the phase space (x1 , x4 ) as k1 → 0 using: k1 = 0.8 (black); k1 = 0.4 (red); k1 = 0.2 (blue); k1 = 0.1 (green); k1 = 0.05 (magenta); and k1 = 0.01 (orange). Insets: (orange) k1 = 0.01; (violet) k1 = 10−4 ; and (green) k1 = 10−5 . Here we have set k2,3,4 = 1. (b) Time series on the attractor for x4 (blue) and x1 (black), x2 (red), x3 (green), with k1 = 0.1, k2 = 0.9, k3 = 0.6, and k4 = 0.8. In all panels we used C = 10.

so that the asymptotic behaviour depends on how {xm 2 } tends to 0. We first prove, as in the previous case, that m {xm 1 /x2 } converges to 0. Next we prove that 1 1 + ((k3 +

εν4 )/C))xj2

≤

1 + k2 xm 1 1 ≤ C + φ(xm ) C + ((k + εν3 )/C))xm 3 2 2

for m ≥ m0 , with ν3 = k1 − 2k2 − k3 and ν4 = k2 . We then check that the number m of iterations to converge from x0 to a ball of radius δ centered at q 3,3 satisfies √ √ C C 0 0 ( 2/δ − 1/xm ( 2/δ − 1/xm 2 ) + m0 < m < 2 ) + m0 . k3 + εν4 k3 + εν3 6.

269

Convergence to q 2,2 (when k1 = 0).

rna

This case is very particular since the map is one dimensional. Written in terms of x1 it has the form xm+1 = 1

For m ≥ m0 we have

C m m ) x1 . C + k2 xm (1 − x 1 1

1 1 m+1 m xm ≤ 1 ≤ x1 m x1 . 1 + (k2 /C)xm 1 + (k (1 − ε)/C)x 2 1 1

Jou

Arguing in a similar way, the number of iterations satisfies √ C √ C 0 0 ( 2/δ − 1/xm ( 2/δ − 1/xm 1 ) + m0 < m < 1 ) + m0 . k2 k2 − ε

270

271 272 273 274 275 276

E.

Invariant curve and study of bifurcations for n = 4

As previously mentioned and, as a difference from time-continuous models (where oscillations appear for n ≥ 5 [43, 44, 54]), the dynamics of the map F defined in (1) for n = 4 and ki > 0 (i = 1, ..., 4) is governed by an invariant curve [55, 56]. The bifurcation diagrams in Fig. 4 display how the local maxima and minima of all the variables, obtained from time series once the invariant curve has been reached, change at decreasing k1 . For 0 < k1 ≤ 1 the dynamics is governed by self-sustained, periodic oscillations (see also Fig. 6). Figure 4 also displays how the invariant curve changes within the range 0 < k1 ≤ 1. The invariant curve shrinks to q 4,4 = (0, 0, 0, 1) as k1 → 0. This change

Journal Pre-proof

≥ k1

lP repro of

17

≥ k1

FIG. 7. The two eigenvalues, λ1,2 , of the fixed point resulting from the study of the invariant curve on a Poincar´e-like section for the case n = 4 (see Section III E). Note that the two eigenvalues are smaller than one within the range 0 < k1 ≤ 1, indicating that the invariant curve is stable.

279

280 281 282 283 284 285 286 287 288 289

290 291 292 293 294 295 296 297 298 299

300 301 302 303 304 305

306 307 308 309 310 311 312

in the size of the invariant curve can be visualised in Fig. 6(a), where projections of the attractor in the phase space (x1 , x4 ) are shown for decreasing values of k1 . Note that the invariant curve changes in size until it collapses at k1 = 0 (see also Fig. 4). Figure 6(b) displays time series setting k1 = 0.1, k2 = 0.9, k3 = 0.6, and k4 = 0.8. So far, the invariant curve when n = 4, described in Ref. [55], was obtained by numerical iteration. Also, the emergence of periodic oscillations for this hypercycle dimension were provided by the presence of a ’Hopf’ bifurcation [56] in the asymptotic limit C → ∞. Nowadays many authors call Neimark-Sacker to the bifurcation of families of maps analogous to the Hopf bifurcation for differential equations. The collapse of the invariant curve at k1 = 0 is through a degenerate transcritical-Neimark-Sacker bifurcation different from the one found in [56] when C → ∞. There are several methods to look for invariant curves (and invariant tori). See [14–17] for description and history of these methods. They are based either on conjugating the map to a rotation (parameterization method), on studying the iterations that fall in a thin region (slices method) or interpolating the map in some way. Our results have been obtained using a method based on interpolation similar, but simpler, to the one proposed in [11]. A further elaboration in a much more sophisticated way is found in [59]. To compute the invariant curve we choose a suitable transversal section M (depending on the parameters) close to the expected invariant curve. We choose it as a hyperplane (intersected with S 3 ) determined by the first variable x1 fixed at x1 = xh1 . Since the invariant curve should be not so far from the inner fixed point p4 , we take xh1 as the first component of p4 . Since the domain of the map is S 3 we will work with the variables x2 , x3 , the variable x4 being recovered from x1 + x2 + x3 + x4 = 1. Now, given a point x0 ∈ M , we iterate it until the iterations cross M in the same sense as x0 goes to F (x0 ). This means the second time they cross M . We consider the previous three iterations before reaching M and the three ones after crossing it. To obtain a point in M we interpolate the six points by a (vector) polyomial p(t) = (p1 (t), p2 (t), p3 (t), p4 (t)), and then look for t∗ such that P1 (t∗ ) = xh1 ; solving the equation using Newton’s method. Then p(t∗ ) ∈ M . We call G : M → M the map that sends x0 to p(t∗ ) obtained by the previous procedure. It can be seen as a pseudo Poincar´e map. We emphasise that it is a two dimensional map.

rna

278

Jou

277

Next we look for a fixed point of G by using Newton’s method, approximating the derivatives numerically by the central difference quotient. In this way we have an approximation of a point on the invariant curve. Iterating this point we recover it. In our example two iterations are sufficiently close so that the polynomial interpolation gives a good local representation of the curve. Moreover, the derivative of G at the fixed point provides a good estimate of the hyperbolicity of the invariant curve. The corresponding eigenvalues, computed as a function of k1 , are displayed in Fig. 7. When k1 → 0, as we have already mentioned, the invariant curve shrinks to q 4,4 and disappears for k1 ≤ 0 in a Neimark-Sacker bifurcation. At the same time p4 collides with q 4,4 undergoing a transcritical bifurcation. All eigenvalues of DF (q 4,4 ) are 1 except 1 + k1 /C which passes from bigger to less than 1 when k1 decreases. As for p4 , for k1 > 0, DF (p4 ) has two eigenvalues bigger than 1 and one less than 1. For k1 < 0, all its eigenvalues are bigger than 1 (note that in this case p4 no longer belongs to S n−1 ). Also, q 4,4 belongs to the line of fixed points {(0, α, 0, 1−α) |, α ∈ R}. Moreover, just at the bifurcation (k1 = 0) a new line of fixed points {(α, 0, 0, 1−α) |, α ∈ R} containing q 4,4 appears, making the bifurcation even more degenerate.

Journal Pre-proof 18 IV.

313

CONCLUSIONS

361

ACKNOWLEDGMENTS

317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359

362 363 364 365

rna

316

Jou

315

lP repro of

360

Hypercycles have been a subject of intensive research within the last 40 years. This theory has become of paramount importance since it suggests a plausible path towards the origins of life from biochemical self-organisation [1–3, 5]. One of the most important properties of hypercycles is their potential to overcome the so-called error threshold, suggested to be a major constraint in the increase of complexity of the first self-replicating systems in prebiotic ages [1–3]. The hypercycle may allow the stable coexistence of all its members, and thus larger information contents could be stored, as a difference from self-replicating, non-cooperative systems, in which the survival of the fittest may limit species’ coexistence and thus genetic diversity [1, 2, 4]. It has been suggested that catalytic RNAs (i.e., ribozymes) could have been the first self-replicating systems in prebiotic evolution [30–33]. RNAs are good candidates since this macromolecules are known to have catalytic activities [26–28, 35–40] as well as the capacity of genetic information storage. The dynamics and stability of catalytic networks is largely determined by its graph structure [4]. For example, several works have investigated the impact of catalytic parasites (i.e., replicators receiving catalytic aid but not providing catalysis to the next members of the cycle) in hypercycles persistence [47–49, 52]. Also, the so-called catalytic short-circuits [50, 51], although less explored, have been studied to determine its impact on hypercycles’ persistence. In this contribution, we have analysed a different scenario in which a functional shift in a given species changes the cooperative interaction to an antagonistic one. Specifically, we have studied small hypercycles in which a heterocatalytic interaction shifts to a density-dependent degradation (trans-cleaving activity). Several experimental studies have described trans-cleaving activities in ribozymes [39, 57, 58]. Despite hypercycle dynamics have been widely investigated, most of the research has been performed using timecontinuous approaches [1, 43–46, 48, 51, 52]. Only few discrete-time hypercycle systems have been explored [55, 56]. We here have considered that discrete-time hypercycle introduced by Hofbauer [55]. However, here, in contrast with [55], we have investigated how functional shifts impact the dynamics of small hypercycles with n = 2, 3, 4 species. Fixed points and stability analyses are developed for these systems. In this discrete-time setting, hypercycles with n = 4 display an oscillatory state allowing the coexistence of all the species via an invariant curve, while smaller hypercycles achieve coexistence via an interior fixed point. We provide a proof for the ω−limit of hypercycles when one replicator undergoes directed degradation, shown to be given by the out-competition of all the cooperative species by the one conducting the degradation. This functional change from cooperation to directed degradation makes the hypercycle become more similar to a catalytic chain. Our results are in agreement with previous research describing the impossibility of replicators’ coexistence in linear catalytic chains [4]. The convergence times to the fixed points have been analytically obtained and the relevant parameters in the n asymptotic expressions Pn identified. Concretely, we have obtained that for the points p , n = 2, 3, the times are proportional to C( i=1 k1i ), n = 2, 3, and for the points q n,n , n = 2, 3, 4, the times are proportional to C/kn . Numerical computations confirm the results and illustrate the behaviour. We have also described the bifurcations tied to the functional shift in one of the replicators. For cases n = 2, 3 a transcritical bifurcation is responsible for the extinction of the hypercycle. When n = 4, the analytical/numerical computations lead us to conclude there is a degenerate transcritical–Neimark-Sacker bifurcation when k1 → 0 as described at the end of Section III E. We emphasize that this bifurcation is different from the one described in [56] that occurs when C → ∞. As mentioned in the Introduction, hypercycle equations have been used to model the dynamics of different nonlinear systems such as cooperativity in ecosystems [5, 18], virus replication [8–10, 12], and, more recently, experimentallybuilt synthetic systems using bacteria [22] and yeast [21]. We want to notice that our contribution, albeit carrying a deep mathematical background, aims to model the changes introduced by functional shifts that can occur in molecular replicators by mutation processes. Indeed, functional shifts are found in ecological systems and are usually caused by behavioural or environmental changes. We are here focusing on changes in ribozymes switching their phenotype from the cooperative to the degradative one (self-cleaving). In terms of complex ecosystems, such functional shifts can be given by transitions between cooperation and competition. These shifts have been described in plants in semiarid ecosystems (the so-called stress-gradient hypothesis), in which facilitation (cooperation) increase as resources (e.g., water availability) decrease [62].

314

This work has been partially funded by the “Mar´ia de Maeztu” Programme for Units of Excellence in R&D (MDM2014-0445) and by the CERCA Programme of the Generalitat de Catalunya. JS has been funded by a “Ram´ on y Cajal” contract (RYC-2017- 22243) and by the Spain’s “Agencia Estatal de Investigaci´ on” grant RTI2018-098322-BI00. The research leading to these results has been also funded by ”la Caixa” Foundation. EF has been supported by

Journal Pre-proof 19

369 370 371 372

[1] [2] [3] [4] [5] [6]

373 374 375

[7] [8]

376 377 378

[9] [10]

379 380

[11]

381 382

[12]

383 384

[13]

385 386 387

[14]

388 389

[15]

390 391

[16]

392 393

[17]

394 395

[18]

396 397 398

[19] [20]

399 400

[21]

401 402

[22]

403 404 405

[23] [24]

406 407

[25]

408 409 410

[26] [27]

411 412

[28]

413 414

[29]

415 416 417 418

[30] [31] [32]

419 420

[33]

421 422 423 424

[34] [35]

M. Eigen, P. Schuster, The Hypercycle. A Principle of Natural Self-Organization. Springer-Verlag, Berlin, 1979. M. Eigen, Selforganization of matter and the evolution of biological macromolecules, Naturwiss. 58(10) (1971) 465-523. S.A. Kauffman, The origins of order: self-organization and selection in Evolution, Oxford University Press, Inc. 1993. B. Stadler, P. Stadler, Molecular Replicator Dynamics, Adv. Complex Syst. 6(01) (2003) 47-77. J.M. Smith, E. Szathm´ ary, The major transitions in evolution, Oxford University Press, 1995. M.A. Cohen, S. Grossberg, Absolute stability and global pattern formation and parallel memory storage by competitive neural networks, IEEE Trans Syst Man Cybernet 13 (1983) 815-826. M.A. Cohen, S. Grossberg, Pattern recognition by self-organizing neural networks, MIT Press, Cambridge, MA, 1991. E. Szathm´ ary, Natural selection and dynamical coexistence of defective and complementing virus segments, J. Theor. Biol. 157 (1992) 383-406. E. Szathm´ ary, Co-operation and defection: playing the field in virus dynamics, J. Theor. Biol. 165 (1993) 341-356. J. Sardany´es, S.F. Elena, Error threshold in RNA quasispecies models with complementation, J. Theor. Biol. 265 (2010) 278-286. C. Sim´ o, Effective computations in celestial mechanics and astrodynamics in Modern methods of analytical mechanics and their applications (Udine, 1997), 55–102, CISM Courses and Lect., 387, Springer, Vienna, 1998. M. Eigen, C.K. Biebricher, M. Gebinoga, The Hypercycle. Coupling of RNA and Protein Biosynthesis in the Infection Cycle of an RNA Bacteriophage, Biochem. 30 (1991) 11005-11018. J.D. Farmer, S.A. Kauffman, N.H. Packard, A.S. Perelson, Adaptive dynamic networks as models for the immune system and autocatalytic sets, in: Perspectives in biological dynamics and theoretical medicine, Annals of the New York Academy of Sciences, 1987. ` Haro, M. Canadell, J.L. Figueras, A. Luque, J.M. Mondelo, The parameterization method for invariant manifolds. From A. rigorous results to effective computations. Applied Mathematical Sciences, 195. Springer, 2016. ` Haro, R. de la Llave, A parameterization method for the computation of invariant tori and their whiskers in quasiA. periodic maps: rigorous results. J. Differential Equations 228(2) (2006), 530–579. Schilder, Frank; Osinga, Hinke M.; Vogt, Werner Continuation of quasi-periodic invariant tori. SIAM J. Appl. Dyn. Syst. 4(3) (2005), 459–488. C. Froeschl´e, Numerical study of dynamical systems with three degrees of freedom. I. Graphical displays of four-dimensional sections. Astronom. and Astrophys. 4 (1970) 115–128. J. Sardany´es, The Hypercycle: from molecular to ecosystems dynamics, in: Landscape ecology research trend, Nova Publishers, 2009. M.A. Nowak, D.C. Krakauer, The evolution of language, Proc. Natl. Acad. Sci. U.S.A. 96(14) (1999) 8028-8033. D.H. Lee, K. Severin, Y. Yokobayashi, M.R. Ghadiri, Emergence of symbiosis in peptide self-replication through a hypercyclic network, Nature 390 (1997) 591-594. W. Shou, S. Ram, J.M.G. Vilar, Synthetic cooperation in engineered yeast populations, Proc. Natl. Acad. Sci. U.S.A. 104(6) (2007) 1877-1882. D.R. Amor, R. Monta˜ nez, S. Duran-Nebreda, R.V. Sol´e, Spatial dynamics of synthetic microbial mutualists and their parasites, PLoS Comput. Biol. 13(8) (2017), e15689. U. Niesert, D. Harnasch, C. Bresch, Origin of life between scylla and charybdis, J. Mol. Evol. 17 (1981) 348-353. D.A. Usher, A.H. McHale, Hydrolyitc stability of helical RNA: A selective advantage for the natural 30 , 50 − bond, Proc. Natl. Acad. Sci. U.S.A. 73(4) (1976) 1149-1153. R. Lohrmann, L.E. Orgel, Self-condensation of activated dinucleotides on polynucleotide templates with alternating sequences, J. Mol. Evol. 14 (1979) 243-250. T.R. Cech, The chemistry of self-replicating RNA and RNA enzymes, Science 236 (1987) 1532-1539. J.A. Dar` os, Eggplant latent viroid: a friendly experimental system in the family Asunviroidae, Mol. Plant Pathol. 17 (2016) 1170-1177. R.M. Jimenez, J.A. Polanco, A. Lupt´ ak, Chemistry and biology of self-cleaving ribozymes, Trends Biochem. Sci. 40(11) (2015) 648-661. K. Kruger, P.J. Grabowski, A.J. Zaug, J. Sands, D.E. Gottschling, T.R. Cech, Self-splicing RNA: autoexcision and autocyclization of the ribosomal RNA intervening sequence of Tetrahymena, Cell 31(1) (1982) 147-157. M. Neveu, H.J. Kim, S.A. Benner, The ”strong” RNA world hypothesis: fifty years old, Astrobiology. 13(4) (2013) 391-403. J. A. Doudna, T.R. Cech, The chemical repertoire of natural ribozymes, Nature 418 (2002) 222-228. M.P. Robertson, G.F. Joyce, The origins of the RNA world, Cold Spring Harbor Perspectives in Biology. 4(5) (2012) a003608. A.V. Vlassov, Mini-ribozymes and freezing environment: a new scenario for the early RNA world, Biogeosci. Discuss. 2 (2005) 1719-1737. F.H.C. Crick, A speculation on the origin of protein synthesis, Orig. Life 7(4) (1976) 389-397. B. Zhang, T.R. Cech, Peptidyl-transferase ribozymes: trans reactions, structural characterization and ribosomal RNA-like features, Chem Biol. 5(10) (1998) 539-53.

lP repro of

368

rna

367

grants MTM2016-80117-P (MINECO/FEDER, UE) and 2017-SGR-1374 (Catalan Autonomous Government).

Jou

366

Journal Pre-proof 20

428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470

lP repro of

427

[36] M.D. Been, E.T. Barford, J.M. Burke, J.V. Price, N.K. Tanner, A.J. Zaug, T.R. Cech, Structures involved in Tetrahymena rRNA self-splicing and RNA enzyme activity, Cold Spring Harb. Symp. Quant. Biol. 52 (1987) 147-157. [37] W.S. Zielinksi, L.E. Orgel, Oligomerization of activated derivatives of 30 −amino−30 −deoxyguanosine on poly(C) and poly(dC) template, Nucl. Acid Res. 13(24), (1985) 8999-9009. [38] G. Von Kiedrowski, A self-replicating hexadeoxynucleotide, Angew. Chem. Int. Ed. 119(10) (1986) 932-935. [39] A. Carbonell, R. Flores, S. Gago, Trans-cleaving hammerhead ribozymes with tertiary stabilizing motifs: in vitro and in vivo activity against a structured viroid RNA, Nucl. Acid Res. 39(6) (2011) 2432-2444. [40] C.E. Weinberg, Z. Weinberg, C. Hammann, Noverl ribozymes: discovery, catalytic mechanisms, and the quest to understand biological functions, Nucl. Acid Res. 47(18) (2019) 9480-9494. [41] C.E. Weinberg, Z. Weinberg, C. Hammann, In-ice evolution of RNA polymerase ribozyme activity, Nature Chem. 5(12) (2013) 1011-1018. [42] N. Vaidya, M.L. Manapat, I.A. Chen, R. Xulvi-Brunet, E.J. Hayden, N. Lehman, Spontaneous network formation among cooperative RNA replicators, Nature 491(18) (2012) 72-77. [43] D.A.M.M. Silvestre, J.F. Fontanari, The information capacity of hypercycles, J. Theor. Biol. 254(4) (2008) 804-806. [44] A. Guillamon, E. Fontich, J. Sardany´es, Bifurcations analysis of oscillating hypercycles, J. Theor. Biol. 387 (2015) 23-30. [45] G. Farr´e, J. Sardany´es, A. Guillamon, E. Fontich, Coexistence stability in a four-member hypercycle with error tail through center manifold analysis, Nonlinear Dyn. 90 (2017) 1873-1883. [46] J Sardany´es, R.V. Sol´e, Bifurcations and phase transitions in spatially extended two-member hypercycles, J. Theor. Biol. 243 (2006) 468-482. [47] M.C. Boerlijst, P. Hogeweg, Spiral wave structure in prebiotic evolution: Hypercycles stable against parasites, Physica D 48 (1991) 17-28. [48] J. Sardany´es, R.V. Sol´e, Spatio-temporal dynamics in simple asymmetric hypercycles under weak parasitic coupling, Physica D 231 (2007) 116-129. [49] M.C. Boerlijst, P. Hogeweg, Spatial gradients enhance persistence of hypercycles, Physica D 88 (1995) 29-39. [50] J. Sardany´es, J.T. L´ azaro, A. Guillamon, E. Fontich. Full analysis of small hypercycles with short-circuits in prebiotic evolution, Physica D 34 (2017) 90-108. [51] P.J. Kim, H. Jeong, Spatio-temporal dynamics in the origin of genetic information, Physica D 203 (2005) 88-99. [52] M.B. Cronhjort, Cluster compartmentalization may provide resistance to parasites for catalytic networks, Physica D 101 (1997) 289-298. [53] J. Sardany´es, R.V. Sol´e, Delayed transitions in non-linear replicator networks: About ghosts and hypercycles, Chaos, Solitons & Fractals 31 (2007) 305-315. [54] J. Hofbauer, J. Mallet-Paret, H. L. Smith, Stable periodic solutions for the hypercycle system. J. Dyn. Diff. Eq. 3(3) (1995) 423-436. [55] J. Hofbauer, A Difference Equation Model for the Hypercycle, SIAM J. Appl. Math. 44 (1984) 762-772. [56] J. Hofbauer, G. Iooss, A Hopf Bifurcation Theorem for Difference Equations. Approximating a Differential Equation, Monatsh. Math. 98 (1984) 99-113. [57] O.C. Uhlenbeck, A small catalytic oligoribonucleotide, Nature 328 (1987) 596-600. [58] J. Haseloff, W.L. Gerlach, Simple RNA enzymes with new and highly specific endoribonuclease activities, Nature 334 (1988) 585-591. [59] V. Gelfreich, A. Vieiro, Interpolating vector fields for near identity maps and averaging, Nonlinearity 31(9) (2018) 42634289. [60] R. M. Gray, Toeplitz and Circulant Matrices: A Review, Foundations and Trends in: Communications and Information Theory 2(3) (2006) 155-239. [61] V.I. Arnol’d, Geometrical methods in the theory of ordinary differential equations, Springer-Verlag, New York-Berlin, 1983. [62] M.J. O’Brien, F.J. Pugnaire, C. Armas, S. Rodriguez-Echeverria, C. Sch¨ ob, The shift from plant-plant facilitation to competition under severe water deficit is spatially explicit. Ecol Evol 31(9) (2018) 4263-4289.

rna

426

Jou

425

Journal Pre-proof

Research Highlights

lP repro of

We explore hypercycles with a functional shift from cooperation to degradation Dynamics, bifurcations and transients are investigated

The functional shift involves the outcompetition of the cooperative species

Jou

rna

Our results have implications in origins of life and systems and synthetic Biology

Journal Pre-proof

Conflict of interest statement

lP repro of

To: the Editor of Physica D: Nonlinear phenomena

Barcelona, December 2nd 2019

Dear Editor,

The authors declare no conflicts of interest for the submitted article

‘Dynamical effects of loss of cooperation in discrete-time hypercycles’,

Dr. Josep Sardanyés (on behalf of all authors)

rna

Nonlinear Dynamics and Evolution Laboratory Centre de Recerca Matemàtica, Edifici C

Jou

Campus de Bellaterra 08193 Bellaterra, Barcelona

Journal Pre-proof

Author statement

Jou

rna

lP repro of

Júlia Perona: Model conceptualization; Mathematical calculations; Numerical simulations; Software; Writing - review and editing. Ernest Fontich: Model conceptualization,; Mathematical calculations; Methodology; Supervision; Writing- review & editing. Josep Sardanyés: Model conceptualization; Supervision; Writing - initial draft; Numerical simulations; Software, Writing - review & editing.