key: cord-0886552-cycc3fvl
authors: Fudolig, Miguel; Howard, Reka
title: The local stability of a modified multi-strain SIR model for emerging viral strains
date: 2020-03-23
journal: nan
DOI: 10.1101/2020.03.19.20039198
sha: ada9208c7c35d9480a4ddeebbd9c9e9421bc407c
doc_id: 886552
cord_uid: cycc3fvl

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 equilibrium. 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 demonstrates 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 for some of its strains, 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.

More recently, the anti-vaccination movement has been gaining traction in different 2 parts of the world. Individuals who do not advocate vaccination commonly cite reasons 3 such as fear of adverse side-effects, perceived low efficacy of vaccines, and perceived low influenza-like illness in Israel accounting for weather and antigenic drift by adding terms 48 that account for weather signals and lost of immunity. Finkenstadt et al. [31] created a 49 predictive stochastic SIRS model for weekly flu incidence accounting for antigenic drift. 50 Roche et al. [32] used an agent-based SIR model to describe the spread of a multi-strain 51 epidemic, while Shi et al. [33] used the same approach and empirical data from Georgia, 52 USA to model an influenza pandemic that incorporates viral mutation and seasonality. 53 However, these approaches have been stochastic in nature, which does not provide 54 information regarding the stability and existence of equilibrium points in an infected 55 population. The abovementioned articles also do not take vaccination and the presence 56 of other strains into account in their models. Consequently, one can model the presence 57 of a mutated virus spreading into a population using a multi-strain model, which was 58 used in the following studies for avian flu [23, 24] . These models study the birds and the 59 humans as one population in an SI-SIR model, where the first SI corresponds to the 60 compartments for the avian species. However, the two infected compartments in this 61 model do not cross since they are separated by species. Casagrandi et al. [25] 62 introduced a non-linear deterministic SIRC epidemic model to model the influenza A 63 virus undergoing an antigenic drift. The SIRC model is a modified SIR model with an 64 additional compartment, C, for individuals that receive partial immunity from being 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 72 population, but only have vaccine for a certain strain that will not be effective for 73 others. The fact that viruses undergo changes regularly indicates that people who have 74 recovered from the virus, as well as individuals who have been vaccinated for a specific 75 strain of the virus, can be susceptible again to a newly-emerged strain. It is important 76 to determine the conditions in which a newly emerged strain and a common strain that 77 has a means of immunity will coexist in a population. 78 Another apt example for emerging diseases is the emergence of the COVID- 19 virus 79 in 2019 [39, 40] . As of March 17, 2020, there have been 167,545 confirmed cases of 80 COVID-19 in 150 different countries that has led to 6,606 deaths since it was declared 81 as an outbreak in January 2020 according to the WHO situation report [41] . At the 82 time that this paper is being written, there are papers that have modeled the dynamics 83 of the virus using different modifications of the SIR model. Zhou et. al. [42] included 84 compartments corresponding to suspected cases, which consists of the individuals that 85 show similar symptoms but are not confirmed cases, and indirectly infected individuals. 86 Pan et. al. [43] used a modified SEIR model which included asymptomatic and given that one of these strains has a vaccine available. This paper introduces a model 97 that approaches the lack of cross-immunity across different viral strains by introducing 98 new compartments to the SIR with vaccination model. This enables us to introduce 99 acquired immunity through vaccination and cross-immunity between strains in a simple 100 compartmental model and investigate the existence and stability of the resulting 101 equilibrium points. We aim to model the dynamics of an epidemic where a new emergent strain of an 103 existing virus affects a closed population. The existing virus will be modeled using a 104 modified SIR model with vaccination, however we assume that the vaccine does not 105 provide immunity to the newer strain. The equilibrium points were determined for the 106 system based on the transition equations and local stability was investigated for each 107 point. Once the stability conditions have been established, the epidemic model was 108 simulated using R [46] to check the steady-state behavior of the surveillance data for next two subsections will explain the dynamics before and after the emergence of the 124 mutated strain.

Before emergence 126 The system starts off as a population exposed to the original strain of the virus. The 127 spread of the virus is described by a modified SIR model that accounts for 128 vaccination [34] . The vaccinated members of the population can be treated as members 129

