key: cord-0921556-1pojgpyt
authors: Fudolig, Miguel; Howard, Reka
title: The local stability of a modified multi-strain SIR model for emerging viral strains
date: 2020-12-09
journal: PLoS One
DOI: 10.1371/journal.pone.0243408
sha: b89d49760553a6346afee3faa5cbdf101098409d
doc_id: 921556
cord_uid: 1pojgpyt

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.

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] [2] [3] . The drop in numbers in vaccination has led to outbreaks of diseases such as mumps [4] [5] [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] [8] [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] [11] [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 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 reassortment which can make it feasible to create a more virulent strain than the original strain [12] [13] [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] [16] [17] , Bayesian inference using Monte Carlo methods to estimate parameters through simulation based on the stochastic model [18] [19] [20] , or agent-based network dynamics to model the spread of the disease through interactions between individuals [21] [22] [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] [25] [26] [27] [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] [32] [33] [34] [35] [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 R 0 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 R 0 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 susceptibleinfected-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] [47] [48] [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] [59] [60] [61] [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 crossimmunity 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.

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.

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 I 1 : 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, I 1 , 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. 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:

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 bSI 1 N or die due to natural causes at a rate μS. In equation form, this translates to

For the infected compartment I 1 , the number of infected individuals increase when susceptible individuals get infected at a rate bSI 1 N . The infected individuals can either die at a rate of μI 1 or get removed and not interact with the system again at a rate γI 1 when they recover. Translating this to an ordinary differential equation yields

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

Since the population is closed, the number of individuals in the removed compartment can be expressed as R = N − S − I 1 − 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,

where (s,

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:

. The latter case corresponds to the endemic equilibrium point

Solving for s � , we get

We can use this result to solve for i � 1 in Eq 6. The resulting endemic equilibrium point

According to Chauhan et. al [46] , the reproduction number of the disease for the vaccinated SIR model is given by

. This means that the endemic equilibrium point will be asymptotically stable if R v > 1, while the DFE will be asymptotically stable if R v < 1. The reproduction number will be discussed in the Section ''Reproduction Number''.

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 β 0 and removal rate coefficient γ 0 . This results to the existence of another compartment I 2 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 (

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:

• R 1 : Individuals who have recovered from the original strain but are now susceptible to the newer strain.

• I 2 : Individuals who are infected by the newer strain.

• R 2 : 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. 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:

and

Similar to the simple SIR with vaccination scenario, the solution for the variables should follow the constraint s(t), i 1 (t), v(t), i 2 (t), r 1 (t), r 2 (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: 

; m½bð1 À pÞ À g À m� bðg þ mÞ ; p; g½bð1 À pÞ À g À m� bðg þ mÞ ; 0 � � 3. New strain equilibrium:

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

where

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:

where R 0 0 ¼ b 0 =ðg 0 þ mÞ is the corresponding reproduction number of the newer strain if modeled using a standard SIR model. Equilibrium point 4 is the endemic equilibrium where

For the endemic equilibrium to exist, the following condition should be satisfied:

The next section will discuss the local stability of the four equilibrium points.

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 C 0 = (C 1 , C 2 , . . .) T be the vector of the population number of each compartment. For a general compartment C i , the components of the Jacobian, J ij can be obtained using the following equation:

For our system, the Jacobian of the system can be obtained by applying Eq 24 to Eqs 11-15. For any equilibrium point, ð� s;

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 À 1lÞ ¼ 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 β 0 and γ 0 that satisfy the conditions for R 0 0 for each equilibrium point to exist.

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:

and the corresponding characteristic equation is (λ − μ) 3 (λ − β(1 − p) + γ + μ)(λ − β 0 + γ 0 + μ) = 0. This means that the eigenvalues are λ = −μ, −μ, −μ,

Recall that for the DFE to be locally asymptotically stable, all eigenvalues should have negative real parts.

Hence, the following conditions should hold:

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 

Modified multi-strain SIR model for emerging viral strains 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.

For this equilibrium scenario, i 1 6 ¼ 0 and since Eq 12 is equal to zero then

where s 1 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,

where s 1 , v 1 , r 1 are the respective equilibrium values for the susceptible, vaccinated, and the initially recovered compartments for the Disease 1 equilibrium. The resulting eigenvalues are

ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi mbð1 À pÞ g þ m � � 2 À 4mðbð1 À pÞ À g À mÞ

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, R 0 (1 − p) > 1. This means that

which ensures that λ ± is negative when R 0 (1 − p) > 1. When R 0 (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,

When these conditions are satisfied, then the Disease 1 equilibrium point is locally asymptotically stable. These conditions also imply that this equilibrium point can only be achieved if the system before the emergence is already in endemic equilibrium with Disease 1. Fig 4 shows the simulation of the system using parameters that satisfy Eq 31. In this case, the emergent strain has a reproduction number of 1.18 while the original strain has a reproduction number of 3.33. Unlike the DFE case, the proportion of individuals infected by the original strain, denoted by the red dashed line, along with the recovered individuals, denoted by the green dashed line, increase up to their respective equilibrium values. This is accompanied by the decrease in susceptible individuals in the population.

At t = 100, the number of individuals infected by the newer strain, denoted by the pink dashed line, was shown to have a small spike, but quickly went to zero while the number of individuals infected by the original strain remained relatively unchanged. Note that the reproduction number of the emergent strain is 1.11, which is greater than one, which means that the emergent strain will be able to survive on its own in this population. This implies that the new strain is not strong enough to infect enough people to dominate over the original strain.

For this equilibrium scenario, i 2 6 ¼ 0 and thus for Eq 15 to be zero, we have to have

where � s 2 , � v 2 , and � r 2 are the respective equilibrium values of the susceptible, vaccinated, and initially recovered compartments corresponding to the Disease 2 equilibrium. The resulting Jacobian, J, for the equilibrium where only the mutated disease exists is given by

The corresponding characteristic equation is given by,

The eigenvalues are À mb 0 g þ m ;

ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi mb 0 g 0 þ m

Similar to the λ ± in the Disease 1 equilibrium, the real part will be negative if the discriminant is negative. For the equilibrium point to exist, R 0 0 > 1, which means that when the discriminant is positive, the following inequality holds:

which means that both λ ± will be negative as long as R 0 0 > 1. As for the remaining eigenvalue, it will be negative if

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 R 0 (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 i 1 is zero before the emergence of the strain, which happens when R 0 0 > 1 > R 0 ð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 (i 2 ), 

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.

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 i 1 is nonzero, which occurs when R 0 0 > R 0 ð1 À pÞ > 1. The reproduction number of the original strain for this simulation is R 0 = 3, while the reproduction number of the emergent strain is R 0 0 ¼ 3:33. Upon emergence of the second strain however, the proportion of the population infected by the original strain, i 1 , and the proportion of the individuals who recovered from the original strain, r 1 , 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. 

For the endemic equilibrium case, both i 1 and i 2 are nonzero and thus both Eqs 28 and 32 hold. Since Eq 11 is zero,

The resulting Jacobian for the endemic equilibrium case is given by

The characteristic equation is given by

where

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 (a 1 , a 2 , a 3 , a 4 ) that all have negative real parts when:

Since the values of the equilibrium points should be positive, then all coefficients (a 1 , a 2 , a 3 , a 4 ) 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

It is safe to say that the endemic equilibrium for the two strains can only occur when the system is initially in endemic equilibrium for the original strain, which can explain why the condition is more restrictive than the one for the equilibrium with Disease 2. This highly restrictive criterion for endemic equilibrium might be challenge for simulating stochastic data for strains that have a relatively low reproduction number and a high vaccination rate. A fluctuating value for the incidence rate might lead to either one of the single-strain equilibria. Fig 7 shows that both i 1 and i 2 are asymptotically nonzero after the emergence of the new strain given that Eq 44 is satisfied. For this simulation, the reproduction number of the original strain is R 0 = 3.33 while the reproduction number of the emergent strain is R 0 0 ¼ 1:33. As the endemic equilibrium between the two strains is reached, the proportion of the vaccinated individuals that are healthy decreased considerably compared to the proportion before the emergence while unvaccinated individuals return to the same proportion after the system has stabilized after the emergence of the new strain. This implies that the newer strain mostly survives on infecting the vaccinated and the initially recovered individuals. The proportion of individuals infected by the original strain is also observed to decrease upon reaching endemic equilibrium after emergence, which implies that some of the susceptible population get infected by the newer strain. 

One of the important parameters to be calculated for an epidemic model is the reproduction number, which quantifies how infectious a certain disease. Formally, the reproduction number is defined as the expected number of secondary infections caused by a single infected individual for the whole duration that they are infectious. A value for the reproduction number that is greater than one indicates that the epidemic persists in the population, while a value less than one means that the disease will die out in the population [37, 38] . Since the values are dependent on how fast the disease is transmitted and how fast infected individuals are removed from the population, estimates for the reproductive number can vary for different countries and locations within countries. The reproduction number of different influenza viruses were compiled by Biggerstaff et. al. [72] , while other studies focus on the estimation of the reproduction number for specific locations in Hong Kong [73] and Japan [74] . Park et. al. [75] provided a range for the estimated reproduction number of COVID-19 from various published and preprint papers that used data from different locations around the world. To provide the reader with relative values for the reproductive number, the 2009 A/H1N1 influenza virus is estimated to have a reproduction number from 1.30-1.70 [72] , while the estimates for the reproduction number of COVID-19 for various locations ranged from 1.9-6.5 [75] .

The reproduction number R 0 of this epidemic model was calculated using the approach formulated by van den Diessche and Watmough [37] . This approach does not account for any measures taken to control the epidemic, but will give us an idea of the conditions needed for the disease to spread on its own.

To obtain the reproduction number for this model, we need to solve for the next generation matrix. The next generation matrix describes the expected number of new infections that an infected individual produces from each susceptible compartment. Eqs 11 to 15 can be written in vector form as

where X ¼ ði 1 ; i 2 ; s; v; r 1 Þ T , F i corresponds to the vector that describes the rate of new infections in compartment i, and V i is the vector that correspond to the transitions from compartment i to the other non-infected compartments such as R 1 and R 2 [76]. Define F and V such that for the disease free equilibrium X 0 ,

where (i, j) corresponds to the index of the infected compartments. The resulting matrices F and V are m × m matrices, where m is the number of infected compartments. F ij describes the rate at which the infected individuals at compartment j contribute to the infection of compartment i, while V ij corresponds to the rate at which the infected individuals are removed from the infected compartments. This means that FV −1 is related to the rate at which individuals are infected by the disease within an average time span that an infected individual remains infected. For the system discussed in this paper, there are two infected compartments after the emergence of the new strain: I 1 and I 2 . Therefore, m = 2 and 1 � i, j � m.

The DFE was calculated to be given by (1 − p, 0, p, 0, 0) . Based on the transition equations (Eqs 11 to 15), F and V are given by

Note that m = 2, so the corresponding F and V matrices are then,

The reproduction number is obtained by taking the maximum eigenvalue of the next generation matrix FV −1 . The next generation matrix is the product of the rate of infection (F) and the average time that an individual remains infected (V −1 ). The next generation matrix is given by,

where F is the matrix that describes the infection rates for the two infections at the DFE and V −1 describes the average time an infected individual stays infectious. It is easy to see that the eigenvalues of the next generation matrix are R 0 (1 − p) and R 0 0 . This means that the reproduction number R 0 for this system is given by the larger of the two. Formally,

Note that the resulting threshold equation for the system is R 0 < 1, which means that the system will only approach the disease free equilibrium when Eq 53 is less than one. For an outbreak to occur, at least one of these two strains should be able to persist in the population on its own, that is, to have an individual reproduction number greater than one. This is consistent with Eqs 26 and 27 which give us the condition of the stability of the DFE.

We started by modeling the emergence of a new strain by adding and modifying compartments to the existing SIR model with vaccination. The emergent strain was assumed to be unaffected by the existing vaccine designed for the original strain in the population. After establishing the possible transitions between compartments, the system was found to have four equilibrium points: the disease free equilibrium, the existing strain equilibrium, the emergent strain equilibrium, and the endemic equilibrium. Fig 8 shows the resulting equilibrium state of the population for reproduction numbers within the estimated ranges for influenza and COVID-19 [72, 75] . For illustrative purposes, the values used for the transmission rate coefficients (β = β 0 = 2) were based on the value used in the simulations performed Chauhan et. al. [46] . Upon examining the conditions for existence and local stability, the disease-free equilibrium was determined to be locally stable when both R v and R 0 0 are less than one. This is consistent with the reproduction number for this multi-strain model. The existing strain equilibrium surprisingly does not impose that R 0 0 < 1, which is the condition for the DFE of a normal SIR model for a single-strain without immunity. This implies that if the original strain has a higher reproduction number than the emergent strain, the emergent strain would still die out eventually even if it would persist in the population if it was on its own. On a similar note, the local stability condition for the emergent equilibrium condition also does not impose that R v < 1. This means that the state of the system prior to the emergence of the new strain does not matter as long as the emergent strain is more contagious than the existing strain. This leaves a highly restrictive condition for endemic equilibrium to exist: the emergent strain must be able to survive by itself and the original strain must be contagious enough to infect enough people, which is given by Eq 44. Eq 44 also implies that the two-strain endemic equilibrium will only be locally asymptotically stable if the original strain is endemic to the population upon the emergence of the newer strain. These restrictions for the stability of endemic equilibrium would highly affect simulation studies about emergent strains especially when stochasticity is added to the model. Knowing the stability conditions for a multi-strain SIR model without cross-immunity will also be able to give us insights about when a newer strain emerges while the original strain still exists. One of the advantages of this model is that there was no assumption of any specific strain of influenza virus as an original strain, as long as the reproduction number and the vaccination assumption is satisfied. This is common for the flu virus which changes every season and the vaccines lose their efficacy after a new mutated strain emerges. This result might also be relevant when a highly contagious strain like the COVID-19 virus, a viral strain that does not exhibit cross-immunity with the influenza virus, emerges in a population. Our results show that COVID-19 will be able to co-exist with the flu in an endemic equilibrium under certain circumstances, which might become a problem in terms of prioritizing patients in health care centers especially since the two diseases show similar symptoms [77] . Our research also suggests that if the spread of the emergent virus like COVID-19 is not contained, the virus could infect more people than the existing influenza virus, which highly influences the development of protocol in receiving patients with flu-like symptoms and allocating resources in healthcare centers as we see in recent times.

In summary, the modified multi-strain SIR model of an emerging disease that affects both susceptible and previously immune individuals was studied. The local stability of the equilibrium points as well as the reproduction number for the model was calculated. Based on the results, it is found out that the original and the emergent strain can coexist in an endemic equilibrium if the emergent strain has a lower reproduction number than the original strain and that the system should already be in endemic equilibrium with the original disease before the emergence. The requirement for the endemic equilibrium to exist is quite strict especially for low values of R 0 and high values of p, which presents a challenge in simulating surveillance data with stochastic incidence rates. This modified SIR model can be improved further by using time-dependent infection rates to account for the environmental factors and the inherent characteristics of the virus. Common measures against COVID-19 could also be incorporated in the model by adding compartments or modifying transition probabilities between the existing compartments. As an example, self-quarantine and self-isolation can be analyzed as special transitions to the R 2 compartment because infected individuals that take these precautions effectively reduce their contact with the other members of the population. Exploring the effect of vaccination on other epidemic models such as epidemics with animal vectors and epidemics with latency and asymptomatic compartments can also be studied in the future.

Supporting information S1 Table. List of variables used in describing the emergent strain model. (PDF) 

Attitudes to vaccination: a critical review

A systematic review of factors affecting vaccine uptake in young children

Review of vaccine hesitancy: Rationale, remit and methods. Vaccine

Mumps outbreak at a university and recommendation for a third dose of measles-mumpsrubella vaccine-Illinois

Recent resurgence of mumps in the United States

Centers for Disease Control and Prevention. Mumps Cases and Outbreaks

Centers for Disease Control and Prevention. Measles Cases and Outbreaks

Measles in a Micronesian Community-King County

Centers for Disease Control and Prevention. How the Flu Virus Can Change

Vaccination and antigenic drift in influenza

Influenza vaccine: the challenge of antigenic drift

The evolution of epidemic influenza

Whole-genome analysis of human influenza A virus reveals multiple persistent lineages and reassortment among recent H3N2 viruses

Phenomenological Modelling of COVID-19 epidemics in Sri Lanka, Italy and Hebei Province of China. medRxiv

Using phenomenological models to characterize transmissibility and forecast patterns and final burden of Zika epidemics. PLoS currents

Using phenomenological models for forecasting the

A tutorial introduction to Bayesian inference for stochastic epidemic models using Markov chain Monte Carlo methods

Modelling and Bayesian analysis of the Abakaliki smallpox data

Statistical inference for multi-pathogen systems

The analysis of an epidemic model on networks

Modelling seasonal influenza: the role of weather and punctuated antigenic drift

Modelling antigenic drift in weekly flu incidence

An agent-based model to study the epidemiological and evolutionary dynamics of influenza viruses

Modelling seasonality and viral mutation to predict the course of an influenza pandemic

Stability analysis of SIR model with vaccination

On pulse vaccination strategy in the SIR epidemic model with vertical transmission

Stability properties of pulse vaccination strategy in SEIR epidemic model

Pulse vaccination strategy in the SIR epidemic model

Predictions of the emergence of vaccine-resistant hepatitis B in The Gambia using a mathematical model

Real-Time Estimation of the Risk of Death from Novel Coronavirus (COVID-19) Infection: Inference Using Exported Cases

Novel coronavirus 2019 (COVID-19): Emergence and implications for emergency care

World Health Organization. WHO Coronavirus Disease (COVID-19) Dashboard; 2020

The Outbreak Evaluation of COVID-19 in Wuhan District of China

Effectiveness of control strategies for Coronavirus Disease 2019: a SEIR dynamic modeling study. medRxiv

Effective containment explains sub-exponential growth in confirmed cases of recent COVID-19 outbreak in Mainland China

SEIR modeling of the COVID-19 and its dynamics

A modified SEIR model to predict the COVID-19 outbreak in Spain and Italy: simulating control scenarios and multi-scale epidemics. Available at SSRN 3576802

Regression Model based COVID-19 outbreak predictions in India

Epidemic analysis of COVID-19 in China by dynamical modeling

CoronaTracker: worldwide COVID-19 outbreak data analysis and prediction

Outbreak dynamics of COVID-19 in China and the United States. Biomechanics and modeling in mechanobiology

R: A Language and Environment for Statistical Computing

Four predator prey models with infectious diseases

Stability analysis and optimal vaccination of an SIR epidemic model

SARS-CoV-2 and influenza virus co-infection

Ordinary differential equations and stability theory: an introduction. Courier Corporation

Centers for Disease Control and Prevention. Flu Vaccination Coverage

Estimates of the reproduction number for seasonal, pandemic, and zoonotic influenza: a systematic review of the literature. BMC infectious diseases

The effective reproduction number of pandemic influenza: prospective estimation

Estimation of the effective reproduction number of influenza based on weekly reports in Miyazaki Prefecture

A systematic review of COVID-19 epidemiology based on current evidence

Reproduction numbers of infectious disease models

We would like to acknowledge Dr. Yu Jin and Dr. Clay Cressler for their respective for their inputs on mathematical derivation and the applications of the model. We would also like to acknowledge the two anonymous reviewers who have provided valuable input that led to the improvement of this paper.