key: cord-322395-xeq4yix4
authors: Nurtay, Anel; Hennessy, Matthew G.; Sardanyés, Josep; Alsedà, Lluís; Elena, Santiago F.
title: Theoretical conditions for the coexistence of viral strains with differences in phenotypic traits: a bifurcation analysis
date: 2019-01-09
journal: R Soc Open Sci
DOI: 10.1098/rsos.181179
sha: 
doc_id: 322395
cord_uid: xeq4yix4

We investigate the dynamics of a wild-type viral strain which generates mutant strains differing in phenotypic properties for infectivity, virulence and mutation rates. We study, by means of a mathematical model and bifurcation analysis, conditions under which the wild-type and mutant viruses, which compete for the same host cells, can coexist. The coexistence conditions are formulated in terms of the basic reproductive numbers of the strains, a maximum value of the mutation rate and the virulence of the pathogens. The analysis reveals that parameter space can be divided into five regions, each with distinct dynamics, that are organized around degenerate Bogdanov–Takens and zero-Hopf bifurcations, the latter of which gives rise to a curve of transcritical bifurcations of periodic orbits. These results provide new insights into the conditions by which viral populations may contain multiple coexisting strains in a stable manner.

The combination of very large population sizes, very short generation times, and lack of proof-reading mechanisms during genome replication confer viral populations with an extremely high evolutionary plasticity that allow them to quickly adapt to system. Nearly all mathematical models in epidemiology detecting various dynamical behaviours with multi-strain infections illustrate the necessity of numerical approaches and the dependency of such models on a large number of parameters [39] [40] [41] [42] [43] . A classical approach used in epidemiology [44] to circumvent this difficulty is to introduce dimensionless parameter groups, such as the basic reproduction number R 0 , as in [45] . However, the number of dimensionless groups can still become large as additional complexity is introduced into the model, as is the case here. Thus, we perform a bifurcation analysis to systematically track how the dynamics of the system change as multiple parameters are varied. In general, the study of parameters in terms of their effect on the stability of certain states of a model is a highly effective way to gain important insights into the investigated system [46] [47] [48] .

The paper is organized as follows. We first introduce a detailed description of the model in §2. Then, in §3, we provide the equilibrium points of the system, study their stability and investigate, in detail, the effect of phenotypic differences in infection rate, virulence and mutation rate. The biological interpretations of the results are present throughout the work, however, concrete statements are found in §4. Some technical details of the bifurcation analyses are discussed in appendix A.

Here we introduce the mathematical model describing the infection dynamics of wt and mutant strains. The model is based on a nonlinear system of five ordinary differential equations. The state variables of the model are: uninfected susceptible cells, x; two different virus strains, given by the wt (z w ) and the mutant (z m ) species; and two types of infected cells, one infected by the wt strain, y w , and another one infected by the mutant strains, y m (figure 1). The time evolution of the interacting populations is described with the following model: . Schematic diagram of the processes included in the mathematical model for the main system investigated, which does not consider backward mutations. The system is composed of two pathogenic strains: wild-type (variable z w ) and mutant (variable z m ) species, that compete for the infection of healthy cells (variable x, green cells). Infection of healthy cells gives rise to two different populations of infected cells with wt (variable y w ) and mutant (variable y m ) strains, displayed by red and blue cells, respectively. Viral strains are assumed to grow and mutate (dashed lines) within the cells, being released again to the system for further infection after cell lysis. The terms in model (2.1) are labelled with the biological processes. The meaning of these processes as well as of the parameters is discussed in the following lines. The uninfected host cells, which are limited by the carrying capacity of their environment, K, proliferate and die proportionally to parameters b and d, respectively.

As previously mentioned, the mutant and the wt strains infect the host cells at rates a m and a w , respectively. Mathematical and statistical models of multiple infections, coinfection and superinfection have been studied in detail [49] [50] [51] [52] and are not considered here. Acknowledging the importance of differences between lytic and lysogenic infection cycles in the production of virions, yet unlike [53] [54] [55] , we consider only lytic infections, i.e. the populations of infected cells do not grow as uninfected cells do [56] . In fact, the model is based on the same assumptions of the classic Lotka-Volterra equations, which were adapted to virus dynamics by Nowak & May [57] and many others [58] [59] [60] [61] [62] . In the study of coexistence of viral populations, the main role of infected cells is viewed in the scope of the processes of the lytic cycle, i.e. a direct impact in change and production of free virions. Viruses replicate and produce their own type of virions via infected cells. Strains can change and mutations occur only during the process of viral replication inside an infected cell (see dashed arrow in figure 1 ). Therefore, cells from population y w infected by the wt strain can mutate at rate m into cells of population y m . That would increase the size of the latter population at the same rate m. We note that our analyses are mainly developed considering mutations from the wt to the mutant strains. This is a standard strategy to keep the model as simple as possible, assuming that the probability of backward mutations is extremely small due to the enormous size of the sequence space. However, some results considering backward mutations will also be presented. While considering mutation in infected cells, we implicitly study the mutation of the viral genomes. It is clear that infection affects the life-span of infected cells in a nontrivial way [63] . In this study, the virulence of the strains is depicted by including new decay mechanisms for the uninfected cells. However, for simplicity, we consider death rates of infected cells to be g m and g w , and as opposed to d þ Dg m and d þ Dg w .

The absence of a mechanism for virus replication makes the multiplication of viral strains entirely dependent on the machinery of the infected cells. Therefore, the overall number of virions produced in one lytic cycle must be proportional to the virulence and the number of virions produced by a single-infected cell, the latter of which we refer to as the burst size. We average the burst sizes of viral strains and consider them as constants k m and k w for mutant-and wt-infected cells, respectively. Au contraire, the infection from the perspective of the virus is associated with an average number of virions spent to ensure a successful infection process. Owing to the absence of co-infection by different strains in our model, this number can be considered as the multiplicity of infection of the viral strain. We take the multiplicity of infections to be constants n m and n w for mutant and wt strains, respectively. In this model, we assume that populations are not being harvested, but do consider a decay of all the populations in the system. In the case of the virus populations, this decrease is described as an outflow or 'death' of free virus particles from the system. The overall 'death' rates of the virions in the system are z w and z m for wt and mutant strains, respectively. Initial conditions for system (2.1) are non-negative values:

