## Abstract

Evolutionary and dynamical investigations into real viral populations indicate that RNA replication can range between the two extremes represented by so-called ‘stamping machine replication’ (SMR) and ‘geometric replication’ (GR). The impact of asymmetries in replication for single-stranded (+) sense RNA viruses has been mainly studied with deterministic models. However, viral replication should be better described by including stochasticity, as the cell infection process is typically initiated with a very small number of RNA macromolecules, and thus largely influenced by intrinsic noise. Under appropriate conditions, deterministic theoretical descriptions of viral RNA replication predict a quasi-neutral coexistence scenario, with a line of fixed points involving different strands' equilibrium ratios depending on the initial conditions. Recent research into the quasi-neutral coexistence in two competing populations reveals that stochastic fluctuations fundamentally alter the mean-field scenario, and one of the two species outcompetes the other. In this article, we study this phenomenon for viral RNA replication modes by means of stochastic simulations and a diffusion approximation. Our results reveal that noise has a strong impact on the amplification of viral RNAs, also causing the emergence of noise-induced bistability. We provide analytical criteria for the dominance of (+) sense strands depending on the initial populations on the line of equilibria, which are in agreement with direct stochastic simulation results. The biological implications of this noise-driven mechanism are discussed within the framework of the evolutionary dynamics of RNA viruses with different modes of replication.

## 1. Introduction

An essential yet poorly understood process during intracellular amplification of viral genomes is the mode of genome replication. Very few theoretical studies have explored the dynamical properties of alternative modes of RNA virus genome replication [1–3], and even fewer experimental studies have collected the necessary data for making a proper description of the temporal increase of sense and antisense viral genomes, the latter being the unavoidable intermediates of replication [4]. Nevertheless, the scarce available data suggest different models of replication for different viruses. Two extreme mechanisms of RNA virus genome replication have been discussed. The first mechanism is referred to as the stamping machine replication (hereafter, SMR) mode. For SMR, and considering an infecting virus with a (+) sense RNA genome, the whole progeny of (+) sense strands will be synthesized from a few (−) antisense strands complementary to the infecting (+) one. In SMR, an asymmetric accumulation of strains of both polarities exists: (+) strands are massively produced while intermediate (−) strands stay at a low concentration. Under this mode of replication, the expected fraction of mutant genomes produced per infected cell follows 1 − e^{−μ}, where *μ* is 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 of mutants has been described for bacteriophage Q*β* [5]. The second possible mode of replication is geometric replication (hereafter, GR). GR is symmetrical and both (+) and (−) sense RNA strands are used as templates for viral amplification with equal efficiency. 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 expression 1 − e^{−nμ}. The resulting distribution of mutants then follows the Luria–Delbrück distribution. Deviations from the Poisson distribution were found for poliovirus [6], thus suggesting that this virus replicates according to the GR model. In principle, these two modes of replication represent the two extremes of a continuum of possible strategies. In this sense, intermediate modes of replication have been described for bacteriophage *ϕ*6 [7] and turnip mosaic potyvirus (TuMV) [4]. In the former 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.

Results from experiments specifically designed to determine the mode of replication for eukaryotic RNA viruses have been published only in recent years. Of particular interest, Martínez *et al.* [4] monitored and quantified the accumulation of both RNA polarities of TuMV infecting *in vitro* protoplasts of the host plant *Nicotiana benthamiana*. In this study, a simple dynamical model describing the production of (−) and (+) sense RNAs and interference of (+) strands on the synthesis of (−) ones was proposed (see §2). The model was then used to fit the experimental data and to estimate replication parameters; most remarkably, the fraction of RNA viral molecules produced by SMR and GR. Inspired by this study, Combe *et al.* [8], working with the (−) sense RNA vesicular stomatitis virus (VSV), and Schulte *et al.* [6], working with the (+) sense poliovirus (as mentioned above), showed that GR is the main mode of replication for these very different viruses.

Stochasticity is an essential component of molecular biosystems, with virus replication being a paradigmatic case. Firstly, infection is generally initiated by one or very few genomes per cell (a variable known as multiplicity of infection or MOI). Secondly, these genomes have to bind to cellular components in a crowded within-cell environment to begin replication. Accordingly, the outcome of an infection is likely to be affected by the variability in the initial molecular interactions between the virus and the host cell's resources, resulting in heterogeneity among infected cells (e.g. [8,9]). Viruses have developed mechanisms to minimize the effect of noise by means of actively controlling MOI [10] or limiting replication to well-defined and space-limited membrane-associated compartments known as viral replication factories [11].

