Quasispecies Spatial Models for RNA Viruses with Different Replication Modes and Infection Strategies

Empirical observations and theoretical studies suggest that viruses may use different replication strategies to amplify their genomes, which impact the dynamics of mutation accumulation in viral populations and therefore, their fitness and virulence. Similarly, during natural infections, viruses replicate and infect cells that are rarely in suspension but spatially organized. Surprisingly, most quasispecies models of virus replication have ignored these two phenomena. In order to study these two viral characteristics, we have developed stochastic cellular automata models that simulate two different modes of replication (geometric vs stamping machine) for quasispecies replicating and spreading on a two-dimensional space. Furthermore, we explored these two replication models considering epistatic fitness landscapes (antagonistic vs synergistic) and different scenarios for cell-to-cell spread, one with free superinfection and another with superinfection inhibition. We found that the master sequences for populations replicating geometrically and with antagonistic fitness effects vanished at low critical mutation rates. By contrast, the highest critical mutation rate was observed for populations replicating geometrically but with a synergistic fitness landscape. Our simulations also showed that for stamping machine replication and antagonistic epistasis, a combination that appears to be common among plant viruses, populations further increased their robustness by inhibiting superinfection. We have also shown that the mode of replication strongly influenced the linkage between viral loci, which rapidly reached linkage equilibrium at increasing mutations for geometric replication. We also found that the strategy that minimized the time required to spread over the whole space was the stamping machine with antagonistic epistasis among mutations. Finally, our simulations revealed that the multiplicity of infection fluctuated but generically increased along time.