The model (2.1) can be simplified by introducing dimensionless variables that are based on characteristic timescales and population sizes. The quantity b 2 d describes the effective growth rate of uninfected cells and its inverse, (b 2 d) 21 , is used to define the characteristic time scale of the system. In a virus-free environment, the maximum size of the uninfected cell population isx max ¼ (1 À d=b)K, which is used to define the characteristic population size for both the uninfected and infected cells. The characteristic population size of the viral strains is chosen to be the product of the mutant burst size k m and the characteristic size of the infected cell population,x max . The former is required due to differences in sizes and measurement units of the viral loads and the cell populations in the system. We therefore non-dimensionalize the variables according to

where bars are used to denote dimensionless quantities. We also introduce dimensionless parameters defined by

where p i (recall that i ¼ m, w) stands for other 'rate' parameters, namely: g m , g w , z m and z w . Notice that all non-dimensional 'rate' parameters, including mutation m, are relative rates, i.e. the rate relative to the effective growth rate of uninfected cells. Meanwhile, non-dimensional parameters related to the virus populations are inevitably linked to the burst size due to the choice of z i . The non-dimensional burst size k is simply the ratio of the wt and mutant burst sizes. Both multiplicity of infections are scaled with respect to the burst size of the mutant strain. The dependence of the dimensionless infection rates a i on the burst size of the mutant-type-infected cells, the carrying capacity, growth rate of uninfected cells and dimensional infection rate shows how each of these quantities affects the overall rate of infection.

Understandably, the non-dimensionalization places a restriction on the parameter values and requires b . d. However, this restriction is biologically justified. Taking b . d, as shown later, forces the trivial equilibrium to be unstable, thus avoiding scenarios where all populations become extinct. Upon ignoring bars for clarity purposes, we obtain the non-dimensionalized system

and _ z w ¼ kg w y w À n w a w z w x À z w z w : (2:4e)

The dimensionless model (2.4) is now analysed to understand how the dynamics change under parameter variation. We calculate the equilibria of (2.4), conduct a linear stability analysis, and identify analytical conditions that lead to transcritical and Hopf bifurcations. We find that curves of these bifurcations can intersect at specific points in parameter space, giving rise to degenerate Bogdanov-Takens (DBT) and zero-Hopf (DZH) bifurcations. The role of DBT and DZH bifurcations is to organize the stability (or phase) diagram into regions with distinct dynamics. The numerical continuation package MATCONT [64] is used to track how the equilibria and bifurcations evolve as the infection rate, virulence, and mutation rate are varied. The MATCONT source code and data are available online [65] .

The non-dimensional model given by equations (2.4) has four equilibria. For mathematical convenience, we define the population vector v(t) ¼ {x(t), y m (t), y w (t), z m (t), z w (t)}. The first equilibrium is the origin

If stable, the trivial solution v 0 corresponds to the extinction of all of the populations. The second equilibrium

describes, whenever stable, a virus-free state whereby only the uninfected cells persist. Owing to our choice of non-dimensionalization, the maximum population of uninfected cells is equal to one. The third equilibrium is given by

and, if stable, corresponds to the persistence of uninfected cells and the mutant strain of the virus. There is no wt strain of the virus and the viral population is only composed of mutant genotypes. In order for the wt-free state v 2 to be biologically meaningful, the condition n m , 1 must hold. This condition corresponds to the burst size of the mutant virus being greater than its multiplicity of infection. Throughout the remainder of the paper, it will be assumed that n m , 1. Finally, the fourth equilibrium point is given by

where A ¼ g w k À (m þ g w )n w , B ¼ À((n w À k)a w þ z w )g w À m(a w n w þ z w ) and

C ¼ a m z w (g w (1 À n m ) À mn m ) À a w z m (g w (k À n w ) À mn w ):

This equilibrium point involves, whenever stable, a state of coexistence for the wt and mutant strains. Similar to the wt-free state v 2 , the coexistence state v 3 can only be biologically meaningful if the burst size of wt strain satisfies k . n w (1 þ m/g w ). This inequality is a generalization of that derived for the mutant virus (n m , 1) which accounts for mutation. The singularity that occurs in the coexistence state v 3 when C ¼ 0 leads to difficulties when using numerical methods to track how this equilibrium evolves under parameter variation. The assumption of uni-directional mutation prevents the existence of an equilibrium point that is analogous to v 2 whereby only the wt virus exists.

A linear stability analysis is carried out to determine the behaviour of the system close to the equilibrium points. The Jacobian matrix for equations (2.4) is given by Straightforward calculations for the trivial solution v 0 yield eigenvalues of the Jacobian given by

Thus, for all parameter values, the trivial solution v 0 , as mentioned before, is always unstable, i.e. it is a saddle point with a one-dimensional unstable manifold.

Although it is possible to find analytical expressions for the eigenvalues of the Jacobian evaluated at v 1 , they are sufficiently complicated that little insight is gained from analysing them directly. In order to conduct a stability analysis for v 1 , it is more useful to consider the characteristic polynomial for the eigenvalues l. The polynomial det (J(v 1 ) À lI) ¼ 0 can be regrouped into the form of three factors

and

