Abstract
Background
Viral quasispecies can be regarded as a swarm of genetically related mutants. A common approach employed to describe viral quasispecies is by means of the quasispecies equation (QE). However, a main criticism of QE is its lack of frequencydependent selection. This can be overcome by an alternative formulation for the evolutionary dynamics: the replicatormutator equation (RME). In turn, a problem with the RME is how to quantify the interaction coefficients between viral variants. Here, this is addressed by adopting an ecological perspective and resorting to the niche theory of competing communities, which assumes that the utilization of resources primarily determines ecological segregation between competing individuals (the different viral variants that constitute the quasispecies). This provides a theoretical framework to estimate quantitatively the fitness landscape.
Results
Using this novel combination of RME plus the ecological concept of niche overlapping for describing a quasispecies we explore the population distributions of viral variants that emerge, as well as the corresponding dynamics. We observe that the population distribution requires very long transients both to A) reach equilibrium and B) to show a clear dominating master sequence. Based on different independent and recent experimental evidence, we find that when some cooperation or facilitation between variants is included in appropriate doses we can solve both A) and B). We show that a useful quantity to calibrate the degree of cooperation is the Shannon entropy.
Conclusions
In order to get a typical quasispecies profile, at least within the considered mathematical approach, it seems that pure competition is not enough. Some dose of cooperation among viral variants is needed. This has several biological implications that might contribute to shed light on the mechanisms operating in quasispecies dynamics and to understand the quasispecies as a whole entity.
Background
The concept of quasispecies [1] refers to the equilibrium spectrum of closely related mutants, dominated by a master sequence, generated by a mutationselection process. It has become an adequate descriptor of RNA viruses at the population level [2,3] and provides natural links between population biology and virology.
Competitive exclusion and displacement when two viral populations with nearly equal starting fitness compete has been observed [4]. Thus, a challenging puzzle is how so many (sometimes very similar) variants can coexist in nature. In the case of ordinary ecosystems niche differentiation is obviously an important aspect to promote biodiversity [5]. However, the often striking similarity in coexisting variants suggests that other mechanisms must be involved.
An infected individual can harbor several genetically related variants of the same species (assuming that this concept is valid for RNA viruses), and thus the host can be regarded as an ecosystem in which distinct viruses interact.
These interactions are in general quite complex and may be direct or indirectfor example, when the host's immune system responds against one particular viral variant, affecting the whole fitness landscape. There is a process called complementation, in which one virus provides in trans a useful product that cannot be made by another variant. On the other hand, it was proposed that in trans acting defective gene products, expressed by mutagenized viruses can interfere with replication and prevent replication of pathogenic viruses [6]. Thus, under such circumstances viral phenotype may not necessarily reflects viral genotype.
Obviously, if the viruses provide each other with useful resources the interaction is of mutual benefit. Indeed this is what it was found in the evolution of competitive interactions among RNA phage viruses [7]. When two mutants coinfect a cell the common resource pool allows the viruses to use each other's protein products. The coinfection rescues the mutants, allowing them to reproduce when they would be otherwise unable to do so [8,9].
Furthermore, it was recently found that for polio viruses the diversity of the quasispecies population, rather than selection of individual adaptative mutations, may determine pathogenesis through cooperative interactions [10,11].
One of the common approaches employed to describe viral quasispecies is by means of the quasispecies equation [12] (QE) which is limited by the fact that it lacks frequencydependent selection. In other words, the fitness of a particular phenotype is set to a constant value independently of the other competing 'players'. This can be fixed by using the replicatormutator equation [13] (RME), which represents a general formulation for the evolutionary dynamics.
A far from trivial problem with the RME is how to get the interaction coefficients among variants. Here we adopt an ecological perspective and resort to the niche theory [14] of competing communities, which assumes that exploitation (or utilization) of 'resources' primarily determines ecological segregation. Thus, we regard the swarm of mutants that constitute a quasispecies as a set of interacting RNA molecules distributed along a hypothetical niche axis [15]. This allows, in a simple way, to estimate the size of the interaction coefficients of the RME. Basically, the degree of overlap between variants is what ultimately determines the intensity of these coefficients (see Method section).
The other novel ingredient, besides the use of niche theory, is to consider the effect of some measure of cooperation or facilitation among viral variants. We start by showing that in general niche theory does not lead by itself to a 'typical' quasispecies profile (see Results section). Then we show that this problem can be solved if some dose of cooperative interactions between variants is taken into account. Therefore, in order to get a 'typical' quasispecies profile both competition and cooperation or facilitation between viral variants is needed.
Results
In the case of pure competition, i.e. all interaction coefficients α_{ij }negative (see Method section), we obtain distributions like the ones shown in Figure 1, which represent in each row the evolution for a different value of σ:σ = 0.2, σ = 0.15 and σ = 0.1 in rows 1, 2 and 3, respectively. The left panels are histograms of the fractions x_{i }of phenotype i vs. their corresponding positions on the niche axis μ_{i }for N_{g }= 500 (see Method section for a detailed description of parameters). In the three cases, a pattern of lumps 2, 3 and 4 respectivelyseparated by gaps along the niche axis emerges. This multilump population profile is in complete agreement with the one recently found by Scheffer and van Nes [16] for LV competition models, dubbed by them as selforganized similarity (SOS). This is not surprising since the LV competition model is mathematically equivalent to a replicator dynamics. In fact, up to now, the main difference between both treatments is in the way mutations are implemented, but this seems not to affect qualitatively the results. Furthermore, the observed inverse dependence in the number of lumps with σ was recently explained [17]. The right panels in Figure 1 show the temporal evolution of each fraction during the N_{g }= 500 generations.
Figure 1. Population fractions for the case of pure competition. Relative abundance of each variant after 500 generations for, respectively, σ = 0.2, σ = 0.15 and σ = 0.1. Corresponding temporal evolution of the fractions of each variant is shown (right panels).
Notice that this multilump pattern resembles several coinfecting quasispecies (one per lump) rather than a population profile of a single quasispecies. Moreover, after 500 generations, a) there are many phenotypes that survived, and b) the viral ecosystem has not reached equilibrium yet (the population fractions of several variants are still considerably changing with time). So the pure competitive niche model does not seem to work properly.
Let us now explore what happens when some degree of cooperation or facilitation between variants is included in this model by turning a random fraction f_{C }of the interactions from competitive to cooperative (see Table 1 in Method section). Figure 2 is completely analogous to Figure 1 but when 1% (f_{C }= 0.01) of all α _{ij }are positive, illustrating the dramatic changes introduced by taking into account even a very small amount of cooperation. First of all, notice that the multilump pattern disappeared and only some variants remain. This is particularly clear for σ = 0.2 and σ = 0.15 with a master sequence that represents around 50 and 40% of the population, respectively. In addition, the population profile does not change very much after a few hundred generations, in agreement with the concept of a quasispecies. In contrast, if σ is too small (Figure 2 row 3) the equilibrium can not be easily attained, and therefore these results offer a criterion to estimate a lower limit for this parameter.
Figure 2. Population fractions for the case f_{C }= 0.01. Analogous to Figure 1 for the case of f_{C }= 0.01. Relative abundance of each variant after 500 generations for σ = 0.2, σ = 0.15 and σ = 0.1. Corresponding temporal evolution of the fractions of each variant is also shown.
Table 1. Parameters involved in the model.
It is also interesting to evaluate the dependency on f_{C}. Consider, for example, increasing it one order of magnitude, from 0.01 to 0.1. In this case there are no qualitative changes. However, if we increase the degree of cooperation even further to f_{C }= 0.5 (50% of cooperative interactions) and the other two parameters are set to σ = 0.15 and P_{m }= 0.1, the whole picture changes again. Figure 3 clearly shows that after a few hundred generations, much more variants survive, the master sequence is below 14% of the total population and still the populations are changing drastically. That is, the profiles for large values of f_{C }become widely distributed, with a master sequence that seems underrepresented. In order to further address this question we have performed, for different values of f_{C }(and always for σ = 0.15), 1000 simulations corresponding to 1000 distinct realizations. Each realization corresponds to a different random distribution of variant populations, and consequently different interaction matrices α_{ij}. Moreover, different growth rates r_{i }were assigned. We have computed two quantities by averaging over these 1000 realizations. First, , the number of variants, among the N = 100 considered, which after N_{g }= 500 generations represent at least 1/2 of the master sequence (of course 1/2 is an arbitrary threshold). Second, another measure of variability given by the normalized Shannon entropy (S) calculated as S_{n }=  [Σ_{i }x_{i }ln(x_{i})]/ln(N)]. The possible values of S_{n }range from zero (when there is just one variant, the master) to one (when all variants differ from one another and are equally represented). For normalized Shannon entropy applied to viral quasispecies see, for instance, reference [18].
Figure 3. Population fractions when half of the interactions are cooperative. Relative abundance of each variant after 500 generations (Left) and temporal evolution of the fractions of each variant (Right).
We found that these average percentages vary from 1.6% for f_{C }= 0.01 to 5.3% for f_{C }= 0.5 (see Table 2). Since we don't have experimental data of quasispecies profiles with such detail, it is difficult to determine a critical value for f_{C}, f_{C}* However, the scant experimental data seem to favor small values of
Table 2. Role of f_{C }in the model.
In the case of S_{n }we have found that the average varies from 0.474 for f_{C }= 0.01 to 0.471 for f_{C }= 0.1 and to 0.592 for f_{C }= 0.5, while in the case without cooperation, = 0.708. Therefore, the case in which f_{C }= 0.5 appears to be closer to the purely competitive situation, indicating again too much variability. So these results also favor smaller values of f_{C}.
Discussion
In this study we have modeled quasispecies by the RME plus the ecological concept of niche overlapping to quantify the intensity of interactions within a quasispecies. In other words, the swarm of viral variants that constitute a quasispecies is regarded as a set of interacting viral variants distributed along a niche axis.
The population distributions that emerge when only competition between variants is taken into account are different from what one would expect for a quasispecies profile: they are very unstable and resemble more to several coexisting quasispecies than to just a single one. Both problems can be solved if some dose of cooperation or facilitation between variants is included. It is worth stressing that this is just one possible solution to get a population profile that resembles that of a quasispecies from a general evolution model. The existence of cooperation among viruses is supported by recent experimental data, both in the case of coinfection [8,9] and among the variants that constitute a quasispecies [10]. In fact, the composition of a quasispecies mutant spectrum may indeed determine the viral population behavior through complementation or interference, i.e., positive and negative interactions, respectively [10]. It was recently demonstrated that interfering potential of replicationcompetent mutants may eventually modulate viral infection and population behavior of highly variable RNA viruses [19].
In addition, we have found two, at first sight, surprising results. First, the cooperation among variants seems to undermine, rather than enhance, biodiversity: much more variants within a quasispecies survive in the case of pure competition between them than when some facilitation occurs. Second, the effect of changing a small percentage of interactions from competitive to cooperative is more drastic than changing a larger fraction. This becomes evident from the comparison between results for f_{C }= 0.01 with those for f_{C }= 0.5 cooperative interactions. These two outcomes can be understood, from a mathematical point of view, as follows. When a small number of positive interactions, completely at random, are incorporated, the fitness of the few 'lucky variants', which receive help from others or become involved in cooperative cycles, soars. This has a strong destabilizing influence and, in consequence, many variants or RNA molecules disappear. So the whole effect is that diversity decreases. As the percentage of cooperative interactions is increased, the balance is gradually restored and so is diversity. This is also consistent with and S_{n }analysis (see Table 2); the higher the value of the proportion of cooperative interactions f_{C}, the higher the value of these two parameters. These data offer a guide to estimate the order of the percentage of cooperation that is needed to reproduce a quasispecies profile. In fact, Figure 2 fits better with the idea one has of a quasispecies than Figure 3.
One might speculate that the profile of randomly assigned growth rates r_{i }can have an important effect in the above results. However, we checked that with a uniform profile i.e. r_{i }= r = 1 for all variants results are qualitative similar. Thus in general the major role in determining the fitness of each variant comes mostly from the interaction with other variants.
The above discussed results have several direct biological implications that may contribute to shed light over the mechanisms operating in population dynamics, particularly in quasispecies mutant spectra behavior. As it is known, in an infected organism during normal infection, multiple viral genomes compete for resource needed to complete life cycle. Due to the high mutation rate of RNA viruses [20], defective mutants or genome defectors coding for nonfunctional proteins able to use these resources arise, likely interfering with efficient replication of viable viruses. In fact, what lethal defection model proposes is that when a slight increase in the level of mutant defectors occurs (e.g., by mutagenic activity), viral extinction may be achieved [21].
Conclusions
In this work we introduce two main novelties to mathematically describe the dynamics of RNA viral populations. First, we resort to niche theory in order to quantitatively estimate the interactions between viral variants constituting a quasispecies. Second, we include cooperative interactions between these viral variants and analyze their effects from an evolutionary point of view.
The model we propose is constituted then by two main parameters: the overlapping parameter σ, which controls the intensity of the interactions, and the percentage of cooperative or positive interactions f_{C}, which determines the sign balance of these interactions (see Method section).
Concerning potential applications, this model could be eventually applied to better understand complexity of viral population behavior in the course of chronic infections (e.g. those caused by Hepatitis C Virus and Human Immunodeficiency Virus).
Finally, it is important to stress that in the RME model diversity by itself cannot promote facilitation or cooperation since the interaction matrix do not evolve. In fact, as in the canonical RME model there is always just competition. Therefore we explicitly introduce some dose of cooperation measured by f_{C}. This is a nuance to ref. [10] in which the whole quasispecies diversity emerges as relevant in evolutionary terms. An evolving interaction matrix from which cooperation might emerge would imply a different and quite more complex model. Indeed this is a very interesting aspect concerning viral evolution, which is beyond the scope of this work but worth analyzing in the future.
Method
Model
We approach the quasispecies as an ecosystem, what implies that the fitness of each viral phenotype or virus variant results from the combination of its intrinsic properties (collected in a maximum growth rate parameter) and those resulting from its interactions (interaction coefficients, negative in the case of competition and positive for cooperation) with the other variants. Thus, in a very simplified way, each variant is represented by a single quantity which can be thought as an aggregate of several relevant properties like the affinity of the virus variant to bind to a cell receptor, its resistance to interferon, etc. In other words we are considering a hypothetical continuous onedimensional niche axis. One may think in the case of ordinary ecosystems, involving plants or animals, that the 'position' of a species on this niche axis is related to its body size.
Replicator Dynamics plus Mutations
We denote by x_{i }the relative abundance or frequency of phenotype i and then we have . The fitness of phenotype i is denoted by f_{i}, a nonnegative real number corresponding to the rate at which the phenotype i replicates. The average fitness of the population is given by . During replication of a genome, mistakes can happen. The probability of replication of a variant i results in variant j is given by q_{ij }and then these mutation coefficients obey . Since q_{ii }is the probability that the variant doesn't change by mutation, the probability that it does is given by 1q_{ii}. A common equation used to model replication of viruses is the quasispecies equation:
A problem with this equation is that the fitness of a particular type is set to a constant value. We want to allow the fitness landscape to incorporate the distribution of the population variants rather than setting the fitness f_{i }of a particular variant constant. Alternatively, one could treat the swarm of virus variants as an ecosystem and use the LotkaVolterra competition (LV) equations:
where r_{i }is the maximum per capita growth rate and the coefficients α_{ij }represent the interaction of variant j over variant i. These N LV equations are in turn mathematically equivalent to a set of N+1 replicator equations [22] given by:
where
The replicator equation succeeds in incorporating to the fitness landscape the population distribution of variants but has the problem that, since it does not incorporate mutations, it is noninnovative. Hence we will resort to a hybrid model, a kind of replicator with mutations or replicatormutator equations, given by
to describe N virus variants evolving by selection and mutations. On the one hand, if there were no mutations, i.e. q_{ij }= 1 if i = j and 0 otherwise, equation 4 would reduce to the well known replicator equations. On the other hand, by replacing by a fixed f_{i}, we would recover the quasispecies equations. Therefore, this hybrid model generalizes the quasispecies and replicator equations joining the advantages of both of them: selection and innovation. Moreover, previously reported RME describing frequency dependant selection and mutation, has been used in population genetics [23] and language evolution [24].
Niche overlap, interaction coefficients and mutation probabilities
To compute the coefficients α_{ij}, adopting an ecological point of view, let's consider that the interaction between species (virus variants in our case) depends on the niche overlap. That is, on the 'proximity' along the niche axis: the closer the stronger. A given variant i is represented by a certain probability distribution position on the niche axis (ξ): P (ξ) around its average position μ_{i}.
A reasonable assumption for the shape of P (ξ) is the normal distribution . To avoid border effects, the niche axis is defined circular, i.e. we impose periodic boundary conditions, so that each variant has equal numbers of competitors on both sides. The intensity of the interaction between variants i and variants j is related to niche overlap, and thus to the probability P that individuals of the two variants are at the same position on the niche axis, which is the product of both probabilities. For two variants with probability distributions centered at μ_{1 }= 0.4 and μ_{2 }= 0.6 and both with σ = 0.1, a schematic picture is given (Figure 4)
Figure 4. RNA molecules niche position (ξ) and overlapping. Scheme showing two normal probability distributions centered at μ_{1 }= 0.4 and μ_{2 }= 0.6 of two RNA molecules, both with σ = 0.1. The region in black corresponds to the overlap between the two normal distributions.
Hence the absolute values of the interaction coefficients α_{ij } can be expressed as the ratio of the probability of matching an individual of variant j and the probability of matching one of the same variant, that is [15]:
(hence α_{ii }= 1). The sign of the α_{ij }is negative or positive depending if the interaction is competitive or cooperative, respectively.
We assume that the mutation probabilities q_{ij }also depend on the proximity along the niche axis, and then they are proportional to the absolute value of the α_{ij }for the case of i different from j. The proportionality factor m is chosen in such a way that the probability of an error in each variant replication P_{m }≡ 1 ⟨q_{ii}⟩ where the brackets ⟨..⟩ denoting the average over the population, is of the order of 10%. In any case, we checked that the model is quite robust against variation of this parameter. Figure 5 shows the evolution of the population profile for σ = 0.15 after 250 and 500 generations. From panels (A) to (C) it becomes clear that the initial phenotype is disappearing (indeed after 1000 generations it becomes 0 for all practical purposes).
Figure 5. The independence from the initial state. Panel (A) shows the evolution of the population profile for σ = 0.15 starting from only one viral variant at an arbitrary position μ_{j }on the niche axis (that is, x_{i }= 0 for all i ≠ j and x_{j }= 1). Panel (B) and (C) corresponds to the population profile after 250 and 500 generations, respectively. Emergence of three different lumps becomes clearer after 500 generations.
Moreover, we checked that the emergence of this pattern is independent from the chosen initial configuration, e.g. the same occurs if one takes an arbitrary distribution of random 'seeds' consisting in many different viral phenotypes distributed along the niche axis. In this case the diagonal probabilities q_{ii }are computed using the normalization condition .
We proceed by partitioning the segment [01] representing the niche axis into N = 100 subsegments each corresponding to a given virus variant, i.e. the viral phenotypes are binned into N = 100 categories. We start with these variants uniformly distributed at positions μ_{i }on the niche axis, each with the same niche width given by the standard deviation σ (typically σ = 0.15). Then we let the system evolve over N_{g }generations according to eq. 4.
Initially we consider the case of pure competition, i.e. all the α_{ij }are negative. Next we allow some dose of cooperation or facilitation between species. We implement this in the simplest way: a fraction f_{C }of the interactions, chosen at random, turns from competitive to cooperative i.e. if some of the α_{ij }change of sign and become positive. As in the case of pure competition the whole interaction matrix remains fixed during evolution. Thus in particular we want to stress that the interaction between two given variants doesn't switch between being cooperative and defective during the evolution of the quasispecies.
See Table 1 for brief description of parameters involved in the model.
Authors' contributions
HF and JA conceived and designed the model and experiments. HF, JA and SM analyzed the data and wrote the paper. All authors read and approved the final manuscript.
Acknowledgements
We would like to thanks Esteban Domingo for useful discussion and critical reading of this manuscript. We also would like to thanks administrative and financial support from PEDECIBA and Agencia Nacional de Investigación e Innovación.
References

