Subthreshold Coexistence of Strains : The Impact of Vaccination and Mutation

We consider a model for a disease with two competing strains and vaccination. The vaccine provides complete protection against one of the strains (strain two) but only partial against the other (strain one). The partial protection leads to existence of subthreshold equilibria of strain one. If the first strain mutates into the second, there are subthreshold coexistence equilibria when both vaccine-dependent reproduction numbers are below one. Thus, a vaccine that is specific towards the second strain, and in absence of other strains should be able to eliminate it by reducing its reproduction number below one, cannot do so because it provides only partial protection to another strain that mutates into the second strain.


Introduction.
Many diseases are caused by more than one antigenically different variant of the causative agent.Such variants are referred to as different strains of a microorganism.The number of strains that give rise to the same disease depends on the mutability of the parasite.For some highly mutable viruses, like HCV (hepatitis C virus) that causes hepatitis C, there are more than one hundred strains of the virus identified so far, classified into six genotypes.Bacterial pneumonia is caused by more than ninety different serotyps of Streptococcus pneumoniae, some of which are much more common than others.Influenza type A viruses mutate continuously.These changes in the virus are called antigenic drifts.Although an infection with one strain of influenza type A leads to life long immunity, the antigenic drift produces new virus strains to which the host has only partial immunity or no immunity at all and can be reinfected with the disease.

Multi-strain disease interactions
The dynamics of the pathogen-host interactions involving multiple strains and the implications for the dynamics of the disease has fascinated researchers for a long time.Bremermann and Thieme [5] justify a competitive exclusion principle for the interactions of multiple strains by considering a multistrain SIR ODE model with possible acquired immunity to all strains and demographic renewal.In particular they show that the strain with the largest reproduction number (R 0 ) will outperform and eliminate the remaining strains in the system, provided that the growth of the host is limited by the carrying capacity of the environment.Castillo-Chavez, Huang and Li also establish that competitive exclusion is the norm in a two-sex SI model of gonorrhea with demographic renewal and two strains [10].
Many diseases, however, are represented by several or a multitude of strains which appear to coexist in nature.Dengue fever has four different serotypes which coexist often in the same geographical region.Infection with one of the serotypes gives permanent immunity but the same host remains vulnerable to infections with the remaining serotypes.Some particular sequences of infections with these four serotypes are believed to lead to the deadly hemorrhagic fever.The competitive exclusion and coexistence of two of the four dengue serotypes is discussed in [15].Feng and Velasco-Hernándes [15] consider the possibility of infecting an already infected individual with the other strain (a process called super-infection).They present numerical evidence that coexistence of the two strains is possible.A model of consecutive infections with two dengue fever serotypes is considered in [12] where it is also established that the two strains can coexist.
In fact, both the super-infection in [15] and the cross-immunity in [12] are known mechanisms that lead to coexistence of strains.Super-infection has been found earlier to lead to coexistence in a two-strain model in [26,35] and to more complex dynamics in multi-strain models in [38].Cross-immunity has been discussed mainly in relation to influenza and several articles report the presence of coexistence equilibria in this case [7,8,2,28].Similarly to super-infection, coinfection, which is defined as the simultaneous infection of the same host with two different parasites, or two different strains of the same parasite, also leads to coexistence of the strains on the population level [36].Under the same conditions as the ones considered in [5] the competitive exclusion principle will no longer be valid if the host population is allowed to grow exponentially in time.Apparently, in an exponentially growing host population there is enough "room" for two strains to coexist [31,1].Parasite polymorphism is obtained also via density dependent host death rate [3].
In this paper we consider another mechanism that is known to generate coexistence, namely mutation [4].By mutation of the strains on population level we understand a process of substitution on personal level of one of the strains by the other.In particular, we assume that if in a host, infected by the first strain, a small amount of the second strain is produced, then it takes over the host as a result of intra-host competition.Consequently, the same host is then infected by the second strain.We call the combined effect of these two processes mutation and denote the rate at which that happens by ρ.Li, Zhiu et al. [27] investigate an epidemic model of mutating pathogens in a recent article and find a unique coexistence equilibrium which may lose stability.In fact, mutation enters several other articles as a mechanism that promotes coexistence.In [6,14] incomplete treatment of individuals infected with tuberculosis (TB) leads to emergence of drug-resistant strain of the disease.The authors conclude that natural TB strains will not coexist under their models but a natural and a drug-resistant strain which is in fact a mutant of the natural strain will coexist under appropriate conditions.In [14] the natural strain can dominate only in the absence of mutation; if mutation is present there are only two possible outcomes: either the two-strains are both present in the population or the drug-resistant strain dominates.We obtain similar results in this article.

