## Abstract

We propose and study nonlinear mathematical models describing the intracellular time dynamics of viral RNA accumulation for positive-sense single-stranded RNA viruses. Our models consider different replication modes ranging between two extremes represented by the geometric replication (GR) and the linear stamping machine replication (SMR). We first analyse a model that quantitatively reproduced experimental data for the accumulation dynamics of both polarities of turnip mosaic potyvirus RNAs. We identify a non-degenerate transcritical bifurcation governing the extinction of both strands depending on three key parameters: the mode of replication (*α*), the replication rate (*r*) and the degradation rate (*δ*) of viral strands. Our results indicate that the bifurcation associated with *α* generically takes place when the replication mode is closer to the SMR, thus suggesting that GR may provide viral strands with an increased robustness against degradation. This transcritical bifurcation, which is responsible for the switching from an active to an absorbing regime, suggests a smooth (i.e. second-order), absorbing-state phase transition. Finally, we also analyse a simplified model that only incorporates asymmetry in replication tied to differential replication modes.

## 1. Introduction

A multitude of theoretical and computational models considering different levels of detail and complexity have been proposed to better understand the dynamics of viral populations. Broadly speaking, such models can be divided into structured and unstructured. On the one hand, the structured models consider a high degree of detail in the interactions as well as in the processes governing replication and infection of viruses and have been applied to cases such as the bacteriophage T7 [1,2], human immunodeficiency virus type 1 [3], subgenomic hepatitis C virus [4], influenza A virus [5] and vesicular stomatitis virus [6]. This type of model often entails high dimensionality and a large number of parameters, which make analytical calculations and the characterization of bifurcations a difficult task. On the other hand, unstructured models can also be used to study viral dynamics in a more abstract way. Contrary to structured models, such models only consider the main processes linked to virus replication and/or infection, ignoring the details of the interactions, the cellular compartments where they take place and the different macromolecules participating in the entire virus infectious cycle. Many examples of unstructured models for viruses are found in the literature (see [7–10] and references therein). This type of model, although carrying more assumptions than the structured models, often enables detailed analytical characterization of the equilibrium points and their stability, thus providing clearer and valuable information about the role of the parameters in the overall dynamics of the system.

RNA viruses are obligate cellular parasites infecting bacteria, fungi, plants and animals. They are characterized by large populations, short generation times and high error rates [11–14]. The reproductive cycle of RNA viruses consists of several processes crucial for the success in replication and spreading of the virus [15]. Viruses can enter into the cell as complete particles, as ribonucleoprotein complexes or even as naked nucleic acid molecules. For the particular case of positive-sense single-stranded RNA (ssRNA) viruses such as those belonging to the picorna-like group (studied in this article), the next step is viral uncoating, by which the genome becomes available for translation, genome replication and the production of new viral particles. After infection, the genomic RNA is translated by the host ribosomes giving rise to the viral polyprotein, which is self-processed giving place to non-structural and structural proteins. The non-structural proteins give rise to the enzymes needed to replicate the genomes, a process conducted by an RNA-dependent RNA polymerase. After some rounds of amplification of the viral genomes, the structural proteins pack the new genomes forming the virions, which eventually will infect other available host cells.

A key but poorly understood parameter during intracellular amplification of viral genomes is the mode of genome replication. Few studies have explored the effect of the mode of replication on the population dynamics of viral genomes from a dynamical point of view [16,17], and even fewer experimental studies have investigated the dynamics of viral amplification quantitatively [18]. Nevertheless, the few available data suggest different models of replication for different viruses such as the stamping machine replication (SMR) mode. For SMR, and considering an infecting virus with a positive-sense RNA genome, the progeny of strands will be synthesized from negative strands complementary to the infecting positive-sense genome, and thus the expected fraction of mutant genomes produced per infected cell follows 1 − e^{−μ}, *μ* being the genomic mutation rate. In this case, the distribution of mutants per infected cell before the action of selection follows a Poisson distribution. Such a distribution was identified for bacteriophage *ϕ*X174 [19]. Another suggested mode of replication, opposed to the SMR, is the geometric replication (GR) mode. For GR, both positive- and negative-sense RNA strands are used as templates for viral amplification. For this mode of replication, the expected fraction of mutant genomes produced per infected cell depends on the number of replication cycles, *n*, according to the expression 1 − e^{−nμ}. The resulting distribution of mutants then follows the Luria–Delbrück distribution. Deviations from the Poisson distribution were found for the phage *T*2 [20], thus suggesting that such a virus replicates according to the GR model. Intermediate modes of replication may also exist, as described for bacteriophage *Φ*6 [21]. In this case, the distribution of mutants slightly deviated from the Poisson distribution, thus suggesting that the replication was mainly achieved by an SMR strategy plus a small contribution of GR.