The first factor, P 1 , provides a constant eigenvalue, l ¼ 21, which reserves the possibility for stability of v 1 . Although the next two factors P 2 and P 3 do not provide obvious eigenvalues, they enable the identification of critical parameter groups for which v 1 undergoes a bifurcation. Based on our knowledge of the bifurcations that can occur in models similar to (2.4), we may expect to find transcritical and Hopf bifurcations. For a varying parameter (or set of parameters), the onset of a Hopf bifurcation leads to the creation of periodic orbits (POs) after the change of stability of the equilibrium point. Furthermore, Hopf bifurcations are characterized by the Jacobian matrix having a single pair of complex conjugate eigenvalues with zero real part. An analysis of the factors P 2 and P 3 given by (3.7b,c) reveals that v 1 cannot undergo a Hopf bifurcation. This is because the coefficients of the linear terms are strictly positive, thus preventing the eigenvalues from ever being purely imaginary.

Transcritical bifurcations occur when two equilibria collide non-destructively, exchanging their stability and resulting in the Jacobian matrix having a single eigenvalue that is equal to zero. Since v 1 exists for all parameter combinations, transcritical bifurcations can be straightforwardly detected by forcing l to be zero in P 2 and P 3 . By solving

It will be shown below that R w 0 and R m 0 are the basic reproductive numbers for the wt and mutant strains, respectively. When R w 0 ¼ 1, the virus-free state v 1 and the coexistence state v 3 intersect. By expanding v 3 around R w 0 ¼ 1, we find that it has negative components when R w 0 , 1 and thus lies outside of the biologically meaningful phase space. However, all of the components of v 3 become positive when R w 0 . 1. Likewise, expanding the equilibrium v 2 around R m 0 ¼ 1 shows that this condition corresponds to points where the virus-free state v 1 and wt-free state v 2 intersect. For R m 0 , 1, some components of v 2 are negative; for R m 0 . 1, all components are positive. In the case when both R w 0 ¼ R m 0 ¼ 1 hold, there is a triple intersection of v 1 , v 2 and v 3 . Furthermore, the Jacobian has a double zero eigenvalue at this point, indicating the onset of a non-degenerate or degenerate Bogdanov-Takens bifurcation. As will be shown in §3.1, the Bogdanov-Takens bifurcations in this model are of degenerate type. More generally, we find that v 2 and v 3 intersect when

Thus, there are three curves of transcritical bifurcations defined by

The dynamics of the system near the DBT bifurcation can be determined by expanding the equilibria and their eigenvalues around R w 0 ¼ 1 and R m 0 ¼ 1. The eigenvalues of all of the equilibria are real, ruling out the possibility of a branch of Hopf bifurcations emanating from the DBT point, which is a generic feature of non-degenerate Bogdanov-Takens bifurcations [66] . When R w 0 , 1 and R m 0 , 1, the virusfree state v 1 is locally asymptotically stable and the wt-free state v 2 and the coexistence state v 3 are unstable, meaning that the virus-free state is achieved. For R m 0 . 1 and R w 0 , 1, the virus-free and wtfree state exchange stability, with v 1 becoming unstable and v 2 locally asymptotically stable, with v 3 remaining unstable. Similarly, for R w 0 . 1 and R m 0 , 1, the virus-free and coexistence state exchange stability: v 1 becomes unstable, v 3 becomes locally asymptotically stable, and v 2 remains unstable. The stability in the region R m 0 . 1 and R w 0 . 1 is governed by the transcritical bifurcation occurring along 7 royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181179 the curve R m 0 ¼ R w 0 , which leads to an exchange of stability between wt-free and coexistence states, v 2 and v 3 . The transcritical bifurcations that occur for

with R w 0 , 1 do not result in an exchange of stability. Thus, the local stability diagram near the DBT bifurcation can be drawn as in figure 2. There are three distinct regions, denoted by I, II and III, where v 1 , v 2 and v 3 are local attractors, respectively.

Based on these results, we can interpret the quantities R m 0 and R w 0 as basic reproduction numbers for the mutant and wt viral strains, respectively. In general, for a strain to persist, its basic reproduction number has to be strictly greater than one. Besides, the value of a basic reproduction number is proportional to the number of new infections arising in a following generation of virions from an infected cell. Identifying basic reproduction numbers enables different models of population dynamics to be compared in a consistent manner. However, as it is not possible to write the dimensionless model (2.4) solely in terms of the basic reproduction numbers R m 0 and R w 0 , we will consider how the dynamics change under the variation of specific individual parameters that are of biological interest.

To examine the stability of the wt-free state v 2 , we obtain the eigenvalues from the characteristic polynomial computed from det (J(v 2 )ÀlI) ¼ 0. This equation can be factorized and rewritten as Q 2 (l) Á Q 3 (l) ¼ 0, where Q 2 and Q 3 are quadratic and cubic polynomials in l, respectively. The exact form of Q 2 is not required here. We write

Following the analysis scheme discussed earlier, to obtain conditions for Hopf bifurcations of v 2 , we search for a purely imaginary pair of eigenvalues. Setting Q 3 (l ¼ +iv)¼ 0, we obtain the critical condition for a Hopf bifurcation,

being the angular frequency of the emerging POs. Interestingly, for this Hopf bifurcation, there is no dependence on parameters associated with the wt strain of virus. For all biologically meaningful solutions of ad ¼ bc, the mutant strain of virus has the potential to gain periodic behaviour through the creation of a stable PO. The other factor, Q 2 (l), provides no possibility for a Hopf bifurcation due to a strictly positive linear coefficient in the quadratic polynomial.

Transcritical bifurcations of v 2 may occur if

which is equivalent to R m 0 ¼ 1 and corresponds to an intersection of v 1 and v 2 . Second, from Q 2 (l ¼ 0) ¼ 0, we obtain an expression which matches with (3.10), that is, the parameter combination yielding an intersection of v 2 and v 3 .