Vaccination in multi-strain diseases.
Vaccination is a wide spread method for disease prevention and control.It produces better results against diseases generated by pathogens of low mutability.One of the most successful vaccination campaigns is the campaign for the eradication of smallpox.The World Health Organization declared the disease eliminated in 1980.Measles has been essentially also eradicated in developed countries and vaccination against mumps and chicken pox gives promising results since the causative agents of these diseases show little tendency to vary antigenically.Although poliomyelitis is caused by three antigenic types, since they do not change significantly, vaccination against each one is necessary but produces promising results.On the other hand, vaccination against highly mutable viruses is either not very efficient or not at all possible at present time.Vaccination against bacterial pneumonia whose causative agent is represented by more than ninety serotypes, is carried out with vaccines containing agents of up to twenty-three of the most common serotypes.
Vaccines against highly mutable viruses such as as HIV and HCV at this time are not available.Providing adequate immunity against influenza has been a particularly challenging task.Because the virus continuously mutates and generates new strains any immunity furnished by a vaccine or by infection is short lived as the next flu season the same host faces a new set of strains.Thus, annual vaccinations are necessary.The vaccines are also updated every season.They are typically trivalent and consist of two type A strains and one type B strain.For the US, CDC estimates early in the year what strains are likely to be most distributed in the following flu season and makes a recommendation for the composition of the vaccine.The decision which strains to be included in the vaccine is based on methods for predicting the evolution of influenza A [25,22] as well as global surveillance.We model this epidemiological situation when the present virus, which is only partially targeted by the vaccine, mutates into a strain to which the vaccine is specific.
The impact of vaccination on the evolution of strain interactions in multistrain diseases has been investigated trough models in several articles [37,17,35,39,40].McLean [37] supports the view that as vaccination provides weaker immunity compared to infection, it gives favorable environment for the emergence of vaccine-resistant strains.T. Porco and S. Blower investigate [39] how the mode of action of potential HIV vaccine will influence the coexistence of two HIV subtypes.In particular, they find that if the vaccine provides full protection against subtype one for a given fraction of the vaccinated individuals and a complete protection against subtype two for a fraction of those protected against subtype one, then coexistence of the subtypes is possible.On the other hand if the vaccine acts only by decreasing the infectivity in vaccinated individuals infected with either subtype, then coexistence is not possible.In [40] the authors assume that the vaccine reduces the susceptibility of the vaccinated individuals and it is established that coexistence is possible.In both articles vaccination is applied before individuals enter the system.
Lipsitch [29] considers the interplay of two serotypes of bacteria subjected to serotype-specific or bivalent vaccine.He applies his theoretical results to shed light on existing data on serotype replacement in Haemophilus influenzae.The ability of the vaccines to target only specific strains of the causative agent has generated significant concern in epidemiology as this could increase the incidence of the disease from other strains not represented in the vaccine.This has not occurred with the use of H. influenzae type b vaccines but has occurred in trials of pneumococcal vaccines.In [30] these different outcomes are investigated with the use of mathematical models.

Multiple and subthreshold coexistence equilibria.
In simple epidemic models typically when the reproduction number is below one only the disease-free equilibrium exists.This equilibrium is locally and globally stable which implies that the disease will disappear from the population.Recommendations for disease control can be made based on that observation.In particular, measures which act to reduce the reproduction number below one will lead to the disease disappearance.
Recently, it has been observed in theoretical considerations that nontrivial equilibria can be present even when the reproduction number of the disease is smaller than one.One way for this situation to occur is through a phenomenon called backward bifurcation.In the case of backward bifurcation the endemic equilibrium which bifurcates from the disease-free equilibrium at the critical value one of the reproduction number exists for values of the reproduction number smaller than one.In fact, for values of the reproduction number between some minimal value R * , called minimal transition value, and one there are two or more endemic equilibria.In [34] it is established that there is typically an even number of equilibria with alternating stability so that the one with the lowest number of infectives is unstable.If backward bifurcation occurs, it is not sufficient to reduce the reproduction number below one to eradicate the disease but it is necessary to reduce it below a much lower value -the minimal transition value.Although this phenomenon is not as readily observed in data as oscillations, it plays a significant role in the dynamics of the disease and our ability to combat it effectively.
The presence of backward bifurcation in epidemic models has led to significant interest in the literature in recent years.In many cases backward bifurcation seems to be caused by the presence of several classes with different susceptibilities to the disease.Thus this phenomenon often occurs in multigroup models [20,11].As a special case, it can also be observed when the population is divided into never infected and previously infected individuals [16] and in educated and uneducated individuals [19].In addition, backward bifurcation appears when super-infection is present [13].Most of these models consider homogeneous populations with respect to age structure and, as a result, the corresponding models consist of ordinary differential equations.However, recently the heterogeneity of the host in age structure both chronological [9,33] and disease-induced [34] have also been found to lead to subthreshold equilibria.
One of the significant consequences of the backward bifurcation is the presence of multiple stable equilibria, which in turn leads to the fact that initial conditions determine to which equilibrium solutions may tend.This gives the opportunity, even in systems with time-independent coefficients, for multiple outcomes in the long-term development of the disease.
Backward bifurcation is also very typical for models involving partially effective vaccination as vaccination creates a class of susceptible individuals (namely the vaccinated individuals) with lower susceptibility to the disease compared to the regular susceptible class.Several articles report the existence of subthreshold equilibria in the presence of vaccination [18,23,24] but in those cases only one strain of the disease is considered.On the other hand two-strain models with or without vaccination are associated with two dominance equilibria, one for each strain present, and a unique coexistence equilibrium.
The impact of vaccination as a mechanism capable of generating multiple subthreshold equilibria on the dynamics of a disease in the context of multiple strain interactions has not been investigated so far.We address that impact in this article.