A direct experimental evaluation of the mode of replication for a eukaryotic RNA virus has only been provided very recently. Martínez *et al.* [18] monitored and quantified the accumulation of both RNA polarities of turnip mosaic virus (TuMV) infecting protoplasts of the plant *Nicotiana benthamiana*. There, we developed a simple dynamical system describing the production of negative- and positive-sense RNA strands including constitutive transcription of positive strands from a plasmid and interference of positive strands on the synthesis of negative ones (see §2). Such a model was used to fit the experimental data (figure 1*a*). In the present study, we first analyse a simplified version of this model that may be useful to better understand the within-cell dynamics of positive-sense RNA viruses. We analyse in detail the dynamics of the model, studying the equilibrium points and their stability depending on the key parameters governing the intracellular dynamics of viral RNA replication, paying special attention to the role of the replication mode. In short, we characterize a non-degenerate transcritical bifurcation governing the shift from an active phase, where both viral strands coexist, to an absorbing state, for which both sequences become extinct. The critical expression responsible for the bifurcation depends on the asymmetry of replication (i.e. the replication mode), on the replication rate as well as on RNA degradation rates. Our results suggest that SMR is more sensitive to the bifurcation and thus GR confers viral genomes with a dynamical advantage. We finally analyse a simpler model that only considers the mode of replication, also studying its equilibria and stability properties.

## 2. Mathematical model

This section introduces a mathematical model that we recently used to fit experimental data on the accumulation for TuMV genomic positive and antigenomic negative RNAs [18]. The model describes the time dynamics of well-mixed populations of positive and negative strands during intracellular viral amplification considering key parameters such as different replication modes and degradation of viral strands in the limit of infinite diffusion (figure 1*b*). The model is given by the next couple of nonlinear differential equations:
2.1and
2.2

The state variables *p* and *m* denote, respectively, the concentration of positive (*p*: plus) and negative (*m*: minus) strands. Note that for simplicity, we obviate the intermediates of replication given by double-stranded RNAs as well as the explicit consideration of mutation. We consider that the replication of the viral RNA is constrained by a logistic-like function, given by
assuming finite resources. Here, *K* is the cellular carrying capacity. The parameter r > 0 corresponds to the replication rate of the strands, which is assumed to be symmetric. However, we introduce the parameter *α* (with 0 < *α* ≤ 1) in equation (2.2) to model all the scenarios of asymmetric replication between the GR mode (with *α* = 1, i.e. both strands replicate at the same rate) and the SMR mode (with *α* ≳ 0). For the latter case, the replication proceeds mainly from negative to positive strands. Finally, the model also considers, respectively, that positive- and negative-sense strands are degraded at rates *δ*_{p} > 0 and *δ*_{m} > 0.

Equations (2.1) and (2.2) provide a simple scheme to study the dynamics of the strands during intracellular viral replication considering asymmetries in the replication rates owing to different replication modes. As we mentioned in §1, this model was successfully applied to reproduce experimental data obtained during the amplification phase of TuMV in *N. benthamiana* protoplasts [18] ( figure 1*a*). In that work, we estimated the key parameters in the replication of this virus, especially the mode of replication. The parameter *τ* was introduced because in the experiments viral infection was initiated by the constitutive transcription of positive strands from an infectious cDNA containing the viral genome (figure 1). Moreover, the term 1/(1+*ψ**p*) was introduced in equation (2.2) to account for possible interference of positive strands in the synthesis of negative strands that could result in differences between the two RNA accumulation curves. The inference of parameter *ψ* (which measures the negative effect of accumulating positive strands on the synthesis of the negative ones) indicated that no major interference existed and thus that the differences in accumulation observed between the two strands emerged mainly from the different replication modes [18].

Given that interference was not significant in the TuMV system, in the following, we will not consider it (i.e. *ψ* = 0). Moreover, the strand initiating the process of viral replication for positive-sense ssRNA viruses in natural infections is the positive genome of the virus. Henceforth, in our analyses, we will set *τ* = 0, and the biologically meaningful initial conditions will be *p*(0) > 0 and *m*(0) = 0.