As will be shown in §3.3, it is possible to simultaneously satisfy Q 3 (l ¼ +iv)¼ 0 and Q 2 (l ¼ 0) ¼ 0, implying that the Jacobian matrix at v 2 ¼ v 3 has a pair of purely imaginary complex conjugate eigenvalues and a zero eigenvalue. This corresponds to the onset of a degenerate zero-Hopf (DZH) bifurcation. Like the DBT bifurcation, the DZH bifurcation will be shown to have an unusual structure that does not coincide with any of the standard normal forms (this point is further discussed in appendix A). In particular, our analysis reveals that a curve of global transcritical bifurcations of POs (TPO bifurcation) emanates from the DZH point rather than a curve of torus bifurcations [66] . A TPO bifurcation occurs when an unstable PO and a stable PO collide with each other in a non-destructive way (differently from what would happen in a saddle-node bifurcation of POs, with destruction of POs), and subsequently exchange stability. The detailed stability and bifurcation analysis of v 2 , v 3 , as well as the DZH point, will be performed numerically. In particular, we will construct one-and two-dimensional bifurcation diagrams using specific pairs of parameters.

We first examine the influence of the infection rates a m and a w on the dynamics. The values of the other parameters are set to

We begin our investigation by constructing one-dimensional bifurcation diagrams using a m as the bifurcation parameter and fixing a w . We first consider two values of the wt strain infections rate given by a w ¼ 0.5 and a w ¼ 2, corresponding to basic reproduction numbers R w 0 ¼ 0:48 and R w 0 ¼ 1:93, respectively. These values of a w therefore lie on opposite sides of the curve of transcritical bifurcations involving v 1 and v 3 defined by R w 0 ¼ 1. The resulting bifurcation diagrams are shown in figure 3 . Solid and dashed lines represent stable and unstable equilibria, respectively. Solid circles represent maxima and minima of stable POs. Blue and orange lines denote equilibria that exist in biologically meaningful and non-meaningful phase space, respectively. The parametric dependence of the population of uninfected cell population x is used as this is the only component that differs between all four equilibria. Figure 3a shows, for the case a w ¼ 0. When the infection rate of wt virus is increased to a w ¼ 2, the virus-free and coexistence states undergo a transcritical bifurcation. Consequently, the bifurcation diagram in figure 3b shows that for mutant-type infection rates given by a m , 0.86, the coexistence state v 3 is stable whereas the virus-free state v 1 is unstable. This implies that Region I (v 1 stable) is replaced by Region III (v 3 stable). Although v 1 and v 2 undergo a transcritical bifurcation at a m ¼ 0.44, there is no exchange of stability. Instead, it is v 2 and v 3 which exchange stability through a transcritical bifurcation at a m ¼ 0.86, corresponding to the case of equal basic reproduction numbers, R m 0 ¼ R w 0 ¼ 1:93. Both Regions II and IV persist for a w ¼ 2, with the supercritical Hopf bifurcation occurring at a m ¼ 1.55, the same location as in figure 3a .

The one-dimensional bifurcation diagrams shown in figure 3 can be rebuilt using the infection rate of the wt virus a w as the bifurcation parameter for fixed values of a m that lie on opposite sides of R m 0 ¼ 1. However, the resulting diagrams are qualitatively similar to those in figure 3 , with the exception that the wt-free state v 2 is switched with the coexistence state v 3 and vice versa. The supercritical Hopf bifurcation now occurs from v 3 and thus gives rise to Region V, characterized by the periodic coexistence of the viral populations.

The five regions identified from the bifurcation analysis can be conveniently visualized by constructing a two-dimensional bifurcation diagram where both infection rates a m and a w are continuously varied. In this diagram, displayed in figure 4a, the locations of the bifurcations that separate the five regions are traced out as the infection rates vary. Not all of the bifurcations shown in figure 4a lead to a biologically meaningful change in the dynamics, i.e. they do not represent a boundary between two distinct regions. Thus, in figure 4b we represent the two-dimensional diagram with only the biologically meaningful bifurcations shown. Figure 4a reveals a complicated bifurcation scenario that is largely centred about DBT and DZH bifurcations occurring at a w ¼ 1.04, a m ¼ 0.44 and a w ¼ 3.62, a m ¼ 1.55, respectively. The DBT bifurcation lies at the intersection of three curves of transcritical bifurcations satisfying R w 0 ¼ 1, R m 0 ¼ 1 and R m 0 ¼ R w , which may be written, respectively, in terms of the infection rates as

(3:13)

The bifurcation structure around the DBT point is identical to that predicted from the local stability analysis. Thus, we can eliminate from figure 4a the biologically irrelevant curves of transcritical bifurcations emerging from the DBT point in order to obtain the boundaries between Regions I, II and III shown in figure 4b. The DZH bifurcation occurs at the simultaneous intersection of two curves of Hopf bifurcations associated with v 2 and v 3 and the curve of transcritical bifurcations involving v 2 and v 3 . As predicted from linear stability analysis, the Hopf bifurcation curve associated with the wt-free state v 2 does not depend on the infection rate of the wt virus a w and thus appears as a straight line given by a m ¼ 1.55. Emanating from the DZH point is a curve of TPO bifurcations. To better understand the dynamics that occur near the DZH point, one-dimensional bifurcation diagrams are obtained by setting a w ¼ 8 and treating a m as a bifurcation parameter and then setting a m ¼ 2 and treating a w as the bifurcation parameter. The resulting diagrams are shown in figure 5 , with open circles denoting unstable POs. Figure 5a shows that as a m is increased from zero when a w ¼ 8, the wt-free state v 2 first undergoes a transcritical bifurcation with the virus-free state v 1 , then a subcritical Hopf bifurcation, and finally a transcritical bifurcation with the coexistent state v 3 . None of these bifurcations change the stability of the equilibria. Thus, the associated curves of transcritical and Hopf bifurcations in figure 4a do not represent boundaries between the five characteristic regions and are not shown in figure 4b. For these infection rates, all of the equilibria are unstable and the dynamics are therefore determined by the stability of POs. For sufficiently small values of a m , the system is in Region V, and there is a stable PO around the coexistence state v 3 . This PO is created by the Hopf bifurcation that defines the boundary between Regions III and V shown in figure 4b. As a m is increased, the unstable POs around the wt-free state v 2 that emerge from the subcritical Hopf bifurcation grow in size. Eventually, this unstable PO collides and exchanges stability with the stable PO around v 3 through a TPO bifurcation at a m ¼ 2.3972. During the TPO bifurcation, the stable and biologically meaningful PO around v 3 becomes unstable and some components enter negative phase space. Afterwards, the PO around the wt-free state is stable and the system is in Region IV.