Organization of this article.
We introduce our two-strain model with vaccination and mutation in Section 2. The model consists of three ordinary differential equations and one partial differential equation structured by the time spent in the corresponding class.The differential equation for the total population size is given by the simplest population model of logistic growth.We also introduce several parameters which appear often in our discussion.
In Section 3 we discuss the existence of steady states.In the first Subsection we consider the case when there is no mutation ρ = 0. We provide the vaccine-dependent reproduction number of strain one in absence of mutation R • 1 (ψ), where ψ is the per capita vaccination rate.We include the vaccine-dependent reproduction number of strain two R 2 (ψ) which is independent of mutation.We establish that if R 2 (ψ) > 1 there is always the equilibrium E * 2 which corresponds to dominance of the second strain.If R • 1 (ψ) > 1 there is always the equilibrium E * 1 which corresponds to dominance of the first strain.If R • 1 (ψ) < 1 there might be none or two equilibria E * 11 , E * 12 .We introduce the invasion reproduction numbers of strain one and two, defined as the number of secondary cases one infectious individual with strain i will produce during the time it is infectious in a population where strain j is at equilibrium.We show that if each strain can invade the stable dominance equilibrium of the other there exists a unique coexistence equilibrium E • * * .In the second Subsection we introduce the reproduction number of strain one in the presence of mutation R 1 (ψ).We establish that if R 1 (ψ) < 1 there might be none or two coexistence equilibria: E * * 21 , E * * 22 and if R 1 (ψ) > 1 there might be up to three coexistence equilibria.
Section 4 is devoted to the local stability of equilibria.In Subsection 4.1 we investigate the stability of the disease-free equilibrium where we obtain the typical result: if both reproduction numbers are below one it is locally asymptotically stable, if any of the reproduction numbers is above one -it is unstable.Subsection 4.2 investigates the local stability of the equilibrium E * 2 and establishes that it is locally stable if R 1 (ψ) < R 2 (ψ) and unstable otherwise.Subsection 4.3 investigates the stability of equilibria which exist when there is no mutation.It establishes that the equilibrium E * 11 is always unstable whenever it exists.Section 5 overviews our observations and their relation to disease control and prevention.Some of the more technical derivations and proofs are presented in the Appendix.

A two-strain model with vaccination.
In this section we introduce a two-strain epidemic model.We consider a population whose total population size at time t is denoted by N (t).The dynamics of the total population in the absence of a disease is described by the simplest demographic model which accounts for a limited population size, namely, we assume constant birth/recruitment of individuals in the population, denoted by Λ, and constant per capita natural death rate, denoted by µ: This equation has the globally stable steady state and we assume that this state has been attained so that the total population stays constant at all times.Now we assume that a disease is spreading in the population.The presence of the disease divides the population into non-intersecting subclasses.The individuals who are healthy but can contract the disease and have not been previously vaccinated, form the class of susceptible individuals whose size at time t is denoted by S(t).Individuals who are healthy but have been vaccinated against the causative agent of the disease form the vaccinated class.The size of the total population in the vaccinated class is denoted by V (t).The per capita vaccination rate is denoted by ψ.It is assumed that the vaccine protection does not wane.
We assume that two genetically distinct forms of the causative agent of the disease are present and can infect the individuals in the population: we call the first form strain one and the second form -strain two.Individuals who are infected with strain one form the class whose size is denoted by I(t).The individuals in this class are stratified according to their time-since-infection θ, and their density is given by i(θ, t).Individuals who are infected with strain two form the class whose size is denoted by J(t).We assume there are no individuals who are simultaneously infected by both strain one and strain two, that is, we assume there is no coinfection.However, strain one can mutate into strain two at a rate ρ(θ).A completely susceptible individual can come into contact with an infective individual and become infected.If the infectious individual is a carrier of strain one the susceptible individual becomes infected also with strain one, at a rate β 1 (θ).Such an individual progresses to the infectious class but enters it with age-since-infection equal to zero.If the infectious individual is a carrier of strain two, the susceptible individual becomes infected with strain two, at a rate β 2 .The newly infected individual moves to the J-class.
The susceptible individuals are not the only ones who are healthy and can become sick after coming into contact.We assume that the vaccine is tailored to protect against strain two and it is completely effective against it.We call strain two vaccine-sensitive strain or simply vaccine strain.However, the vaccine is only partially effective against strain one.For its ability to elude the immune response promoted by the vaccine, we call strain one vaccine-evasive strain.A vaccinated individual can become infected after being in contact with an individual infected with strain one at a rate β 1 (θ)δ, where δ reflects vaccine imperfection with respect to strain one.The vaccine is perfect if δ = 0 and no vaccinated individual can be infected.On the other hand if δ = 1 the vaccine plays no protective role.We assume 0 ≤ δ ≤ 1.
The removal rate from the class i(θ, t) is given by the function γ(θ).The total rate at which individuals recover from the class i(θ, t) alive is given by the quantity A proportion χ of those recover to the vaccinated class, to account for those who entered the class i(θ, t) already vaccinated.We assume that 0 ≤ χ ≤ 1.A proportion (1 − χ) recover to the susceptible class.The recovery rate from the class J is α.Since only susceptible individuals can get infected with the second strain, they recover to the susceptible class only.
The model takes the form (see Figure 1): The parameters of the model and their meanings are listed in Table 1.We remark that this model with strain two and the age-since-infection structure removed is exactly the model considered in [18].
We assume that all parameters are nonnegative and µ > 0. We also assume that To ensure well-posedness we endow the system (2.2) with initial conditions: where S 0 , J 0 , and V 0 are given nonnegative constants, while i 0 (θ) is a given nonnegative, integrable function.The initial conditions must satisfy the relation (see (2.1)) so that the system is consistent with the assumption on the total population and (2.3) The system (2.2) with the initial conditions is well-posed, that is, independently of what nonnegative initial conditions are taken, the model has a unique nonnegative solution which depends continuously on the initial data.This result can be established with standard techniques.
The following notations will be used throughout the paper.The quantities give the proportion of individuals remaining in the infectious class of strain one until progression age θ in the cases when ρ(θ) = 0 and when ρ(θ) = 0 respectively, given that the individuals have survived till that age.Next we define It is easy to see that 0 ≤ Γ ≤ Γ • < 1 (see identity (2.7)).The parameters Γ and Γ • give the proportion of individuals leaving the infectious period of strain one through recovery in the cases without and with mutation respectively.The proportion of individuals leaving the infectious class of strain one through mutation of the strain with which they are infected is given by The proportion of individuals who die while infectious with strain one is given by µ∆ • in the case of no mutation, and by µ∆ in the case with mutation where ∆ • and ∆ denote the integrals Naturally, the sum of the proportions of surviving the infectious with strain one period and dying while in it is equal to one: Γ + φ + µ∆ = 1. (2.7) In fact, this equality can be justified rigorously by integrating by parts the integral in Γ.Finally, we introduce the following notation In the next section we discuss the equilibria of the model.