## 3. Results

### 3.1. Analyses with degradation of strands

In this section, we provide analytical and numerical results to characterize the dynamics of equations (2.1) and (2.2) (with *ψ* = 0, *τ* = 0 and *δ*_{m} = *δ*_{p} ≡ *δ*). We note that setting *δ*_{p} = *δ*_{m}≡ *δ* enables clearer identification of the effect of the mode of replication on the population dynamics of the two strands, which is the central question we sought to address in this study. Hence, the dynamical system we will explore is given by
3.1and
3.2

The concentration variables or population numbers span the two-dimensional open space: and only part of which is physically meaningful:

According to the logistic-like constraint (hereafter assuming *K* = 1), the biologically meaningful equilibrium points are in the triangular phase plane, *Γ*_{K} ∈ *Γ*, with

It can be shown that the system of equations (3.1) and (3.2) has three fixed points or equilibria, which are calculated from d*p*/d*t* = 0 and d*m*/d*t* = 0. The first equilibrium is trivial, given by *P**_{1} = (0,0), which involves, if stable, the extinction of both types of strands. The other two fixed points, denoted as *P**_{±} = (*p**_{±}, *m**_{±}), are given by
with
and

The next step is to study the stability of these equilibrium points, which will depend on the model parameters. To do so, we perform linear stability analysis, evaluating the Jacobian matrix at each of the equilibrium points of the system. The general form of the Jacobian matrix for the dynamical system under study is given by

The equilibria can be categorized as stable or unstable depending on the sign of the resulting eigenvalues. For a two-dimensional dynamical systems such as the present system, two negative eigenvalues mean that the fixed point is stable. However, if one or both eigenvalues are positive, then such an equilibrium will be unstable, and initial conditions in the immediate vicinity of this fixed point will move away from this equilibrium. The general expressions of the eigenvalues, obtained from det(*J* − *λ**I*) = 0, take the form
with

The corresponding eigenvectors are

We can study the qualitative behaviour of equations (3.1) and (3.2), especially for those scenarios of extinction and persistence of both strands. To do so, and for simplicity, we first study the stability of the trivial fixed point, *P*_{1}^{*} (0, 0), with linear stability analysis, while the stability of the other equilibria is numerically studied (all numerical results shown in this study are obtained using the fourth-order Runge–Kutta method with a time stepsize *Δ**t* = 10^{−2}).

It can be shown that the two eigenvalues evaluated at the trivial equilibrium point (i.e. computed from det(*J*(0) − *λ**I*) = 0) are , and the eigenvectors are . Note that *λ*_{−} is always negative because of the positivity of all parameters. Hence, the stability of the equilibrium point *P*_{1}^{*} entirely depends on *λ*_{+}. From this second eigenvalue, we can compute the critical value of the mode of replication, *α*_{c}, associated with the changes in the stability of the extinction equilibrium given by *P*_{1}^{*}. Such a critical value is given by
3.3

From the previous calculations, it follows that if *α* < *α*_{c} then *λ*_{+} < 0, and the origin will be stable because *λ*_{±} < 0, and both strands will become extinct. However, if *α* > *α*_{c} then *λ*_{+} > 0, and the origin will be a saddle and the flows will be attracted by the equilibrium point *P*_{+}^{*}, which is a stable node (figure 2). Figure 2 shows that when the replication mode gets closer to the SMR (*α* → 0), the asymptotic concentration of positive-sense strands increases much more pronouncedly than the concentration of negative-sense strands. Moreover, when the system is analysed for values of *α* closer to 1 (i.e. closer to a purely GR), the equilibria of both strands get closer. The bifurcation diagram of figure 2 indicates that the decrease in *α* below the critical value given by equation (3.3) results in the extinction of the two strands (see thin lines in the bifurcation diagram of figure 2). The presence of a critical value of the mode of replication responsible for the switch between coexistence to extinction allows such a change to be interpreted as a standard *absorbing-state phase transition*. This is a class of non-equilibrium transition in which the system crosses from an active to an absorbing phase, by varying the control parameter. Once the absorbing phase is achieved, the system remains in such a phase forever, with no possibility to escape [22] (see the description of the bifurcation below). For the selected parameters (*r* = 0.1249 and *δ* = 0.05), *α*_{c} = 0.1602 …. If *α* < *α*_{c}, the populations of strands cannot self-maintain and becomes extinct. That is, a decrease in *α* entailing the approach to the SMR mode causes the extinction of the whole population of viral strands.