Eigen M, Schuster P: The Hypercycle: a principle of natural selforganization. SpringerVerlag New York; 1979.

Domingo E, Escarmis C, Sevilla N, Moya A, Elena SF, Quer J, Novella IS, Holland J: Basic concepts in RNA virus evolution.

Eigen M, Biebricher CK: Sequence space and quasispecies distribution. In RNA Genetics. Volume 3. Edited by Domingo E, Holland J, Ahiquist P. CRC Press. Inc., Boca Raton, Florida; 1998::211245.

Clarke D, Duarte E, Elena S, Moya A, Domingo E, Holland J: The red queen reigns in the kingdom of RNA viruses.
PNAS 1994, 91:48214824. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

May RM: Stability and complexity in model ecosystems. Princeton; Princeton University Press; 1974.

GonzálezLópez C, Arias A, Pariente N, GomezMariano G, Domingo E: Preextinction viral RNA can interfere with infectivity.
J Virol 2004, 78(7):33193324. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Turner PE, Chao L: Sex and the evolution of intrahost competition in RNA virus ϕ6.
Genetics 1998, 150:523532. PubMed Abstract  PubMed Central Full Text

Turner PE, Chao L: Escape from prisoner's dilemma in RNA phage ϕ6.
Am Nat 2003, 161(3):497505. PubMed Abstract  Publisher Full Text