Steady states.
Using (2.1) and the notation where S * , i * (θ), J * , V * is a time-independent solution of (2.2), we obtain the following system for the equilibria which consists of four algebraic equations and one ordinary differential equation: From the equation for the total population size (2.3) we get also the following algebraic condition: that is also a consequence of (3.1).
The equilibrium which is always present is the disease-free equilibrium, that is, an equilibrium in which there are no infected individuals: where is the proportion of susceptible individuals in the disease-free population and is the proportion of vaccinated individuals in the disease-free population.
To find endemic equilibria, that is, equilibria in which the number of infectives is not zero, we solve the second equation in (3.1), with i(0) given by the third equation, getting then substitute the solution in the remaining equations to obtain a system of nonlinear algebraic equations: where i * denotes i(0).From the last equation in the system (3.5)we express v * in terms of i * and s * as Solving for s * in the first equation we have Observe that s * thus defined is positive and smaller than one if B 1 > (1 − χ)Γ and β 2 > α.
In this section we analyze the case in which there is no mutation (ρ(θ) = 0) so that the two strains are not directly connected but, nevertheless, they compete through susceptibles.In this case system (3.5)becomes The existence and stability of endemic equilibria depends on two parameters -the reproduction number of each strain.The reproduction number of strain 1 in absence of mutation is given by (see the proof of Proposition 4.1 for derivation as well as the interpretation below) and in the absence of vaccination it is Concerning the interpretation of the reproduction number, we notice that the quantity R • 1 = B • 1 gives the number of secondary infections that strain one will generate in a completely susceptible population.However, in the absence of the disease our population consists of both susceptible and vaccinated individuals.The proportion of individuals who are susceptible is µ µ+ψ (see (3.3)).Thus, the first term in R • 1 (ψ), given by , gives the number of secondary infections of susceptibles an infected individual can produce in a disease-free population.The number δB •  1 gives the number of secondary infections an infected individual can produce in a vaccinated population; ψ µ+ψ is the proportion of vaccinated individuals in a disease-free population.Thus, the second term in R • 1 (ψ), given by (µ+ψ) , gives the number of secondary infections of vaccinated individuals an infected individual can produce in a disease-free population.
The reproduction number of the second strain is given by and in the absence of vaccination it is R 2 = R 2 (0) = β2 µ+α .We note that the reproduction number of the vaccine-sensitive strain is not influenced by the presence or absence of mutation because mutation does not lead to infections of healthy individuals with the second strain, and hence, has no impact on the number of secondary infections one infected with strain two individual can produce in a population of susceptible and vaccinated individuals.Since the second strain does not infect vaccinated individuals, its reproduction number consists of a term that corresponds to the first term in R • 1 (ψ) only.Since 1 µ+α is the mean time spent as infected with strain two, the expression β2 µ+α gives the number of secondary infections the vaccine-sensitive strain can produce in a entirely susceptible population, and µ µ+ψ is the proportion of susceptibles in a disease-free population.It can be seen that both R • 1 (ψ) and R 2 (ψ) are decreasing functions of ψ which is the expected impact of a vaccination campaign.Because the vaccine protects completely against the vaccine-sensitive strain, its reproduction number can be made very small if the vaccination rate is sufficiently large.On the other hand the vaccine is only partially protective against the vaccine-elusive strain.Thus, even very large levels of vaccination cannot reduce the reproduction number of the first strain below some threshold value We recall that B Vaccination acts to decrease the reproduction numbers of the strains but it often does so by a different amount.As a result it has the ability to switch the relationship between them.This observation is the source of the concern expressed in the literature that although vaccination might lead to the effective control of some strains, it might also lead to the proliferation of others.We notice this phenomenon occurring in this model too.In particular, if in the absence of vaccination the second strain dominates the first one, that is, R • 1 < R 2 then there is a vaccination level ψ * given by The relation between the reproduction numbers is not changed by vaccination if In this case we also have R 2 (ψ) < R • 1 (ψ) for all ψ.First we investigate existence of dominance equilibria.It is easy to see that besides the disease-free equilibrium E 0 , a unique equilibrium corresponding to nonzero levels of infected with the vaccine-sensitive strain J is also feasible: exists.
In addition one or more nontrivial equilibria corresponding to nonzero values of the vaccine-evasive strain may exist.We consider the system (3.8) with j * = 0.
From the first equation we get for s * (see (3.7)): In the second equation we can first cancel i * , substitute the value for v * from (3.6) to get We determine i * so that the two expressions for s * are equal.Thus, the equilibrium value for i is a solution for the equation δB In particular, we have the following result similar to the one above.
exists where i * 0 (θ) is determined from formula (3.4) with i(0) given by the unique positive solution of the equation (3.12), s * 0 is determined from formula (3.11), and v * 0 is obtained from (3.6).
Proof.To see the existence and uniqueness denote by ) is actually equivalent to a quadratic equation and so it has at most two solutions.Since Hence, equation (3.12) has at least one solution and the number of the intersections of the two functions is odd.Therefore, the equation has exactly one positive solution i * .♦ If R • 1 (ψ) < 1 there may be no equilibria with the vaccine-evasive strain present or, under some additional condition on the parameters of the model, there may be two equilibria.The existence of these two equilibria depends on the occurrence of backward (subcritical) bifurcation at the critical value i = 0 (R • 1 (ψ) = 1).Since the reproduction number R • 1 = B • 1 , we choose for a bifurcation parameter B • 1 .We rewrite equation (3.12) in a more convenient form: The equation above defines B • 1 as a function of i, that is, The bifurcation in B • 1 at the critical value i = 0 is backward if and only if dB • 1 di (0) < 0. Using implicit differentiation in the equation above we obtain dB where B • 1 (0) is the value of B • 1 at the critical value i = 0: 0) < 0, existence of the two subthreshold equilibria occurs for values of B • 1 in some non-empty interval with right end-point B • 1 (0).We have the following proposition then there exists a constant R • inf ∈ [0, 1) such that for Backward bifurcation will be feasible if the inequality (3.14) is satisfied for some values of the parameters.It is not hard to see that the parameter space for which this inequality is true is not empty.If, in particular, we take Γ • = 0.9 (or γ(θ) = 0.9 (time) −1 ), µ = 0.1 (time) −1 , χ = 0.0, δ = 0.1, ψ = 1 (time) −1 then the inequality will be valid with the left-hand side having value 0.21 while the right-hand side having value 0.235.Inequality (3.14) is a necessary and sufficient condition for existence of backward bifurcation.It reveals the mechanisms in the model (2.2) which promote this phenomenon and those which obstruct it.Inspection of (3.14) shows that if χ = 1, that is if there is no recovery to the susceptible class there will be no backward bifurcation.Vaccination itself is the chief mechanism responsible for the presence of subcritical bifurcation in this model ψ = 0. Indeed, in the case when ψ = 0 inequality (3.14) becomes which after cancelling µ and rearranging the terms becomes which is clearly impossible since the right-hand side is smaller than one.Some specific components of the vaccination process are particularly involved.In particular, if δ = 0 the inequality (3.14) cannot be satisfied.Consequently, if the vaccine is perfect with respect to both strains, the disease will not be able to establish itself for R • 1 < 1.To summarize, vaccination, and particularly, vaccine imperfection is the main mechanism for supporting the presence of the disease even when the reproduction number is below one.
Finally, we consider the coexistence endemic states.As it turns out, under some conditions on the parameters there exists a unique state with i * > 0 and j * > 0. To address the question of coexistence equilibrium we define several invasion reproduction numbers.First, define the invasion reproduction number of strain one in the case of no mutation ρ = 0 (see proof of Proposition 4.2 for derivation): The invasion reproduction number of strain one gives the number of secondary cases one infected individual with strain one can infect in a population where strain two is at equilibrium E * 2 .We note that the condition R• 1 > 1 is a condition that strain one can invade the equilibrium of strain two.Strain two might be invading one of three equilibria of strain one E * 1k for k=0,1,2.To each of these equilibria of strain one corresponds one invasion reproduction number of strain two.Define the invasion reproduction numbers of strain two in the case of no mutation ρ = 0: where i * k is i(0) from equation (3.4) which gives the proportion of infected with strain one in the corresponding equilibrium E * 1k .We note that the condition R• 2 (E * 1k ) > 1 says that strain two can invade the corresponding equilibrium of strain one.
The following result establishes the existence of coexistence equilibrium in the case of no mutation.We note that the conditions in the proposition below are only sufficient and there may be coexistence even if some of them are not satisfied.Furthermore, as the proof shows there is no coexistence without vaccination ψ = 0 in the region R At the same time positive vaccination levels lead to coexistence if each strain can invade the stable equilibrium of the other and R 2 > 1.In Figure 2 we illustrate the coexistence in the absence of mutation.The simulation in Figure 2 suggests that the coexistence equilibrium is stable.
Proposition 3.4 Assume and R 2 > 1.There are two cases: Then there is a unique coexistence equilibrium Then there is a unique coexistence equilibrium Proof.We denote the coexistence equilibrium by E • * * = (ŝ, î(θ), ĵ, v) where î(θ) = i(0)π 0 e −µθ with i(0) = î.We show that î exists so that system (3.8) is satisfied.From the third equation in system (3.8)we have ŝ = 1 R2 < 1.From the second equation in (3.8), after cancelling î we have The expression for v from (3.6) becomes v = ψŝ + χ îΓ • δB • 1 î + µ Substituting v above we have that ŝ is given by the expression in (3.11) with î in place of i * .Thus î is a solution of the following equation δB The left hand side of this equation is a function of i, denoted as before with f 1 (i).One can see that f 1 (i) is a monotone function of i.Thus, this equation has at most one solution.Consequently, the coexistence equilibrium, if it exists, is unique.Case 1: To see the existence in this case, notice that R• This, in particular, leads to the fact that equation (3.17) has no solution if f 1 (i) is decreasing.We note here that in case there is no vaccination ψ = 0 it follows that f 1 (i) is decreasing and there is no coexistence in the region R Consequently, there exists î in the interval (0, i * 0 ) such that equality (3.17) is satisfied.To finish the proof for Case 1 notice that since î < i * 0 we have f 1 ( î) < g 1 ( î), where g 1 is the function in the proof of Proposition 3.2.Then from the first equation in (3.8) we have That establishes the existence of coexistence equilibrium in Case 1, given that β 2 ŝ − α = µ > 0. Case 2: On the other hand, assumption R• That implies that there exist a solution î of equation (3.17) that lies in the interval (i * 1 , i * 2 ) where i * 1 gives the equilibrium E * 11 and i * 2 gives the equilibrium E * 12 .Consequently in this case f 1 ( î) < g 1 ( î) and ĵ > 0. If R• 2 (E * 11 ) > 1 and R• 1 > 1 then î is in the interval (0, i * 1 ) and f 1 ( î) > g 1 ( î).Consequently, ĵ < 0. This concludes the proof.♦ 3.2 Nontrivial equilibria in presence of mutation: the case ρ(θ) = 0.
In this case the vaccine-elusive strain mutates into the vaccine-sensitive strain and there are coexistence equilibria such that the possible ultimate outcomes are either coexistence of the two strains, or competitive dominance of the second strain.It is interesting that genetic changes alone can give the competitive advantage to the vaccine-sensitive strain.The existence and stability of nontrivial equilibria again depends on two reproduction numbers of the strains.The reproduction number of strain 1 in presence of mutation is given by and in the absence of vaccination it is The interpretation of the reproductive number is similar to the one before.The quantity R 1 = B 1 gives the number secondary infections that strain 1 will generate in a completely susceptible population in the absence of vaccination.
As before, R 1 (ψ) is a decreasing function of ψ whose minimal value is We note that the reproduction number of the second strain remains unchanged and is given by (3.9).Finally, if in the absence of vaccination the second strain has larger reproduction number than the first, that is R 1 < R 2 , then there is a vaccination level If in the absence of vaccination the first strain has a larger reproduction number than the second, that is R 2 < R 1 , then this relation is preserved for all vaccination levels: R 2 (ψ) < R 1 (ψ).
Besides the disease-free equilibrium E 0 , a unique equilibrium corresponding to nonzero levels of infected with the vaccine-sensitive strain J is also feasible: We note that there is no endemic state corresponding to the absence of strain two, but one or more coexistence equilibria may exist.The existence of coexistence equilibria depends on the invasion reproduction number of the first strain when ρ(θ) = 0: Since in the case ρ(θ) = 0 there is only one dominance equilibrium E * 2 , there is also only one invasion reproduction number.To find the coexistence equilibria, we consider the system (3.5).From the second equation where we use the value of v * from (3.6) we express s * as a function of i * : From the third equation in (3.5) we express j * in terms of i * : where the function ω(i) is defined as We notice that ξ(i) is a monotone function of i but it could be increasing or decreasing.We have The function ξ(i) is defined and positive at i * (and therefore j * is defined and nonnegative) if and only if ω(i * ) < 1, that is, if and only if the following inequality is satisfied: for the corresponding i * .From here a condition for the absence of coexistence equilibria can be derived: Proposition 3.6 If δψ > χΓ(µ + δψ) and R1 ≤ 1 then there are no coexistence equilibria.
♦ On the other hand we note that if R 2 ≤ R 1 the inequality (3.21) is satisfied for all i * (we recall that B 1 = R 1 ) so that the function ξ(i) is defined and positive at any i * .Another sufficient condition for inequality (3.21) to be satisfied is given by the following hypothesis In fact, by a similar argument as in Proposition 3.6, if this condition is satisfied then inequality (3.21) is valid for every i.
From (3.7) we get a second expression for s * in terms of i * , using the fact that j * = ξ(i * )i * : We determine i * so that the two expressions for s * are equal.Thus, the equilibrium value for i is a solution for the equation Establishing the existence of coexistence equilibria can be done under weaker and more natural conditions than (3.22).In the proposition below we assume that the invasion reproduction number of the first strain R1 > 1.
We notice that if we know that (3.22) is satisfied that implies that R1 > 1 and the proposition is still valid.
Proposition 3.7 If R 1 (ψ) > 1 and R1 > 1 are satisfied then there exists at least one and up to three coexistence equilibria where each i k , a positive solution of the equation (3.24), gives i(0 , then there is an odd number of them.
Case 2: R 2 (1 − χΓ) ≤ R 1 .This inequality implies that Consequently, inequality (3.21) is satisfied for all i since f 2 (i) is a monotone function and f 2 (0) < 1 R2 .This implies that both functions f 2 and g 2 are defined and continuous on the interval (0, ∞).We have again that f 2 (0) < g 2 (0).So we show that as i → ∞ the limits of the two functions satisfy the opposite inequality.First, we note that lim = ξ(∞).
This gives the following limit for g 2 : Consequently, in order to show that lim i→∞ f 2 (i) > lim i→∞ g 2 (i) we have to show the inequality: Replacing ξ(∞), cancelling B 1 from both sides and rearranging terms we arrive at the following inequality, that has to be established: Adding µφR 1 to left-hand side we get a stricter inequality Cancelling a common term from both sides Moving the αφ-term inside the brackets and using identity (2.7) in the left-hand side, it becomes Cancelling (α + µ)R 1 term and moving the remaining term to the left-hand side while moving the µ∆-term to the right-hand side we obtain Dividing by (α + µ) and moving the first term from the right-hand side to the left-hand side we have Dividing by µ∆ we arrive at a true inequality (compare with the inequality in Case 2) Therefore, the initial inequality (3.25) is also true.
In both cases the number of the intersections of the two functions is odd.This completes the proof.♦ If R 1 (ψ) < 1 there may be no coexistence equilibria or, under some additional condition on the parameters of the model, there may be two coexistence equilibria.The existence of these two coexistence equilibria depends on the occurrence of backward (subcritical) bifurcation at the critical value i = 0 (R 1 (ψ) = 1).We again choose for a bifurcation parameter B 1 .As before, we cross-multiply in equation (3.24) to obtain: This equation defines B 1 as a function of i, that is, B 1 = B 1 (i).The bifurcation in B 1 at the critical value i = 0 is backward if and only if B ′ 1 (0) < 0. Using implicit differentiation in equation above, setting i = 0 and solving for B ′ 1 (0) we obtain where B • 1 is the value of B 1 at the critical value i = 0: .
We have the following proposition for the existence of the two subthreshold equilibria: then there exists a parameter R inf ∈ [0, 1) such that, if then the equilibria Inequality (3.26) is nontrivial.In particular, it is satisfied for the parameters in Figure 3.The presence of stable subthreshold coexistence equilibria is illustrated in Figure 3 where coexistence is possible when both reproduction numbers are below one.
In the next section we investigate the local stability of equilibria.

Local stability of equilibria
To investigate the local stability behavior of equilibria we look at the linearized right-hand side of system (2.2).This operation is the analogue of taking the Jacobian in ordinary differential equation models.In particular, we consider an equilibrium (S * , i * (θ), J * , V * ) of system (2.2) and we set where we denote by (x, ȳ(θ), z, w) the deviations from such an equilibrium, and by n the deviation for the total population size.Then we can linearize the nonlinear terms and determine the eigenvalues of the linearized problem by looking for solutions of the form x = e λt x, ȳ = e λt y(θ), z = e λt z, w = e λt w, n = e λt n.
We note that (2.1) implies that We obtain the following linear eigenvalue problem: where i * above is i(0) from (3.4) in the corresponding equilibrium.We introduce the following notation which will be useful in this section: These two quantities satisfy a relation similar to the relation between φ, Γ and ∆: This equality can be established through integration by parts.Then, solving the ordinary differential equation y(θ) = y(0)e −(λ+µ)θ π(θ) and substituting it in the remaining equations we obtain the following linear eigenvalue system for the real variables x, y, z, and w, where y is a shorthand notation for y(0): Our aim in the following subsections is to examine the eigenvalues of this problem in correspondence with each equilibrium that has been proved to exist in the previous sections.

Local stability of the disease-free equilibrium.
When the infection-free equilibrium is concerned, we have i * (θ) = 0, j * = 0. Since the existence of this equilibrium is not influenced by the presence or the absence of mutation we consider both cases simultaneously, that is for a general ρ.In this case, the linear eigenvalue problem (4.3) above takes the form Thus we have Proof.We can rewrite the first characteristic equation (4.5) as We first note that G In addition, G 1 (λ) is a decreasing function of λ for λ real with G 1 (λ) → 0 as λ → ∞.Hence there is a positive real λ * which solves equation (4.7) and the disease-free equilibrium E 0 is unstable.If R 1 (ψ) < 1, then for λ with ℜλ ≥ 0 we have Thus equation (4.7) has no solutions with nonnegative real part.The second characteristic equation can be explicitly solved for λ: λ The corresponding eigenvalue is clearly negative if and only if R 2 (ψ) < 1.Consequently, if R 1 (ψ) < 1 and R 2 (ψ) < 1 then the eigenvalue problem (4.4) has eigenvalues with only negative real part and the disease-free equilibrium is locally asymptotically stable.If R 1 (ψ) > 1 the first characteristic equation (4.5) has a positive real solution.If R 2 (ψ) > 1 the corresponding solution of the second characteristic equation is positive.In these two cases the disease-free equilibrium is unstable.♦ 4.2 Local stability of the vaccine-sensitive strain equilibrium E * 2 .
Now we consider the stability of the vaccine-sensitive strain equilibrium E * 2 .In this case also the existence of the equilibrium is not depending on the mutation parameter and we can treat the cases ρ = 0 and ρ = 0 simultaneously.We recall that in either case the equilibrium with only the vaccine-sensitive strain present E * Proposition 4.2 Let R 2 (ψ) > 1.Then the equilibrium E * 2 is locally asymptotically stable if R1 < 1 and unstable if R1 > 1.
Proof.We consider the linear eigenvalue problem (4.3) with i * = 0 and with s * , j * , v * given by the coordinates of E *  In order for this system to have a nonzero solution we need the determinant to be zero.From that condition we get that one eigenvalue is λ = −µ and all remaining eigenvalues of the problem are provided by the following two characteristic equations (λ + β 2 j * + µ + ψ)(λ + µ + α) = β 2 s * (λ + µ + ψ) + αβ 2 j * (4.9) (s * + δv * ) B1 (λ) = 1 (4.10) We consider first (4.9).We notice that we can cancel αβ 2 j * from both sides of this equation thus obtaining Furthermore, observing that β 2 s * = µ + α we can simplify this equation to the following quadratic equation in λ: λ 2 + (µ + ψ + β 2 j * )λ + µβ 2 j * = 0 whose solutions have negative real parts or are negative numbers.Consequently, the equation (4.9) has no solutions with ℜλ ≥ 0. Now we turn our attention to equation (4.10).We rewrite it in the form where G 3 (λ) = (s * + δv * ) B1 (λ) First, we notice that Second, G 3 (λ) is a decreasing function of λ for λ real and positive.In addition, G 3 (λ) → 0 as λ → ∞.Consequently, if there exists a real positive solution of the equation and the equilibrium E * 2 is unstable.Third, if R 1 (ψ) < R 2 (ψ), for λ with ℜλ ≥ 0 we have and, therefore the equation G 3 (λ) = 1 has no solutions with ℜλ ≥ 0. It follows that the vaccine-sensitive strain equilibrium is locally stable.♦ coexistence equilibria, some of which occur when both reproduction numbers are below one: R 1 (ψ) < 1 and even if R 2 (ψ) < 1.We call such equilibria strongly subthreshold equilibria.Thus, although the vaccine is designed to provide 100% protection with the respect to the second strain, the second strain will persist if R 1 (ψ) < 1 and even if R 2 (ψ) < 1 (see Figure 3).The presence of strongly subthreshold coexistence equilibria has significant implications for the disease control as reducing both reproduction numbers below one will not eradicate either strain or the disease.Thus, a vaccine that is specific and perfect to the vaccine-strain, and has the potential to eliminate it by reducing its reproduction number below one if the vaccine-elusive strain is not present, may not necessarily do so if the vaccine strain is a mutant of another strain to which the vaccine is only partially effective.Thus, in summary, the partial effectiveness of the vaccine enables (in some cases) the backward bifurcation, which in turn enables the vaccine-evasive strain to persist when R 1 (ψ) < 1, which in turn enables the vaccine target strain to persist when it shouldn't.Actually, we suspect that mutation is only one of a whole range of mechanisms generating coexistence which leads to reducing the effectiveness of the vaccine with respect to the second strain as a results of its partial effectiveness with respect to the first strain.Many of the well known coexistence mechanisms, e.g.cross-immunity, coinfection, super-infection and others, may have similar consequences.The main implication of this work is that through mutation (ρ(θ) = 0) the subthreshold existence of strain one, generated by the vaccine's partial effectiveness to this strain, translates into subthreshold existence of strain two to which the vaccine is fully effective.Another effect of vaccination that we observe in this paper is its ability to lead and enhance pathogen polymorphysm in multistrain diseases.In particular we establish that in the absence of mutation ρ(θ) = 0, nonzero vaccination levels ψ > 0 lead to coexistence of the two strains under certain conditions while coexistence is ruled out for that region of the parameter space when ψ = 0. Thus effectively vaccination leads to coexistence.In fact the role of vaccination as a coexistence mechanism can be derived from the considerations of the first model in [37].

Figure 2 :
Figure 2: The left figure illustrates that the number infected with strain one I(t) and the number infected with strain two J(t) may tend toward a coexistence equilibrium when ρ = 0, and ψ = 0.5.The right figure illustrates that if ψ = 0 strain two eliminates strain one.The remaining parameters used for these figures are β 1 = 6, β 2 = 4.5, γ = 0.8, µ = 0.1, χ = 0.0, δ = 0.04, Λ = 5.The units of Λ are (number of people)/(unit of time), χ and δ are dimensionless.The remaining parameters have units given by (unit of time) −1 .The corresponding reproduction numbers are dimensionless and are given by R • 1 (ψ) = 1.333 and R 2 (ψ) = 1.25.The reproduction numbers in the absence of vaccination are R • 1 = 6.66667 and R 2 = 7.5.

Figure 4 : 2 = 9 .
Figure 4:  This figure shows that the number infected with strain one I(t) and the number infected with strain two J(t) tend toward a coexistence equilibrium even though vaccine-dependent reproduction number of the first strain is above one while the vaccine-dependent reproduction number of the second strain is below one.The parameters used for this figure are ρ = 0.1, γ = 0.8, µ = 0.1, χ = 0.0, δ = 0.1, ψ = 1.0,Λ = 5, β 1 = 6.5 and β 2 = 9.The unites for Λ are (number of people)/(unit of time), δ and χ are dimensionless, and the units for the remaining parameters are (unit of time) −1 .The corresponding reproduction numbers are R 1 (ψ) = 1.1818 and R 2 (ψ) = 0.7438.

Table 1 :
List of Parameters per capita recovery rate from the class i(θ, t) α per capita recovery rate from the class J χ proportion of recoveries to the vaccinated class, given recovery from strain 1 1 − δ efficiency of vaccine ρ(θ) per capita mutation rate of strain 1 into strain 2 is not satisfied then there are no such equilibria for R 1 (ψ) < 1.