Introduction
The dynamics and evolution of RNA virus populations is a current and important topic of research because RNA viruses are the most abundant pathogens of bacteria, humans and plants [1]. The role of these pathogens as a source of new emerging infectious diseases is also a very important subject of research in Virology and Epidemiology. RNA viruses present high population diversities, that with more or less precision are described in the virological literature as quasispecies [2][3][4][5]. A quasispecies can be roughly defined as a master sequence surrounded by a cloud of mutant genomes at the mutation-selection balance. Such a complex and polymorphic population structure may arise because of the large number of replication rounds that take place during intracellular amplification associated with the high mutation rates of the viral RNA-dependent RNA-polymerase owed to their lack of proof-reading activity [6][7][8]. Due to these peculiarities, RNA viruses have also served as excellent models for experimentally addressing important questions in evolutionary biology [9][10][11]. Several works on theoretical quasispecies [12][13][14][15][16][17][18][19][20] have been developed to understand key phenomena in virus dynamics and evolution. The convergence between theoretical and experimental results about virus dynamics and evolution is pivotal for the advance and success of future antiviral strategies [1,21].
Although new insights for theoretical quasispecies can be extracted from nonlinear dynamical models, bifurcation theory or statistical physics [12,20,22], models usually take assumptions or simplifications that jeopardize experimental validation. In this sense, a very common assumption of viral quasispecies models has been that replication follows a geometric scheme. However, empirical data suggest that viral replication may strongly depart from this model (see below). The main goal of the present work is to study differential replication modes for RNA viruses incorporating other relevant features of viral infections, such as spatial structuring of host cells, epistasis among mutations and different mechanisms of infection. The consideration of all these features into a single model framework, and especially, the consideration of differential modes of replication, may cover the gap of previously existing models. A second common assumption of theoretical quasispecies involves oversimplified and unrealistic fitness landscapes. Similarly, the consideration of determinism, or the analysis of mean field models, which do not incorporate the effect of spatial correlations, has been of common practice. The latter assumption may pose serious constraints to the interpretation of results about real viral populations replicating in spatially structured host cells, like occur in plant or animal tissues. Empirical observations suggest strong spatial structuring of different genotypes in different areas of a leaf [23][24][25] or different parts of the plant [26,27]. Similarly, the analysis of multiple samples from different tissues suggest that many animal viruses differentiate in tissue-specific subpopulations [28][29][30]. Broadly speaking, it is known that spatial correlations can influence the dynamics of nonlinear dynamical systems [20,31,32]. The effect of space on quasispecies dynamics has been investigated in several works. For example, limited diffusion was shown to provide mutant classes with a competitive advantage, also decreasing the critical mutation rate, m c (i.e., the mutation rate beyond which the mutant genomes outcompete the master sequence) at which the error threshold phase transition occurs [14,17]. The effect of space in the competition dynamics of two quasipecies has been also studied in the context of the survival of the flattest effect [20]. More recently, the effect of spatial competition on the diversity for structured quasispecies has been investigated by Aguirre et al. [19].
During the earlier stages of infection by plant viruses, the spreading of the viral population within a host starts from the initially infected cells to the nearest neighbors through the plasmodesmata, in a process known as cell-to-cell movement. Although some studies on different viruses infecting their hosts show that systemic movement can cause strong population bottlenecks and highly heterogeneous viral subpopulations in different organs [27,[33][34][35][36], the effects of the population bottlenecks during cell-to-cell movement have not been deeply studied. In this context, a key parameter in virus evolution is the number of virus genomes infecting a given cell, a parameter known as the multiplicity of infection (MOI). MOI is important as it determines processes such as the rate of genetic exchange among genomes, selection intensity on viral genes, epistatic interactions, and the evolution of multipartite viruses [37,38]. Several models of virus evolution have explored the role of MOI in host-pathogen interactions [38][39][40][41][42], but experimental estimations of MOI along infection are still scarce, and only a handful of studies have estimated MOI in bacteriophages [43][44][45] and in the larvae of insects [46]. In the recent years, plant virologists have turned their attention to this problem. In a seminal study, González-Jara et al. [47] have obtained estimates of MOI for the Tobacco mosaic virus (TMV) infecting Nicotiana benthamiana plants. They followed the process of infection and characterized the temporal variation of MOI for two TMV genotypes, finding that MOI decreased as infection progressed. These authors suggested that such a reduction in MOI could be explained by mechanisms limiting superinfection and/or by genotype competition. More recently, Gutiérrez et al. have provided a spatio-temporal monitoring of the cellular MOI for the Cauliflower mosaic virus (CaMV) [48]. This second study revealed the presence of dynamic changes of MOI throughout the infectious cycle in the plant, with a maximum MOI reached at intermediated times post infection.
Theoretical and computational quasispecies models have mainly considered RNA populations replicating exponentially, more generally, geometrically (hereafter GR). In the case of singlestranded RNA viruses, GR implies that both the genomic and antigenomic viral strands are used as templates for replication, and thus the accumulation of mutations is large because mutant genomes also serve as templates for replication. From this replication mode, the distribution of the number of mutants per infected cell follows the Luria-Delbrück distribution [49]. Experimental studies carried out with bacteriophage T2 supported such strategy [50]. Alternatively, viruses may replicate according to the stamping machine replication mode (SMR). Under this scenario, the initially infecting genomic strand is used for the production of one or few antigenomic ones, which are then used as templates for the generation of all the progeny of positive-sense strands that will then be encapsidated to continue the infection process. In this case, the number of mutant genomes per infected cell follows a Poisson distribution. Such a distribution of mutants was found for phage wX 174 [52]. Intermediate modes of replication, where some fraction of positive-sense strands may be also replicated, have been described for phage w6, whose distribution of mutants slightly differed from the Poisson distribution [53]. Recently, the mode of replication was inferred for Turnip mosaic virus (TuMV) [51], which was largely dominated by the SMR. The role of the replication mode in the accumulation of deleterious mutations as well as in the mutational robustness of well-mixed quasispecies populations was recently investigated by Sardanyés et al. [54]. These authors developed theoretical and computational models to characterize the effect of the replication mode on the accumulation of mutations for positive-sense, singlestranded RNA viruses under different fitness landscapes, paying especial attention to the epistatic fitness landscape, which has been confirmed in several examples of RNA vuses (see [55] for a review). In short, the main conclusion of this study was that the SMR was less sensitive to the effect of mutations and compatible with higher critical mutation rates.
The aim of the present work is to extend the results of [54] by incorporating the effect of space. Theoretical or computational works exploring the effects of the mode of replication on RNA virus dynamics are scarce, and previous attempts to tackle this question have not considered spatially-distributed viral populations [54,56]. The question we are addressing in this study is precisely what is the effect of space for viral quasispecies replicating under GR and SMR. To do so, we have developed stochastic cellular automata (CA) simulation models that consider replication, mutation and cell-to-cell infection in a two-dimensional environment. As we did in [54], here we also take a very general modeling approach simulating single-stranded RNA viral populations in silico, using digital quasispecies. Hence, the results of our study might serve as an approach to the dynamics of RNA viruses like arteriviruses, picornaviruses, flaviviruses or togaviruses, among others. Finally, our modeling approach also allows us to explore two different infection strategies that have been widely observed in experiments with viruses, namely, free superinfection (i.e., already infected cells are susceptible to additional infections) and superinfection exclusion (i.e., viruses of infected cells block the entrance of new viruses) [57][58][59]. Although our model is still a simplified picture of real viral infections, it represents a major step forward from previous models since it incorporates key features of real viral populations (e.g., different replication modes and different infection strategies).

Simulation Model
The effect of the replication mode (geometric replication, GR; and stamping machine replication, SMR) in the spatial dynamics of replication and infection of a quasispecies is studied by using stochastic cellular automata (CA) models. The CA works on a square state space C(i,j) [ Z 2 , with L|L cells (we use L~21) and zero-flux boundary conditions simulating the bounded system, for examle, of plant leaves ( Figure 1A). Following the approach of Leuthä usser [60,61], we use a bit-string description of the quasispecies population structure [20,54,62,63]. Hence, we do a mapping between RNA sequence, defined as a chain of nucleotides involving a four-letter alphabet V, and a binary sequence, given by: F : V~fU,C,G,Ag?
P~f 0,1g. In this case, the strings contain a sequence of purines or pyrimidines that only incorporate the linear information encoded in the genotype.
With such an approach we can analyze the spatial dynamics of RNA viruses using digital genomes and taking into account the mode of replication, mutation and cell-to-cell infection. Each lattice cell in C(i,j) has the potential to contain a maximum population of N~N max sequences (we use N max~5 00). That is, each cell has N max sites which can be occupied by newly produced strings or by strings entering from neighbor cells during cell-to-cell movement. Each one of these sequences, S k (i,j), with k~1,:::,N max , is a small digital genome of length l~1,:::,n (we use strings of length n~16 bits) i.e., S k (i,j)~(S k,1 ,S k,2 ,:::,S k,n ), representing a vertex S k (i,j), of a discrete, n-dimensional sequence space (hypercube, H n ): S k (i,j) [ H n , living in the (i,j) lattice cell ( Figure 1A). Hence the total population of strings once the lattice is full is N max |L 2~2 20,500 strings, being &3:36-fold the number of different strings of the entire sequence space. In order to model the effect of mutations in the fitness associated to the replication process as well as to infection, we consider that the digital genomes contain two different loci, each of them with a length of 8 bits. These two loci will be used to assign the fitness of each genome tied to the processes of replication and cell-to-cell infection. By doing so, we decouple replication from transmission: a genome may be able of replicating but not of transmiting itself, and viceversa. As a first approach and for the sake of simplicity, we make a direct mapping between the genotype and both infection and replication success of a given string, without explicitly modeling the production of replicase or movement proteins. Moreover, our model also obviates recombination and/or complementation between different genotypes. For both replication modes we study two different deleterious fitness landscapes with epistatic interactions ( Figure 1B). We specifically study the antagonistic fitness landscape, as being the more commonly described in RNA viruses [64][65][66][67], in which the deleterious effect of multiple mutations together is lower than expected from their individual effects. For the sake of completeness, we also implemented a synergistic landscape, in which an increasing number of mutations has a stronger deleterious effect than expected from the effects of each individual mutation. The fitness function associated to replication success, w k , for the k-th sequence is given by: where d w (S k (i,j),S m )~P n=2 l~1 jS k,l (i,j){S m,l j. The fitness associated to infection success, c k , for the k-th sequence is given by: is the Hamming distance between each locus of the k-th sequence and the corresponding locus of the master sequence (i.e., the all-ones string labeled S m~1 :::1). Note that d computes the number of mutations of the k-th string (i.e., number of bits 0). The parameter j denotes the sign and the strength of epistasis (see [54,55]). For the antagonistic landscape we use jv1, while for the synergistic one we will use jw1. The model does not explicitly include beneficial mutations, but, according to the studied landscapes, backward mutations will always involve a fitness increase.
The CA works as follows: we first choose a random cell of the lattice whenever the lattice is not full. If the chosen cell is empty, we continue with the same process. However, if the chosen cell is not empty, then we consider that a generation has taken place. Then, for such a cell and generation time, we first apply N max times the rule of intracellular replication in order to ensure that, on average, all the strings inside the chosen cell are updated once per generation. After the N max rounds of replication we apply the rule of cell-to-cell infection. These two rules are applied until the whole lattice is full of strings. Next, we describe the state-transition rules of the CA.
1. Intracellular replication. We choose at random a replicating string of the cell which is copied with a probability proportional to its fitness w k into another randomly chosen empty site. Although we do not explicitly consider the polarity of the strings, we can differentiate between GR and SMR. On the one hand, GR is implemented by considering that all the strings will replicate proportionally to their replicative fitness, w k , giving rise to more replicating genomes. On the other hand, SMR is implemented as follows: the strings belonging to the initial conditions or to the newly infecting strings entering into the cell (assumed to be positive-sense strands) will replicate only once, giving place to the template genomes that will be used for further replication. Such initial strings will not continue replication in the following generations. Instead, their offspring (acting as antigenomic templates) will continue replicating, giving rise to non-replicating genomes. By doing this, the progeny of strands will only be generated from the templates synthesized from the first infecting genomes. Replication mechanism presents per-bit mutation probability, m b [ ½0,1, per replication cycle. 2. Cell-to-cell infection. We will assume that the strings within a cell can infect neighboring cells when N~N max . When this condition is fulfilled, we choose a string at random inside that cell, which moves with probability c k [ ½0,1 to a randomly chosen empty site of a neighboring cell (we use a von Newmann neighborhood i.e., 4 nearest neighbours). We will also consider that all the strings infecting neighboring cells will become replicators independently if they were previously replicators or not. However, in order to simulate a more realistic scenario for positive-sense RNA viruses (where the genomic strands are encapsidated to infect neighboring cells), we will assume that for GR a given string will infect the neighboring cell with probability c k =2, because, on average, GR is producing the same number of genomic and antigenomic strands (i.e., *50% of each type of string). On the contrary, for SMR, we will take into account that all the strings can infect neighboring cells with probability c k because the majority of the offspring are positive-sense strands. For both replication modes and fitness landscapes, we will study two possible mechanisms of infection (see Figure 1B): (i) superinfection exclusion (SE) and (ii) free superinfection (FS). In case (i) we consider that when a string of a given cell, e.g., S k (i,j), infects an empty neighbor or a neighbor that has not reached yet N max , i. e., S empty (s,p) with s[fiz1,i{1g, p~j; and s~i, p[fjz1,j{1g, the cell S(i,j) cannot be infected again by strings from neighboring cells. In (ii) we will consider that, independently of their current infection status, a given cell can be infected by a new string from a neighboring cell.
The algorithm starts with the central cell of the lattice inoculated with an initial population of N 0~1 master sequences (assuming that are positive-sense strands). As previously mentioned, for GR these strings will always replicate producing the offspring that will also replicate. For SMR, they will replicate once, giving rise to templates that will be the responsible of producing the entire progeny of genomes that will not further contribute to the generation of more strands. We note that in our computational model there is only one real free parameter given by mutation rate.

Results
The lowest critical mutation rate, m c b , is found for GR and antagonistic epistasis We first investigate the value of critical mutation rate per bit, m c b , at which the population experiences the transition to the error threshold. To do so, we study how the concentration of master sequences changes at increasing mutation rates [ Figure 2(a)]. The per-bit critical mutation rate, m c b (i.e., the mutation rate beyond which the population of strings is dominated by mutants) is considered as the lowest value of mutation rate at which the concentration of master sequences is lower than 10 {4 . Although our model does not incorporate degradation of strings, the error threshold we are characterizing corresponds to extremely low population numbers of master sequences due to mutation processes. The numerical value attributed to the critical mutation rate involves an upper bound with a small population of master sequences (i. e., ƒ22 over 220,500). The increase in mutation rate involves a decrease of the master sequences for all the studied combinations. However, the magnitude of such a decrease strongly depends on the replication mode ( Figure 2 is, on average, 269:70% larger for GR, this difference drops to 49:08% for SMR. Although previous results for well-mixed populations suggest that for SMR the population of master sequences might be less sensitive to mutation [54], the inclusion of spatial correlations allows the stable existence of master sequences at very high mutation rates provided the combination of GR and a synergistic fitness landscape. This phenomenon may be due to the effect of purifying selection, which under GR and the synergistic fitness landscape is more efficient because the production of strings with a very low or no fitness increases, and the competition for space is not so strong as for the SMR mode, where the master sequences are competing with higher fitness sequences and then are suffering the error threshold at lower values of mutation.
Next, we evaluate the effect that the infection strategy (FS vs SE) may have on m c b . First, we find that, on average, SE is compatible with m c b values that are 12:58% larger than if FS is allowed, being this difference highly significant (Figure 2 associated to SE is only 0:85% larger than for FS. By contrast, for SMR m c b is 28:30% larger if superinfection is not allowed in the system. This differential effect can be rationalized as follows: if infection is limited and no strings can enter into a cell once the quasispecies has started infection to a neighboring cell, the population of master strings is more robust and can persist under larger mutation rates. On the other hand, if superinfection takes place, the critical mutation rate diminishes and the quasispecies enters into error catastrophe at lower values of m c b . Third, in a lesser extent the effect of the infection strategy also depends on the topography of the fitness landscape ( Figure 2

Dynamics and spatial distribution of digital quasispecies
The space-time dynamics of the digital quasispecies has also been investigated for all combinations of fitness landscapes, replication modes and infection strategies. For illustrative purposes, Figure 3 shows the results obtained for the antagonistic fitness landscape considering SE. The other cases are shown as supplementary material ( Figure S1 antagonistic landscape with FS, Figure S2 synergistic landscape with SE, and Figure S3 synergistic landscape with FS). As expected, for both replication modes, an increase in mutation rates involves a decrease of the concentration of the master genomes in all the studied cases once all lattice cells are totally full. Such decrease is much more accentuated for the GR mode. Generically, the master sequences can persist with SMR. However, for the GR mode, the master sequences maintain very low numbers even for small values of mutation rate (i.e., m b~0 :0025). The most important differences between fitness landscapes correspond to the spatial distribution of the fitness of the quasispecies in the lattice. For both fitness landscapes we show that the mean fitness of both loci per cell is lower for the GR mode, while for the SMR mode the mean fitness per cell displays darker gray colours, being nearer to the maximum value (~1, in black). Due to the nature of the fitness landscape we see that for the synergistic landscape, the mean fitness drastically reduces as depicted for the clearer spatial patterns shown in Figures S2 and  S3, being much more pronounced for the GR mode. For the antagonistic landscape (with j~0:6) the value of minimum fitness per locus is w i &0,5647.
Increasing mutation for GR rapidly stabilizes digital quasispecies loci at linkage equilibrium Next we evaluate whether mutations at the two loci (i.e., replication and movement) associate randomly in the resulting viral population or linkage disequilibrium is created. Three forces may create linkage disequilibrium, namely mutation, selection and the sampling events that take place to initiate new cell infections. We are not intended to disentangling the contribution of these three mechanisms to the disequilibrium but just to determine wheter it may exist. To this end, we compute the linkage disequilibrium coefficient, D, between two alleles of the two loci [68]. The alleles are differentiated for each locus, given by master (i.e., all-ones locus, indicated with 1) and mutant (i.e., a locus with one or more zeros, indicated with 0) loci. Hence, D is computed from D~x 11 {p 11 q 11 , with p 11~x11 zx 10 and q 11~x11 zx 01 . Here, x 11 is the relative frequency of master sequences in the whole lattice, once the quasispecies has infected the whole space. The value of x 10 corresponds to the relative frequency of strings in the whole lattice with master replication locus and mutant infection locus and x 01 is the relative frequency of strings once the lattice is filled with mutant replication locus and master infection locus. We compute the mean linkage disequilibrium, D D, over 50 independent runs (after the lattice was completely full of strings), at increasing mutation rates for both replication modes also considering the antagonistic and synergistic landscapes for both infection strategies. The results, displayed in Figure 4 considering SE, clearly show that the main force determining the linkage disequilibrium is the mode of replication. For both replication modes, D D first increases reaching a maximum value and then declines at increasing mutation. However, in the entire range of mutation rates analyzed (i.e., 0ƒm b ƒ0:04), D D remains significantly larger than zero only for the SMR mode and regardless of whether superinfection was allowed or not. By contrast, for relatively low values of mutation rate, GR quickly reaches random association of alleles at both loci. This effect of the mode of replication can be explained by the fact that SMR produces genomes with a lower number of mutations than GR. By using already mutated templates, GR generates molecules carrying multiple mutations and thus, breaking any association between alleles at both loci. Figure 4 also shows that the topography of the fitness landscape has a minor effect on D D, specially for the case of GR. For the case of SMR the effect depends on whether superinfection is allowed in the system. If SE is the norm, then the maximum linkage disequilibrium is reached at mutation rates larger for synergistic landscapes (m b &0:0075) than for antagonistic ones (m b &0:0025). By contrast, if superinfection is allowed, the maximum is reached at the same mutation rate regardless the landscape topography (m b &0:0075). On average, the maximum value of D D is 1.5-fold larger if SE exists than if superinfection occurs in a free manner.

The multiplicity of infection generically fluctuates and increases along time as infection progresses
Another interesting phenomenon that can be explored with our computational models is how MOI changes among replication modes for both fitness landscapes, as well as the spatial patterns of infections inside each cell of the lattice. As we mentioned in the Introduction, MOI is the number of viral particles infecting a host cell. Figure 5 shows the distribution of MOIs per cell for simulations with and without SE and for the case of antagonistic fitness landscapes (the results for the synergistic fitness landscape are shown in Figures S5 and S6 and are not discussed because are qualitatively similar to those shown in Figure 5). The frequency distribution of MOIs is computed once the whole lattice is completely full of strings. That is, after whole infection, we compute how many lattice cells are infected by n~1,2,3,:::, strings. In both cases the spatial distribution of infections is disordered, and no clear spatial patterns are found ( Figure S4). Overall, for both fitness landscapes, the average number MOI per cell is higher for the SMR than for the GR.
To analyze in a more quantitative way the data shown in Figure 5, we fit MOI values to a generalized linear model. MOI is assumed to be Poisson distributed and the mode of replication and whether superinfection is allowed or not are treated as main factors whereas the mutation rate (m b ) is incorporated in the model as a covariable. All interactions between factors and the covariable are included in the model. Overall, the mode of infection has a highly significant effect on MOI (Pv0:0001), being it 67:20% larger for SMR than for GR. As expected, allowing for FS makes MOI to increase up to 68:42% when compared to the case of SE (Pv0:0001). However, the effect of the mode of replication is not independent on the superinfection status. These two main factors interact in a very significant way (Pv0:0001), being the effect of the replication mode on MOI larger when superinfection is free (95:44%) than when it is excluded (43:04%). The covariable, m b , has a significant negative effect on MOI (Pv0:0001). Increasing it in the range 0ƒm b ƒ0:07 implies a decline in MOI of 19:70%. However, the magnitude of this decline depends on whether superinfection was free or limited (test of interaction term, P~0:0010), being smaller in the former situation (14:26%) than in the latter (24:79%). The effects of mutation rate and replication mode are independent (P~0:0689). Similarly, the interaction between the two factors is not affected by m b (P~0:2830).
Next, we explore the temporal dynamics of MOI. To do so, we compute, at each time generation, the mean number of strings that have entered into the cells of the lattice (averaged over all infected cells). The results are shown in Figure 6 for the antagonistic fitness landscape and m b~0 :0025 (no qualitative difference exists for the synergistic landscape). We specifically show two plots considering SE (left panel) and FS (right panel). For each case we display the time evolution of the MOI(t) with three trajectories corresponding to three independent replicas. The quantity MOI(t) is a measure of how the lattice is filled along time, so it is a kind of cumulative measure which fluctuates because the data is normalized at each time generation by the total number of infected cells. The results of Figure 6 show that, independently of the replication mode and of the infection strategy, MOI fluctuates, but significantly increases as time progresses. Indeed, linear regression analyses confirm that the slope is significant in all four cases (Pv0:0001 in all cases). An ANCOVA using mode of replication and superinfection status as factors and time as covariable shows that both factors have a significant effect on the average MOI reached in the lattice (in both cases Pv0:0001). Interestingly, both factors show a significant interaction with the covariable (in both cases Pv 0:0001), suggesting that the rate at which MOI increases with time depends on them. For instance, when superinfection is excluded, the slope of the regression line obtained for the SMR is 222:11% larger than when the virus replicative strategy is GR. This difference in the rates of MOI change among replicative strategies is even larger for the case of FS (750:60%), a difference supported by a significant three-ways interaction term in the ANCOVA (Pv0:0001). Therefore, we conclude that MOI increases with time but that the increase is faster if SMR is the replication strategy followed by the virus and if no SE mechanisms are at play. Qualitatively similar results have been found at increasing mutation rates and for the two fitness landscapes analyzed (data not shown).
The mean values of genome entries per cell once the lattice is completely full are represented, for each studied case, in Table 1 and 2. These results indicate that larger MOIs are found for the SMR mode under the antagonistic fitness landscape with superinfection. Generically, MOI decreases for the synergistic fitness landscape, for both SMR and GR modes. This phenomenon occurs when considering both SE and FS. We note that for both infection strategies and fitness landscapes, the values of MOI are always higher for the SMR mode. This result may reflect the implicit consideration of the sense of the strings (recall the assumption that for the GR the infection probability was halved because approximately the 50% of the progeny might correspond to positive-sense strands). However, this is a consequence of the nature of the infections, since a single-stranded RNA virus might be able to encapsidate more genomic strands when replicating under SMR.

SMR together with antagonistic epistasis involves the fastest genomes colonization time
Finally, as a measure of the performance of the different strategies in spreading the infection, we study the mean infection time, T T i , computed as the time it takes to complete the infection of all cells in the lattice (i.e., all cells contain N max strings), for the two replicative strategies and fitness landscapes for increasing mutation rates but considering only the case of superinfection inhibition ( Figure 7). As we did before, the data are fitted to an ANCOVA model in which landscape topography and replication mode are treated as fixed factors and m b as covariable. A first result is that T T i significantly increases with m b (test of covariable, Pv0:0001). The mean infection times for the quasispecies replicating via the SMR are, for the whole range of mutation rates analyzed, systematically lower (68:03%, on average) than for quasispecies replicating via GR mode (test of replication mode main effect, Pv0:0001).
Hence, the spread of the strands is faster for the SMR mode, and such result is consistent independently of the fitness landscape assumed (test of the landscape topography main effect, P~0:2285; test of the interaction between landscape topography and replication mode, P~0:1074). However, if one compares the time of infection between the two fitness landscapes for a given mode of replication we find that when mutations interact in a synergistic manner in determining fitness, the time required to complete an infection increases. This phenomenon is observed for both replication modes, although it is much more accentuated for the GR (test of the three-ways interaction term, Pv0:0001). For  this case, and for large mutation rates (e.g., m b &0:04), the mean time of infection for the synergistic landscape is *1:98-fold the time of infection for the antagonistic one. This may occur because for the synergistic landscape, mutations have a stronger deleterious effect and the quasispecies is producing less efficient mutants, who are not able to replicate and infect optimally, then needing a Table 1. Mean values (+1 SEM) of the multiplicity of infection for the model with SE, computed over 50 independent replicas once the lattice is completely infected by the quasispecies at increasing mutation rates, m b , for both replication modes and both fitness landscapes (antagonistic with j~0:6; synergistic with j~1:4).  Table 2. Mean values (+1 SEM) of the multiplicity of infection for the simulations considering FS, computed over 50 independent replicas once the lattice is completely infected by the quasispecies at increasing mutation rates, m b , for both replication modes and both fitness landscapes (antagonistic with j~0:6; synergistic with j~1:4).
Stamping machine replication (SMR) Geometric replication (GR) longer time to complete the infection of the whole lattice. For the combination of SMR and synergistic fitness landscape, the quasispecies also undergoes longer infection times as mutation rates grow, but such times remain always below the time needed for both landscapes under the GR mode. Under the SMR, the mean time of infection with m b~0 :04 for the synergistic landscape is *1:25-fold the one for the antagonistic landscape also under SMR. Finally, we also want to notice that for small mutation rates the values of T T i are quite similar for the SMR mode, but for the GR mode they rapidly diverge as m b becomes larger (see. Figure 7 for details). Qualitatively similar results are found for the simulations run with FS (results not shown).

Discussion
Different replication modes have been suggested for different viruses. Depending on whether genomes are replicated according to a stamping machine (SMR) or geometrically (GR), one may observe different compositions in the mutant spectrum of the quasispecies. For the SMR model, the initial antigenomic templates are the ones used for further replication, as observed for the phage wX 174 [52], where the number of mutants per infected cell followed a Poisson distribution. For the GR mode, such a distribution of mutants deviates from a Poisson, and follows a more complex distribution, termed Luria-Delbrück distribution. Experimental analysis showed that the bacteriophage T2 did not fit well a Poisson distribution, suggesting a GR strategy [50]. Intermediate modes of replication may also exist, as illustrated by experiments with phage w6 showing a distribution of mutants slightly different from the Poisson expectation [53], and with TuMV showing assymmetric accumulation of strands of both polarities [51]. Although the replication strategy for RNA viruses might have a deep impact in the accumulation of mutations and therefore in the fitness of the sequences, few works have explored the effects of different replication strategies in the population dynamics of viral quasispecies [54,56]. Not to say that no work at all has been published exploring the interplay between the mode of replication and spatial correlations nor the existence of mecha-nisms of controlling superinfection. In this study we intend to cover the gap existing in biologically unrealistic models of virus replication and spread proposing new models that incorporate some of the most basic features of viral genomes. We analyzed the dynamics of replication and infection of quasispecies on a twodimensional space by means of stochastic cellular automata models, especially focusing on the effect of different replication modes, topography of the fitness landscapes and existence of mechanisms inhibiting superinfection. Among the theoretical and computational approaches to study spatially-extended biological systems (e.g., metapopulation models, partial differential equations, etc.) we have chosen to simulate single-stranded RNA quasispecies by means of digital genomes using a cellular automaton approach. Such a strategy allows us to use spatial individual-based modeling taking into account the heterogeneous population structure of the quasispecies together with stochasticity. The digital genomes were constituted by two independent loci, one determining replication and the other determining the efficiency of cell-to-cell movement. Hence, our results may be useful to understand the dynamics and evolution of widely different viruses as far as they fulfill these basic assumptions. To the extend of our knowledge, our study is the first one simulating the spatio-temporal dynamics of single-stranded RNA viruses under different modes of replication. Moreover we modeled the fitness effects of mutations on each of these loci assuming two different epistatic fitness landscapes, one antagonistic and one synergistic. Supporting these choices, epistasis has been widely found in real RNA viruses [55,69]. Together with the incorporation of two different replication modes and fitness landscapes in our simulations, we also investigated two possible infection strategies given by limitation of infection and superinfection, strategies known to occur in real viral population.
Non-spatial computational models with digital quasispecies showed that when replication proceeds via SMR, the population of master sequences is less sensitive to mutation and the critical mutation rates involved in the error threshold were always lower under GR, independently of the fitness landscape assumed [20]. In agreement with the previous work, our results show that the extinction of the master sequences occurs at a larger mutation rate with GR for the synergistic fitness landscape. However, in [20], the critical mutation rate of the quasispecies replicating under GR was always lower than the critical values under the SMR. In the simulations developed in this work we show that the largest critical mutation rate corresponds to the GR mode under synergistic interactions between mutations. This may be due to an enhanced synergy between space and purifying selection, where mutants with extremely deleterious mutations in both replication and infection loci could not spread over the lattice, favoring the selection of master sequences because reduced competition with other mutant sequences during replication and infection.
The analyses of two different infection strategies, given by superinfection exclusion (SE) and free superinfection (FS), revealed that when SE is considered, the critical mutation rates for SMR become larger, and thus the robustness of the quasispecies increased because the master sequence was able to persist for larger mutation rates. This result was in agreement with several works suggesting that when a cell is coinfected by different viral genomes, the fitness of individual genotypes may decrease in comparison with their fitness in a single infection due to competition processes [41,70]. This phenomenon, however, seemed not to be important for the GR mode, probably because the accumulation of mutations was so large that the quasispecies in a given cell was dominated by the pool of mutants, and thus the entry of new mutants was irrelevant. Nevertheless, the entrance of new mutants in a cell for the SMR could drastically change the mutant spectrum of the progeny, especially if the infecting mutants carry many mutations.
Several examples of viruses that have evolved mechanisms to avoid superinfection have been reported and studied. For example, the wild-type sie z P22 prophage is able to prevent growth of superinfecting phage P22 as well as of other phages like the L, the MG178 or the MG40, by means of exclusion mechanisms [71]. Other examples of viruses with mechanisms to avoid coinfection are the phage W6 [43] as well as the Vesicular stomatitis virus (VSV) [57]. Whitaker-Dowling et al. [57], showed that the presence of mechanisms to limit superinfection is virusdependent. They showed that infection of BHK 21 cells with either influenza viruses, Encephalomyocarditis virus or Newcastle disease virus did not inhibit superinfection by VSV, whereas cells initially infected with VSV were not susceptible for VSV superinfection. Our work suggests that a possible answer to the complex question of why some viruses limit superinfection but some others do not present such a property could rely on the mode of replication of each virus type. As previously mentioned, an important result we obtained in our simulations was that when replication proceeds via SMR, the limitation of superinfection largely affects the survival and maintenance of the master sequence, as a difference from the GR, where no significant differences were found between the two studied infection strategies.
The multiplicity of infection (MOI) (i.e., number of viral genomes infecting a host cell) is a key parameter in virus evolution because it can determine selection intensity on viral genomes, exchange among genomes or epistatic interactions. MOI has been mainly studied in different DNA and RNA bacteriophages [43][44][45]. However, very few experimental studies reporting estimates of MOI are found for virus infecting eukaryotic hosts. For example, MOI was studied for Autographa californica nuclear polyhedrosis virus during the infection of larvae of the lepidoptera Trichoplusia ni [46]. Only very recently, González-Jara et al. [47] and Gutiérrez et al. [48] have taken the task of studying the evolution of MOI during infection of plants by RNA viruses. In the former work, the authors carried out experiments of host colonization by two genotypes of the TMV infecting N. benthamiana plants. They found that MOI decreased during infection, and suggested that two nonexclusive processes could cause such a decrease: by mechanisms limiting superinfection and/or by genotype competition. Interestingly, our analyses of MOI dynamics showed that MOI fluctuated and that, independently of the replication mode, the fitness landscape and despite the existence of SE mechanisms, MOI increased in time. This finding suggests that the results found by González-Jara et al. [47] could not be explained by genotype competition or by SE. Hence, other mechanisms lowering MOI might be operating during infection. In the work of Gutiérrez et al. [48], the spatio-temporal dynamics of two competing variants of CaMV was monitored during the infection of turnip plants. They reported great changes of MOI at different infection phases during plant development [48]. Actually, our results, which might be interpreted as a quasispecies infecting a single tissue, reflected these fluctuations in MOI during the process of infection. A direct consequence of high MOI is the recombination and complementation between genotypes [72][73][74]. For the sake of simplicity, our models do not take into consideration both recombination and complementation processes between genotypes. These important phenomena, together with differential replication modes, should be considered in future research. Moreover, the consideration of a full model considering intracellular amplification and both cell-tocell and systemic movement under the previous scenarios also remains a theoretical and a computational challenge.
Finally, a take-home message of our work is that important differences between non-spatial and spatially-structured models exist for quasispecies replicating under different replication modes. For instance, previous results indicated that the critical mutation rate of quasispecies was lower if replication was GR than if it was SMR independently of the unferlying fitness landscape [54]. However, our simulations have shown that by considering spatial correlations, the outcome would be the opposite for the synergistic landscape: GR would result in a more robust replication strategy. Moreover, we have also shown that mechanisms of superinfection exclusion during cell-to-cell movement might play an important role in virus robustness to mutations. Our findings also suggest that other mechanisms beyond limiting superinfection and/or genotype competition should be considered to explain the decrease in MOI reported in [47].