Similarly, as we did for the parameter tied to the mode of replication, we can also obtain the critical values of replication rate as well as of degradation rate responsible for the extinction of the strands. The critical values for these two parameters are given, respectively, by 3.4and 3.5

If *r* > *r*_{c} : *λ*_{+} > 0, and the trivial equilibrium point will be a saddle point (recall that *λ*_{−} < 0). If *r* < *r*_{c} : *λ*_{+} < 0, and then the origin will be stable because both eigenvalues are negative. For the degradation parameter, the trivial fixed point will be stable if *δ* > *δ*_{c}. However, if *δ* < *δ*_{c}, *λ*_{+} > 0 and then the trivial equilibrium will be a saddle. Under this scenario, the two RNA strands will achieve a non-trivial steady state. Such steady state is given by the fixed point *P*_{+}^{*}. Actually, the third fixed point, *P*_{−}^{*}, is outside the phase plane *Γ*: the second coordinate of the equilibrium point *P*_{−}^{*}, which corresponds to the equilibrium concentration of the negative-sense strands, is always negative under the biologically meaningful ranges of the parameters. Note that under the ranges 0 < *α* < 1 and , the denominator of the coordinate *m**_{−} will be always positive. Hence *m**_{−} will be negative always.

The dynamics of equations (3.1) and (3.2) are illustrated in figures 2 and 3 by means of the solutions obtained numerically. Some illustrative examples of the dynamics arising using differential values of *α*, parametrizing differential replication modes, are shown in figure 2. As previously mentioned, the equilibrium concentration between both positive and negative strands largely differs when *α* is near 0, and thus replication proceeds via SMR. As *α* grows towards 1, the concentration of both strands become more similar, being equal for *α* = 1, which actually correspond to the pure GR. In the analyses shown in figure 2, the degradation rate of the strands is below the critical value, and the equilibrium points of *Γ*_{K} are an unstable fixed point (the origin) and the stable node *P*_{+}^{*}. The phase portraits shown in figure 2*a*–*d* display the effect of increasing *α*: the stable node travels towards symmetric equilibrium values for both strands. We note that the dynamics associated with the meaningful equilibria for equations (3.1) and (3.2) (given by the origin and *P*_{+}^{*}) are independent of the initial conditions. It means that from any arbitrary initial condition, the flows will achieve a single attractor depending on whether the system is in the survival or in the extinction scenario.

Figure 3*a* shows the equilibrium concentration of both types of strands in the parameter space (*α*, *δ*). In agreement with the results reported above, two different scenarios are found: (i) survival and (ii) extinction of strands. Scenario (i) occurs for those values of *δ* < *δ*_{c}, and the dynamics asymptotes towards the non-trivial steady state are given by *P*_{+}^{*}. As we discussed above, for a given value of *δ* and *r*, the system will become extinct as *α* is decreased and the replication gets closer to the SMR mode. An interesting result shown in figure 3*a* is that GR confers a higher resistance to degradation. This is due to the fact that all produced strands are used as templates for further replication. On the contrary, for the SMR mode, only a small fraction of negative-sense strands are used as templates and thus the population is much more sensitive to degradation.

The results from linear stability analysis and from the bifurcation diagram shown in figures 2 and 3*b* suggest that the transition from survival to extinction of strands associated with the critical parameters previously characterized is governed by a transcritical bifurcation. The transcritical bifurcation occurs when two equilibrium points collide and interchange their stability [23]. This type of bifurcation suggests the presence of a smooth, absorbing-state phase transition because the order parameter (i.e. equilibrium values of both types of strands) decreases in a continuous way as the control parameter is driven towards its bifurcation value, and no sharp transitions occur as one might found for saddle-node bifurcations tied to first-order phase transitions. For our system, it can be shown that
3.6

Hence, as the critical value of the parameter tied to the replication mode, both fixed points *P*_{1}^{*} and *P*_{+}^{*} collide. The stability of the fixed point *P*_{+}^{*} was numerically investigated by evaluating the sign of the two eigenvalues from
with