of an additional compartment that do not interact with the infected individuals. This 130 means that the modified SIR model will have four compartments instead of three, which 131 March 17, 2020 5/29

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted March 23, 2020. For this model, µ be the natural birth rate of the population, and consequently the 144 natural death rate of the population to keep the population size constant. It is assumed 145 that the individuals are vaccinated at birth with a vaccination rate p. β is the standard 146 incidence transmission coefficient, which assumes that the infection occurs based on how 147 many susceptible individuals interact with the infected [47] . For standard incidence, the 148 contact rate between infected and susceptible individuals is constant over all infected 149 individuals regardless of the population size [48] . The removal rate coefficient for the 150 infected individuals is denoted by γ.

The dynamics of the system is described by ordinary differential equations that 

March 17, 2020 6/29 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted March 23, 2020. 

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

Unlike the members of the susceptible compartment, the individuals in the . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted March 23, 2020. .

where ( 

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

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

. This means that 188 the endemic equilibrium point will be asymptotically stable if R v > 1, while the DFE 189 will be asymptotically stable if R v < 1. The reproduction number will be discussed in . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. to the existence of another compartment I 2 for those who are infected with the new 196 strain, which we will refer to as Disease 2.

The existence of the newer strain will be constrained by the following assumptions: 198

Since the vaccine is assumed to only work on the original strain, the vaccinated 199 and the previously removed individuals are susceptible to the newer strain. 

This means that the number of compartments that need to be monitored will increase 207 from four to six, with the addition or modification of the following compartments: . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. . As in Chauhan et. al's [34] work, the standard incidence was used to model infection 221 of the susceptible individuals for the newer strain. Based on Eq 1 and the compartment 222 diagram in Fig 2, the dynamics of the system can be expressed in terms of the following 223 ordinary differential equations:

and 

To solve for the equilibrium points of the system, Eqs 11-15 should be equal to zero. 228

Wolfram Mathematica [49] was used to obtain solutions for the system of equations, 229 which are the following:

3. New strain equilibrium:

is the reproduction number of the original 240 strain for the SIR model with vaccination [34] .

The third equilibrium corresponds to the scenario where only the new strain survives. 242

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

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 250 these points are stable. These conditions dictate which equilibrium will describe the 251 steady state behavior of the system. The local stability of the equilibrium points will be 252 March 17, 2020 11/29

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. . https://doi.org/10.1101/2020.03. 19.20039198 doi: medRxiv preprint determined based on the eigenvalues of its Jacobian evaluated at a specific equilibrium 253 point [17] . Let C 0 = (C 1 , C 2 , ...) T be the vector of the population number of each 254 compartment. For a general compartment C i , the components of the Jacobian, J ij can 255 be obtained using the following equation:

For our system, the Jacobian of the system can be obtained by applying Eq 25 to 257 Eqs 11-15. For any equilibrium point, (s,ī 1 ,v,r,ī 2 ) yields

Local stability is attained when the eigenvalues of the Jacobian, λ, are negative or 259 have negative real parts. In other words, the solutions for λ such that det(J − 1λ) = 0 260

should be negative or have negative real parts if the solution is complex [50] .

Simulations will then be used to check if the system approaches the equilibrium The Jacobian for the DFE can be obtained by substituting the respective values to 279 (s,ī 1 ,v,r,ī 2 ) in the expression for the Jacobian. This yields:

and the corresponding characteristic equation is

This means that the eigenvalues are 282

Recall that for the DFE to be locally 283 asymptotically stable, all eigenvalues should have negative real parts. Hence, the 284 following conditions should hold:

This is consistent with the local stability of the disease-free equilibrium for the 286 regular standard incidence SIR and the SIR with vaccination models [34] . The . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. . https://doi.org/10.1101/2020.03. 19.20039198 doi: medRxiv preprint reproduction number of the emergent disease is less than that of the original disease, 292 then the system will remain in DFE in the long run. Fig 3 shows Jacobian is then given by,

where s 1 , v 1 , r 1 are the respective equilibrium values for the susceptible, vaccinated, and 308

March 17, 2020 14/29 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. 

The discriminant of λ ± can dictate whether the eigenvalues will have a negative real 311

part. If the discriminant is negative or zero, then the eigenvalues will be negative. If the 312 discriminant is positive, recall that for the system to not go to the DFE, R 0 (1 − p) > 1. 313

This means that

which ensures that λ ± is negative when

would have a positive real part which makes the equilibrium point unstable. This 316 suggests that this equilibrium point will not be stable if the system was not already in 317 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 320 asymptotically stable. These conditions also imply that this equilibrium point can only 321 be achieved if the system before the emergence is already in endemic equilibrium with 322

Disease 1. Fig 4 shows the simulation of the system using parameters that satisfy Eq 31. 323

Unlike the DFE case, the proportion of individuals infected by the original strain, . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. . https://doi.org/10.1101/2020.03. 19.20039198 doi: medRxiv preprint pink dashed line, was shown to have a small spike, but quickly went to zero while the 329 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 331 one, meaning the emergent strain will be able to survive on its own in this population. 332

This implies that the new strain is not strong enough to infect enough people to 333 dominate over the original strain. 

The corresponding characteristic equation is given by,

March 17, 2020 16/29 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. . https://doi.org/10.1101/2020.03. 19.20039198 doi: medRxiv preprint Similar to the λ ± in the Disease 1 equilibrium, the real part will be negative if the 345 discriminant is negative. For the equilibrium point to exist, R 0 > 1, which means that 346 when the discriminant is positive, the following inequality holds

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

This indicates that the second disease will be locally stable if the mutated disease . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. . https://doi.org/10.1101/2020.03. 19.20039198 doi: medRxiv preprint of the second strain, the proportion of the population infected by the original strain, i 1 , 368 and the proportion of the individuals who recovered from the original strain, r 1 , 369 decrease and go to zero asymptotically. The behavior of the second strain is similar to 370 that in Fig 5. This implies that the new emergent strain is much more infectious than 371 the older strain that the new strain infects more susceptible individuals compared to the 372 original strain. This causes the original strain to die down at steady state. 373 Fig 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 R 0 = 5.71, while the reproduction number of the emergent strain is R 0 = 6.67. The vaccination rate used is 0.7.

For the endemic equilibrium case, both i 1 and i 2 are nonzero and thus both Eq 28 and 375 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 379 March 17, 2020 18/29 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. .

Recall that for the endemic equilibrium to be locally stable, all eigenvalues should 380 have a negative real part. Eq 39 shows that one of the eigenvalues is λ = −pµ/v * which 381 is negative. For all the roots of the quartic term to have negative real parts, the 382 Routh-Hurwitz criteria for stability should be applied [16, 17] . According to the 383 Routh-Hurwitz criterion, a polynomial with degree 4 will have roots (a 1 , a 2 , a 3 , a 4 ) that 384

all have negative real parts when:

Since the values of the equilibrium points should be positive, all coefficients 386 (a 1 , a 2 , a 3 , a 4 ) are positive. Based on the stability of the first three equilibrium points 387 and the existence criterion, the endemic equilibrium point is expected to be stable when 388

This implies that the endemic equilibrium for the two strains can only occur when 389 the system is initially in endemic equilibrium for the original strain, which can explain 390 why the condition is more restrictive than the one for the equilibrium with Disease 2. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. . the new strain given that Eq 44 is satisfied. Reproduction number [28] 406 One of the important parameters to be calculated for an epidemic model is the 407 reproduction number, which quantifies how infectious a certain disease is. Formally, the 408 reproduction number is defined as the expected number of secondary infections caused 409 by a single infected individual for the whole duration that they are infectious. A value 410 for the reproduction number that is greater than one indicates that the epidemic 411 persists in the population, while a value less than one means that the disease will die 412 out in the population [28, 29] .

The reproduction number R 0 of this epidemic model was calculated using the 414 approach formulated by van den Diessche and Watmough [28] . This approach does not 415 account for any measures taken to control the epidemic, but will give us an idea of the 416 conditions needed for the disease to spread on its own. from compartment i to the other non-infected compartments such as R 1 and R 2 [51] . 425 We define F and V such that for the disease free equilibrium X 0 ,

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

where α = γ + µ. Note that m = 2, so the corresponding F and V matrices are then 438 March 17, 2020 21/29 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. 

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

where F is the matrix that describes the infection rates for the two infections at the 444 DFE, and V −1 describes the average time an infected individual stays infectious. It is 445 easy to see that the eigenvalues of the next generation matrix are R 0 (1 − p) and R 0 .

This means that the reproduction number R 0 for this system is given by the larger of 447 the two. Formally,

Note that the resulting threshold equation for the system is R 0 less than one, which 449 means that the system will only approach the disease free equilibrium when Eq 53 is 450 less than one. For an outbreak to occur, at least one of these two strains should be able 451

to persist in the population on its own, that is, to have an individual reproduction . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 23, 2020. . https://doi.org/10.1101/2020.03. 19.20039198 doi: medRxiv preprint population. After establishing the possible transitions between compartments, the 459 system was found to have four equilibrium points: the disease free equilibrium, the 460 existing strain equilibrium, the emergent strain equilibrium, and the endemic 461 equilibrium. Upon examining the conditions for existence and local stability, the 462 disease-free equilibrium was determined to be locally stable when both R v and R 0 are 463 less than one. This is consistent with the reproduction number for this multi-strain 464 model. The existing strain equilibrium surprisingly did not impose that R 0 < 1, which 465 is the condition for the DFE of a normal SIR model for a single-strain without 466 immunity. This implies that if the original strain is much more contagious than the 467 emergent strain, the emergent strain would still die out eventually regardless of the fact 468 that it would persist in the population if it was on its own. On a similar note, the local 469 stability condition for the emergent equilibrium condition also does not impose that 470 R v < 1. This means that the state of the system prior to the emergence of the new 471 strain does not matter as long as the emergent strain is more contagious than the 472 existing strain. This leaves a highly restrictive condition for endemic equilibrium to 473 exist: the emergent strain must be able to survive by itself and the original strain must 474 be contagious enough to infect enough people, which is given by Eq 44. Eq 44 also 475 implies that the two-strain endemic equilibrium will only be locally asymptotically . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted March 23, 2020. . the existing influenza virus, which can highly influence the development of protocol in 491 receiving patients with flu-like symptoms and allocating resources in health care centers. 492

In summary, the modified multi-strain SIR model of an emerging disease that affects 493 both susceptible and previously immune individuals was studied. The local stability of 494 the equilibrium points as well as the reproduction number for the model were calculated. 495

Based on the results, we found that the original and the emergent strain can coexist in 496 an endemic equilibrium if the emergent strain has a lower reproduction number than 497 the original strain and that the system should already be in endemic equilibrium with 498 the original disease before the emergence. The requirement for the endemic equilibrium 499 to exist is quite strict especially for low values of R 0 and high values of p, which 500 presents a challenge in simulating surveillance data with stochastic incidence rates. 

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

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

Recent resurgence of mumps in the United States

Center for Disease Control and Prevention. Mumps Cases and Outbreaks

Center for Disease Control and Prevention. Measles Cases and Outbreaks

Measles in a Micronesian Community-King County

Measles Cases and Outbreaks

Center 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

Introduction to mathematical biology. Pearson/Prentice Hall

Compartmental analysis in biology and medicine

Containing papers of a mathematical and physical character

Viral infections leave a signature on human immune system

Evaluation of preexisting anti-hemagglutinin stalk antibody as a correlate of protection in a healthy volunteer challenge with influenza A/H1N1pdm virus

Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1)

Avian-human influenza epidemic model

A modeling study of human infections with avian influenza A H7N9 virus in mainland China

The SIRC model and influenza A

Emergence of drug resistance during an influenza epidemic: insights from a mathematical model

Methods of modelling viral disease dynamics across the within-and between-host scales: the impact of virus dose on host population immunity

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

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. Novel Coronavirus (2019-nCoV) situation reports

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

Center for Disease Control and Prevention. The Flu Season

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

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

Reproduction numbers of infectious disease models. Infectious Disease Modelling

Center of Disease Control and Prevention. Symptoms of Coronavirus Disease

We would like to acknowledge Dr. Yu Jin and Dr. Clay Cressler for their input on the 507 mathematical derivation and applications of the model.