The one-dimensional bifurcation diagram created with a m ¼ 2 and shown in figure 5b shares the same qualitative features as figure 5a, with v 2 and its POs exchanging roles with v 3 and its POs. None of the bifurcation of equilibria have biological significance. The main difference in this case is that the subcritical Hopf bifurcation and the POs it creates exist outside the biologically meaningful phase space. As these POs are involved in a biologically meaningful TPO bifurcation, this figure demonstrates how it can still be useful to understand and resolve features that lie outside of the biologically meaningful space.

Upon removing the biologically redundant bifurcation curves from figure 4a, the stability diagram in figure 4b is obtained. The transcritical bifurcations at the boundaries between Regions I and II, and I and III, capture the transition from virus-free to virus-persistent states that tend to a stable equilibrium. Supercritical Hopf bifurcations separate Regions II from IV and III from V, and mark the transition between stationary and periodic population dynamics. The curve of TPO bifurcations separates Regions IV and V, both of which are characterized by periodic viral populations, in the same way that a curve of transcritical bifurcations separates Regions II and III, which describe stationary viral populations. Thus, despite the system exhibiting a range of complex dynamics, they can be elegantly classified and organized in terms of the regions shown in figure 4b .

Overall, in figure 4b , we notice quite a large area of coexistence, Regions III and V. In population genetics, this corresponds to the most simple case of the emergence and persistence of a polymorphism in a population and maintenance of biodiversity. An example of an experiment that illustrates a simplified competition was conducted for Escherichia coli in [67] . These experiments focused on studying the behaviour of a newly emerging mutant in a population of a few strains of bacteria that compete with each other in a spatially homogeneous environment for the same type of nutrients. The results of the competition experiments demonstrate that in the vast majority of cases, the competitors stably coexist in steady or periodic states, which align with the outcome of our theoretical model. When even a slight change is possible in the genome, new mutations will appear. However, if all the other characteristics of an initial and a mutant population are the same, then, for both populations to persist, the initial population must have larger fitness than a mutant population. Otherwise, the initial population will be out-competed because of the constant 'leak' into the mutant population. A more detailed analysis of the impact of mutation rates in the dynamics can be found in §3.5.

We now study the effect of the strains' virulence ( parametrized by g w and g m ) on the dynamics of the system by building bifurcation diagrams. The infection rates are chosen to be from the different regions of the stability diagram shown in figure 4b , which are based on the reference values g m ¼ g w ¼ 0.25. The other parameters are fixed to

We focus on Regions I, III and IV of figure 4b as a starting point, covering all the other regions from there.

We begin with the case where the infection rates are chosen to coincide with Region I at the virulence reference values. These infection rates therefore correspond to basic reproduction numbers that are less than one. From the definition of R m 0 given by (3.9) , we see that the basic reproduction number of the mutant is independent of the virulence. Thus, changes in g m or g w cannot increase R m 0 beyond one. Therefore, Regions II and IV, where the mutant-type virus persists at the expense of the wt virus becoming extinct, cannot be entered. The basic reproduction number for the wt virus given by (3.8) is an increasing function of the virulence of the wt virus. In the limit of very large virulence, g w ! 1, we find that R w 0 ! a w (k À n w )=z w . Thus, if the infection rate a w is so small that the limit of R w 0 is less than one, then it will not be possible to leave Region I and both strains of the virus will always become extinct. However, if a w is sufficiently large and the basic reproduction number increases beyond one, then a change in dynamics will be observed as g w is increased. By solving R w 0 ¼ 1, a critical value of the virulence of the wt virus is obtained

The above equation is the condition for a transcritical bifurcation between the virus-free state v 1 and the coexistence state v 3 and defines the boundary between Regions I and III. Thus, only for values of g w . g crit w will the virus persist in the system for this choice of infection rates. The stability diagram for a w ¼ 0.5 and a m ¼ 0.1 has been numerically computed and is shown in figure 6a . These values of the infection rate correspond to Region I at the reference values of the virulence (figure 4b). As predicted, Regions II and IV are absent from the stability diagram. However, Region V is also missing. Thus, for this choice of infection rates, only Regions I and III can be entered by changing the values of the virulence. Regions I and III are separated by a straight dash line given by the critical condition (3.15) .

We now consider infection rates given by a w ¼ 3, a m ¼ 1, which correspond to Region III at the reference virulence (figure 4b). The resulting stability diagram in terms of the virulence is shown in figure 6b . The choice of a m ¼ 1 along with the values of the parameters in (3.14) leads to R m 0 ¼ 2:25. Thus, variations in the virulence cannot bring the system to Region I and at least one type of virus will always persist. A transition between Regions II and III can occur via the transcritical bifurcation between the wt-free and coexistence states v 2 and v 3 . The critical condition for this transcritical bifurcation is given by equation (3.10) , and depends only on the virulence of the wt viral strain. Thus, the transition between Regions II and III appears as the vertical dashed line in figure 6b . Interestingly, this figure shows that as the virulence of the wt virus is increased, the system transitions from Region III to Region V and then back to Region III. These transitions occur via supercritical Hopf bifurcations, denoted by solid lines. This scenario corresponds to the so-called bubble bifurcation (found in other epidemiological systems, e.g. [68] ), in which the system is in a stationary state, then enters into an oscillating one, and finally goes back to a stationary state as the control parameter is changed.