Numerical computations of the eigenvalues for the fixed point *P*_{+}^{*} using the same parameter range shown in the bifurcation diagram of figure 2 were performed, obtaining (results not shown)
and *λ*_{+}(*P*_{+}^{*})|_{α >αc} < 0. Moreover, the study of the eigenvalues in the parameter range studied in figure 3*b* (i.e. 0 ≤ *δ* ≤ 0.06) also indicated that *λ*_{−}(*P*_{+}^{*}) was always negative and that
(results not shown). The two insets in figure 3*b* show, respectively, the phase portraits for the scenarios of survival and extinction of the strands (we show the flows in the phase plane *Γ*). The origin with *δ* < *δ*_{c} is an unstable equilibrium (indicated with a white circle), and the non-trivial equilibrium *P*_{+}^{*} is a stable node (indicated with a black circle). Once the degradation rate overcomes the critical value, the equilibrium *P*_{+}^{*} goes outside the phase plane *Γ*, and *P*_{1}^{*} becomes stable.

All the previous analyses indicate that the fixed points *P*_{1}^{*} and *P*_{+}^{*} collide at the critical parameter values and that one eigenvalue (i.e. *λ*_{+}) interchanges its stability. Under such conditions, the transcritical bifurcation is called non-degenerate. The same type of transition occurs as the other two parameters (i.e. *r* and *δ*) cross their critical values given by equations (3.4) and (3.5). From equation (3.6), we also obtain

### 3.2. Analyses without degradation of strands

Since we are mainly interested in the dependence between the viral RNA dynamics and the asymmetry in the mode of replication, we further simplify the model by only keeping parameter *α*, and assuming *r* = 1 and *δ* = 0. The parameters estimated in Martínez *et al.* [18] revealed that the degradation rates of the viral RNAs were approximately between one and two orders of magnitude lower than the replication rate. Hence, the assumption *δ* = 0 is a first good approach to study the simplest model of viral strand amplification under differential replication modes. The model is then simplified to the following couple of dynamic equations (also setting *K* = 1):
3.7and
3.8

The dynamical system above has four fixed points given by and

Note that all the points of the form (1 − *m**, *m**), that is, *p** = 1 − *m** (which is the diagonal border in *Γ*_{K}), are equilibrium points. The Jacobian matrix for this reduced model is given by
3.9

The eigenvalues of matrix (3.9) have the general form
(here also with *ω* = *m*^{2} + *α* [4 + *p*(*p**α* + 8*p* + 18*m* − 12) − 4*m*(3 − 2*m*)]). The general expressions for the respective eigenvectors are

The eigenvalues evaluated at the trivial fixed point, obtained from det(*J*(0) − *λ**I*) = 0, are ; and the eigenvectors are . Then, this equilibrium point will always be a saddle because, under the biologically meaningful range (i.e. 0 < *α*≤ 1), *λ*_{+}^{(a)} > 0 and *λ*_{−}^{(a)} < 0. As a difference from the model given by equations (3.1) and (3.2), the trivial equilibrium is always unstable because no degradation is considered. We note that the fixed points *P*_{b}^{*} = (0, 1) and *P*_{c}^{*} = (1, 0) are the extremes of the equilibria contained in . The eigenvalues of *P*_{b}^{*}, computed from det(*J*(*P*_{b}^{*}) − *λ**I*) = 0, are *λ*_{1}^{(b)} = − 1 and *λ*_{2}^{(b)} = 0, with eigenvectors *ν*_{1}^{(b)} = {1, 0} and *ν*_{2}^{(b)} = { − 1, 1}. Moreover, the eigenvalues and the corresponding eigenvectors for the fixed point *P*_{c}^{*} = (1, 0) and the line of equilibria are equal, and are given by *λ*_{1}^{(c)} = 0 and *λ*_{2}^{(c)} = − *α*, respectively, and with eigenvectors *ν*_{1}^{(c)} = { − 1, 1} and *ν*_{2}^{(b)} = { 0, 1}. Note that for the equilibria, *P*_{b,c}^{*} and , the flows are attracted in one subspace.

