==== Front PLoS One PLoS One plos plosone PLoS ONE 1932-6203 Public Library of Science San Francisco, CA USA PONE-D-20-07823 10.1371/journal.pone.0243408 Research Article Biology and Life Sciences Immunology Vaccination and Immunization Medicine and Health Sciences Immunology Vaccination and Immunization Medicine and Health Sciences Public and Occupational Health Preventive Medicine Vaccination and Immunization Computer and Information Sciences Systems Science System Stability Physical Sciences Mathematics Systems Science System Stability Medicine and Health Sciences Medical Conditions Infectious Diseases Infectious Disease Control Vaccines Viral Vaccines Biology and Life Sciences Microbiology Virology Viral Vaccines Biology and Life Sciences Evolutionary Biology Evolutionary Processes Evolutionary Emergence Biology and life sciences Organisms Viruses RNA viruses Orthomyxoviruses Influenza Viruses Biology and Life Sciences Microbiology Medical Microbiology Microbial Pathogens Viral Pathogens Orthomyxoviruses Influenza Viruses Medicine and Health Sciences Pathology and Laboratory Medicine Pathogens Microbial Pathogens Viral Pathogens Orthomyxoviruses Influenza Viruses Biology and Life Sciences Organisms Viruses Viral Pathogens Orthomyxoviruses Influenza Viruses Medicine and Health Sciences Epidemiology Infectious Disease Epidemiology Medicine and Health Sciences Medical Conditions Infectious Diseases Infectious Disease Epidemiology Medicine and Health Sciences Medical Conditions Infectious Diseases Viral Diseases Covid 19 Physical Sciences Mathematics Algebra Linear Algebra Eigenvalues The local stability of a modified multi-strain SIR model for emerging viral strains Modified multi-strain SIR model for emerging viral strainshttps://orcid.org/0000-0002-0148-8090Fudolig Miguel ConceptualizationData curationFormal analysisMethodologySoftwareValidationVisualizationWriting – original draftWriting – review & editing* Howard Reka ConceptualizationFunding acquisitionSupervisionWriting – original draftWriting – review & editing Department of Statistics, University of Nebraska-Lincoln, Lincoln, NE, United States of America Khudyakov Yury E Editor Centers for Disease Control and Prevention, UNITED STATES Competing Interests: The authors have declared that no competing interests exist. * E-mail: mfudolig@huskers.unl.edu 2020 9 12 2020 9 12 2020 15 12 e024340818 3 2020 22 11 2020 © 2020 Fudolig, Howard2020Fudolig, HowardThis is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.We study a novel multi-strain SIR epidemic model with selective immunity by vaccination. A newer strain is made to emerge in the population when a preexisting strain has reached equilbrium. We assume that this newer strain does not exhibit cross-immunity with the original strain, hence those who are vaccinated and recovered from the original strain become susceptible to the newer strain. Recent events involving the COVID-19 virus shows that it is possible for a viral strain to emerge from a population at a time when the influenza virus, a well-known virus with a vaccine readily available, is active in a population. We solved for four different equilibrium points and investigated the conditions for existence and local stability. The reproduction number was also determined for the epidemiological model and found to be consistent with the local stability condition for the disease-free equilibrium. The authors received no specific funding for this work. Data AvailabilityThis study was mostly theoretical and used simulated data from R for illustration. The R code used for the simulation can be accessed through https://github.com/MiguelFudolig/multistrainSIR.OutbreaksCOVID-19Data Availability This study was mostly theoretical and used simulated data from R for illustration. The R code used for the simulation can be accessed through https://github.com/MiguelFudolig/multistrainSIR. ==== Body Introduction In recent times, the anti-vaccination movement has been gaining traction in different parts of the world. Individuals who do not advocate vaccination commonly cite reasons of fear of adverse side effects, perceived low efficacy of vaccines, and perceived low susceptibility to diseases amongst others [1–3]. The drop in numbers in vaccination has led to outbreaks of diseases such as mumps [4–6] that could have been prevented by vaccination. Another example would be measles, which was declared eliminated in the United States back in 2000, has had outbreaks reported in the country since 2008 [7]. According to the Centers for Disease Control and Prevention (CDC), the reemergence is due to the presence of unvaccinated individuals and their interaction with other people who got the disease from other countries such as Israel, Philippines, and Micronesia [7–9]. According to the CDC, 880 individual cases of measles have been confirmed in 24 states as of May 2019, the highest since 1994. As of May 2019, there are 10 active measles outbreaks in ten jurisdictions in the US. Another way for a disease to reemerge is through change in its antigenic properties, which is the case for the influenza virus. The influenza virus can mutate in two ways: through antigenic shift or antigenic drift [10–12]. Antigenic drift is defined as the result of frequent mutations of the virus, which happens every 2-8 years. On the other hand, the antigenic shift occurs around three times every one hundred years and only happens with influenza A viruses [12]. Although more unlikely to happen than the antigenic drift, the antigenic shift involves genetic reassortment which can make it feasible to create a more virulent strain than the original strain [12–14]. Infectious diseases such as the influenza virus can be modeled in a variety of ways. One can use phenomenological methods on available empirical data [15–17], Bayesian inference using Monte Carlo methods to estimate parameters through simulation based on the stochastic model [18–20], or agent-based network dynamics to model the spread of the disease through interactions between individuals [21–23]. This paper focuses on modeling an epidemic using compartmental systems, which involves separating the population to multiple components and describing infection and recovery as transitions between the set components. The simplest compartmental model to describe a viral infection is called the standard Susceptible-Infected-Removed (SIR) model, the dynamics of which has been studied in different references [24–28]. The SIR model separates the population into three compartments: the susceptible (S), infected (I), and removed (R) compartments. The susceptible compartment is comprised of individuals that are healthy but can contract the disease. The infected compartment is comprised of individuals who have already contracted the disease. Lastly, the removed compartment is comprised of individuals who have recovered from the disease. Individuals who have recovered from a certain strain of a viral infection are likely to be immune to infection of the same strain [29, 30], which is why the SIR model is used to model viral infections. SIR models can provide insight on the dynamics of the system and has been used to model different influenza virus strains such as the swine and avian flu focusing on the spatio-temporal evolution and equilibrium dynamics of the system for both disease free and endemic equilibrium cases [31–36]. One important parameter resulting from the SIR models is the reproduction number. The reproduction number of an infectious disease is defined as the expected number of secondary infections caused by a single infected individual for the whole duration that they are infectious [37, 38]. The reproduction number R0 describes how infectious a disease can be, and can also be used as a threshold parameter to determine whether a disease would survive in a healthy population. A value of R0 greater than one indicates the epidemic persists in the population [38]. The reproduction number of a virus is related to how fast the infection spreads in the population due to the contact between susceptible and infected individuals, which is described by the transmission rate coefficient β, and how fast the infected recover or are removed from the population, described by the removal rate coefficient γ. Measures to contain the disease, such as social distancing and isolation, quarantine, and closing of establishments to prevent interaction between individuals can be factored into β and γ [39]. Further details on the calculation of the reproduction will be discussed later in a separate section. The SIR model is preferred by some infectious disease modelers because of the low number of parameters that need to be estimated for the full model to be defined. However, this advantage comes from oversimplifying the model through relatively unrealistic assumptions such as the population being closed and homogeneous and having only three compartments [39, 40]. Even with these limitations, the SIR model can still provide basic estimates on whether the disease is expected to persist in the population or the proportion of the population expected to be infected by the disease, which would be helpful in development of public policy about medical response to the epidemic [39]. Its inherent simplicity also makes it easier for modelers to modify the models through addition of new compartments and stochasticity in the system to reflect aspects of epidemics such as immunity, incubation, and variation in individual movements [41]. Modifications of the SIR model have been used to describe mutations and changes in an infectious virus such as influenza. Yaari et al. [42] used a discrete time stochastic susceptible-infected-removed-susceptible (SIRS) model to describe influenza-like illnesses in Israel accounting for weather and antigenic drift by adding terms that account for weather signals and loss of immunity. Finkenstadt et al. [43] created a predictive stochastic SIRS model for weekly flu incidence accounting for antigenic drift. Roche et al. [44] used an agent-based approach based on the SIR compartment model to model the spread of a multi-strain epidemic, while Shi et al. [45] used the same approach and empirical data from Georgia, USA to model an influenza pandemic that incorporates viral mutation and seasonality. However, these approaches have been stochastic in nature, which does not provide information regarding the stability and existence of equilibrium points in an infected population. The aforementioned articles also do not take vaccination and the presence of other strains into account in their models. Consequently, one can model the presence of a mutated virus spreading into a population using a multi-strain model, which was used in the following studies for avian flu [32, 33]. These models study the birds and the humans as one population in an SI-SIR model. However, the two infected compartments in this model do not cross since they are separated by species. Casagrandi et al. [34] introduced a non-linear deterministic SIRC epidemic model to represent the antigenic drift for the influenza A virus. The SIRC model is a modified SIR model with an additional compartment, C, for individuals that receive partial immunity from being infected by one of the present strains. Although able to account for cross-immunity between strains, the model does not include the effect of vaccination into the system. Papers which have considered vaccination only consider one strain propagating within the population [46–49]. There has been very few studies that investigate the effect of vaccination in the presence of multiple strains like Wilson et al. did for Hepatitis B [50], which did not investigate the equilibrium model in detail. In the case of the influenza virus, it is possible to have multiple strains exist in a population, but only have vaccine for a certain strain that will not be effective for others. The fact that viruses undergo changes regularly indicates that people who have recovered from the virus, as well as individuals who have been vaccinated for a specific strain of the virus, can be susceptible again to a newly-emerged strain. It is important to determine the conditions in which a newly emerged strain and a common strain that has a means of immunity will coexist in a population provided that the two strains have a common subset for their susceptible pools. From a modeling standpoint, a highly infectious emergent strain can infect the susceptible population before the original strain which can impede the spread of the original strain or the two strains can coexist in an endemic equilibrium. An apt example for emerging disease that fits this description is the emergence of the COVID-19 virus in 2019 [51, 52]. As of August 31, 2020, there have been approximately 25 million confirmed cases of COVID-19 worldwide that has led to approximately 844,000 deaths since it was declared as an outbreak in January 2020 according to the WHO situation report [53]. At the time that this paper is being written, there are papers that have modeled the dynamics of the virus using different modifications of the SIR model. Zhou et. al. [54] included compartments corresponding to suspected cases, which consists of the individuals that show similar symptoms but are not confirmed cases, and indirectly infected individuals. Pan et. al. [55] used a modified SEIR model which included asymptomatic and treatment compartments for occurrences in Wuhan, China, the city where the outbreak started, and outside of Wuhan. Maier and Brockmann [56] included a separate compartment for quarantined individuals in the SIR model to account for the containment measures applied by the public for the virus. They then estimated the reproduction number of COVID-19 in different locations in China. He, Peng, and Sun [57] used the particle swarm optimization to approximate the parameters of the SEIR model with additional compartments for quarantined and hospitalized individuals. Similar models that are specific for each country/region have been formulated since the pandemic has spread worldwide [58–62]. It is notable that this virus emerged during the flu season [63] and had managed to infect a large number of individuals around the world in such a short time even when the threat of the influenza virus still exists. The CDC has recommended getting the flu vaccine for the incoming flu season to reduce the risk of getting the flu even if the vaccine will be ineffective against COVID-19 [64]. Although both these viruses have been modeled individually, there has been little to no work done in modeling the dynamics of COVID-19 and the influenza virus coexisting within a population where vaccination for the influenza virus is an option. This paper introduces a model that approaches the lack of cross-immunity across different viral strains by introducing new compartments to the SIR with vaccination model. This paper will give researchers insight about the conditions in which one strain can dominate another or if two different strains can coexist in a population, given that one of these strains has a vaccine available. This enables us to introduce acquired immunity through vaccination and cross-immunity between strains in a simple compartmental model and investigate the existence and stability of the resulting equilibrium points. We aim to model the equilibrium dynamics of an epidemic at the population level where a new emergent strain of an existing virus affects a closed population. The existing virus will be modeled using a modified SIR model with vaccination, however we assume that the vaccine does not provide immunity to the newer strain. The equilibrium points were determined for the system based on the transition equations and local stability was investigated for each point. Once the stability conditions have been established, the epidemic model was simulated using R [65] to investigate the steady-state behavior of the surveillance data for each compartment of the population. The values for the transmission and removal coefficients were dictated by the existence and stability conditions for each equilibrium point during the simulation. The reproduction number for the epidemic was also determined for this modified SIR epidemic model and compared to existing SIR models. Modeling the emergence of the new strain This section describes how the emergence of the new strain of the virus will be incorporated into the model. This emergence can either be due to mutation, antigenic drift/shift, or an introduction of a different strain from an external source. Let us assume that initially, there is only one strain of the virus that exists in the population. Immunity can be achieved either by recovering from the infection or getting vaccinated. After equilibrium has been established with the original strain, the new strain is introduced to the population. In addition to the individuals in the susceptible compartment, the new strain can affect individuals previously infected by the original strain and those who are vaccinated against the original strain; the only way to be immune to the mutated strain is to recover from the infection of the new strain. This model focuses on how the disease persists in the population in the long run without additional intervention aside from the preexisting vaccination for the original strain. The analysis will be focused on the macroscopic behavior of each compartment, which means that the population will not be studied at the individual level. The next two subsections will explain the dynamics before and after the emergence of the mutated strain. Before emergence The system begins as a population exposed to the original strain of the virus. The spread of the virus is described by a modified SIR model that accounts for vaccination [46]. The vaccinated members of the population can be treated as members of an additional compartment that do not interact with the infected individuals. This means that the modified SIR model will have four compartments instead of three, which are given by: Susceptible S: Individuals in this compartment are healthy, but are susceptible to be infected by the disease since they are not vaccinated. Vaccinated V: Individuals that were given a vaccine, making them immune to the disease. This also includes individuals with natural immunity to the disease. Infected I1: Individuals that are infected by the disease Removed R: Individuals that were infected but are now immune to the disease upon recovery. Because of their immunity, the members of this compartment do not interact with the remaining compartments. Let S, V, I1, and R be the respective number of individuals in the susceptible, vaccinated, infected, and removed compartments. The transition between the compartments is summarized by the compartmental diagram shown in Fig 1. A list of variables that explains each parameter used in describing the transitions between compartment can be found in S1 Table. 10.1371/journal.pone.0243408.g001Fig 1 Compartment diagram with transitions for the SIR with vaccination model. The arrows show the transitions between the compartments, as well as the exits due to natural death. The transition rates are shown next to the arrows. For this model, let μ be the natural birth rate of the population, and consequently the natural death rate of the population to keep the population size constant. It is assumed that the individuals are vaccinated at birth with a vaccination rate p. β is the standard incidence transmission coefficient, which assumes that the infection occurs based on how many susceptible individuals interact with the infected [66]. For standard incidence, the rate at which infected and susceptible individuals interact (also known as contact rate) is constant over all infected individuals regardless of the population size [67]. The removal rate coefficient for the infected individuals is denoted by γ. β and γ serve the same purpose as a rate constant in chemical kinetics. The dynamics of the system is described by ordinary differential equations that describe the rate of change in individuals belonging to a specific compartment. For any number of individuals in a general compartment C the rate of change of membership in the compartment can be expressed by the following equation: dCdt=(rateofinput)-(rateofoutput),(1) where the rates are denoted in the compartment diagram shown in Fig 1. For the susceptible compartment S the rate of increase comes from the birth of new members of the population who are not vaccinated, which is given by (1 − p)μN. Meanwhile, susceptible individuals can either get infected at a rate of βSI1N or die due to natural causes at a rate μS. In equation form, this translates to dSdt=(1-p)μN-βSI1N-μS(2) For the infected compartment I1, the number of infected individuals increase when susceptible individuals get infected at a rate βSI1N. The infected individuals can either die at a rate of μI1 or get removed and not interact with the system again at a rate γI1 when they recover. Translating this to an ordinary differential equation yields dIdt=βSI1N-γI1-μI1(3) Unlike the members of the susceptible compartment, the individuals in the vaccinated compartment will not get infected by the virus. This implies that the changes in the number of vaccinated individuals can only be due to the rate of vaccination in the population and death due to natural causes. Using a similar approach, the rate equation for the vaccinated compartment V is given by dVdt=pμN-μV(4) Since the population is closed, the number of individuals in the removed compartment can be expressed as R = N − S − I1 − V. This implies that solving Eqs 2–4 is enough to describe the system completely at any time t. Without loss of generality, Eqs 2–4 can be normalized with respect to the total population, N, so that the equations would be invariant to population scaling. This yields the following equations, dsdt=(1-p)μ-βsi1-μs(5) di1dt=βsi1-(γ+μ)i1(6) dvdt=(p)μ-μv(7) where (s, i1, v) = (S/N, I1/N, V/N) and r = R/N = N(1 − s − i1 − v) are functions of time t. Note that plausible solutions only exist when s(t), i(t), v(t), r(t) ≥ 0. To achieve equilibrium, there should not be any changes in the proportion for each compartment, which implies that Eqs 5–7 should be zero. If Eq 6 is zero, then two conditions emerge: i1 = 0 or i1 ≠ 0. The first case corresponds to the disease-free equilibrium (DFE) point (s(t), i1(t), v(t)) = (1 − p, 0, p). The latter case corresponds to the endemic equilibrium point (s*,i1*,v*). If i1 ≠ 0, then the following condition should be satisfied, βs*-(γ+μ)=0(8) Solving for s*, we get s*=(γ+μ)β(9) We can use this result to solve for i1* in Eq 6. The resulting endemic equilibrium point (s*,i1*,v*) is given by, (s*,i1*,v*)=(γ+μβ,μ[β(1-p)-γ-μ]β(γ+μ),p)(10) According to Chauhan et. al [46], the reproduction number of the disease for the vaccinated SIR model is given by Rv = R0(1 − p) = β(1 − p)/(γ + μ). This means that the endemic equilibrium point will be asymptotically stable if Rv > 1, while the DFE will be asymptotically stable if Rv < 1. The reproduction number will be discussed in the Section ‘‘Reproduction Number’’. Emergence of new strain Suppose that at time T when equilibrium has been achieved in the SIR with vaccination model, a new strain of the disease is introduced to the population. This new strain will have a different transmission coefficient β′ and removal rate coefficient γ′. This results to the existence of another compartment I2 for those who are infected with the new strain, which we will refer to as Disease 2. The existence of the newer strain will be constrained by the following assumptions: Since the vaccine is assumed to only work on the original strain, the vaccinated and the previously removed individuals are susceptible to the newer strain. Once infected by the newer strain, the individual cannot be infected by the original strain. The individuals infected by the newer strain will be removed from the population or die. For the case of COVID-19, we can interpret this removal from the population as isolation from the compartments susceptible to the second infection. Individuals infected by the original disease have to be removed first before being susceptible to the newer strain; meaning that there is no chance of super-infection (I1 → I2) [37]. Although there are cases of co-infection, these are rare cases and are usually diagnosed after they were initially removed from the population [68]. This means that the number of compartments that need to be monitored will increase from four to six, with the addition or modification of the following compartments: R1: Individuals who have recovered from the original strain but are now susceptible to the newer strain. I2: Individuals who are infected by the newer strain. R2: Individuals who were previously infected by the newer strain but have now been removed due to recovery or treatment. The members of the vaccinated compartment, which was initially an isolated compartment, can now be infected by the new strain. The same can be said for the individuals who have recovered from the original strain. For mathematical simplicity, the infection coefficients for the new strain are assumed to be the same for the susceptible, vaccinated, and initially recovered compartments. These assumptions and the increase in number of compartments also introduces the possibility of new transitions between compartments as shown in Fig 2. A list of variables that explains each parameter used in describing the transitions between compartment can be found in S1 Table. 10.1371/journal.pone.0243408.g002Fig 2 Compartment diagram for the emerging disease model. The transitions between compartments, together with the corresponding rates, are described by the arrows directed in and out of each compartment. As in Chauhan et. al’s [46] work, the standard incidence was used to model infection of the susceptible individuals for the newer strain. Based on Eq 1 and the compartment diagram in Fig 2, the dynamics of the system can be expressed in terms of the following ordinary differential equations: dsdt=(1-p)μ-βsi1-β′si2-μs(11) di1dt=βsi1-(γ+μ)i1(12) dvdt=(p)μ-β′vi2-μv(13) dr1dt=γi1-β′ri2-μr(14) di2dt=β′(s+v+r)i2-(γ′+μ)i2(15) and r2 = 1 − s − i1 − v − r1 − i2. Similar to the simple SIR with vaccination scenario, the solution for the variables should follow the constraint s(t), i1(t), v(t), i2(t), r1(t), r2(t) ≥ 0 for any time t. To solve for the equilibrium points of the system, Eqs 11–15 should be equal to zero. Wolfram Mathematica [69] was used to obtain solutions for the system of equations, which are the following: DFE: (s, i1, v, r1, i2) = (1 − p, 0, p, 0, 0) Original strain equilibrium: (s,i1,v,r1,i2)=(γ+μβ,μ[β(1-p)-γ-μ]β(γ+μ),p,γ[β(1-p)-γ-μ]β(γ+μ),0) New strain equilibrium: (s,i1,v,r1,i2)=((1-p)(γ′+μ)β′,0,p(μ+γ′)β′,0,μ[β′-(γ′+μ)]β′(γ′+μ)) Endemic equilibrium: (s, i1, v, r1, i2) = (s*, i1*, v*, r1*, i2*) The second equilibrium point corresponds to the scenario where only the original strain is present. Applying the constraint for the plausible solution, the original strain equilibrium exists if β(1-p)-γ-μ>0→Rv>1(16) where Rv = R0(1 − p) = β(1 − p)/(γ + μ) is the reproduction number of the original strain for the SIR model with vaccination [46]. The third equilibrium corresponds to the scenario where only the new strain survives. For this equilibrium point to exist, the following condition should be satisfied: β′-γ′-μ>0→R0′>1(17) where R0′=β′/(γ′+μ) is the corresponding reproduction number of the newer strain if modeled using a standard SIR model. Equilibrium point 4 is the endemic equilibrium where s*=γ+μβ(18) i1*=μ[β(1-p)(γ′+μ)-β′(γ+μ)](γ+μ)[β(γ′+μ)-μβ′](19) v*=p(γ+μ)[β(γ′+μ)-μβ′]ββ′(γ+pμ)(20) r*=γ[β(1-p)(γ′+μ)-β′(γ+μ)]ββ′(γ+pμ)(21) i2*=μ[(γ+μ)(μβ′-β(γ′+μ))+ββ′(γ+pμ)](γ+μ)β′[β(γ′+μ)-μβ′](22) For the endemic equilibrium to exist, the following condition should be satisfied: R0(1-p)>R0′(23) The next section will discuss the local stability of the four equilibrium points. Stability analysis and simulations After solving for the equilibrium points, we need to determine the conditions in which these points are stable. These conditions dictate which equilibrium will describe the steady state behavior of the system. The local stability of the equilibrium points will be determined based on the eigenvalues of its Jacobian evaluated at a specific equilibrium point [25]. Let C0 = (C1, C2, …)T be the vector of the population number of each compartment. For a general compartment Ci, the components of the Jacobian, Jij can be obtained using the following equation: (Jij)|C=C0=∂∂Ci(dCjdt)|C=C0(24) For our system, the Jacobian of the system can be obtained by applying Eq 24 to Eqs 11–15. For any equilibrium point, (s¯,i¯1,v¯,r¯,i¯2) yields J=(-βi¯1-β′i¯2-μ-βs¯00-β′s¯βi¯1βs¯-(γ+μ)00000-β′i¯1-μ0-β′v¯0γ0-βi¯2-μ-β′r¯β′i¯20β′i¯2β′i¯2β′(s¯+v¯+r¯)-(γ′+μ)) Local stability is attained when the eigenvalues of the Jacobian, λ, are negative or have negative real parts. In other words, the solutions for λ such that det(J-1λ)=0 should be negative or have negative real parts if the solution is complex [70]. The system is simulated to reach equilibrium. As described in the previous sections, the system starts as a one-strain SIR model with vaccination as discussed in the section “Modeling the emergence of the new strain” with the following values for the parameters: μ = 0.5 (birth/death rate) and p = 0.5 (vaccination rate). The values for μ and p were based on the simulations performed by Chauhan et. al. [46] for their stability analysis of the SIR with vaccination model. The vaccination rate of 0.5 is a close estimate to the overall vaccination coverage for the influenza virus in the United States during the 2018-2019 influenza season [71]. For this simulation, the time is discretized in units of the average time between compartment interactions, i.e. the average time it takes for individuals to transition from one compartment to another. At time t = 0, we allow 1% of the population to be infected by the original strain and the system is made to evolve in time using the values of infection coefficients β and removal rate γ that satisfies the respective requirements for the reproduction number for each equilibrium point to exist. The new (mutated) strain was made to emerge at time t = 100, when we expect the system to be in equilibrium. The new strain is introduced to the population by infecting 1% of the susceptible population with the newer strain. The evolution will then be dictated by the modified multi-strain SIR model developed in section “Modeling the emergence of the new strain” using values for β′ and γ′ that satisfy the conditions for R0′ for each equilibrium point to exist. Disease free equilibrium (DFE) The Jacobian for the DFE can be obtained by substituting the respective values to (s¯,i¯1,v¯,r¯,i¯2) in the expression for the Jacobian. This yields: J=(-μ-β(1-p)00-β′(1-p)0β(1-p)-(γ+μ)00000-μ0-β′p0γ0-μ00000β′-(γ′+μ))(25) and the corresponding characteristic equation is (λ − μ)3(λ − β(1 − p) + γ + μ)(λ − β′ + γ′ + μ) = 0. This means that the eigenvalues are λ = −μ, −μ, −μ, β(1 − p) − γ − μ, β′ − γ′ − μ. Recall that for the DFE to be locally asymptotically stable, all eigenvalues should have negative real parts. Hence, the following conditions should hold: β(1-p)γ+μ=R0(1-p)<1(26) β′γ′+μ=R0′<1(27) This is consistent with the local stability of the disease free equilibrium for the regular standard incidence SIR and the SIR with vaccination models [46]. The conditions stated in Eqs 26 and 27 give the threshold conditions of the transmission and removal coefficients for this model, which is related to the basic reproduction number of the model to be discussed in “Reproduction Number” section. Eqs 26 and 27 also imply that if the system is in DFE before the emergence and the reproduction number of the emergent disease is less than that of the original disease, then the system will remain in DFE in the long run. Fig 3 shows the simulation of the DFE using parameters that satisfy Eqs 26 and 27. The original strain was simulated to have a reproductive number of 0.80, while the new strain was simulated to have a reproductive number of 0.57. Note that the plot legends will be the same for the succeeding surveillance plots for the other equilibrium plots. The plot shows that the proportion of vaccinated individuals (v), denoted by the solid blue line, and the proportion of susceptible individuals (s), denoted by the solid black line, remained relatively constant at long times. Since the vaccination rate is set to be 0.7, we expect more individuals to be vaccinated than susceptible. Due to the emergence of the new strain at time point t = 100, there appears to be a slight dip in s but it quickly stabilized to the disease free equilibrium. 10.1371/journal.pone.0243408.g003Fig 3 Surveillance data of the compartments for the disease free equilibrium. The reproduction number of the original strain is R0 = 0.8, while the reproduction number of the emergent strain is R0′=0.57. The vaccination rate used is 0.5. Disease 1 equilibrium (original strain) For this equilibrium scenario, i1 ≠ 0 and since Eq 12 is equal to zero then βs1-(γ+μ)=0(28) where s1 is the equilibrium value corresponding to the susceptible compartment for the Disease 1 equilibrium. We can calculate the resulting Jacobian for Equilibrium Point 2 by substituting the corresponding values to the Jacobian equation given in Eq 24. The Jacobian is then given by, (-μ[β(1-p)(γ+μ)-(γ+μ)00-β′(γ+μ)βμ[β(1-p)-γ-μ](γ+μ)000000-μ0-β′p0γ0-μγ[β(1-p)-γ-μ](γ+μ)0000β′(s1+v1+r1)-γ′-μ) where s1, v1, r1 are the respective equilibrium values for the susceptible, vaccinated, and the initially recovered compartments for the Disease 1 equilibrium. The resulting eigenvalues are (-μ,-μ,β′[μβ+γ+μpγ+μ]-(γ′+μ),λ±), where λ±=-(12)[μβ(1-p)γ+μ±(μβ(1-p)γ+μ)2-4μ(β(1-p)-γ-μ)](29) The discriminant of λ± can dictate whether the eigenvalues will have a negative real part. If the discriminant is negative or zero, then the eigenvalues will be negative. If the discriminant is positive, recall that for the system to not go to the DFE, R0(1 − p) > 1. This means that [μβ(1-p)γ+μ]2>[μβ(1-p)γ+μ]2-4μ(β(1-p)-γ-μ)(30) which ensures that λ± is negative when R0(1 − p) > 1. When R0(1 − p) < 1, λ± would have a positive real part which makes the equilibrium point unstable. This suggests that this equilibrium point will not be stable if the system was not already in endemic equilibrium with Disease 1. For the third eigenvalue to be negative, R0′1, which means that when the discriminant is positive, the following inequality holds: [μβ′γ+μ]2>[μβ′γ+μ]2-4μ(β′-γ′-μ)(35) which means that both λ± will be negative as long as R0′>1. As for the remaining eigenvalue, it will be negative if R0′>R0(1-p)(36) This indicates that the second disease will be locally stable if the mutated disease has a higher reproduction number compared to the original disease. Once this happens the endemic equilibrium can not be achieved. Note that there are no conditions for the value of R0(1 − p), which means that this equilibrium point can occur whether the system was initially in DFE or endemic equilibrium with the original strain. The emergence from a system in DFE is shown in the simulation in Fig 5 where i1 is zero before the emergence of the strain, which happens when R0′>1>R0(1-p). For this simulation, the reproduction number of the original strain is 0.8 while the emergent strain has a reproduction number of 3.33. The proportion of vaccinated and the susceptible individuals remained constant before the emergence. Upon the emergence of the newer strain, the system behaves like a regular SIR model with a susceptible compartment comprising of the S, V, and R compartments as shown in Fig 5. The proportion of both susceptible and vaccinated individuals decreased drastically shortly after the emergence at t = 100, while the proportion of individuals infected by the emergent strain (i2), denoted by the pink dashed line, had an upward spike before settling into its equilibrium value. The original strain showed no sign of reemergence after it has settled into DFE, which is what is expected. 10.1371/journal.pone.0243408.g005Fig 5 The reproduction number of the original strain is R0 = 0.8, while the reproduction number of the emergent strain is R0′=3.33. This scenario corresponds to the emergence of the new strain under DFE conditions. The vaccination rate used is 0.5. Another possible case is when the system is initially in endemic equilibrium but the original strain dies because of the introduction of the new strain. Fig 6 shows that before time t = 100 and i1 is nonzero, which occurs when R0′>R0(1-p)>1. The reproduction number of the original strain for this simulation is R0 = 3, while the reproduction number of the emergent strain is R0′=3.33. Upon emergence of the second strain however, the proportion of the population infected by the original strain, i1, and the proportion of the individuals who recovered from the original strain, r1, decrease and go to zero asymptotically. The behavior of the second strain is similar to that in Fig 5. This implies that the new emergent strain is much more infectious than the older strain and the new strain infects more susceptible individuals compared to the original strain. This causes the original strain to die down at steady state. 10.1371/journal.pone.0243408.g006Fig 6 Surveillance data of the compartments for the new strain equilibrium where the system is originally in endemic equilibrium with Disease 1. The reproduction number of the original strain is R0 = 3, while the reproduction number of the emergent strain is R0′=3.33. The vaccination rate used is 0.5. Endemic equilibrium For the endemic equilibrium case, both i1 and i2 are nonzero and thus both Eqs 28 and 32 hold. Since Eq 11 is zero, βi1*+β′i2*-μ=-(1-p)/s*(37) The resulting Jacobian for the endemic equilibrium case is given by J=(-(1-p)μ/s*-βs*00-β′s*βi1*000000-β′i2*-μ0-β′v*0γ0-βi2*-μ-β′r1*β′i20β′i2β′i20)(38) The characteristic equation is given by (λ+pμ/v*)(λ4+a1λ3+a2λ2+a3λ+a4)=0(39) where a1=ps*μ+v*(1-p)μs*v*(40) a2=i1*s*v*β(γ+μ)+μ2p(1-p)+i2*s*v*β′2(s*+r*+v*)s*v*(41) a3=i1*ps*βμ(γ+μ)+i2p(s*)2μβ′2+(1-p)i2*vβ′2μ(r*+v*)s*v*(42) a4=i2*i1*ββ′2[(r*s*v*+s*(v*)2)(γ+μ)+(s*)2v*γ]s*v*(43) Recall that for the endemic equilibrium to be locally stable, all eigenvalues should have a negative real part. Eq 39 shows that one of the eigenvalues is λ = −pμ/v* which is negative. For all the roots of the quartic term to have negative real parts, the Routh-Hurwitz criteria for stability should be applied [24, 25]. According to the Routh-Hurwitz criterion, a polynomial with degree 4 will have roots (a1, a2, a3, a4) that all have negative real parts when: a1,a3,a4>0a1a2a3>a32+a12a4. Since the values of the equilibrium points should be positive, then all coefficients (a1, a2, a3, a4) are positive. Based on the stability of the first three equilibrium points and the existence criterion, the endemic equilibrium point is expected to be stable when R0[γ+μμ+R0(γ+pμ)]