The impact of stochasticity on the fate of replicator systems is a current subject of interest, especially, the effects of noise on the so-called ‘quasi-neutral’ or ‘degenerate scenarios’, found in the deterministic limit [12–14]. Such scenarios can be characterized by a line of steady states, the long time behaviour of the system (i.e. the final composition of the population) being determined by the initial conditions. It is known that, when noise is being considered, the dynamics is dominated by the noise-driven drift on such a quasi-neutral line, whereby the composition of the population slowly and randomly changes until eventual outcompetition of one of the species occurs [12,13]. In this article, we analyse a simplified version of the model introduced in Martínez *et al.* [4], used to quantify the within-cell replication dynamics of (+) sense RNA viruses. This model displays a quasi-neutral coexistence scenario between (+) and (−) RNA strands when only considering the strands' replication and competition (see [3] for details). We provide a detailed analysis of the quasi-neutral dynamics of this model, focusing on the impact of stochastic drift on the invariant manifold which determines the fate of viral RNAs. Based on a diffusion approximation of the stochastic dynamics, we present analytical results regarding the ability of the system to achieve a given asymptotic state as a function of the initial conditions on the quasi-neutral line; that is, the fate of the viral RNA genomes under strong stochastic fluctuations and under differential modes of replication. The analytical results are complemented with stochastic simulations. Finally, we describe a novel mechanism for noise-induced bistability.

## 2. Mathematical models

### 2.1. Deterministic model

A simple mathematical model considering asymmetric RNA replication modes has been recently introduced in [3] (see also [4]). This model is given by the differential equations 2.1and 2.2

Equations (2.1) and (2.2) are deterministic rate equations for two strands of viral RNA with different polarities (*m*: antisense or (−) and *p*: sense or (+)). Strands compete for the available cellular resources and differ only in the time scales of their evolution (birth and death rates), defined by the density amplification rates *αr* and *r*. The constant *C* is the carrying capacity. Since we are interested in the impact of the replication mode, we will focus our study on parameter *α*. When *α* = 1 replication proceeds via GR as both strands are replicated at the same rate. The cases *α* → 0 follow the SMR mode, because replication becomes asymmetric and there is much more production of (+) sense RNAs from the (−) one(s), which act as template(s). The two strands' common carrying capacity *C* indicates the total number of individuals in the non-empty steady state.