To understand why Region IV does not appear in the stability diagrams on figure 6a ,b, we first recall that Region IV is separated from Region II by a Hopf bifurcation curve of the wt-free state v 2 . This Hopf curve has been calculated analytically from the equality ad ¼ bc, where a, b, c and d are given in (3.11) , and is a quadratic expression for g m independent of g w and a w . Solving ad ¼ bc to find g crit m (not shown due to complexity of the expression), at the values of infection rates chosen for figure 6a,b, we notice they are outside of the positive parameters space. In fact, only for values of the mutant-type infection rate that satisfy

equivalent to a m . 1.3332 for the chosen parameters, does the Hopf curve of v 2 appear on a g w versus g m stability diagram ( figure 6c) . Greater values of a m lead to increases in the area of Region IV under Region II. The stability diagrams of figure 6 indicate the impact of virulence on the complexity of the dynamics. Clearly, the virulence of the mutant strain does not affect the survival of the wt strain. That is, increases in g m do not lead to an appearance or disappearance of the wt strain. However, the survival of the mutant strain does not appear to depend on g m either. Importantly, the virulence of the wt strain, g w , strongly controls the dynamics and the persistence of the wt strain. As can be seen for small values of g w , this strain can become extinct as shown in figure 6b. In other words, even for a superior infection rate of the wt viral strain, there is a threshold of g w that must be surpassed in order for this strain to exist. For inferior values that are below this threshold, the survival of the wt strain is impossible because the rate of wt virion production, i.e. the death of rate infected cells, is insufficient. The qualitative behaviour of solutions, i.e. whether the populations of wt strains undergo oscillations or stabilize to a certain value, depends on the virulence in a non-trivial way. As shown in figure 6b, there is an interval of values for g w for which the system has a stable PO governing the coexistence of strains, marked as Region V. Outside of this interval, the populations reach stable stationary states.

In the previous analyses, the mutation rate has been fixed to m ¼ 0.1. We now consider the effect of m on the dynamics. This is a key parameter that has been largely investigated within the framework of the error threshold [69, 70] and lethal mutagenesis [71] in quasi-species theory. In general, the coexistence of two strains requires the basic reproduction number of the wt virus to be greater than one, R w 0 . 1, which can be interpreted as a condition on the mutation rate: m , g w k z w a w À1 þ n w À 1

:

Thus, for two strains to coexist, the mutation rate must be sufficiently small. This result supports a conjecture that excess mutation exhausts the population of the wt strain, thereby leading to a process similar to the well-known error catastrophe [70] . We would expect that a gradual increase of the mutation rate contributes to a better success of the mutant strain as the frequency of mutants generated de novo increases. The right-hand side of (3.17) is an increasing function of a w , reflecting the fact that a greater rate of infection by the wt virus will offset a greater mutation rate. A stability diagram in terms of the mutation rate m and the wt infection rate a w is shown in figure 7a for the Figure 6 . Two-dimensional bifurcation diagrams for g w versus g m at different sets of a i , with i ¼ m, w. Parameter values are given by (3.14) and (a) (3.12) . The dashed line marking the boundary between Regions I and III, and also defining the region of coexisting populations, has been obtained by replacing the inequality with equality in (3.17) . As the infection rate a w increases, the boundary between Regions I and III reaches a horizontal asymptote given by

For mutation rates that satisfy m , m c , Regions III and V can be entered from Region I by increases in the infection rate, promoting coexistence. However, for m . m c , only the wt-free state can occur. Hence, equation (3.18) defines a critical, finite mutation rate at which coexistence no longer becomes possible due to extinction of the wt virus due to the outcompetition by the mutant strains. To explore this in more detail, we have repeated the stability diagram shown in figure 4b using a value of m ¼ 0.5 . m c . The stability diagram changes drastically, becoming that shown in figure 7b , and contains only three regions: I, II and IV. The stability diagrams in figure 4b and 7b are linked through the fact that as the mutation rate increases, the DBT bifurcation shifts to the right, eventually tending to infinity as the critical value is approached. Finally, figure 7c,d illustrates how the equilibrium populations of infected cells and viral strains change with increasing mutation rate m. Specifically, we have used a w ¼ 3.0 and a m ¼ 1.0, along with the parameters in (3.12), in both panels. As m increases beyond the point where (3.17) is satisfied, the population of mutants (virions and infected cells) outcompetes the wt populations ( figure 7c,d) . This phenomenon is similar to the error threshold defined in quasi-species theory. Here we show that mutation is not only involved in this shift, but also depends (at the infection-cell level) on virulence, burst size and multiplicity of infection.

Throughout the previous analyses, we have assumed that the wt strain of the virus mutates and the new strain has no chance of mutating exactly into the original strain. Although it is highly unlikely, the full picture of the proposed model would require a consideration of a possibility for backwards mutation, especially when modelling phenotypic traits. We therefore replace (2.4b,c) with