The previous results on linear stability analysis are illustrated in figure 4. The dependence of the equilibrium concentrations of the positive- and negative-sense strands is shown at increasing values of *α*. Here, as a difference of the previous model analysed, such concentrations are dependent on the initial conditions, because of the invariant line of infinite point attractors given by (see the phase portraits displayed in figure 4*a*–*d*). The highest difference between the equilibrium concentration of both strands is found for values of *α* close to 0, where replication proceeds closer to the stamping machine mode and there is greater production of positive-sense strands. As *α* grows towards 1, approaching the geometric mode of replication, the two equilibria approach each other. If the initial condition of positive-sense strands is low (as one might expect during viral infection) and *m*(0) = 0, the equilibrium concentrations become very close for GR. However, if the initial condition for positive-sense strands is increased, the equilibrium concentrations for both strands are not so symmetric as one might expect considering *α* → 1 (figure 4).

## 4. Conclusions

Whether the choice of the mode of replication might confer viruses with increased robustness against the accumulation of deleterious mutations has been a recent subject of debate and some studies have theoretically handled questions related to this hypothesis [17,18]. Experimental studies giving clues about the model of replication for viruses, mainly obtained from the analysis of mutant distributions, were performed many years ago. For instance, the studies of Luria [20] and Denhardt *et al.* [19] already suggested different models of replication for different viruses. Some other investigations have shown that the replication strategy in viruses as different as TuMV [18], *Φ*6 [21] and *ϕ*X174 [19] is closer to the SMR model. These findings suggest that selection may have favoured this replicative scheme operating on independent viral lineages, perhaps as a way of reducing the population mutational load, thus increasing mutational robustness. Together with the mode of replication, several mechanisms have been proposed to contribute to the robustness of RNA virus populations [24]. Robustness, defined as the constancy of the phenotype under mutations or perturbations, has been proposed to arise owing to mechanisms like *trans*-complementation or neutrality, intrinsic to virus replication, as well as due to other extrinsic mechanisms exploiting cellular buffering mechanisms like the heat-shock chaperones [25].

In this present study, we have analysed the dynamics of simple models describing the intracellular amplification dynamics of viral RNA taking into account differences in the mode of replication. We first studied a simple model obtained from a dynamical system and successfully used this to reproduce experimental data on the accumulation of positive- and negative-sense strands during the amplification phase of TuMV [18]. The conclusions of our previous study were that, for TuMV, the model of replication occurred through a mixed strategy but approximately 90 per cent of the genomes were produced via SMR. Here, we studied the equilibria and stability properties for such a model describing intracellular dynamics of positive-sense ssRNA viruses. We reported a non-degenerate transcritical bifurcation responsible for the extinction of the viral genomes. Such a bifurcation separates two different regimes, given by the coexistence of both types of strands (i.e. condition found in the experimental data) and by the extinction of strands corresponding to an absorbing state. The transition between such phases is governed by a smooth change, suggesting the presence of a second-order phase transition. Other studies have reported this type of transition for quasi-species dynamical systems modelling RNA viruses considering mutation [26,27]. Moreover, recent investigations have also characterized transcritical bifurcations for quasi-species models considering both mutation and complementation phenomena [28].

A second model only considering the mode of replication was also investigated. This model indicated that when degradation rates of the viral RNAs are not considered, there exists an invariant line of attractors resulting in the coexistence of both strands. Under these dynamics, the equilibrium values of the viral RNAs depend on the initial conditions, as a difference from the first analysed model, for which the asymptotic dynamics was independent of the initial conditions.

The results reported with the model considering viral RNA degradation reveal that as replication gets closer to the SMR model, the viral strands can suffer the bifurcation and become extinct because of increased degradation rates or decreased replication rates. Although SMR might confer RNA viruses with mutational robustness [17], our study suggests that GR could provide RNA viruses with some other dynamical advantages. For instance, with an increased resistance to stop the replication process owing to the degradation of the negative strands in an environment dominated by positive-sense strands, or to the degradation of double-stranded intermediates of replication by the RNA-silencing machinery [29]. In this sense, one might expect that RNA viruses may have evolved towards replication strategies optimizing the interplay between both mutational and dynamical robustness.

## Acknowledgements

This work was funded by the Human Frontier Science Program Organization grant RGP12/2008, by the Spanish Ministerio de Ciencia e Innovación grants BIO2008-01986 (J.A.D.) and BFU2009-06993 (S.F.E.) and by the Santa Fe Institute. F.M. is the recipient of a predoctoral fellowship from Universitat Politècnica de València. We also thank the hospitality and support of the Kavli Institute for Theoretical Physics (University of California at Santa Barbara), where part of this work was developed (grant NSF PHY05-51164).

- Received July 15, 2011.
- Accepted August 15, 2011.

- This journal is © 2011 The Royal Society