Theoretical conditions for the royalsocietypublishing.org/journal/rsos coexistence of viral strains with differences in Research phenotypic traits: a Cite this article: Nurtay A, Hennessy MG, bifurcation analysis Sardanyés J, Alsedà L, Elena SF. 2019 Theoretical conditions for the coexistence of viral strains with differences in phenotypic traits: a bifurcation Anel Nurtay1,2,3,4, Matthew G. Hennessy1,2, analysis. R. Soc. open sci. 6: 181179. Josep Sardanyés1,2, Lluı́s Alsedà3,1,2 http://dx.doi.org/10.1098/rsos.181179 and Santiago F. Elena4,5 1 Centre de Recerca Matemàtica, 2Barcelona Graduate School of Mathematics (BGSMath), and 3 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Campus de Bellaterra, Received: 23 July 2018 Edifici C, 08193 Bellaterra, Spain Accepted: 26 November 2018 4 Instituto de Biologı́a Integrativa de Sistemas, CSIC-Universitat de València, Parc Cientı́fic UV, Paterna, València 46980, Spain 5 Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA AN, 0000-0001-7107-1656; MGH, 0000-0002-5928-6256; Subject Category: SFE, 0000-0001-8249-5593 Mathematics We investigate the dynamics of a wild-type viral strain which generates mutant strains differing in phenotypic properties for Subject Areas: infectivity, virulence and mutation rates. We study, by means of mathematical modelling/health and disease and a mathematical model and bifurcation analysis, conditions under epidemiology/applied mathematics which the wild-type and mutant viruses, which compete for the same host cells, can coexist. The coexistence conditions are Keywords: formulated in terms of the basic reproductive numbers of the bifurcations, epidemiology, infection dynamics, strains, a maximum value of the mutation rate and the virulence of the pathogens. The analysis reveals that parameter space can mathematical biology, virus evolution 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 Authors for correspondence: transcritical bifurcations of periodic orbits. These results provide Anel Nurtay new insights into the conditions by which viral populations may contain multiple coexisting strains in a stable manner. e-mail: anurtay@crm.cat Matthew G. Hennessy e-mail: mhennessy@crm.cat 1. Introduction 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 & 2019 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. environmental changes such as new host species, the presence of antiviral drugs, new transmission 2 routes or to new vectors [1]. This tremendous evolvability, coupled with densely populated animal royalsocietypublishing.org/journal/rsos and plant susceptible hosts (in many cases lacking genetic variability for resistance to infection), are the reasons for the persistence and emergence of new viral diseases or the re-emergence of new strains with novel properties of already known diseases. The continuous emergence of new mutants leads to an overlap in existence of wild-type (hereafter referred to as wt) and mutant genotypes of the same virus within individual infected hosts [2,3]. This cloud of mutants is usually known as a viral quasi-species [1] and it constitutes the target of selection, instead of the individual viral genomes. This reservoir of coexisting genetic variants may lead to the emergence of new genotypes with different host ranges, pathologies and epidemiological properties that may result in outbreaks [4,5]. With the development of high-coverage, ultra-deep sequencing techniques, it is now possible to characterize in great detail virus genetic diversity along the course of infection of individual hosts, demonstrating the coexistence of multiple mutant genotypes within individual hosts, some even during long periods of time [6–8]. Furthermore, some of these studies have also shown that dynamics R. Soc. open sci. 6: 181179 are highly complex and do not only depend on the differences in replicative fitness among individual genotypes, but on other parameters such as the size and frequency of within-host bottlenecks, complementation of strains, fixation of additional mutations on the same genotype, epistasis, the availability of beneficial mutations (which indeed depends on the degree of adaptation to the host), the load of deleterious mutations or clonal interference among coexisting beneficial mutations [9–12]. What evolutionary mechanisms determine the long-term coexistence of different genetic variants and strains that, in principle, shall be competing for the same resources (e.g. target cells)? Understanding the evolutionary forces of such intraspecific competition or strain coexistence are essential for understanding the long-term fate and composition of viral populations and for a thoughtful design of more robust control strategies for known and future outbreaks [13]. Consequently, the coexistence of evolving pathogens has been the target of extensive research [14 –17]. In the mathematical theory of population genetics, mutation, which is understood as any change in the genome of an organism, is often modelled as a ‘flow’ between populations of initial wt individuals and emerging mutant individuals [18 –21]. In the present article, we introduce mutation into a mathematical model in a similar way although avoiding forceful restrictions put upon mutant strains and allowing the mutant virus to have similar characteristics to the wt virus. The difference between mutant and wt strains occurs when focusing the study on specific phenotypic characteristics. A number of studies employ fitness to investigate the survival of a population in dynamical systems that present competition or coexistence [22 –24]. As an accumulative property resulting from various phenotypic traits, fitness is a convenient measure of the overall success of the population. In the case of virus evolution, fitness can be considered proportional to the infection rate of the virus [25 –27], and measured by the infection rate when other parameters are fixed. It is known that infection rates differ between strains, a fact which has significant implications for the evolution of virulence and strain coexistence in nature [28]. Even the balance between genetic diversity and competition is believed to be achieved due to the possibility of coexistence among strains with differences in infection rates [29]. In other words, a direct competition for infecting available cells mediates the stable coexistence only when competitive abilities in viral clones satisfy certain pairwise asymmetries [30]. A secondary phenotypic characteristic that can differ between viral strains is their strategy for exploiting the host cell, i.e. their virulence [31]. The evolution of virulence has received great attention from theoreticians, particularly on the coevolution between resistance and virulence traits and their combined effect on host and virus dynamics [32,33]. However, understanding the evolution of virulence for coexisting viral strains still requires attention due to the complexity of the underlying evolutionary and dynamical processes, being inherently nonlinear. Many of the models brought forward to explain the evolution of virulence take into consideration the processes of coinfection and superinfection [34 –37], where the host or the host cell is infected simultaneously by more than one pathogen particle (coinfection) or sequentially by different pathogens (superinfection). Here, infection of a cell will be modelled with virions of a single strain and superinfection will be neglected. The benefit of this approach is that analytical insights into coexistence of viral strains can be obtained. The purpose of this paper is to illustrate, by means of a dynamical mathematical model, the conditions for coexistence of viral strains that considers both a wt viral strain and its mutants. We present analytical and numerical results focusing on the parameters related to the differential phenotypic traits of the wt and mutant strains. Conditions for coexistence and invasion have previously been studied using a mathematical model of one host shared between two competing parasites [38]. However, this model did not incorporate the mutation of parasites as a factor, and thus neglected the input of new strains into the infected cells with infected cells with 3 the wild-type strain the mutant strain royalsocietypublishing.org/journal/rsos yw ym within-cell viral replication and mutation lysis lysis R. Soc. open sci. 6: 181179 infection infection wild-type strain zw mutant strain proliferation uninfected zm host cells x Figure 1. 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 zw) and mutant (variable zm) 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 yw) and mutant (variable ym) 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. 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–43]. A classical approach used in epidemiology [44] to circumvent this difficulty is to introduce dimensionless parameter groups, such as the basic reproduction number R0, 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–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. 2. Mathematical model 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 (zw) and the mutant (zm) species; and two types of infected cells, one infected by the wt strain, yw, and another one infected by the mutant strains, ym (figure 1). The time evolution of the interacting populations is described with the following model: x_ ¼ bxð1  x=KÞ am zm x  aw zw x dx; (2:1a) |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |ffl{zffl} bounded replication infection decay of of healthy cells of healthy cells healthy cells 4 mutation from wt to mutant royalsocietypublishing.org/journal/rsos zfflfflffl}|fflfflffl{ y_ m ¼ am zm x þ myw gm ym , (2:1b) |fflfflffl{zfflfflffl} |fflfflfflffl{zfflfflfflffl} growth of mutant decay of mutant infected cells infected cells mutation from wt to mutant zfflfflffl}|fflfflffl{ y_ w ¼ aw zw x  myw gw yw , (2:1c) |fflfflffl{zfflfflffl} |fflfflfflffl{zfflfflfflffl} growth of wt decay of wt infected cells infected cells and R. Soc. open sci. 6: 181179 internalization of type i virions by infection zfflfflfflfflffl}|fflfflfflfflffl{ z_ i ¼ ki gi yi ni ai zi x zi zi , with i ¼ m, w: (2:1d) |fflffl{zfflffl} |fflffl{zfflffl} lytic release of degradation of type i virions type i virions 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 am and aw, respectively. Mathematical and statistical models of multiple infections, coinfection and superinfection have been studied in detail [49–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–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–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 yw infected by the wt strain can mutate at rate m into cells of population ym. 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 gm and gw, and as opposed to d þ Dgm and d þ Dgw. 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 km and kw 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 nm and nw 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’ 5 rates of the virions in the system are zw and zm for wt and mutant strains, respectively. royalsocietypublishing.org/journal/rsos Initial conditions for system (2.1) are non-negative values: x(0) ¼ X, ym (0) ¼ Ym , yw (0) ¼ Yw , zm (0) ¼ Zm , zw (0) ¼ Zw : 2.1. Non-dimensionalization 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 is ~xmax ¼ (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 km and the characteristic size of the infected cell population, ~xmax . The former is required due to differences in R. Soc. open sci. 6: 181179 sizes and measurement units of the viral loads and the cell populations in the system. We therefore non-dimensionalize the variables according to t ¼ (b  d)1t, x ¼ ~xmax x, yi ¼ ~xmax yi , zi ¼ ~xmax kmzi , (2:2) where bars are used to denote dimensionless quantities. We also introduce dimensionless parameters defined by km Kai ni kw m pi i ¼ a , ni ¼ , ¼ k , ¼ m , pi ¼ , (2:3) b km km bd bd where pi (recall that i ¼ m, w) stands for other ‘rate’ parameters, namely: gm, gw, zm and zw. 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 zi . 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 x_ ¼ x(1  x)  am zm x  aw zw x, (2:4a) y_ m ¼ am zm x þ myw  gm ym , (2:4b) y_ w ¼ aw zw x  myw  gw yw , (2:4c) z_ m ¼ gm ym  nm am zm x  zm zm (2:4d) and z_ w ¼ kgw yw  nw aw zw x  zw zw : (2:4e) 3. Results and discussion 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]. 3.1. Equilibrium points 6 The non-dimensional model given by equations (2.4) has four equilibria. For mathematical convenience, royalsocietypublishing.org/journal/rsos we define the population vector v(t) ¼ {x(t), ym (t), yw (t), zm (t), zw (t)}. The first equilibrium is the origin v0 :¼ {x ¼ 0, ym ¼ 0, yw ¼ 0, zm ¼ 0, zw ¼ 0}: (3:1) If stable, the trivial solution v0 corresponds to the extinction of all of the populations. The second equilibrium v1 :¼ {x ¼ 1, ym ¼ 0, yw ¼ 0, zm ¼ 0, zw ¼ 0}, (3:2) 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 ( zm ðam ð1  nm Þ  zm Þzm R. Soc. open sci. 6: 181179 v2 :¼ x ¼ , ym ¼ , yw ¼ 0, am ð1  nm Þ a2m ð1  nm Þ2 gm  am ð1  nm Þ  zm zm ¼ , zw ¼ 0 , (3:3) a2m ð1  nm Þ 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 v2 to be biologically meaningful, the condition nm , 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 nm , 1. Finally, the fourth equilibrium point is given by  z (g þ m) Bmzw (C þ am gw zw ) Bz (C  am mz ) v3 :¼ x ¼ w w , ym ¼ , yw ¼ w 2 2 w , Aa w A2 Ca2w gm A Ca w  Bmzw B(C  am mzw ) zm ¼ , zw ¼ , (3:4) ACaw ACa2w where 9 A ¼ gw k  (m þ gw )nw , > = B ¼ ((nw  k)aw þ zw )gw  m(aw nw þ zw ) (3:5) > ; and C ¼ am zw (gw (1  nm )  mnm )  aw zm (gw (k  nw )  mnw ): This equilibrium point involves, whenever stable, a state of coexistence for the wt and mutant strains. Similar to the wt-free state v2, the coexistence state v3 can only be biologically meaningful if the burst size of wt strain satisfies k . nw(1 þ m/gw). This inequality is a generalization of that derived for the mutant virus (nm , 1) which accounts for mutation. The singularity that occurs in the coexistence state v3 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 v2 whereby only the wt virus exists. 3.2. Linear stability analysis and bifurcations 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 2 3 am zm  aw zw  2x þ 1 0 0 aw x am x 6 aw zw gw  m 0 aw x 0 7 6 7 J¼6 6 am z m m g m 0 a m x 7: 7 4 aw nw zw kgw 0 aw nw x  zw 0 5 am nm zm 0 gm 0 am nm x  zm Straightforward calculations for the trivial solution v0 yield eigenvalues of the Jacobian given by L0 ¼ (1,  zw ,  zm ,  gm ,  gw  m)T : (3:6) Thus, for all parameter values, the trivial solution v0, 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 7 v1, they are sufficiently complicated that little insight is gained from analysing them directly. In order to royalsocietypublishing.org/journal/rsos conduct a stability analysis for v1, it is more useful to consider the characteristic polynomial for the eigenvalues l. The polynomial det (J(v1 )  lI) ¼ 0 can be regrouped into the form of three factors P1 (l)  P2 (l)  P3 (l) ¼ 0, where P1 (l) ¼ l þ 1, (3:7a) P2 (l) ¼ l2 þ (aw nw þ gw þ m þ zw )l  aw (gw (k  nw )  mnw ) þ zw (gw þ m) (3:7b) and P3 (l) ¼ l2 þ (am nm þ gm þ zm )l þ am gm (nm  1) þ gm zm : (3:7c) The first factor, P1, provides a constant eigenvalue, l ¼ 21, which reserves the possibility for stability of v1. Although the next two factors P2 and P3 do not provide obvious eigenvalues, they enable the identification of critical parameter groups for which v1 undergoes a bifurcation. R. Soc. open sci. 6: 181179 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 P2 and P3 given by (3.7b,c) reveals that v1 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 v1 exists for all parameter combinations, transcritical bifurcations can be straightforwardly detected by forcing l to be zero in P2 and P3. By solving P2(l ¼ 0) ¼ 0, we find Rw 0 ¼ 1, where   aw k Rw0 :¼  nw : (3:8) zw 1 þ mg1w Similarly, solving P3(l ¼ 0) ¼ 0 results in Rm 0 ¼ 1, where am Rm 0 :¼ (1  nm ): (3:9) zm It will be shown below that Rw m 0 and R0 are the basic reproductive numbers for the wt and mutant strains, w respectively. When R0 ¼ 1, the virus-free state v1 and the coexistence state v3 intersect. By expanding v3 around Rw w 0 ¼ 1, we find that it has negative components when R0 , 1 and thus lies outside of the biologically meaningful phase space. However, all of the components of v3 become positive when Rw m 0 . 1. Likewise, expanding the equilibrium v2 around R0 ¼ 1 shows that this condition corresponds to points where the virus-free state v1 and wt-free state v2 intersect. For Rm 0 , 1, some components of v2 are negative; for Rm 0 . 1, all components are positive. In the case when both Rw m 0 ¼ R0 ¼ 1 hold, there is a triple intersection of v1, v2 and v3. 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 v2 and v3 intersect when Rm w 0 ¼ R0 , that is, when   aw k am  nw ¼ (1  nm ): (3:10) zw 1 þ mg1w zm Thus, there are three curves of transcritical bifurcations defined by Rw m m w 0 ¼ 1, R0 ¼ 1 and R0 ¼ R0 , all of which simultaneously intersect at the DBT point. The dynamics of the system near the DBT bifurcation can be determined by expanding the equilibria and their eigenvalues around Rw m 0 ¼ 1 and R0 ¼ 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 Rw m 0 , 1 and R0 , 1, the virus- free state v1 is locally asymptotically stable and the wt-free state v2 and the coexistence state v3 are unstable, meaning that the virus-free state is achieved. For Rm w 0 . 1 and R0 , 1, the virus-free and wt- free state exchange stability, with v1 becoming unstable and v2 locally asymptotically stable, with v3 remaining unstable. Similarly, for Rw m 0 . 1 and R0 , 1, the virus-free and coexistence state exchange stability: v1 becomes unstable, v3 becomes locally asymptotically stable, and v2 remains unstable. The stability in the region Rm w 0 . 1 and R0 . 1 is governed by the transcritical bifurcation occurring along w 8 R0 m = II R0 royalsocietypublishing.org/journal/rsos Rm 0=1 Rm 0 DBT III Rw0 = 1 Rw0 Figure 2. Co-dimension two bifurcation diagram displaying the local stability near the degenerate Bogdanov– Takens (DBT) bifurcation in terms of the basic reproduction numbers defined in (3.8) and (3.9). Three regions with different qualitative R. Soc. open sci. 6: 181179 behaviours are identified. Region I: virus-free state without infected cells (v1 stable); Region II: wt-free state (v2 stable); and Region III: coexistence of all strains (v3 stable). the curve Rm w 0 ¼ R0 , which leads to an exchange of stability between wt-free and coexistence states, v2 and v3. The transcritical bifurcations that occur for Rm w w m m 0 ¼ 1 with R0 . 1, R0 ¼ 1 with R0 . 1 and R0 ¼ R0 w w with R0 , 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 v1, v2 and v3 are local attractors, respectively. Based on these results, we can interpret the quantities Rm w 0 and R0 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 Rm w 0 and R0 , 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 v2, we obtain the eigenvalues from the characteristic polynomial computed from det (J(v2 ) lI) ¼ 0. This equation can be factorized and rewritten as Q2 (l)  Q3 (l) ¼ 0, where Q2 and Q3 are quadratic and cubic polynomials in l, respectively. The exact form of Q2 is not required here. We write Q3(l) ¼ al3 þ bl2 þ cl þ d ¼ 0, where 9 a ¼ am (1  nm )2 , > > > > > = b ¼ (1  nm )(am gm (1  nm ) þ zm (am þ 1)), (3:11) c ¼ zm (n2m am þ (zm  am  gm )nm þ zm þ gm ) >> > > > ; and d ¼ (1  nm )zm gm (am (1  nm )  zm ): Following the analysis scheme discussed earlier, to obtain conditions for Hopf bifurcations of v2, we search for a purely imaginary pair of eigenvalues. Setting Q3(l ¼ +iv)¼ 0, we obtain the critical pffiffiffiffiffiffiffi condition for a Hopf bifurcation, ad ¼ bc, with v ¼ d=b 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, Q2(l), provides no possibility for a Hopf bifurcation due to a strictly positive linear coefficient in the quadratic polynomial. Transcritical bifurcations of v2 may occur if Q2(l ¼ 0) ¼ 0 or Q3(l ¼ 0) ¼ 0. First, setting Q3(l ¼ 0) ¼ 0 yields d ¼ 0, which is equivalent to Rm0 ¼ 1 and corresponds to an intersection of v1 and v2. Second, from Q2(l ¼ 0) ¼ 0, we obtain an expression which matches with (3.10), that is, the parameter combination yielding an intersection of v2 and v3. As will be shown in §3.3, it is possible to simultaneously satisfy Q3(l ¼ +iv)¼ 0 and Q2(l ¼ 0) ¼ 0, implying that the Jacobian matrix at v2 ¼ v3 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 (a) (b) 9 1.0 T 1.0 T royalsocietypublishing.org/journal/rsos 0.8 0.8 0.6 0.6 T x 0.4 0.4 0.2 H 0.2 H 0 0 0 1 2 3 4 5 0 1 2 3 4 5 am am Figure 3. One-dimensional bifurcation diagrams showing the dependence of the equilibrium population of healthy cells x on the R. Soc. open sci. 6: 181179 infection rate of mutant-type virus am. We set (a) aw ¼ 0.5 and Rw0 ¼ 0:48 and (b) aw ¼ 2 and Rw0 ¼ 1:93. Solid and dashed curves illustrate, respectively, stable and unstable equilibria. Filled markers stand for maxima and minima of a stable periodic orbit (PO). The letters T and H denote points of transcritical and Hopf bifurcations of equilibria, respectively. Blue denotes biologically relevant solutions ( positive variables) while orange describes solutions with negative components. The minimum values of the POs are above zero as am is increased after the Hopf bifurcation. All the other parameter values are given in (3.12). 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 v2, v3, 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. 3.3. Phenotypic differences in infection rates We first examine the influence of the infection rates am and aw on the dynamics. The values of the other parameters are set to m ¼ 0:1, k ¼ 1, gm ¼ gw ¼ 0:25, nm ¼ nw ¼ 0:5, zm ¼ zw ¼ 0:2222: (3:12) We begin our investigation by constructing one-dimensional bifurcation diagrams using am as the bifurcation parameter and fixing aw. We first consider two values of the wt strain infections rate given by aw ¼ 0.5 and aw ¼ 2, corresponding to basic reproduction numbers Rw w 0 ¼ 0:48 and R0 ¼ 1:93, respectively. These values of aw therefore lie on opposite sides of the curve of transcritical bifurcations involving v1 and v3 defined by Rw 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 aw ¼ 0.5, two critical values of the mutant-virus infection rate am that lead to qualitative changes in the dynamics of the system. The point am ¼ 0.44 corresponds to Rm 0 ¼ 1, and marks the position of the transcritical bifurcation involving the virus-free state v1 and the wt-free state v2. As predicted by the linear stability analysis in §3.2, the wt-free state undergoes a supercritical Hopf bifurcation at am ¼ 1.55. The two critical values of am divide the parameter space into three distinct regions. For am , 0.44, the virus-free state v1 is locally asymptotically stable, corresponding to Region I in figure 2. Similarly, for 0.44 , am , 1.55, the wt-free state v2 is stable, corresponding to Region II of figure 2. Finally, for am . 1.55, there is a new region, termed Region IV, where stable POs exist about the equilibrium v2, which is unstable. These POs describe oscillatory populations of uninfected cells and mutant-type virus with extinct populations of the wt virus. (a) 4 (b) 4 10 royalsocietypublishing.org/journal/rsos 3 3 IV am 2 2 DZH DZH 1 1 II V III DBT I DBT 0 2 4 6 8 10 0 2 4 6 8 10 aw aw Figure 4. (a) Two-dimensional bifurcation diagram, illustrating the effect of infection rates on the dynamics of the system. Blue R. Soc. open sci. 6: 181179 denotes biologically meaningful solutions (positive values for all components). Orange describes solutions with negative components, which are not biologically meaningful. The abbreviations DBT and DZH mark the locations of degenerate Bogdanov– Takens and degenerate zero-Hopf bifurcations, respectively. (b) Two-dimensional stability diagram with the biologically meaningless bifurcations removed, showing five regions of distinct dynamics. Solid, dashed and starred lines represent curves of Hopf bifurcations, transcritical bifurcations of fixed points and transcritical bifurcations of periodic orbits (POs). Region I: stability of virus- free solution; Region II: stability of wt-free state; Region III: stability of coexistence state; Region IV: stable PO around the wt-free state; Region V: stable PO around the coexistence state. In both plots, the parameter values are given by (3.12). When the infection rate of wt virus is increased to aw ¼ 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 am , 0.86, the coexistence state v3 is stable whereas the virus-free state v1 is unstable. This implies that Region I (v1 stable) is replaced by Region III (v3 stable). Although v1 and v2 undergo a transcritical bifurcation at am ¼ 0.44, there is no exchange of stability. Instead, it is v2 and v3 which exchange stability through a transcritical bifurcation at am ¼ 0.86, corresponding to the case of equal basic reproduction numbers, Rm w 0 ¼ R0 ¼ 1:93. Both Regions II and IV persist for aw ¼ 2, with the supercritical Hopf bifurcation occurring at am ¼ 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 aw as the bifurcation parameter for fixed values of am that lie on opposite sides of Rm 0 ¼ 1. However, the resulting diagrams are qualitatively similar to those in figure 3, with the exception that the wt-free state v2 is switched with the coexistence state v3 and vice versa. The supercritical Hopf bifurcation now occurs from v3 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 am and aw 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 aw ¼ 1.04, am ¼ 0.44 and aw ¼ 3.62, am ¼ 1.55, respectively. The DBT bifurcation lies at the intersection of three curves of transcritical bifurcations satisfying Rw m 0 ¼ 1, R0 ¼ 1 m w and R0 ¼ R , which may be written, respectively, in terms of the infection rates as zw zm a w zw 1  nm aw ¼ 1 , am ¼ , ¼ 1 : (3:13) k(1 þ mg1 w )  nw 1  nm am zm k(1 þ mg1 w )  nw 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 v2 and v3 and the curve of transcritical bifurcations involving v2 and v3. As predicted from linear stability analysis, the Hopf bifurcation curve associated with the wt-free state v2 does not depend on the infection rate of the wt virus aw and thus appears as a straight line (a) (b) 11 1.0 T 1.0 T royalsocietypublishing.org/journal/rsos 0.8 TPO 0.8 TPO 0.6 0.6 x 0.4 0.4 H H T 0.2 T 0.2 0 0 0 1 2 3 4 5 0 2 4 6 8 10 am aw Figure 5. One-dimensional bifurcation diagrams for fixed values of (a) aw ¼ 8 and (b) am ¼ 2, with only the x-component of R. Soc. open sci. 6: 181179 equilibria ( population of healthy cells) plotted for clarity. Dashed curves illustrate unstable equilibria. Filled and open markers stand for stable and unstable periodic orbits (POs), respectively. Only the maxima and minima of each PO is shown. Blue denotes biologically meaningful solutions and POs (strictly positive values for all components). Orange describes biologically meaningless solutions with negative components. Parameter values are given in (3.12). The letter T represents a transcritical bifurcation of equilibria, H a Hopf bifurcation, and TPO a transcritical bifurcation of POs. given by am ¼ 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 aw ¼ 8 and treating am as a bifurcation parameter and then setting am ¼ 2 and treating aw as the bifurcation parameter. The resulting diagrams are shown in figure 5, with open circles denoting unstable POs. Figure 5a shows that as am is increased from zero when aw ¼ 8, the wt-free state v2 first undergoes a transcritical bifurcation with the virus-free state v1, then a subcritical Hopf bifurcation, and finally a transcritical bifurcation with the coexistent state v3. 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 am, the system is in Region V, and there is a stable PO around the coexistence state v3. This PO is created by the Hopf bifurcation that defines the boundary between Regions III and V shown in figure 4b. As am is increased, the unstable POs around the wt-free state v2 that emerge from the subcritical Hopf bifurcation grow in size. Eventually, this unstable PO collides and exchanges stability with the stable PO around v3 through a TPO bifurcation at am ¼ 2.3972. During the TPO bifurcation, the stable and biologically meaningful PO around v3 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 am ¼ 2 and shown in figure 5b shares the same qualitative features as figure 5a, with v2 and its POs exchanging roles with v3 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 12 on studying the behaviour of a newly emerging mutant in a population of a few strains of bacteria that royalsocietypublishing.org/journal/rsos 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. 3.4. Phenotypic differences in virulence R. Soc. open sci. 6: 181179 We now study the effect of the strains’ virulence ( parametrized by gw and gm) 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 gm ¼ gw ¼ 0.25. The other parameters are fixed to m ¼ 0:1, k ¼ 1, nm ¼ nw ¼ 0:5, zm ¼ zw ¼ 0:2222: (3:14) 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 Rm 0 given by (3.9), we see that the basic reproduction number of the mutant is independent of the virulence. Thus, changes in gm or gw cannot increase Rm 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, gw ! 1, we find that Rw w 0 ! aw (k  nw )=zw . Thus, if the infection rate aw is so small that the limit of R0 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 aw is sufficiently large and the basic reproduction number increases beyond one, then a change in dynamics will be observed as gw is increased. By solving Rw 0 ¼ 1, a critical value of the virulence of the wt virus is obtained m gcrit w ¼ : (3:15) k(zw a1 w þ nw )1  1 The above equation is the condition for a transcritical bifurcation between the virus-free state v1 and the coexistence state v3 and defines the boundary between Regions I and III. Thus, only for values of gw . gcrit w will the virus persist in the system for this choice of infection rates. The stability diagram for aw ¼ 0.5 and am ¼ 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 aw ¼ 3, am ¼ 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 am ¼ 1 along with the values of the parameters in (3.14) leads to Rm 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 v2 and v3. 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. (a) (b) (c) 13 4 4 4 royalsocietypublishing.org/journal/rsos 3 3 3 II III V III II gm 2 I III 2 2 1 1 1 IV 0 1 2 3 4 0 0.5 1.0 1.5 2.0 0 1 2 3 4 gw gw gw Figure 6. Two-dimensional bifurcation diagrams for gw versus gm at different sets of ai, with i ¼ m, w. Parameter values are given by (3.14) and (a) aw ¼ 0.5, am ¼ 0.1; (b) aw ¼ 3, am ¼ 1; and (c) aw ¼ 0.5, am ¼ 2. Dashed and solid lines correspond to transcritical and Hopf bifurcations, respectively. Region I: stability of trivial solution; Region II: stability of wt-free state; Region III: stability of coexistence state; Region IV: stable periodic orbit (PO) around the wt-free state; Region V: stable R. Soc. open sci. 6: 181179 PO around the coexistence state. 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 v2. 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 gm independent of gw and aw. Solving ad ¼ bc to find gcritm (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   zm 1 am . 1þ , (3:16) 1  nm nm equivalent to am . 1.3332 for the chosen parameters, does the Hopf curve of v2 appear on a gw versus gm stability diagram (figure 6c). Greater values of am 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 gm 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 gm either. Importantly, the virulence of the wt strain, gw, strongly controls the dynamics and the persistence of the wt strain. As can be seen for small values of gw, 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 gw 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 gw 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. 3.5. The effect of heterogeneity in mutation rates 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, Rw 0 . 1, which can be interpreted as a condition on the mutation rate:   k m , gw 1 : (3:17) zw aw 1 þ nw 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 aw, 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 aw is shown in figure 7a for the (a) (b) 14 0.4 4 royalsocietypublishing.org/journal/rsos 0.3 3 I IV m 0.2 am 2 III 0.1 1 II V I 0 0 0 2 4 6 8 10 0 2 4 6 8 10 aw aw R. Soc. open sci. 6: 181179 (c) (d) 1.0 0.6 zm ym 0.4 y 0.5 z 0.2 yw zw 0 0 0 0.05 0.10 0.15 0.20 0.25 0.30 0 0.05 0.10 0.15 0.20 0.25 0.30 m m Figure 7. (a,b) Two-dimensional stability diagrams with a pair of bifurcation curves: transcritical bifurcations (dashed) and Hopf bifurcations (solid). Parameters, except m, are given by (3.12). (a) We use am ¼ 0.1 and consider stability diagram over parameters aw versus m. (b) We set m ¼ 0.5 . mc and rebuild the aw versus am stability diagram. As before, Region I: stability of virus-free solution; Region II: stability of wt-free state; Region III: stability of coexistence state; Region IV: stable PO about wt-free state; Region V: stable PO about coexistence state. In the lower panels we display one-dimensional bifurcation diagrams for (c) infected cell and (d) viral strain populations versus m; note, that only stable equilibria (solid curves) and stable POs (filled markers) are shown. We set aw ¼ 3.0 and am ¼ 1.0 and track the outcompetition of the wt strains by the mutant ones at increasing m, resulting in a type of error-threshold found after the oscillatory and static coexistence scenarios for m satisfying (3.17). case of am ¼ 0.1 and using parameter values in (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 aw increases, the boundary between Regions I and III reaches a horizontal asymptote given by   k mc ¼ gw 1 : (3:18) nw For mutation rates that satisfy m , mc, Regions III and V can be entered from Region I by increases in the infection rate, promoting coexistence. However, for m . mc, 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 . mc. 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 aw ¼ 3.0 and am ¼ 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 (a) (b) 15 4 10 royalsocietypublishing.org/journal/rsos 8 3 6 am 2 4 1 2 0 0 0 2 4 6 8 10 0 2 4 6 8 10 aw aw R. Soc. open sci. 6: 181179 Figure 8. Co-dimension two bifurcation diagrams in terms of the wt and mutant virus infection rates aw and am when backwards mutation is possible. The mutation rates are (a) mm ¼ 1023 and mw ¼ 1021 and (b) mm ¼ 1021 and mw ¼ 1021. All other parameters are given by (3.12). Dashed and solid lines correspond to transcritical and Hopf bifurcations, respectively. Blue denotes biologically feasible equilibria undergoing the bifurcations (all components non-negative), while colour orange denotes bifurcations of equilibria with negative components. The non-patterned, polka-dot and checkerboard regions correspond to virus-free, stationary coexisting, and oscillatory coexisting states, respectively. 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 y_ m ¼ am zm x þ mw yw  mm ym  gm ym (3:19a) and y_ w ¼ aw zw x  mw yw þ mm ym  gw yw : (3:19b) 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 mw ¼ 0.1 and we consider two values of the mutation rate of the mutant virus. First, we take mm ¼ 1023  mw, 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 mm ¼ 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 Rw m 0 ¼ 1 and R0 ¼ 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 v2 and v3 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. 4. Conclusion 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 16 parameters space, we identified five different regions, each characterized by distinct dynamics: royalsocietypublishing.org/journal/rsos (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 R. Soc. open sci. 6: 181179 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]. Authors’ contributions. A.N. and M.H. carried out the analytical and numerical computations. A.N., M.H., L.A., J.S. and S.F.E. conceived the mathematical model. A.N. drafted the initial manuscript. All authors revised the article and gave final approval for publication. 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 partially funded by the ‘Marı́a de Maeztu’ Programme for Units of Excellence in R&D (MDM-2014-0445), as well as from projects MTM2014-52209-C2-1-P and MTM2017-86795-C3-1-P from the Spanish MINECO, and from the CERCA Programme of the Generalitat de Catalunya. M.H. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 707658. Work in València was supported by Spain’s Agencia Estatal de Investigación - FEDER grant no. BFU2015-65037-P to S.F.E. J.S. has been also funded by a ‘Ramón y Cajal’ Fellowship (RYC-2017-22243). 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. Appendix A. Further study of degenerate Bogdanov – Takens and zero- Hopf bifurcations Normal 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]) j_ 0 ¼ j1 (A 1a) and 17 j_ 1 ¼ b1 þ b2 j0 þ a2 j20 þ b2 j0 j1 , (A 1b) royalsocietypublishing.org/journal/rsos provided that the coefficients a2 and b2 are not equal to zero, a2b2 = 0. The quantities b0 and b1 in equations (A 1) play the role of bifurcation parameters, with b0 ¼ b1 ¼ 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 a2 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 a2 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 R. Soc. open sci. 6: 181179 at the DBT point b0 ¼ b1 ¼ 0 if a2 = 0; however, the model under study has three. Thus, a2 ¼ 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] 1 1 v_ ¼ g(b) þ G200 (b)v2 þ G011 (b)jwj2 þ G300 (b)v3 þ G111 (b)vjwj2 (A 2a) 2 6 and 1 1 w_ ¼ L(b)w þ H110 (b)vw þ H210 (b)v2 w H021 (b)wjwj2 , (A 2b) 2 2 where b ¼ (b1, b2) 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 Hijk and Gijk. 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 G011(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 G011(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. References 1. Andino R, Domingo E. 2015 Viral quasispecies. (HBs Ag) and anti-HBs antibodies in chronic evolution of species jump. Trends Ecol. Evol. 20, Virology 479 – 480, 46 –51. (doi:10.1016/j.virol. hepatitis B virus carriers: influence of “a” 238–244. (doi:10.1016/j.tree.2005.02.009) 2015.03.022) determinant variants. J. Virol. 80, 2968 –2975. 6. Beerenwinkel B, Zagordi O. 2011 Ultra-deep 2. Tu H, Xiong SD, Trepo C, Wen YM. 1997 (doi:10.1128/JVI.80.6.2968-2975.2006) sequencing for the analysis of viral populations. Frequency of hepatitis B virus e-minus mutants 4. Elena S, Fraile A, Garcı́a-Arenal F. 2014 Evolution Curr. Opin. Virol. 1, 413– 418. (doi:10.1016/j. varies among patients from different areas of and emergence of plant viruses. Adv. Virus Res. coviro.2011.07.008) China. J. Med. Virol. 51, 85– 89. (doi:10.1002/ 88, 161– 191. (doi:10.1016/B978-0-12-800098- 7. Isakov O, Borderı́a A, Golan D, Hamenahem A, (ISSN)1096-9071) 4.00003-9) Celniker G, Yoff L, Blanc H, Vignuzzi M, 3. Lada O, Benhamou Y, Poynard T, Thibault V. 5. Woolhouse M, Haydon D, Antia R. 2005 Shomron N. 2015 Deep sequencing analysis of 2006 Coexistence of hepatitis B surface antigen Emerging pathogens: the epidemiology and viral infection and evolution allows rapid and detailed characterization of viral mutant 23. Nelson R. 2012 The evolution of parasitic J. Appl. Math. 43, 1005 –1018. (doi:10.1137/ 18 spectrum. Bioinformatics 31, 2141 –2150. fitness. In Plant disease: an advanced treatise, 0143066) (doi:10.1093/bioinformatics/btv101) ch. 2, vol. 4, (eds JG Horsfall, EB Cowling), 40. Muratori S, Rinaldi S. 1989 Remarks on royalsocietypublishing.org/journal/rsos 8. Krishnamurthy S, Wang D. 2017 Origins and pp. 23–46. New York, NY: Academic Press. competitive coexistence. SIAM J. Appl. Math. 49, challenges of viral dark matter. Virus Res. 239, 24. Lehmann L, Mullon C, Akcay E, Cleve J. 2016 1462 –1472. (doi:10.1137/0149088) 136–142. (doi:10.1016/j.virusres.2017.02.002) Invasion fitness, inclusive fitness, and 41. Bazykin A, Berezovskaya F, Zudin S. 1993 9. Acevedo A, Brosdsky L, Andino R. 2014 reproductive numbers in heterogeneous Bifurcation approach to the predator-prey Mutational and fitness landscapes of an RNA populations. Evolution 70, 1689 – 1702. (doi:10. population models (Version of the computer virus revealed through population sequencing. 1111/evo.2016.70.issue-8) book). Technical report International Centre for Nature 505, 686–690. (doi:10.1038/ 25. Levy DN, Aldrovandi GM, Kutsch O, Shaw GM. Theoretical Physics. nature12861) 2004 Dynamics of HIV-1 recombination in its 42. Abu-Raddad LJ, Ferguson NM. 2004 The impact 10. Cuevas J, Willemsen A, Hillung J, Zwart M, natural target cells. Proc. Natl. Acad. Sci. of cross-immunity, mutation and stochastic Elena S. 2015 Temporal dynamics of intrahost USA 101, 4204– 4209. (doi:10.1073/pnas. extinction on pathogen diversity. Proc. R. Soc. molecular evolution for a plant RNA virus. Mol. 0306764101) Lond. B 271, 2431 – 2438. (doi:10.1098/rspb. Biol. Evol. 32, 1132 – 1147. (doi:10.1093/ 26. Wodarz D, Levy DN. 2009 Multiple HIV-1 2004.2877) molbev/msv028) infection of cells and the evolutionary dynamics 43. Blyuss KB. 2014 Analysis of symmetries in models 11. Leonard A et al. 2016 Deep sequencing of of cytotoxic T lymphocyte escape mutants. of multi-strain infections. J. Math. Biol. 69, R. Soc. open sci. 6: 181179 Influenza A virus from a human challenge study Evolution 63, 2326 – 2339. (doi:10.1111/evo. 1431–1459. (doi:10.1007/s00285-013-0740-0) reveals a selective bottleneck and only limited 2009.63.issue-9) 44. Macdonald G. 1952 The analysis of equilibrium intrahost genetic diversification. J. Virol. 90, 27. Weitz JS, Li G, Gulbudak H, Cortez MH, Whitaker in malaria. Trop. Dis. Bulletin 49, 813–829. 11 247–11 258. (doi:10.1128/JVI.01657-16) RJ. 2018 Viral fitness across a continuum from 45. Wang X, Wang W. 2012 An HIV infection model 12. Kutnjak D, Elena S, Ravnikar M. 2017 Time- lysis to latency. bioRxiv p. 296897. based on a vectored immunoprophylaxis sampled population sequencing reveals the 28. White SM, Burden JP, Maini PK, Hails RS. 2012 experiment. J. Theor. Biol. 313, 127– 135. interplay of selection and genetic drift in Modelling the within-host growth of viral (doi:10.1016/j.jtbi.2012.08.023) experimental evolution of Potato virus Y. J. Virol. infections in insects. J. Theor. Biol. 312, 34– 43. 46. Liu Wm, Hethcote HW, Levin SA. 1987 91, e00690-1. (doi:10.1128/JVI.00690-17) (doi:10.1016/j.jtbi.2012.07.022) Dynamical behavior of epidemiological models 13. Spellberg B, Guidos R, Gilbert D, Bradley J, 29. Bastolla U, Lässig M, Manrubia SC, Valleriani A. with nonlinear incidence rates. J. Math. Biol. 25, Boucher HW, Scheld WM, Bartlett JG, Edwards J 2005 Biodiversity in model ecosystems, I: 359–380. (doi:10.1007/BF00277162) Jr, of America IDS. 2008 The epidemic of coexistence conditions for competing species. 47. Hethcote HW. 2000 The mathematics of antibiotic-resistant infections: a call to action for J. Theor. Biol. 235, 521–530. (doi:10.1016/j. infectious diseases. SIAM Rev. 42, 599–653. the medical community from the Infectious jtbi.2005.02.005) (doi:10.1137/S0036144500371907) Diseases Society of America. Clinical Infect. Dis. 30. Gjini E, Valente C, Sa-Leao R, Gomes MGM. 48. Madden LV, Hughes G, Bosch F. 2007 The study 46, 155– 164. (doi:10.1086/524891) 2016 How direct competition shapes coexistence of plant disease epidemics. St. Paul, MN: 14. Levin S, Pimentel D. 1981 Selection of and vaccine effects in multi-strain pathogen American Phytopathological Society (APS Press). intermediate rates of increase in parasite-host systems. J. Theor. Biol. 388, 50 –60. (doi:10. 49. Gross KL, Porco TC, Grant RM. 2004 HIV-1 systems. Am. Nat. 117, 308–315. (doi:10.1086/ 1016/j.jtbi.2015.09.031) superinfection and viral diversity. AIDS 18, 283708) 31. Read AF. 1994 The evolution of virulence. Trends 1513 –1520. (doi:10.1097/01.aids.0000131361. 15. Ewald PW. 1994 Evolution of infectious disease. Microbiol. 2, 73 –76. (doi:10.1016/0966- 75328.47) Oxford, UK: Oxford University Press. 842X(94)90537-1) 50. Cohen T, Lipsitch M, Walensky RP, Murray M. 16. Pugliese A. 2002 On the evolutionary 32. Best A, White A, Boots M. 2009 The implications 2006 Beneficial and perverse effects of isoniazid coexistence of parasite strains. Math. Biosci. of coevolutionary dynamics to host-parasite preventive therapy for latent tuberculosis 177, 355– 375. (doi:10.1016/S0025-5564(02) interactions. Am. Nat. 173, 779–791. (doi:10. infection in HIV –tuberculosis coinfected 00083-4) 1086/598494) populations. Proc. Natl. Acad. Sci. USA 103, 17. Seabloom EW, Borer ET, Gross K, Kendig AE, 33. Wild G, Gardner A, West SA. 2009 Adaptation 7042 –7047. (doi:10.1073/pnas.0600349103) Lacroix C, Mitchell CE, Mordecai EA, Power AG. and the evolution of parasite virulence in a 51. Colijn C, Cohen T, Murray M. 2007 Mathematical 2015 The community ecology of pathogens: connected world. Nature 459, 983– 986. models of tuberculosis: accomplishments and coinfection, coexistence and community (doi:10.1038/nature08071) future challenges. In BIOMAT 2006 (eds RP composition. Ecol. Lett. 18, 401– 415. (doi:10. 34. Nowak MA, May RM. 1994 Superinfection Mondaini, R Dilão), pp. 123 –148. Singapore: 1111/ele.2015.18.issue-4) and the evolution of parasite virulence. World Scientific. 18. Bagnoli F, Bezzi M. 1997 Speciation as pattern Proc. R. Soc. Lond. B 255, 81 –89. (doi:10.1098/ 52. Trauer JM, Denholm JT, McBryde ES. 2014 formation by competition in a smooth fitness rspb.1994.0012) Construction of a mathematical model for landscape. Phys. Rev. Lett. 79, 3302 –3305. 35. Castillo-Chavez C, Velasco-Hernandez JX. 1998 tuberculosis transmission in highly endemic (doi:10.1103/PhysRevLett.79.3302) On the relationship between evolution of regions of the Asia-Pacific. J. Theor. Biol. 358, 19. Bonhoeffer S, Nowak MA. 1994 Mutation and virulence and host demography. J. Theor. Biol. 74– 84. (doi:10.1016/j.jtbi.2014.05.023) the evolution of virulence. Proc. R. Soc. Lond. B 192, 437– 444. (doi:10.1006/jtbi.1998.0661) 53. Wilson WH, Mann NH. 1997 Lysogenic and lytic 258, 133– 140. (doi:10.1098/rspb.1994.0153) 36. Dye C, Williams BG, Espinal MA, Raviglione MC. viral production in marine microbial 20. Roberts M, Andreasen V, Lloyd A, Pellis L. 2015 2002 Erasing the world’s slow stain: strategies to communities. Aquat. Microb. Ecol. 13, 95 –100. Nine challenges for deterministic epidemic beat multidrug-resistant tuberculosis. Science (doi:10.3354/ame013095) models. Epidemics 10, 49 –53. (doi:10.1016/j. 295, 2042–2046. (doi:10.1126/science.1063814) 54. Volkova VV, Lu Z, Besser T, Gröhn YT. 2014 epidem.2014.09.006) 37. Boldin B, Geritz SA, Kisdi É. 2009 Modeling the infection dynamics of 21. Doekes HM, Fraser C, Lythgoe KA. 2017 Effect of Superinfections and adaptive dynamics of bacteriophages in enteric Escherichia coli: the latent reservoir on the evolution of HIV at pathogen virulence revisited: a critical function estimating the contribution of transduction to the within-and between-host levels. PLoS analysis. Evol. Ecol. Res. 11, 153–175. antimicrobial gene spread. Appl. Environ. Comput. Biol. 13, e1005228. (doi:10.1371/ 38. Hochberg ME, Holt RD. 1990 The coexistence of Microbiol. 80, 4350 – 4362. (doi:10.1128/AEM. journal.pcbi.1005228) competing parasites. I. The role of cross-species 00446-14) 22. Noble AE, Temme NM, Fagan WF, Keitt TH. infection. Am. Nat. 136, 517– 541. (doi:10. 55. Tian T, Burrage K. 2004 Bistability and switching 2011 A sampling theory for asymmetric 1086/285111) in the lysis/lysogeny genetic regulatory network communities. J. Theor. Bio. 273, 1 – 14. (doi:10. 39. Keener JP. 1983 Oscillatory coexistence in the of bacteriophage l. J. Theor. Biol. 227, 1016/j.jtbi.2010.12.021) chemostat: a codimension two unfolding. SIAM 229–237. (doi:10.1016/j.jtbi.2003.11.003) 56. Eleouet JF, Chilmonczyk S, Besnardeau L, Laude 63. Roulston A, Marcellus RC, Branton PE. 1999 69. Eigen M. 1971 Selforganization of matter and 19 H. 1998 Transmissible gastroenteritis coronavirus Viruses and apoptosis. Annu. Rev. Microbiol. evolution of biological macromolecules. induces programmed cell death in infected cells 53, 577– 628. (doi:10.1146/annurev.micro. Naturwissenschaften 58, 465– 523. (doi:10. royalsocietypublishing.org/journal/rsos through a caspase-dependent pathway. J. Virol. 53.1.577) 1007/BF00623322) 72, 4918 –4924. 64. Dhooge A, Govaerts W, Kuznetsov YA. 2003 70. Bull JJ, Meyers LA, Lachmann M. 2005 57. Nowak M, May RM. 2000 Virus dynamics: MATCONT: a MATLAB package for numerical Quasispecies made simple. PLoS Comput. Biol. 1, mathematical principles of immunology and bifurcation analysis of ODEs. ACM Trans. Math. e61. (doi:10.1371/journal.pcbi.0010061) virology. Oxford, UK: Oxford University Press. Software 29, 141 –164. (doi:10.1145/779359) 71. Bull J, Sanjuán R, Wilke C. 2007 Theory of lethal 58. Nowak MA, Bangham CR. 1996 Population 65. Nurtay A, Hennessy MG, Sardanyés J, Alsedà L, mutagenesis for viruses. J. Virol. 81, dynamics of immune responses to persistent Elena SF. 2018 Data from: Theoretical conditions 2930 –2939. (doi:10.1128/JVI.01624-06) viruses. Science 272, 74 –79. (doi:10.1126/ for the coexistence of viral strains with 72. Kuznetsov YA. 2005 Practical computation of science.272.5258.74) differences in phenotypic traits: a bifurcation normal forms on center manifolds at 59. Bonhoeffer S, May RM, Shaw GM, Nowak MA. analysis. Dryad Digital Repository. (doi:10.5061/ degenerate Bogdanov – Takens bifurcations. 1997 Virus dynamics and drug therapy. Proc. dryad.56g7v08) Int. J. Bifurcation Chaos 15, 3535 –3546. Natl. Acad. Sci. USA 94, 6971 –6976. (doi:10. 66. Kuznetsov YA. 2013 Elements of applied (doi:10.1142/S0218127405014209) 1073/pnas.94.13.6971) bifurcation theory, vol. 112. Berlin, Germany: 73. Tigan G. 2017 Analysis of degenerate fold-Hopf 60. Anderson RM, May RM. 1979 Population Springer. bifurcation in a three-dimensional differential R. Soc. open sci. 6: 181179 biology of infectious diseases: part I. Nature 67. Xu WZ, Kashiwagi A, Yomo T, Urabe I. 1996 Fate system. Qual. Theory Dyn. Syst. 17, 387 –402. 280, 361– 367. (doi:10.1038/280361a0) of a mutant emerging at the initial stage of (doi:10.1007/s12346-017-0241-4) 61. May RM, Anderson RM. 1979 Population evolution. Res. Popul. Ecol. 38, 231– 237. 74. Tigan G. 2017 Degenerate with respect to biology of infectious diseases: part II. Nature (doi:10.1007/BF02515732) parameters fold-Hopf bifurcations. Discrete Cont. 280, 455– 461. (doi:10.1038/280455a0) 68. LeBlanc V. 2016 A degenerate Hopf bifurcation Dyn. A 37, 2115 –2140. (doi:10.3934/dcdsa) 62. Wodarz D, Nowak MA, Bangham CR. 1999 The in retarded functional differential equations, 75. Tigan G, Llibre J, Ciurdariu L. 2017 Degenerate dynamics of HTLV-I and the CTL response. and applications to endemic bubbles. fold-Hopf bifurcations in a Rössler-type system. Immunol. Today 20, 220 –227. (doi:10.1016/ J. Nonlinear Sci. 26, 1 – 25. (doi:10.1007/ Int. J. Bifurcation Chaos 27, 1750068. (doi:10. S0167-5699(99)01446-2) s00332-015-9266-5) 1142/S0218127417500687)