We begin our investigation of the role of backwards mutation by re-building the two-dimensional bifurcation diagram and stability diagrams shown in figure 4 using the parameters in (3.12) . The mutation rate of the wt virus is set to m w ¼ 0.1 and we consider two values of the mutation rate of the mutant virus. First, we take m m ¼ 10 23 ( m w , so that the new form of mutation can be considered as a small perturbation to the original system of equations (2.1). Then, we equate the mutation rates of both strains and set m m ¼ 0.1. The two-dimensional bifurcation diagrams in figure 8a ,b reveal that backwards mutation results in several important changes to the dynamics. Importantly, the behaviour of the system can be described using three main states: the virus-free state (analogous to Region I and displayed in white in figure 8 ), stationary coexistence all populations ( polka dot pattern in figure 8, analogous to Region III), and oscillatory coexistence of all the populations (checkerboard pattern in figure 8, analogous to Region V). Furthermore, the DBT and the DZH bifurcations have vanished, the latter of which implies that the curve of TPO bifurcations no longer exists either. The vertical and horizontal lines of transcritical bifurcations defined by R w 0 ¼ 1 and R m 0 ¼ 1, which intersected at right angles in figure 4a , have now merged into two separate branches that do not intersect. The curves of Hopf bifurcations, which also intersected in the case of uni-directional mutation, have also merged into two distinct non-intersecting branches. Finally, the transcritical bifurcation between the wt-free and coexistent states v 2 and v 3 has vanished as well. The resulting stability diagram, shown in figure 8 by patterns, is now considerably simpler and involves only analogues of Regions I, III and V. Thus, the only virus-persistent states are those in which both virus strains coexist. The wt-free state can no longer occur due to the creation of the wt virus through backwards mutation.

In this paper, we studied a mathematical model of population dynamics of two viral strains infecting a population of the same cell type. Unlike previous models of parasite-host-like interaction, we consider the emergence of a second strain by mutation from a single strain introduced into the environment. The mechanism of virus replication forces the consideration of two types of infected cells, one for each viral strain, in the model. We analysed the system of differential equations and its solutions. By studying parameters space, we identified five different regions, each characterized by distinct dynamics: (Region I) the virus-free state which maximizes the population of host cells, (Region II) the stationary existence of mutant virus with nonzero cell populations supporting it, (Region III) the stationary persistence of both viral strains with nonzero cell populations, (Region IV) the oscillating existence of the mutant virus and corresponding cell populations but extinct wt virus population, and finally, (Region V) the oscillating coexistence of all populations. The population of the mutant-type virus exists as long as the mutation rate of the wt strain is positive. The broken symmetry between the wt and mutant strains is clear from the stability diagrams based on infection rates. From our results, we observe that survival of the wt virus is essential for coexistence. However, the growth of the wt virus population jeopardizes its persistence in the system by creating its own competitor: the mutant-type virus. Moreover, the larger the mutation rate, the greater the infection rate of the original wt strain should be in order to remain in the system. Interestingly, we have found a maximum-critical value-of the mutation rate in the context of the model which allows for the persistence of the wt strain and hence coexistence. Values of the mutation rate that are higher than the critical value make coexistence impossible even for the most infectious wt strain. Furthermore, our model shows that the concept of the error threshold may not be considered as a one-parameter-driven effect, i.e. based solely on the mutation rate [69, 70] . The critical value of mutation is shown here to be proportional to the virulence and burst size of the wt virus and proportionally inverse to its multiplicity of infection. Hence, our results extend the phenomenon of the error threshold at the infected cell population level.

Data accessibility. MATCONT source files, output data and MATLAB files used to create the figures in the manuscript are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.56g7v08 [65] . and _ j 1 ¼ b 1 þ b 2 j 0 þ a 2 j 2 0 þ b 2 j 0 j 1 , ( A 1 b) provided that the coefficients a 2 and b 2 are not equal to zero, a 2 b 2 = 0. The quantities b 0 and b 1 in equations (A 1) play the role of bifurcation parameters, with b 0 ¼ b 1 ¼ 0 corresponding to the BT point. The bifurcation diagram of system (A 1) consists of curves of saddle-node and Hopf bifurcations, neither of which were detected in the local analysis of the model analysed in this article (2.4) near the BT point, as well as curves of global bifurcations involving homoclinic orbits. The procedure outlined by Kuzetnsov [72] enables the normal form coefficients to be related to the model studied here. Remarkably, the coefficient a 2 can be calculated analytically and is found to be equal to zero for all parameter combinations. Thus, the BT bifurcation occurring in our model is always degenerate and the local dynamics will be different from those of equation (A 1). The vanishing of a 2 can be rationalized in terms of the number of equilibria at the degenerate Bogdanov-Takens (DBT) point. An equilibrium analysis shows that (A 1) will only have two equilibria at the DBT point b 0 ¼ b 1 ¼ 0 if a 2 = 0; however, the model under study has three. Thus, a 2 ¼ 0 is needed to capture the triple equilibrium at the DBT point in our model. This also suggests that higher-order terms must be included in (A 1) in order for it to capture the DBT bifurcation in the model investigated here. While an extended system of normal form equations for degenerate BT bifurcations is proposed in Kuznetsov [72] , it cannot capture the transcritical bifurcations found in the original model, suggesting that an alternative form is required. However, due to the simple nature of the local dynamics near the degenerate BT point shown in figure 2 , which can be obtained analytically from the full model (2.4), we do not pursue this point further.

The Poincaré normal form of a ZH bifurcation can be written as [66, §8.5]

is a vector of bifurcation parameters. Analyses of the normal form equations for the ZH bifurcation show that local dynamics depend on the values of the normal form coefficients H ijk and G ijk . A common feature between the various cases is that the ZH bifurcation occurs at the tangential intersection of curves of saddle-node and Hopf bifurcations. However, figure 4a shows that the zero-Hopf bifurcation in our model lies at a transversal intersection of curves of transcritical and Hopf bifurcations. Furthermore, none of the normal forms predict that a curve of TPO bifurcations should emanate from a ZH point. The calculation of the three normal form coefficients for the DZH bifurcation is rather involved and must be performed numerically. We find that the normal form coefficient G 011 (0) ¼ 0 across a range of parameter values, indicating the ZH bifurcation is degenerate. The normal form equations for degenerate ZH bifurcations are not well known and previous studies have focused on different cases [73, 74] . Tigan et al. [75] showed that a degenerate ZH bifurcation with G 011 (0) ¼ 0 occurs in a Rö ssler-type system and used averaging theory to detect periodic orbits, leaving the structure of global bifurcations unresolved. We believe, as far as we know, that this is the first time a curve of TPO bifurcations has been observed to emanate from a degenerate ZH bifurcation of this type. We leave determining the corresponding normal form as an interesting area of future work.

Viral quasispecies. Virology 479 -480

Frequency of hepatitis B virus e-minus mutants varies among patients from different areas of China

Coexistence of hepatitis B surface antigen (HBs Ag) and anti-HBs antibodies in chronic hepatitis B virus carriers: influence of "a" determinant variants

Evolution and emergence of plant viruses

Emerging pathogens: the epidemiology and evolution of species jump

Ultra-deep sequencing for the analysis of viral populations

Deep sequencing analysis of viral infection and evolution allows rapid and 17 royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181179 detailed characterization of viral mutant spectrum

Origins and challenges of viral dark matter

Mutational and fitness landscapes of an RNA virus revealed through population sequencing

Temporal dynamics of intrahost molecular evolution for a plant RNA virus

Deep sequencing of Influenza A virus from a human challenge study reveals a selective bottleneck and only limited intrahost genetic diversification

Timesampled population sequencing reveals the interplay of selection and genetic drift in experimental evolution of Potato virus Y

The epidemic of antibiotic-resistant infections: a call to action for the medical community from the Infectious Diseases Society of America

Selection of intermediate rates of increase in parasite-host systems

Evolution of infectious disease

On the evolutionary coexistence of parasite strains

The community ecology of pathogens: coinfection, coexistence and community composition

Speciation as pattern formation by competition in a smooth fitness landscape

Mutation and the evolution of virulence

Nine challenges for deterministic epidemic models

Effect of the latent reservoir on the evolution of HIV at the within-and between-host levels

A sampling theory for asymmetric communities

The evolution of parasitic fitness

Invasion fitness, inclusive fitness, and reproductive numbers in heterogeneous populations

Dynamics of HIV-1 recombination in its natural target cells

Multiple HIV-1 infection of cells and the evolutionary dynamics of cytotoxic T lymphocyte escape mutants

Viral fitness across a continuum from lysis to latency

Modelling the within-host growth of viral infections in insects

Biodiversity in model ecosystems, I: coexistence conditions for competing species

How direct competition shapes coexistence and vaccine effects in multi-strain pathogen systems

The evolution of virulence

The implications of coevolutionary dynamics to host-parasite interactions

Adaptation and the evolution of parasite virulence in a connected world

Superinfection and the evolution of parasite virulence

On the relationship between evolution of virulence and host demography

Erasing the world's slow stain: strategies to beat multidrug-resistant tuberculosis

Superinfections and adaptive dynamics of pathogen virulence revisited: a critical function analysis

The coexistence of competing parasites. I. The role of cross-species infection

Oscillatory coexistence in the chemostat: a codimension two unfolding

Remarks on competitive coexistence

Bifurcation approach to the predator-prey population models (Version of the computer book)

The impact of cross-immunity, mutation and stochastic extinction on pathogen diversity

Analysis of symmetries in models of multi-strain infections

The analysis of equilibrium in malaria

An HIV infection model based on a vectored immunoprophylaxis experiment

Dynamical behavior of epidemiological models with nonlinear incidence rates

The mathematics of infectious diseases

The study of plant disease epidemics

HIV-1 superinfection and viral diversity

Beneficial and perverse effects of isoniazid preventive therapy for latent tuberculosis infection in HIV -tuberculosis coinfected populations

Mathematical models of tuberculosis: accomplishments and future challenges

Construction of a mathematical model for tuberculosis transmission in highly endemic regions of the Asia-Pacific

Lysogenic and lytic viral production in marine microbial communities

Modeling the infection dynamics of bacteriophages in enteric Escherichia coli: estimating the contribution of transduction to antimicrobial gene spread

Bistability and switching in the lysis/lysogeny genetic regulatory network of bacteriophage l

Transmissible gastroenteritis coronavirus induces programmed cell death in infected cells through a caspase-dependent pathway

Virus dynamics: mathematical principles of immunology and virology

Population dynamics of immune responses to persistent viruses

Virus dynamics and drug therapy

Population biology of infectious diseases: part I

Population biology of infectious diseases: part II

The dynamics of HTLV-I and the CTL response

Viruses and apoptosis

MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs

Data from: Theoretical conditions for the coexistence of viral strains with differences in phenotypic traits: a bifurcation analysis

Elements of applied bifurcation theory

Fate of a mutant emerging at the initial stage of evolution

A degenerate Hopf bifurcation in retarded functional differential equations, and applications to endemic bubbles

Selforganization of matter and evolution of biological macromolecules

Quasispecies made simple

Theory of lethal mutagenesis for viruses

Practical computation of normal forms on center manifolds at degenerate Bogdanov -Takens bifurcations

Analysis of degenerate fold-Hopf bifurcation in a three-dimensional differential system

Degenerate with respect to parameters fold-Hopf bifurcations

Degenerate fold-Hopf bifurcations in a Rössler-type system

Acknowledgements. We thank Antoni Guillamon for interesting discussions about the zero-Hopf bifurcation. M.H., J.S. and L.A. acknowledge the hospitality of the Instituto de Biología Molecular y Celular de Plantas which led to fruitful meetings.

Competing interests. We declare we have no competing interests. Funding. The research leading to these results has received funding from 'la Caixa' Foundation. This work has been also Appendix A. Further study of degenerate Bogdanov -Takens and zero-Hopf bifurcationsNormal form theory provides a means of predicting the local dynamics that occur near a given bifurcation. This is possible because the local structure of a bifurcation is model independent and thus can be immediately determined from pre-existing analyses of simpler systems of equations. The simplest system of equations that completely captures the dynamics of a particular bifurcation is called the normal form. The conducted numerical studies on the system (2.4) revealed that the Bogdanov-Takens and zero-Hopf bifurcations are non-standard in the sense that the local dynamics do not agree with those predicted by their normal forms. Such a mismatch can occur when one of the terms in the normal form equations vanishes due to its coefficient being equal to zero, implying that an extended system of equations needs to be considered. We now examine the normal forms of the Bogdanov-Takens (BT) and zero-Hopf (ZH) bifurcations in order to rationalize the discrepancies between the two-dimensional bifurcation diagram shown in figure 4 and those predicted from normal form theory.The local bifurcation diagram near a BT point can be determined from the normal form equations given by (e.g. Kuznetsov [66, §8.4])