Here we focus on the case *δ* = 0, which relates to a particular topology of the phase space, given by a neutrally stable invariant line (see below). The points on such a line correspond to the steady-state composition of the viral population, (*m**, *p**), which is determined by the initial conditions. The parameters estimated in [4] from experimental data using a similar model to equations (2.1) and (2.2) revealed that the degradation rates of the viral RNAs were between approximately one and two orders of magnitude lower than the replication rates. Hence, the assumption *δ* = 0 is a good first approach to studying the simplest model of viral strand amplification under differential replication modes. The deterministic continuum description is presumably applicable when *m* and *p* are and *C* ≫ 1. Let us, for convenience, introduce the scaled population variables *x* = *m*/*C* and *y* = *p*/*C* and rescale the time variables by the growth rate *r*, to write the system
2.3and
2.4where the ratio of time scales is *γ* = (*αr*)/*r* = *α*.

### 2.2. Stochastic model

The first level of description involves considering the stochastic evolution of the integer-valued random processes *X*_{t} and *Y*_{t}. Let us denote *n* and *p* as the number of (−) and (+) sense viral RNA genomes, and *P*(*n*, *p*, *t*) as the probability of having *n* and *p* strands at time *t*. We can analyse the time change of such probabilities considering a Markov model, given by the master equation (ME),
2.5In the above equation, *R* is the number of reactions, and *W*_{i} and *r*_{i} are the transition rates (propensities) and the stoichiometry of the system, respectively. The ME (2.5) describes the following reactions (after rescaling parameters).

(1) Reactions for negative-sense strands: 2.6 2.7 2.8

(2) Reactions for positive-sense strands: 2.9 2.10 2.11

Reaction (2.6) denotes the synthesis of (−) sense strands from (+) sense ones. Reactions (2.7) and (2.8) involve the decrease of the population (*r*_{2} = *r*_{3} = −1), as they represent, respectively, competition for cellular resources between strands of the same polarity (intra-specific competition in ecological terms) and of different polarity (interspecific competition). Note that density-dependent competition as introduced here is completely different from spontaneous RNA decay (see the last terms of equations (2.1) and (2.2)), which is not accounted for in this model. The same reactions are represented for (+) sense strands in reactions (2.9)–(2.11).

The propensities for reactions (2.6)–(2.8) are given by *W*_{1} = *γp*, *W*_{2} = (*γ*/*C*)*n*(*n* − 1) and *W*_{3} = (*γ*/*C*)*np*; while the propensities for reactions (2.9)–(2.11) are given by *W*_{1} = *n*, *W*_{2} = *C*^{−1}*pn* and *W*_{3} = *C*^{−1}*p*(*p* − 1). As mentioned, we consider a Markov model where the birth and death rates are, respectively, *β*_{Y} and *β*_{X}, and and when *X*_{t} = *n* and *Y*_{t} = *p*, where with *β*_{Y}/*δ*_{X} = *β*_{X}/*δ*_{Y} = *ρ* > 1. evolves according to the ME:
2.12The low-density replication rates found in the deterministic differential equations are *αr* = *β*_{Y} − *δ*_{X} and *r* = *β*_{X} − *δ*_{Y}, and the ratio of evolution time scales is *γ* = *δ*_{X}/*δ*_{Y} = *β*_{Y}/*β*_{X}. Under suitable conditions [15–17], a diffusion approximation can be applied to ME (2.12), so that the variables *x*_{t} = *C*^{−1}*X*_{t/γ} and *y*_{t} = *C*^{−1}*Y*_{t/γ} evolve according to the stochastic differential equations (SDEs)
2.13where *B*^{x}_{t} and *B*^{y}_{t} are independent Wiener processes and the noise amplitude is *ε* = *C*^{−1/2}.

## 3. Results and discussion

### 3.1. Deterministic dynamics

The dynamics of equations (2.1) and (2.2) with *δ* = 0 was investigated in [3]. In the following, we summarize these results from the scaled equations (2.3) and (2.4). The fixed points for this system are given by the point *P*^{*}_{0} = (0, 0) and the line of fixed points . The Jacobian matrix for this model is given by
The local stability of *P*^{*}_{0} and can be studied from the eigenvalues obtained from the characteristic equation , with *f*^{*} being a fixed point and the identity matrix. The eigenvalues of the fixed point *P*^{*}_{0} are . Thus, the point (0, 0) is a saddle as *γ* > 0. The eigenvalues tied to the equilibria of the line can be computed by linearizing around (*x* * = 1 − *y* *, *y* *). It can be shown that the eigenvalues obtained from are given by
Note that for *y* * ∈ [0, 1] and *γ* ∈ (0, 1], *λ*_{2} < 0, there exists a continuous line of marginally stable fixed points. Hence, the equilibria contained in the line have a neutral eigenvalue (*λ*_{1}, which is not locally attracting) and a negative eigenvalue (*λ*_{2}, which is locally attracting). These features mean that the initial conditions outside the line are attracted towards it, and once they achieve the line they stop and remain there in forward time.

These dynamical characteristics can be visualized in figure 1. Figure 1*b* shows a phase portrait (obtained numerically using the fourth-order Runge–Kutta method) for a case close to the SMR, using *γ* = 0.06. Note that different initial conditions end up in different places on line . The changes in the mode of replication are also important as shown in figure 1*c*, where we demonstrate that different initial conditions for (+) sense strands hit within the range 0.8 ≲ *p*(*t*→ +∞) < 1 for the SMR (*γ* = 0.1). By contrast, the equilibrium values of GR systems (*γ* = 1) are found within the range 0.5 ≲ *p*(*t*→ +∞) < 1.

### 3.2. Stochastic dynamics

As shown in the previous section, is a neutrally attracting line. It is known that this neutrality is broken under the presence of noise [12,13], where the paths are rapidly attracted towards , whereby a noise-driven dynamics ensues such that the system moves along , ultimately hitting one of the two states (1, 0) or (0, 1).

#### 3.2.1. Stochastic simulations and asymptotic dynamics

The dynamics of the master equation (2.12) is simulated using the Gillespie algorithm [18,19] applied to reactions (2.6)–(2.11). Specifically, we will use two values of the carrying capacity (*C* = 500 and *C* = 10^{3}) to determine the role of stochastic fluctuations in the time evolution on . Stochastic simulations show that, away from , the dynamics of both strands typically follows the deterministic trajectories. Once the dynamics has settled onto the quasi-neutral line, stochastic fluctuations can drive the populations towards one of the corners of the phase space (*x*, *y*), involving the dominance of one of the strands over the other.

First, we provide some analyses starting from initial conditions far away from the quasi-neutral line, and, in some cases, above the carrying capacity of the system. These analyses, which may not correspond to biologically meaningful initial populations, are displayed to compare the deterministic dynamics with the stochastic simulations, as well as to analyse the dynamics in the full-phase space. However, the main goal of this article is to focus on the stochastic dynamics along the quasi-neutral line, which is always achieved under biologically meaningful initial populations (i.e. with an initial small amount of viral RNAs), depending on the mode of replication. Figure S1*a*,*b* in the electronic supplementary material shows the stochastic dynamics of the strands under two different values of *γ* and different carrying capacities (*C* = 500 in (*a*) and *C* = 10^{3} in (*b*)). For *C* = 500 and *γ* = 0.1 (black trajectory in figure S1*a* in the electronic supplementary material), the stochastic path follows the deterministic orbit (black dashed line in figure S1*a* in the electronic supplementary material); once it reaches , it evolves towards the corner (1, 0) in the phase space, which involves the dominance of the (−) sense strands over the population of (+) strands. A similar result is observed for *γ* = 1 (GR) also with *C* = 500, but the stochastic trajectory ends up in the opposite corner, where (+) strands dominate over the (−) sense strands (red trajectories in figure S1*a* in the electronic supplementary material). Similar results are obtained for *C* = 10^{3}, where the fluctuations become more narrow (electronic supplementary material, figure S1*b*).

The probability of reaching a given state (*x*(*t* ≫ 1) = 1 and *y*(*t* ≫ 1) = 0 or *x*(*t* ≫ 1) = 0 and *y*(*t* ≫ 1) = 1) has been computed as a function of the initial condition within the whole phase space *x* ∈ [0, 1] × *y* ∈ [0, 1]. Here, for each initial condition sampled regularly in the phase space, we have computed the probability *P*_{y} from a set of different replicates, by counting the number of realizations attracted to the state *x*(*t* ≫ 1) = 0 and *y*(*t* ≫ 1) = 1. The results are shown in figures S1*a*.1–*a*.3 (using *C* = 500) and S1*b*.1–*b*.3 in the electronic supplementary material, setting *C* = 10^{3}. For replication modes close to the SMR (electronic supplementary material, figure S1*a*.1,*b*.1), the values of *P*_{y} are close to 1 in a wide region of the phase space (orange-red gradient) for *x*(0) ≈ 0.5, with maximum values for the initial conditions situated in the left-upper corner of the phase space, where the initial number of (+) sense strands is large and of (−) sense ones is low. As the replication mode approaches the GR, the region where *P*_{y} ≈ 1 remains confined to the corner with initial conditions larger than *y*(0) ≈ 0.6 and *x*(0) ≈ 0.6. These results indicate that, under strong stochasticity, the SMR model would ensure the production of (+) sense strands, assuming that the initial conditions in the infection of a cell would involve *x*(0) = 0 and *y*(0) ≪ 1.

In the previous analyses, we have mainly focused on the behaviour of trajectories starting outside the quasi-neutral line. A way to determine the dynamics on the line is to define a joint variable given by *z*_{t} = *x*_{t} − *y*_{t} [12]. By doing so, one can simplify the problem and find analytical approximations for the probabilities of reaching one of the possible asymptotic states as a function of *z*; for example, approximations for the probability of dominance of (+) sense strands as a function of *z*, labelled *P*_{y}(*z*). We first present simulation results regarding *P*_{y}(*z*) to determine how this probability behaves for different replication modes.

The probability of dominance of (+) sense strands, *P*_{y}(*z*), has been computed for several values of *γ* along *z*_{0} over a total number of 25 × 10^{3} independent runs. For replication modes close to the SMR, *P*_{y}(*z*) remains above the antidiagonal in the space (*z*, *P*_{y}(*z*)), meaning that (+) sense strands experience highest fitness (see figure 2*a* with *γ* = 0.1). This effect is observed also for *γ* = 0.2 and *γ* = 0.3. For *γ* ≳ 0.4 the curve for *P*_{y}(*z*) starts crossing the antidiagonal, the intersection being at *z* < 0. The increase in *γ* towards the GR model makes the intersection happen at *z* = 0. We note that, for some values of *γ*, the curve *P*_{y}(*z*) also intersects the antidiagonal for values close to *z* = 1 or to *z* = −1. For instance, for *γ* = 0.1, *P*_{y}(*z*) intersects the antidiagonal at *z* ≈ 0.7. This phenomenon occurs on the line *z* at the corners of the phase space. In the next section, we discuss and interpret this change in the shape of *P*_{y}(*z*). The vertical red lines correspond to the theoretical predictions for the crossing values of the antidiagonal by *P*_{y}(*z*), given by equation (3.11) (see next section for further details).

Finally, following the approach in [12], we computed the mean extinction times (METs) as a function of *z* for two different values of the carrying capacity, *C* = 500 and *C* = 10^{3} (figure 3). The same qualitative results have been obtained for both values of *C*. Figure 3*a*, with *C* = 500, shows that the MET is symmetric regarding its *z*-dependence for *γ* = 1. This means that there is a maximum at *z* = 0 and that, for increasing and decreasing *z* values from *z* = 0, the METs decrease similarly (giving rise to a parabola-like shape). These times decrease as *z* grows above or below 0 because the initial conditions approach the vertices of . As we have argued previously, the switch from GR to SMR involves longer extinction times, since the amplification kinetics shifts from purely exponential to subexponential. This can be clearly seen in figure 3, where a decrease in *γ* involves longer METs. Interestingly, the maximum of the parabola displaces towards *z* > 0. This indicates that the METs decrease faster from the maximum of the parabola as *z* → 1, i.e. with larger *x*_{0} values. This phenomenon, arising from the replication model close to the SMR model, denotes the fitness advantage of SMR models for (+) sense strands.

#### 3.2.2. Random dynamics on the quasi-neutral line

The stochastic dynamics on a quasi-neutral line of fixed points has been considered previously by Lin *et al.* [12] and Kogan *et al.* [13]. Both approaches are predicated upon the assumption that, once the system quickly settles onto the quasi-neutral line, fluctuations away from the quasi-neutral line quickly relax back to it. By comparison, the dynamics on the quasi-neutral line, characterized by the variable
3.1is much slower, and it thus contains the relevant information regarding the long time behaviour of the system. In particular, we here consider the method proposed by Lin *et al.* [12], who assumed that the relaxation of the fluctuations away from the quasi-neutral line occurs along the deterministic trajectories. Such an assumption allows us to write a reduced dynamics only in terms of the variable *z*_{t} (see [12] for a detailed account). The resulting random drift-and-diffusion dynamics proceeds until the system hits one of the two boundaries, *z* = ± 1.

Consider an initial condition (*x*_{0}, *y*_{0}) on the line of fixed points . Under the action of random noise the state of the system is altered so that we can write the new state of the system as (*x*′, *y*′) = (*x*_{0} + *ϕa*, *y*_{0} + *ηb*), where *ϕ* and *η* are random variables such that 〈*ϕ*〉 = 〈*η*〉 = 0 and 〈*ϕ*^{2}〉 = 〈*η*^{2}〉 = 1 [12]. In order to be consistent with the noise terms in the SDEs given by equation (2.13), *a* and *b* are defined as
3.2and
3.3After such perturbation, the system relaxes back to the quasi-neutral line following the mean-field flow: *x*^{2} − *γy*^{2} = *cnt*, to a position given by (*x*_{0} − *ξ*, *y*_{0} + *ξ*). Since *z*_{t} = *x*_{t} − *y*_{t}, the total displacement along the coexistence line, Δ*z*_{t}, is given by *x* =−2*ξ*. Therefore, the drift *v*(*z*) and diffusion *D*(*z*) are *v*(*z*) = −2〈*ξ*〉/Δ*t* and *D*(*z*) = 4〈*ξ*^{2}〉/Δ*t*, where 〈 · 〉 represents averaging over the distributions of *η* and *ϕ*. In order to find explicit expressions of 〈*ξ*〉 and 〈*ξ*^{2}〉 as a function of the coordinate *z* along the coexistence line, we consider that and we use the equation of the mean-field trajectory: (*x*_{0} − *ξ*)^{2} − *γ*(*y*_{0} + *ξ*)^{2} = (*x*_{0} − *ϕa*)^{2} − *γ*(*y*_{0} + *ηb*)^{2}. After averaging over the probability distribution functions of *ϕ* and *η*, we obtain that
3.4and
3.5

The above equations allow us to write the steady-state backward Kolmogorov equation,
3.6whose solution is the probability of absorption at the boundary *z* = −1 starting from initial condition *z*.

Given the functional form of *v*(*z*) and *D*(*z*), a full analytical solution of equation (3.6) in closed form becomes cumbersome. However, it is still possible to define quantities that are informative regarding fitness that can actually be calculated analytically. If the two strands were neutral (i.e. if both strands were equally fit), the absorption probability would be *u*(*z*) = 1 − *z*. On the one hand, if the actual curve *u*(*z*) is such that *u*(*z*) > 1 − *z*, i.e. the probability of absorption exceeds the neutral value, then the fitness of *z* = −1 is larger than that of *z* = 1. On the other hand, for *u*(*z*) < 1 − *z*, the fitness of *z* = −1 is smaller than that of *z* = 1. Since the boundary conditions of equation (3.6) fix the value of *u*(±1), it is straightforward that *u*(*z*) > 1 − *z* if d^{2}*u*/d*z*^{2} < 0 and *u*(*z*) < 1 − *z* if d^{2}*u*/d*z*^{2} > 0. The value of the coordinate *z*_{c} where d^{2}*u*/d*z*^{2}|_{zc} = 0 marks the boundary between the regions of higher and lower fitness (see below).

We consider the cases *γ* = 0 and 0 < *γ* ≤ 1 separately.

(a) *Case* *γ* = 0. This case is not biologically meaningful but will help to understand the dynamics of the system under stochasticity. When *γ* = 0, *x*(*t*) = *cnt* ≡ *x*_{0}. The Itô differential equation is then given by
with *y* ∈ [0, 1 − *x*_{0}]. The drift and diffusion for the reduced process *y*_{t} determine the associated statistical features of the competitive exclusion dynamics. Indeed, let *τ*(*y*) = inf{*t*: |*z*_{t}| = 1 | *z*_{0} = *z*}. Then the probability that *x*_{t} reaches 0 starting from *x*_{0} before *y*_{t} starting from *z*, i.e. the probability of domination of the *Y*-species over the *X*-species, is
This probability satisfies the boundary-value problem
3.7It is easy to see that *D*(*y*) = *ε*^{2}*x*_{0}((*ρ* + 1)/(*ρ* − 1) + *x*_{0} + *y*) and *V* (*y*) = *x*_{0}(1 − *x*_{0} − *y*). Then, we obtain
3.8

(b) *Case* 0 < *γ* ≤ 1. Here we recover the reduced process obtained with *z*_{t} in equation (3.1). Since
3.9cannot be solved analytically, we restrict our calculations to solve d^{2}*u*/d*z*^{2} = 0. These calculations will allow us to determine whether the curve *P*_{y}(*z*) is above or below the antidiagonal in the space (*z*, *P*_{y}(*z*)). It can be shown that
3.10

From equation (3.10), we can obtain the value of *z* that makes d^{2}*u*/d*z*^{2} = 0, meaning that *P*_{y}(*z*) crosses the antidiagonal in the space (*z*, *P*_{y}(*z*)). The crossing value, *z*_{c}, is given by
3.11with *γ*_{0} = *γ* − 1, *γ*_{1} = 1 + *γ* and *κ* = ((1 − *γ*)(*ρ* − 1))/2*γρ*, such that

The results for equations (3.10) and (3.11) are displayed in figure 2, overlapped by the outcome of the Gillespie simulations. The analytical prediction provided by *z*_{c} on the crossing point compares well with the simulation results. However, equation (3.11) predicts only one crossing value within the relevant interval *z*∈[1, −1]. This is in contrast to our numerical results, which exhibit more than one crossing point. We posit that two possible explanations exist for this discrepancy; namely, (i) a statistical sampling effect in the computation of *P*_{y}(*z*) when *z* is close to the corners or (ii) the relaxation towards the quasi-neutral line fails to follow the mean-field trajectories.

To check whether the former explains such a discrepancy, we computed again *P*_{y}(*z*) in the corner 0.5 ≤ *z* ≤ 1 by increasing the size of the sample to 5 × 10^{6}. These simulations were carried out for the same values of *γ* as in figure 2. The same shape of *P*_{y}(*z*) is found using more replicates. Figure S2*a* in the electronic supplementary material displays the same plot as figure 2*c* (*P*_{y}(*z*) computed from 25 × 10^{3} replicates for *γ* = 0.4) and figure S2*b* in the electronic supplementary material displays the value of *P*_{y}(*z*) in the framed region computed from 5 × 10^{6} replicates. Note that the shape does not change and *P*_{y}(*z*) crosses the antidiagonal. The same results have been obtained using *γ* = 0.6 (electronic supplementary material, figure S2*b*.1); *γ* = 0.8 (triangles) and *γ* = 1 (open circles; electronic supplementary material, figure S2*b*.2); and *γ* = 0.1 and *γ* = 0.2 (results not shown). These results suggest that the lack of statistics is not responsible for the discrepancy.

According to hypothesis (ii), deviations from the mean-field model trajectories in the corners could explain the discrepancies between the numerical and the analytical results. To check this hypothesis, we have compared the stochastic and mean-field dynamics for initial conditions away from and close to the corner of the phase space. Examples are shown in the electronic supplementary material, figure S2. Simulation results with initial conditions *x*_{0} = *y*_{0} = 0.75 reveal that the stochastic trajectories follow the deterministic dynamics. Figure S2*c* in the electronic supplementary material shows four stochastic trajectories in the phase space (*x*_{t}, *y*_{t}) overlapped by the deterministic orbit. The trajectories plotted in figure S2*d*,*e* in the electronic supplementary material exhibit good agreement with the deterministic dynamics. By contrast, the same procedure with initial conditions *x*_{0} = 0.98 and *y*_{0} = 0.45 shows that the stochastic paths clearly depart from the deterministic dynamics (see electronic supplementary material, figure S2*f*–*h*). This result strongly hints that, although successfully applied to models with self-replication [12], the theory applied here does not hold for extreme values of *z* in models with complementary replicators.

#### 3.2.3. Stochastic dynamics on the quasi-neutral line drives to noise-induced bistability

Some of our previous analyses indicate the presence of noise-induced bistability in our system. A dynamical view of this phenomenon is displayed in figure 4, where the fate of a large number of replicates is displayed as a function of time. To quantify this behaviour, we have monitored the time evolution of *n* = 10^{3} sample paths, all with the same initial condition *x*_{0} = *y*_{0} = 0.5 (*z*_{0} = 0). We then calculate *P*(*z*; *t*) by computing the frequency with which each *z*-bin has been visited. Figure 4*a* shows the results for a replication strategy close to the SMR, setting *γ* = 0.1. Initially, the population is concentrated at the centre of line *z*, and, as time progresses, the distribution of *P*(*z*; *t*) spreads out. This first phase corresponds to the early stages of noise-driven dynamics on the quasi-neutral line, prior to the sample paths reaching *z* = ±1. Eventually, as the sample paths reach either end of the quasi-neutral line, *P*(*z*; *t*) becomes bimodal, with modes corresponding to the corners of the simplex (red distributions with large *t*). Hence, noise drives the population to two possible asymptotic states.

Figure 4*d* shows six trajectories on the line *z*, three of them that achieve *z*_{t} = 1 and three that achieve *z*_{t} = −1. Note that for *γ* = 0.1 there is a dominance of *P*(*z*; *t*) towards the region of *z* → −1, indicating that a larger fraction of replicates involve the dominance of (+) sense RNAs, as expected for a replication strategy close to the SMR. This effect is reversed at increasing *γ* (figure 4*b*,*c*; see also figure 4*f*, which considers GR). In agreement with figure 4, the transients for *γ* = 0.1 (SMR) are much longer than the ones for the GR model.

## 4. Conclusion

In this article, we have analysed a simple model considering differences in the mode of RNA replication. Direct [4] and indirect [5,7,20] evidence from real RNA viruses indicates that they might amplify their genomes following different strategies. On the one hand, for single-stranded (+) sense RNA viruses, the so-called stamping machine replication (SMR) means that the entire progeny of RNAs is mainly synthesized from one or a few (−) sense templates produced after the initial infection with the (+) sense strand. On the other hand, the so-called geometric replication (GR) means that all synthesized templates have the same chance of becoming templates for further replication. The evolutionary [1,2] and dynamical [3,21] implications tied to different replication modes have mainly been investigated using deterministic approaches, although some works have used Monte Carlo simulations to explore the impact of the mode of replication on important evolutionary aspects such as error thresholds [1,22] or co-infection dynamics [22].

Viral infections are usually initiated by a tiny amount of RNA molecules, and thus demographic fluctuations may have a strong impact on the dynamics of RNA amplification. In this article, we have investigated the impact of stochasticity during replication under differential replication modes. Particularly, we have focused on an interesting phenomenon given by the so-called quasi-neutral coexistence. This phenomenon has been described in two-species competition mathematical models without decay [12] as well as in the competition of strains in disease dynamics [13]. The deterministic dynamics tied to this quasi-neutral line involves a rapid approach of the trajectories to this line and, once the carrying capacity is achieved, the trajectories stop. However, in the presence of intrinsic fluctuations, once the trajectories have reached the line, a slow, noise-driven dynamics ensues.

Following the theory developed by Lin and co-workers [12], we provide analytical approximations for the probability of dominance of (+) sense strands. Scenarios favouring the accumulation or the dominance of these strands might be advantageous for (+) sense RNA viruses, which need to package the (+) sense genomes for further infection. Interestingly, we found that such a theory, developed for a model with non-complementary replication, was not accurate enough to describe the stochastic dynamics on the quasi-neutral line near the boundaries of the phase space in the system studied here.

Quasi-neutral coexistence has not been identified in real virus dynamics so far. The observation of this dynamics might depend on several, virus-specific properties, such as the reproduction cycle of the virus. For lytic viruses, the virions must lyse and destroy the cell in order to go on with further infections. This cycle may mean that, once the virus approaches the cell's carrying capacity (i.e. the maximum viral load the cell is capable of producing), the cell is lysed and the viral particles are released. Since fluctuations on the quasi-neutral line produce a slow dynamics, this dynamics might be easier to observe for non-lytic viruses (e.g. lysogenic or persistent viruses), which can remain within the host cell for a long time before it eventually exhausts the cell's resources. These viruses rely on vertical transmission and on the fact that infected cells replicate. Also, as our results show, the stochastic dynamics on the quasi-neutral line largely depends on the mode of replication. That is, this dynamics is faster under GR (see, for example, figure 3), and therefore is more likely to be observed for this model of replication.

Our results also reveal a novel type of noise-induced bistability. It is known that random fluctuations can generate novel behaviour in dynamical systems; for instance, stochastic resonance [23], noise-induced transitions [24,25] or the so-called noise-enhanced stabilization [26,27]. More recent discoveries have provided different mechanisms giving rise to the so-called noise-induced bistability. Noise-induced bistability is typically found in dynamical systems which display a single asymptotically stable state in their deterministic limit, while being able to achieve different alternative states due to stochasticity. Different examples of noise-induced stability, which can also be understood in terms of the noise-enhanced stability typical of metastable systems (see [28,29] for examples of cancer stochastic systems), have been described in many different systems [30–33].

Summarizing, we have analysed a particular scenario found in a simple model of complementary replication of viral RNAs giving rise to the so-called quasi-neutral coexistence [12,13]. Under this dynamics, stochasticity plays a major role in the asymptotic outcompetition dynamics. We have studied the probability of dominance of one strand over the other depending on the initial populations of both strands as well as on the mode of replication followed during RNA virus amplification. We have also described a novel phenomenon of noise-induced bistability in the quasi-neutral coexistence of viral genomes.

## Data accessibility

This article has no additional data.

## Authors' contributions

J.S. and S.F.E. conceived the model; T.A. and J.S performed the calculations; J.S. and A.A. performed the stochastic simulations. All authors wrote and revised the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

The research leading to these results has received funding from ‘la Caixa’ Foundation. J.S. and T.A. have been partially funded by the CERCA Program of the Generalitat de Catalunya, MINECO grant no. MTM2015-71509-C2-1-R and by a MINECO grant awarded to the Barcelona Graduate School of Mathematics under the ‘María de Maeztu’ Program (grant no. MDM-2014-0445). T.A. is also supported by AGAUR (grant no. 2014SGR1307). S.F.E. has been supported by MINECO-FEDER grant no. BFU2015-65037-P and by Generalitat Valenciana grant no. PROMETEOII/2014/021.

## Acknowledgements

We thank the members of the Computational and Mathematical Research Group at the CRM for useful discussions. We also want to thank Blai Vidiella for technical help.

## Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.4103849.

- Received February 21, 2018.
- Accepted May 8, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.