Turner PE: Prisoner's dilemma in an RNA virus.
Nature 1999, 398:441443. PubMed Abstract  Publisher Full Text

Vignuzzi M, Stone JK, Arnold J, Cameron C, Andino R: Quasispecies diversity determines pathogenesis through cooperative interactions in a viral population.
Nature 2006, 439:344348. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pfeiffer J, Kierkegaard K: Increased fidelity reduces poliovirus fitness and virulence under selective pressure in mice.
PLoS Pathog 2005, 1(2):e11.
102110
PubMed Abstract  Publisher Full Text  PubMed Central Full Text 
Nowak MA: Evolutionary Dynamics: Exploring the Equations of Life. Cambridge, USA: Harvard University Press; 2006.

Page KM, Nowak MA: Unifying evolutionary dynamics.
J Theor Biol 2002, 219:9398. PubMed Abstract  Publisher Full Text

Levins R: Evolution in Changing Environments. Princeton: Princeton University Press; 1968.

MacArthur RH, Levins R: The limiting similarity, converge, and divergence of coexisting species.
Am Nat 1967, 101:377385. Publisher Full Text

Scheffer M, van Nes E: Selforganized similarity, the evolutionary emergence of groups of similar species.
PNAS USA 2006, 103:62306235. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fort H, Scheffer M, van Nes E: The paradox of the clumps mathematically explained.
Theor Ecol 2009, 2(3):171176. Publisher Full Text

Sierra S, Dávila M, Lowenstein P, Domingo E: Response of FootandMouth disease virus to increased mutagenesis: influence of viral load and fitness in loss of infectivity.
J Virol 2000, 74(18):83168323. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Perales C, Mateo R, Mateu M, Domingo E: Insights into RNA virus mutant spectrum and lethal mutagenesis events: replicative interference and complementation by multiple point mutants.
J Mol Biol 2007, 369:9851000. PubMed Abstract  Publisher Full Text

Drake JW, Holland JJ: Mutation rates among RNA viruses.
PNAS 1999, 96(24):1391013913. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

GrandePérez A, Lázaro E, Lowenstein P, Domingo E, Manrubia S: Suppression of viral infectivity through lethal defection.
PNAS 2005, 102:44484452. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hofbauer M, Sigmund K: Evolutionary Games and Population Dynamics. Cambridge, UK: Cambridge University Press; 1998.

Hadeler KP: Stable polymorphisms in a selection model with mutation.
SIAM J Appl Math 1981, 41:17. Publisher Full Text

Nowak MA, Komarova N, Niyogi P: Evolution of universal grammar.
Science 2001, 291:114118. PubMed Abstract  Publisher Full Text