Geometric Singular Perturbation Theory Analysis of an Epidemic Model with Spontaneous Human Behavioral Change
aa r X i v : . [ q - b i o . P E ] J un GEOMETRIC SINGULAR PERTURBATION THEORYANALYSIS OF AN EPIDEMIC MODEL WITHSPONTANEOUS HUMAN BEHAVIORAL CHANGE
STEPHEN SCHECTER
Abstract.
We consider a model due to Piero Poletti and col-laborators that adds spontaneous human behavioral change to thestandard SIR epidemic model. In its simplest form, the Polettimodel adds one differential equation, motivated by evolutionarygame theory, to the SIR model. The new equation describes theevolution of a variable x that represents the fraction of the popu-lation using normal behavior. The remaining fraction 1 − x usesaltered behavior such as staying home, social isolation, mask wear-ing, etc. Normal behavior offers a higher payoff when the numberof infectives is low; altered behavior offers a higher payoff when thenumber is high. We show that the entry-exit function of geometricsingular perturbation theory can be used to analyze the model inthe limit in which behavior changes on a much faster time scalethan that of the epidemic. In particular, behavior does not changeas soon as a different behavior has a higher payoff; current behav-ior is sticky. The delay until behavior changes in predicted by theentry-exit function. Introduction
A disease epidemic in a human population, such as measles, in-fluenza, or covid-19, spreads due to a combination of pathogen char-acteristics and human behavior. Pathogen characteristics determinethe circumstances under which an infected person can readily infectanother. Human behavior determines how frequently those circum-stances occur.Baseline human behavior varies with the society. In a city, crowdedconditions in housing, public transportation, schools and workplaces
Date : June 21, 2020.1991
Mathematics Subject Classification.
Key words and phrases. epidemiological modeling, entry-exit function, geometricsingular perturbation theory, imitation dynamics, evolutionary game theory.I thank Dan Marchesin and Marlon Michael Lopes Flores of IMPA for introduc-ing me to epidemiological modeling. Their work on epidemiology is supported byFAPERJ and Instituto Serrapilheira. may lead to frequent close human interactions; in a rural area this maybe less true. In East Asia mask-wearing in public is fairly common innormal circumstances; in other parts of the world it is rare.During an epidemic, human behavior may change due to govern-ment policies closing schools and businesses, requiring people to stayat home, and encouraging social distancing and mask-wearing.Spontaneous changes in human behavior also affect the course of anepidemic. People may react to an epidemic, or to information presentedto them, by spontaneously reducing social contacts, staying home to theextent possible, adopting more stringent hygiene or social distancing,or wearing a mask. People may adopt these behaviors independentof government policies; and, to the extent that they feel motivated toadopt such behaviors, they are more likely to comply with governmentorders and encouragement to do so.Similarly, when an epidemic wanes, or when people are presentedwith information that an epidemic is waning or that the disease is lessdangerous than originally feared, people may spontaneously return tonormal behavior. If restrictive government policies are still in place,compliance may decline.In the simplest epidemic models, SIR models, the transmissibility ofa disease in captured in a single parameter, β , defined as the numberof “adequate contacts” per unit time that an infected person has withother people [4]. If these other people are susceptible to the disease(not immune due to previous infection and not currently infected),an adequate contact results in a new infected individual. The basicreproduction number of the disease, R , is β times the typical lengthof time that an infected person remains infective. If R >
1, theninitially, when the susceptible fraction of the population is close to 1,the number of infected individuals will grow.Epidemic control measures aim to reduce β by enforcing or encour-aging changes in behavior. To determine what measures to institute,governments rely on epidemic models that estimate β under normalcircumstances and under various restrictive policies.A weakness of all epidemic models in current use, as far as I know, isthat they ignore spontaneous behavioral change. For example, the Im-perial College covid-19 model [5], which influenced the United Kingdomand United States government to institute social distancing measures[1], was based on a very detailed 2006 influenza epidemic model bythe same group [3]. According to the 2006 paper, “We do not assumeany spontaneous change in the behaviour of uninfected individuals as NALYSIS OF AN EPIDEMIC MODEL 3 the pandemic progresses, but note that behavioural changes that in-creased social distance together with some school and workplace clo-sure occurred in past pandemics . . . and might be likely to occur in afuture pandemic even if not part of official policy. . . . Such spontaneouschanges in population behaviour might more easily reduce peak dailycase incidence.”Epidemiologists appear to be well aware that spontaneous behavioralchange should be incorporated in models. There is a fairly extensiveliterature on ways to do it; a review article is [16]. There does notappear to be agreement on what modeling approach is best. Thisprobably should not be regarded as a serious problem; a variety ofdifferent models are commonly used in epidemiology. A more seriousissue is that there has been little work on how to determine the values ofthe parameters in the models [16]. Without approximate values for theparameters, models of spontaneous behavioral change can only yieldqualitative predictions.The goal of this paper is not to deal with the various issues of howbest to account for behavioral change in epidemic models. Instead wewant to direct attention to a particular approach to behavioral change,due to Piero Poletti and collaborators, in its simplest form [14, 13, 15].This model adds one equation, motivated by evolutionary game theory[12, 6], to the standard SIR model. Our goal is to show how the entry-exit function [2] of geometric singular perturbation theory [8, 9] canbe used to analyze this model. Given values for the parameters, theentry-exit function enables one to make precise predictions in the limitwhere behavioral change occurs on a much faster time scale than theepidemic itself.To my knowledge, there have been two previous uses of the entry-exitfunction in epidemiological models, [10] and [7].Figure 1.1 shows a typical simulation of the Poletti model. Thereare three variables. Two, S and I , are the familiar susceptible andinfective population fractions from the SIR model. The third variable, x , represents the fraction of the population using normal behavior.When x = 1, in this simulation, the model reduces to an SIR modelwith R = 3. When x = 0, the entire population uses altered behavior;in this simulation, the model reduces to an SIR model with R = .
6. Inthe simulation, behavior changes on a time scale 200 times faster thanthat of the epidemic itself. Thus, if the time scale for the epidemic isdays, 1000 time units in the simulation equals five days. The simulationshows 20,000 fast time units, or 100 days.Altered behavior yields a negative payoff due to loss of income, lossof social interactions, and so on. However, altered behavior reduces the
SCHECTER time f r a c t i on o f popu l a t i on susceptibleinfectednormal behavior Figure 1.1.
A simulation of the Poletti model. At thestart (
S, I, x ) = ( . , . , . I <
1, almost allthe population quickly adopts normal behavior. After I rises to about .18 (showing behavior stickiness), thepopulation switches to altered behavior. I falls to about.05 (again showing behavior stickiness); the populationreturns to normal behavior; and I rises to about .13 (sec-ond wave). After two more behavioral switches, the epi-demic dies out.chance of getting the disease. In this simulation, normal behavior yieldsa higher payoff to the individual when I < .
1. When
I > .
1, alteredbehavior yields a higher payoff. There is therefore a tendency to adoptaltered behavior, which moderates the epidemic, when I passes . I falls below .
1, there is a tendency to resume normal behavior.Resuming normal behavior can result in a “second wave” of infections,as seen in the simulation.In the Poletti model, behavior changes because of encounters withother people whose behavior offers a higher payoff than one’s own.Thus, in the simulation, behavior does not change immediately when I passes .1; the current behavior is “sticky.” The delay until behaviorchanges can be calculated in the limit from the entry-exit function.The rather fast evolution of the epidemic in the simulation is due tothe choices R = 3 and R = . NALYSIS OF AN EPIDEMIC MODEL 5
The Poletti model, in my view, plays a role similar to the SIR model:it gives the essence of the situation, stripped of complications, andcan form the basis for more realistic models. I expect that geometricsingular perturbation theory will also prove useful in analyzing morerealistic extensions of the model.In the next few sections of the paper we review the SIR model (sec-tion 2), introduce the Poletti model (section 3), and describe and ex-ploit the model’s slow-fast structure (section 4). The main result of thepaper, Theorem 1, is stated at the end of section 4. We then provideexamples (section 5) and proofs (section 6), and conclude with a briefdiscussion (section 7). 2.
The SIR model
The Poletti model is based on the standard SIR model for an epi-demic, ˙ S = − βSI, (2.1)˙ I = βSI − γI, (2.2)˙ R = γI, (2.3)with ˙ = ddt . The variables S , I , and R are population fractions; theysum to 1. (Since ˙ S + ˙ I + ˙ R = 0, the sum S + I + R is constant.) S is thefraction of the population that is susceptible to acquiring the disease; I is the fraction that is currently infected; R is the fraction that hasrecovered. ( R is sometimes called the fraction removed. If there aredeaths due to the disease, they are included in R without change to themodel.) Thus the equation for R can be ignored; R can be recoveredfrom R = 1 − S − I . The system reduces to˙ S = − βSI, (2.4)˙ I = βSI − γI (2.5)on the triangle T = { ( S, I ) : S ≥ , I ≥ , S + I ≤ } .Let ˆ T = { ( S, I ) ∈ T : S >
I >
0. In ˆ T the orbits of (2.4)–(2.5)satisfy the differential equation dIdS = − γβS , so they are curves I + S − γβ ln S = C. (2.6)The parameter β was discussed in the introduction. The averagelength of time an individual is infected is γ . Thus the basic reproduc-tion number of the disease R mentioned in the introduction is βγ . SCHECTER
Phase portraits on T in the cases 0 < R < R > I = 0, 0 ≤ S ≤
1. Each solution approaches one of the equilibria( S f , R = 1 − S f is the fraction of the population that contracted thedisease in the course of the epidemic.In ˆ T , ˙ S <
0, so the number of susceptibles steadily falls. If 0 < R <
1, ˙
I < T as well, so the number of infectives also steadily falls. If R >
1, ˙
I <
I >
0) for 0 < S < γβ (resp. γβ < S < γβ < S <
1, then I increases until S has fallen to βγ ; after that I decreases. NALYSIS OF AN EPIDEMIC MODEL 7
Susceptibles I n f e c t i v e s (a) Susceptibles I n f e c t i v e s (b) Figure 2.1.
Phase portraits of SIR models. (a) β =1 / γ = 1 /
6, so R = 6 /
10. (b) β = 1 / γ = 1 /
6, so R = 3. In case (b), the vertical line S = γ/β = 1 / I is alsoshown. All solutions move to the left as time increases.3. The Poletti Model
In the Poletti model, susceptible individuals have two available be-haviors, normal, for which β = β n with basic reproduction number R = β n γ >
1, and altered, for which β = β a with basic reproductionnumber R = β a γ <
1. Altered behavior may include staying home tothe extent possible, practicing social distancing, mask wearing, etc.
SCHECTER
Each behavior has a payoff to a susceptible who adopts it. Thepayoffs are p n = − m n I and p a = − k − m a I with m n , m a and k positive and m n > m a . The negative payoff − m n I is due to the possibility that a susceptible with normal behavior willcontract the disease; it is proportional to I , the fraction of infectivesin the population. The negative payoff − m a I is due to the possibilitythat a susceptible with altered behavior will contract the disease; it isalso proportional to I , but the proportionality constant is less negative.In addition, altered behavior has a negative payoff − k independent of I that represents loss of income, loss of valued social interactions, etc.The payoff from altered behavior is higher if and only if I > km n − m a , i.e.,if and only if the fraction of infectives in the population is sufficientlyhigh. We assume km n − m a < . This assumption allows altered behavior to sometimes have a higherpayoff.Susceptibles are assumed to change their behavior from normal toaltered, or vice-versa, due to imitation of other susceptibles they en-counter who are using the opposite behavior and experiencing a higherpayoff. The mathematical formulation of this notion is called imitationdynamics and comes from evolutionary game theory [6].Let x denote the fraction of the susceptibles using normal behavior,so that 1 − x is the fraction using altered behavior. We continue to let S and I denote the susceptible and infected fractions of the population.Then the complete model is˙ S = − (cid:0) β n x + β a (1 − x ) (cid:1) SI, (3.1)˙ I = (cid:0) β n x + β a (1 − x ) (cid:1) SI − γI, (3.2)˙ x = x (1 − x )( β a − β n ) I + 1 ǫ x (1 − x ) (cid:0) k − ( m n − m a ) I (cid:1) , (3.3)with ˙ = ddt . The state space is the prism P = { ( S, I, x ) : S ≥ , I ≥ , S + I ≤ , ≤ x ≤ } . There is also an equation for the recovered fraction of the population R , ˙ R = γI ; we ignore it since R can be recovered from S = 1 − R − I .For the derivation of the model, see [14, 13]. It can be intuitivelyunderstood as follows. NALYSIS OF AN EPIDEMIC MODEL 9
The equations for ˙ S and ˙ I come from assuming that both susceptibleswith normal behavior and susceptibles with altered behavior satisfy SIRmodels.The first summand in the equation for ˙ x is negative; it expresses thefact that susceptibles with normal behavior acquire the disease moreeasily than susceptibles with altered behavior, and hence more readilyleave the susceptible group. Thus the fraction of susceptibles usingnormal behavior tends to decrease.The second summand represents the the rate of change of x due toimitiation dynamics. The rate at which susceptibles using different be-haviors encounter each other is proportional to x (1 − x ). The differencein payoffs of the two behaviors, given the current level of I , is p n − p a = k − ( m n − m a ) I. When this number is positive, normal behavior yields a larger payoff, so x increases at a rate proportional to the difference between the payoffs;when this number is negative, x decreases in the same manner.The rate constant multiplying this summand is written as ǫ with ǫ >
0. We will assume that this constant is large, so that ǫ is small.In other words, we assume that behavior can change on a much fastertime scale than that of the epidemic itself.4. Slow-fast structure
Slow-fast structure.
System (3.1)–(3.3), in which we recall that˙ = ddt , is a slow-fast system [8, 9] with two slow variables, S and I ,and one fast variable, x ; t is the slow time. Such systems are morecommonly written with the last equation multiplied on both sides by ǫ : ˙ S = − (cid:0) β n x + β a (1 − x ) (cid:1) SI, (4.1)˙ I = (cid:0) β n x + β a (1 − x ) (cid:1) SI − γI, (4.2) ǫ ˙ x = ǫx (1 − x )( β a − β n ) I + x (1 − x ) (cid:0) k − ( m n − m a ) I (cid:1) . (4.3)The fast time τ satisfies t = ǫτ . With ′ = ddτ , system (4.1)–(4.3)becomes S ′ = − ǫ (cid:0) β n x + β a (1 − x ) (cid:1) SI, (4.4) I ′ = ǫ (cid:0) β n x + β a (1 − x ) (cid:1) SI − ǫγI, (4.5) x ′ = ǫx (1 − x )( β a − β n ) I + x (1 − x ) (cid:0) k − ( m n − m a ) I (cid:1) . (4.6)The slow system (4.1)–(4.3) and the fast system (4.4)–(4.6) have thesame phase portraits for ǫ >
0, but they have different limits at ǫ = 0. For ǫ = 0, the slow system (4.1)–(4.3) becomes the slow limit system˙ S = − (cid:0) β n x + β a (1 − x ) (cid:1) SI, (4.7)˙ I = (cid:0) β n x + β a (1 − x ) (cid:1) SI − γI, (4.8)0 = x (1 − x ) (cid:0) k − ( m n − m a ) I (cid:1) , (4.9)and the fast system (4.4)–(4.6) becomes the fast limit system S ′ = 0 , (4.10) I ′ = 0 , (4.11) x ′ = x (1 − x ) (cid:0) k − ( m n − m a ) I (cid:1) . (4.12)Singular solutions are constructed by combining solutions of the slowlimit system (4.7)–(4.9) and the fast limit system (4.10)–(4.12). Inmany situations, solutions for small ǫ > Fast limit system.
For the fast limit system (4.10)–(4.12), eachline segment (
S, I ) = ( S , I ) is invariant, and the triangles x = 0 and x = 1 consist of equilibria. (The plane I = km n − m a also consists ofequilibria, but we will not make direct use of them.) On line segments( S, I ) = ( S , I ) with 0 ≤ I < km n − m a , ˙ x >
0, so the solution x ( t )of (4.12) satisfies lim t →−∞ x ( t ) = 0 and lim t →∞ x ( t ) = 1. On linesegments ( S, I ) = ( S , I ) with km n − m a < I ≤
1, ˙ x <
0, so the solution x ( t ) of (4.12) satisfies lim t →−∞ x ( t ) = 1 and lim t →∞ x ( t ) = 0. The fastdynamics just reflect the fact that normal behavior gives a higher payoffif I < km n − m a , and altered behavior gives a higher payoff if I > km n − m a .Equilibria of the fast limit system (4.10)–(4.12) are normally attract-ing if ∂ ˙ x∂x < ∂ ˙ x∂x >
0. One can check thatequilibria with x = 0 are normally repelling for I < km n − m a and nor-mally attracting for I > km n − m a . Equilibria with x = 1 are the reverse.4.3. Slow limit system.
The slow limit system (4.7)–(4.9) makessense on the triangles x = 0 and x = 1.On the triangle x = 0, the slow limit system reduces to˙ S = − β a SI, (4.13)˙ I = β a SI − γI. (4.14)This is just an SIR model with β = β a and basic transmission number R < NALYSIS OF AN EPIDEMIC MODEL 11
Similarly, on the triangle x = 1, the slow limit system reduces to˙ S = − β n SI, (4.15)˙ I = β n SI − γI. (4.16)This is just an SIR model with β = β n and basic transmission number R > ǫ >
0, the triangles x = 0 and x = 1 remain invariant. Theslow system (4.1)–(4.3), restricted to x = 0, remains (4.13)– (4.14).Restricted to x = 1 it remains (4.15)– (4.16). Thus the line segments { ( S, I, x ) : 0 ≤ S ≤ , I = 0 , x = 0 } and { ( S, I, x ) : 0 ≤ S ≤ , I =0 , x = 1 } remain equilibria.We will use the following notation where convenient: • φ ǫ (cid:0) ( S , I , x ) , t (cid:1) = solution of (4.4)–(4.6) with φ ǫ (cid:0) ( S , I , x ) , (cid:1) =( S , I , x ). • ψ (cid:0) ( S , I ) , t (cid:1) = solution of (4.13)–(4.14) with ψ (cid:0) ( S , I ) , (cid:1) =( S , I ). • ψ (cid:0) ( S , I ) , t (cid:1) = solution of (4.15)–(4.16) with ψ (cid:0) ( S , I ) , (cid:1) =( S , I ).4.4. Entry-exit function for the triangle x = 0 . In the triangle x = 0, let ( S , I ) ∈ ˆ T with I > km n − m a , so ( S , I ) lies in the attractingportion of the triangle. Let ( S ( t ) , I ( t )) = ψ (cid:0) ( S , I ) , t (cid:1) , let t >
0, andlet ( S , I ) = (cid:0) S ( t ) , I ( t ) (cid:1) . The solution (cid:0) S ( t ) , I ( t ) (cid:1) traces out a curveΓ, which from (2.6) has the equation I + S − γβ a ln S = v , v = I + S − γβ a ln S = I + S − γβ a ln S . (4.17)We define the entry-exit integral I (cid:0) ( S , I ) , ( S , I ) (cid:1) = Z t k − ( m n − m a ) I ( t ) dt = Z S S − k − ( m n − m a )( v − S + γβ a ln S ) β a S ( v − S + γβ a ln S ) dS. (4.18)The second integral follows from the first by making the substitutions S = S ( t ) , dS = − β a S ( t ) I ( t ) dt , and I = v − S + γβ a ln S , whichfollows from (4.19). It cannot be evaluated analytically, but is readilyevaluated numerically.The integrand of the first integral is negative when I > km n − m a andpositive when I < km n − m a . The integral represents accumulated attrac-tion to (resp. repulsion from) the plane x = 0 where the integrand isnegative (resp. positive). Proposition 1.
For each point ( S , I ) in ˆ T with I > km n − m a , thereis exactly one t > such that ( S , I ) = ( S ( t ) , I ( t )) satisfies I (cid:0) ( S , I ) , ( S , I ) (cid:1) = 0 . Of course, ( S , I ) lies in the region I < km n − m a . Intuitively, at( S , I ) the accumulated repulsion from the plane x = 0 balances theaccumulated attraction to the plane. We shall see that for small ǫ > x = 0near ( S , I ) will track a solution of (4.13)–(4.14) near ( S ( t ) , I ( t )) untilit leaves the neighborhood near ( S , I ). See Figure 4.1 and Subsection6.2. Figure 4.1.
A solution of (4.4)–(4.6) approachesthe triangle x = 0 near a point ( S , I ), followsthe solution of (4.13)–(4.14) through ( S , I ) until I (cid:0) ( S , I ) , ( S , I ) (cid:1) = 0, then leaves the triangle.4.5. Entry-exit function for the triangle x = 1 . In the triangle x = 1, let ( S , I ) ∈ ˆ T with I < km n − m a , so ( S , I ) lies in the attractingportion of the triangle. Let ( S ( t ) , I ( t )) = ψ (cid:0) ( S , I ) , t (cid:1) , let t >
0, andlet ( S , I ) = (cid:0) S ( t ) , I ( t ) (cid:1) . The solution (cid:0) S ( t ) , I ( t ) (cid:1) traces out a curveΓ, which from (2.6) has the equation I + S − γβ n ln S = v , v = I + S − γβ n ln S = I + S − γβ n ln S . (4.19) NALYSIS OF AN EPIDEMIC MODEL 13
We define the entry-exit integral I (cid:0) ( S , I ) , ( S , I ) (cid:1) = Z t k − ( m n − m a ) I ( t ) dt = Z S S − k − ( m n − m a )( v − S + γβ n ln S ) β n S ( v − S + γβ n ln S ) dS. (4.20)The second integral follows from the first as in the previous subsection.The integrand of the first integral is negative when I < km n − m a andpositive when I > km n − m a . The integral represents accumulated attrac-tion to (resp. repulsion from) the plane x = 1 where the integrand isnegative (resp. positive).The system (4.15)–(4.16) on T has a unique orbit that is tangent tothe line I = km n − m a . The point of tangency is ( S ∗ , I ∗ ), I ∗ = km n − m a .Let Γ ∗ denote the part of this orbit with S ∗ < S <
1. Let Γ ∗ have theequation S = S ∗ ( I ), 0 < I < km n − m a . Let V − = { ( S, I ) ∈ ˆ T : 0 < S
Figure 4.2.
Phase portrait of the fast limit system(4.15)–(4.16) in the triangle x = 1 with β n = 1 / γ = 1 /
6. The vertical line S = γ/β = 1 / I = km n − m a = are shown, as are the sets V − and V + bounded above by this line. Solutions thatstart in V − approach equilibria without crossing the line I = km n − m a ; solutions that start in V + cross the line.Let ( S , I ) ∈ V − and let ( S ( t ) , I ( t )) = ψ (cid:0) ( S , I ) , t (cid:1) . Then (cid:0) S ( t ) , I ( t ) (cid:1) ∈ V − for all t ≥
0. Thus I (cid:0) ( S , I ) , ( S , I ) (cid:1) is never 0. As t → ∞ , (cid:0) S ( t ) , I ( t ) (cid:1) approaches an equilibrium ( S f ,
0) of (4.15)–(4.16).In this case, for small ǫ >
0, a solution of (4.4)–(4.6) that enters aneighborhood of the plane x = 1 near ( S , I ) will track a solution of(4.15)–(4.16) near ( S ( t ) , I ( t )) and approach an equilibrium ( S ǫf , ,
1) of(4.4)–(4.6) with S ǫf near S f . See Subsection 6.4.Let ( S , I ) ∈ V + and let ( S ( t ) , I ( t )) = ψ (cid:0) ( S , I ) , t (cid:1) . Then (cid:0) S ( t ) , I ( t ) (cid:1) enters the region I ≥ km n − m a at t = t in > t = t out ≥ t = t in .If R t out k − ( m n − m a ) I ( t ) dt <
0, then there is no point ( S , I )where I (cid:0) ( S , I ) , ( S , I ) (cid:1) . As in the previous paragraph, let ( S f ,
0) =lim t →∞ (cid:0) S ( t ) , I ( t ) (cid:1) . In this case also, for small ǫ >
0, a solution of(4.4)–(4.6) that enters a neighborhood of the plane x = 1 near ( S , I )will track a solution of (4.15)–(4.16) near ( S ( t ) , I ( t )) and approach anequilibrium ( S ǫf , ,
1) of (4.4)–(4.6) with S ǫf near S f .If R t out k − ( m n − m a ) I ( t ) dt >
0, then there is a unique point ( S , I ),with I > km n − m a , where I (cid:0) ( S , I ) , ( S , I ) (cid:1) = 0. For small ǫ >
0, asolution of (4.4)–(4.6) that enters a neighborhood of the plane x = 1near ( S , I ) will track a solution of (4.15)–(4.16) near ( S ( t ) , I ( t )) untilit leaves the neighborhood near ( S , I ). See Subsection 6.3.4.6. Singular orbits.
Motivated by the previous subsections, we con-struct singular orbits of the system (4.4)–(4.6) (or equivalently (4.1)–(4.3)).Consider a starting point ( S , I , x ) with ( S , I ) ∈ ˆ T , I < km n − m a ,and 0 < x <
1. At this point, ˙ x > S is a fast orbit: the portion of the line ( S, I ) =( S , I ) with x ≤ x < S , there are three cases. Let( S ( t ) , I ( t )) = ψ (cid:0) ( S , I ) , t (cid:1) , and, given t >
0, let ( S , I ) = (cid:0) S ( t ) , I ( t ) (cid:1) .2a. Suppose there exists t > I (cid:0) ( S , I ) , ( S , I ) (cid:1) = 0 and I > km n − m a . The next orbit of S is { ( S ( t ) , I ( t )) : 0 ≤ t ≤ t } , a sloworbit.2b. Suppose there exists t > I (cid:0) ( S , I ) , ( S , I ) (cid:1) = 0,and I = km n − m a . In this case the construction of the singular orbitfails. (Notice t = t out from the previous subsection.)2c. Suppose there is no t > I (cid:0) ( S , I ) , ( S , I ) (cid:1) = 0.Then the next orbit of S is { ( S ( t ) , I ( t )) : t ≥ } , a slow orbit. Thisorbit approaches an equilibrium of (4.15)–(4.16). The construction of S terminates. NALYSIS OF AN EPIDEMIC MODEL 15
3. We continue the construction of S in case 2a. The next orbit in S is a fast orbit: the portion of the line ( S, I ) = ( S , I ) with 0 < x < S ( t ) , I ( t )) = ψ (cid:0) ( S , I ) , t (cid:1) . By Proposition 1 there is exactlyone t > S , I ) = ( S ( t ) , I ( t )) satisfies I (cid:0) ( S , I ) , ( S , I ) (cid:1) =0. We have I < km n − m a . The next orbit of S is { ( S ( t ) , I ( t )) : 0 ≤ t ≤ t } , a slow orbit.5. The next orbit in S is a fast orbit: the portion of the line ( S, I ) =( S , I ) with 0 < x < S , I ) equal to( S , I ).Next we consider a starting point ( S , I , x ) with ( S , I ) ∈ ˆ T , I > km n − m a , and 0 < x <
1. In this case the first orbit in S is again a fastorbit: the portion of the line ( S, I ) = ( S , I ) with 0 < x ≤ x . Wecontinue the construction of S at step 4 above, setting ( S , I ) equal to( S , I ).In both cases the singular orbit S is an alternating sequence of fastand slow orbits, with the first orbit fast. The slow orbits alternatebetween orbits in x = 0 and orbits in x = 1. The last orbit is a sloworbit in x = 1 that approaches an equilibrium, for which I = 0. Theorem 1.
Let ( S , I , x ) satisfy ( S , I ) ∈ ˆ T , I = km n − m a , and < x < . Suppose the construction of the singular orbit S thatstarts at ( S , I , x ) never fails at step 2 and terminates after a finitenumber of steps at ( S f , , . Let Γ ǫ denote the orbit of (4.4) – (4.6) that starts at ( S , I , x ) . Then as ǫ → , Γ ǫ → S . The terminal point ( S ǫf , , of Γ ǫ converges to ( S f , , . Roughly speaking, the fast jumps between x = 0 and x = 1 occurbecause the predominant behavior among the susceptibles has becomeless rewarding than the alternative. When normal behavior predom-inates ( x near 1), if the fraction of infectives becomes high, behaviormay switch to the altered form ( x near 0). When altered behaviorpredominates ( x near 0), the fraction of infectives becomes low, andbehavior swiches to the normal form ( x near 1).However, the switch does not occur immediately when the number ofinfectives crosses the threshhold value I = km n − m a . As was mentionedin the introduction, behavior changes because of encounters with otherpeople whose behavior offers a higher payoff than one’s own, so thecurrent behavior is “sticky.” The delay until behavior changes can becalculated in the limit from the entry-exit function. Examples
We consider the system (4.1)–(4.3) with the parameter values β n = 1 / , β a = 1 / , γ = 1 / , k = 3 / , m n = 5 , m a = 2 . From the values of β n , β a , and γ , we see R = 3 for normal behaviorand .6 for altered behavior. The phase portrait of (4.13)–(4.14) in thetriangle x = 0 (resp. (4.15)–(4.16) in the triangle x = 1) is given byFigure 1(a) (resp. Figure 1(b)). The plane I = km n − m a is I = 1 / P start = ( S , I , . I < /
10. Such a singular orbit starts with a fast solution from P start to ( S , I , S , I ,
1) to a point P end = ( S f , , P start , P end ). Otherwise we represent the singular orbit by asequence ( P start , P , P , . . . , P k , P end ) , where • the first fast orbit goes from P start = ( S , I , .
98) to ( S , I , • the first slow orbit goes from ( S , I ,
1) to P = ( S , I , • the second fast orbit goes from P = ( S , I ,
1) to ( S , I , • the second slow orbit goes from ( S , I ,
0) to P = ( S , I , • the third fast orbit goes from P = ( S , I ,
0) to ( S , I , • the last fast orbit goes from P k = ( S k , I k ,
0) to ( S k , I k , • the last slow orbit goes from ( S k , I k ,
1) to P end = ( S f , , P , . . . , P k are the starting points of fast jumps; P i with i odd is in x = 1, and P i with i even is in x = 0.Using the Matlab routines in the appendix, one can compute singularorbits for this system. We give three examples. Corresponding to eachexample we show the solution of the fast system (4.4)–(4.6) with thesame starting point and ǫ = . ≤ t ≤ , ǫ = . I at x = 1 / S at t = 30 , NALYSIS OF AN EPIDEMIC MODEL 17
Example 1.
A singular solution with two jumps. P start = ( . , . , . P = ( . , . , P = ( . , . , P end = ( . , , x occur successively at I = . I = . S = . t = 30 , Example 2.
A singular solution with four jumps. This example wasshown in the introduction. P start = ( . , . , . P = ( . , . , P = ( . , . , P = ( . , . , P = ( . , . , P end = ( . , , x occur successively at I = . I = . I = . I = . S = . t = 30 , Example 3.
A singular solution with six jumps. P start = ( . , . , . P = ( . , . , P = ( . , . , P = ( . , . , P = ( . , . , P = ( . , . , P = ( . , . , P end = ( . , , See Figure 5.3. For the computed solution, jumps in x occur succes-sively at I = . I = . I = . I = . I = . I = . S = . t = 30 , NALYSIS OF AN EPIDEMIC MODEL 19
Susceptibles I n f e c t i v e s PstartP1P2Pend (a) time f r a c t i on o f popu l a t i on susceptibleinfectednormal behavior (b) Figure 5.1.
Example 1, P start = ( . , . , . x = 1, with the vertical line S = β n γ = and the hori-zontal line I = km n − m a = shown. The slow orbits from P start to P and from P to P end are shown in this phaseportrait. The slow orbit from P to P lies in the triangle x = 0; see Figure 1(a). (b) Solution of the fast system(4.4)–(4.6) with the same starting point and ǫ = . Susceptibles I n f e c t i v e s PstartP1P2P3P4Pend (a) time f r a c t i on o f popu l a t i on susceptibleinfectednormal behavior (b) Figure 5.2.
Example 2, P start = ( . , . , . x = 1, with the vertical line S = β n γ = and the hori-zontal line I = km n − m a = shown. The slow orbits from P start to P , from P to P , and from P to P end are shownin this phase portrait. The slow orbits from P to P andfrom from P to P lie in the triangle x = 0; see Figure1(a). (b) Solution of the fast system (4.4)–(4.6) with thesame starting point and ǫ = . NALYSIS OF AN EPIDEMIC MODEL 21
Susceptibles I n f e c t i v e s PstartP1P2P3P4P5P6Pend (a) time f r a c t i on o f popu l a t i on susceptibleinfectednormal behavior (b) Figure 5.3.
Example 3, P start = ( . , . , . x = 1, with the vertical line S = β n γ = and the hori-zontal line I = km n − m a = shown. The slow orbits from P start to P , from P to P , from P to P , and from P to P end are shown in this phase portrait. The slow orbitsfrom P to P , from P to P , and from P to P lie inthe triangle x = 0; see Figure 1(a). (b) Solution of thefast system (4.4)–(4.6) with the same starting point and ǫ = . Proofs
Entry-exit function.
Let U be an open subset of R n , n ≥ c ′ = ǫp ( c, z, ǫ ) , (6.1) z ′ = zq ( c, z, ǫ ) , (6.2)with ( c, z, ǫ ) ∈ U × [0 , z ) × [0 , ǫ ) and ′ = ddτ . We assume(E1) p and q are of class C r , r ≥ q ( c, ,
0) = 0, then Dq ( c, , p ( c, , > q ( c, ,
0) = 0 defines a C r codimension-one submanifold S of U .Let t = ǫτ and let ˙ = ddt . For c ∈ U with q ( c , , <
0, let φ ( c , t )denote the solution of ˙ c = p ( c, ,
0) with φ ( c ,
0) = c . Assumption(E2) implies that φ ( c , t ) crosses the manifold S at most once.Given t >
0, let c = ψ ( c , t ). Define the entry-exit integral I ( c , c ) = Z t q ( ψ ( c , t ) , , dt. (6.3) Theorem 2.
For system (6.1) – (6.2) satisfying (E1)–(E2), assume I (¯ c , ¯ c ) = 0 . For a small neighborhood ˜ U of ¯ c in U , define the entry-exit function P : ˜ U → U by P ( c ) = c , where c is defined implicitlyby I ( c , c ) = 0 . Fix δ > sufficiently small. For a given ǫ > ,consider the solution of (6.1) – (6.2) that starts at ( c, z ) = ( c , δ ) , with c ∈ ˜ U . Then: (1) For ǫ > sufficiently small, the solution reintersects the section z = δ at a point ( c, z ) = ( P ǫ ( c ) , δ ) . (2) P ǫ and P are C r functions, and P ǫ → P in the C r sense as ǫ → . (3) Let Γ ǫ denote the orbit of (6.1) – (6.2) from ( c , δ ) to ( P ǫ ( c ) , δ ) .As ǫ → , Γ ǫ approaches the singular orbit of (6.1) – (6.2) con-sisting of (a) the line segment [( c , δ ) , ( c , ; (b) (Γ , , where Γ is the orbit of ˙ c = p ( c, , from c to c = P ( c ) ; (c) the line segment (cid:0) ( c , , ( c , δ ] (cid:1) . Theorem 2 is proved in [2] under the assumption that for ǫ = 0, thesystem (6.1)–(6.2) has been written in a standard form. The relationof the standard form to the form (6.1)–(6.2) is explained in [11]. NALYSIS OF AN EPIDEMIC MODEL 23
Entry-exit function for the Poletti model in the triangle x = 0 . The Poletti model (4.4)–(4.6) satisfies the hypotheses of The-orem 2 for n = 1 and any r ≥
2, with (
S, I ) corresponding to p and x corresponding to z . The set U is ˆ T , and S is the line I = km n − m a . ( ˆ T is not open, since it includes a segment of the line S + T = 1, but thisdoes not cause any difficulty.)Let ( S , I ) ∈ ˆ T with I > km n − m a . Let ( S ( t ) , I ( t )) = ψ (cid:0) ( S , I ) , t (cid:1) ,let t >
0, and let ( S , I ) = ( S ( t ) , I ( t )). The formula (4.18) for theentry-exit integral I follows immediately from (6.3).Proof of Proposition 1: We consider the solution ( S ( t ) , I ( t )) of (4.13)–(4.14) defined above. Since I ( t ) is decreasing, there is a unique t ∗ > I ( t ∗ ) = km n − m a . The integral R t k − ( m n − m a ) I ( t ) dt is nega-tive and decreasing for 0 < t ≤ t ∗ and is increasing for t > t ∗ . To proveProposition 1 it suffices to show that R ∞ k − ( m n − m a ) I ( t ) dt = ∞ .Since R ∞ k dt = ∞ , it suffices to show that R ∞ ( m n − m a ) I ( t ) dt isfinite. To see this, just note that ( S ( t ) , I ( t )) approaches a normallyattracting equilibrium ( S f , I ( t ) → Entry-exit function for the Poletti model in the triangle x = 1 . To treat the Poletti model (4.4)–(4.6) near x = 1, we first makethe change of variables y = 1 − x . We obtain the system S ′ = − ǫ ( β n (1 − y ) + β a y ) SI, (6.4) I ′ = ǫ ( β n (1 − y ) + β a y ) SI − ǫγI, (6.5) y ′ = ǫ (1 − y ) y ( β n − β a ) I − (1 − y ) y ( k − ( m n − m a ) I ) . (6.6)Define the curve C to be the union of the line I = km n − m a , 0 < S ≤ S ∗ , and Γ ∗ defined in Subsection 4.5. Let U be the part of T above C .The system (4.4)–(4.6) satisfies the hypotheses of Theorem 2 for n =1 and any r ≥
2, with (
S, I ) corresponding to p and y correspondingto z . The set U is defined above, and S is the line segment I = km n − m a , S ∗ < S <
1. (Again the set U is not open because it includes a segmentof the line S + T = 1, but this does not cause any difficulty.)As in the previous subsection, the formula for the entry-exit integral I (cid:0) ( S , I ) , ( S , I ) (cid:1) follows immediately from (6.3).6.4. Solutions that approach equilibria.
Recall the sets V − and V + defined in Subsection 4.5. Proposition 2.
Let K be a compact subset of V − . For each ( S , I ) ∈ K , let ( S f ,
0) = lim t →∞ ψ (cid:0) ( S , I ) , t (cid:1) . Define Q : K → R by Q ( S , I ) = S f . Let δ > be small. Then: (1) For small ǫ > and for each ( S , I ) ∈ K , there exists S ǫf ∈ R such that lim t →∞ φ ǫ (cid:0) ( S , I , δ ) , t (cid:1) = ( S ǫf , , . (2) Define Q ǫ : K → R by Q ǫ ( S , I ) = S ǫf . Then Q ǫ and Q are C r − functions, and Q ǫ → Q in the C r − sense as ǫ → .Proof. Let ˆ K be a compact subset of V − that contains K in its interior.Let ˜ K denote the union of ˆ K , solutions of (4.15)–(4.16) that start in K , and the limits of these solutions.For the system (4.4)–(4.6) with ǫ = 0, ˜ K is a union of equilbria thatis compact, normally hyperbolic, and normally attracting. The point( S , I , δ ) is in the stable fiber ( S , I , ǫ >
0, the set ˜ K remains normally hyperbolic and normallyattracting. ( S , I , δ ) is in the stable fiber of a point ( S ǫ , I ǫ ,
1) near( S , I , x = 1, is still(4.15)–(4.16). The solution of (4.15)–(4.16) through ( S ǫ , I ǫ ) lies nearthe solution of (4.15)–(4.16) through ( S , I ).Given these observations, the proposition follows from the theory ofnormally hyperbolic invariant manifolds. (cid:3) For each ( S , I ) ∈ V + , there exists t in ( S , I ) > t out ( S , I ) ≥ t in ( S , I ) such that ψ (cid:0) ( S , I ) , t (cid:1) enters the region I ≥ km n − m a at t = t in ( S , I ) and leaves that region at t = t out ( S , I ). Proposition 3.
Let K be a compact subset of V + . Assume that for each ( S , I ) ∈ K , R t out ( S ,I )0 k − ( m n − m a ) I ( t ) dt < , where I ( t ) is definedby ( S ( t ) , I ( t )) = ψ (cid:0) ( S , I ) , t (cid:1) . Choose t > sup K t out (S , I ) . Define Q : K → V − by Q ( S , I ) = ψ (cid:0) ( S , I ) , t (cid:1) . Let δ > be small. Forsmall ǫ > , define Q ǫ ( S , I ) and Z ǫ ( S , I ) by ( Q ǫ ( S , I ) , Z ǫ ( S , I )) = φ ǫ (cid:0) ( S , I , δ ) , t (cid:1) . Then: (1) Q ǫ , Z ǫ and Q are C r functions. (2) Q ǫ → Q in the C r sense. (3) There exists
A > such that Z ǫ ≤ δe − At .Proof. See Proposition 3 of [2] and the remark that follows it. Thekey assumption needed is that for ( S , I ) ∈ K and 0 < t ≤ t , R t k − ( m n − m a ) I ( t ) dt < (cid:3) Proposition 4.
Proposition 2 also holds for a compact subset K of V + that satisfies the assumption of Proposition 3.Proof. Roughly speaking, we want to apply Proposition 2 to the com-pact set φ ( K, t ) in V − . However, corresponding to the point ( S , I ) = φ (cid:0) ( S , I ) , t (cid:1) , we want to look not at the solution of (4.4)–(4.6) that NALYSIS OF AN EPIDEMIC MODEL 25 starts at ( S , I , δ ), but at the solution that starts at φ ǫ (cid:0) ( S , I , δ ) , t (cid:1) .This requires minor changes to Proposition 2. (cid:3) In the situation of Proposition 2 or 3, we can also describe the lim-iting position of orbits.
Proposition 5.
Let K be a compact subset of V − that satisfies theassumption of Proposition 2, or a compact subset of V + that satisfiesthe assumption of Proposition 3. Let ( S , I ) ∈ K . Then there is anequilibrium ( S f , of (4.15) – (4.16) such that φ ( S , I ) , t ) → ( S f , as t → ∞ . Let Γ ǫ denote the orbit of (4.4) – (4.6) that starts at ( S , I , δ ) .As ǫ → , Γ ǫ approaches the singular orbit of (4.4) – (4.6) consisting of (1) the line segment [( S , I , δ ) , ( S , I , ; (2) { φ ( S , I ) , t ) : t ≥ } . Proof of Theorem 1.
We will only consider one type of singu-lar orbit; the proof for other types is similar. Let ( S , I , x ) satisfy( S , I ) ∈ ˆ T , I > km n − m a , and 0 < x <
1. We will consider thefollowing singular orbit S :(1) Fast orbit from ( S , I , x ) to ( S , I , from ( S , I ,
0) to ( S , I ,
0) with I < km n − m a and I (cid:0) ( S , I ) , ( S , I ) (cid:1) = 0.(3) Fast orbit from ( S , I ,
0) to ( S , I , from ( S , I ,
1) to ( S f , , δ >
0, let E = { ( S, I, x ) ∈ P : x = δ } and E = { ( S, I, x ) ∈ P : x = 1 − δ } . For a small ǫ >
0, let Γ ǫ be the orbit of(4.4)–(4.6) that starts at ( S , I , x ). We break Γ ǫ into parts.(1) Γ ǫ from ( S , I , x ) to ( ˜ S , ˜ I , δ ) ∈ E .(2) Γ ǫ from ( ˜ S , ˜ I , δ ) to the next intersection with E at ( ˜ S , ˜ I , δ ).(3) Γ ǫ from ( ˜ S , ˜ I , δ ) to ( ˜˜ S , ˜˜ I , − δ ) ∈ E .(4) Γ ǫ from ( ˜˜ S , ˜˜ I , − δ ) ∈ E to an equilibrium ( ˜ S f , , ǫ → δ ≤ x ≤ x , wesee that Γ ǫ converges to the line segment [ S , I , x ) , ( S , I , δ )].(2) From Theorem 2, Γ ǫ converges to the union of the line seg-ment [( S , I , δ ) , ( S , I , , and the line segment[( S , I , , ( S , I , δ )].(3) By considering the fast system (4.4)–(4.6) on δ ≤ x ≤ − δ , wesee Γ ǫ converges to the line segment [( S , I , δ ) , ( S , I , − δ )]. (4) From Proposition 5, Γ ǫ converges to the union of the line seg-ment [( S , I , − δ ) , ( S , I , . Moreover, thelimiting equilibrium ( S ǫf , ,
1) converges to ( S f , , Discussion
One mathematical issue has been left hanging. Theorem 1 appliesonly to singular orbits of finite length. I suspect that all singular orbitshave finite length, but have not been able to prove it.The model discussed in this paper could be generalized in several tan-talizing directions.. One is to replace the susceptible group by severalsubgroups with different payoff functions. The groups could represent,for example, those with sufficient resources to survive staying home,or with the ability to work from home, and those who need to workoutside the home. A second direction, suggested by the covid-19 pan-demic, is to replace the infective group by subgroups. There could bea group that is infected, and infective, but so far asymptomatic, so un-aware of being infective. Those in this group would continue to use thebehavior they used when susceptible. Some in this group would laterbecome symptomatic; they would presumably change their behavior atthis point.
References [1] W. Booth, A chilling scientific paper helped upend U.S. and U.K. coronavirusstrategies, Washington Post, 17 March 2020, https://tinyurl.com/t8nfc34
NALYSIS OF AN EPIDEMIC MODEL 27 [8] Jones, C. K. R. T., Geometric singular perturbation theory, in DynamicalSystems (Montecatini Terme, 1994), Lecture Notes in Math. , Springer,New York, 1995, pp. 44–118.[9] C. Kuehn, Multiple Time Scale Dynamics, Applied Mathematical Sciences 191,Springer, New York, 2015.[10] M. Li, W. Liu, C. Shan, Y. Yi, Turning points and relaxation oscillation cyclesin simple epidemic models, SIAM J. Appl. Math. 76 (2016), 663–687.[11] W. Liu, Exchange lemmas for singular perturbation problems with certainturning points, J. Differential Equations 167 (2000),13–180.[12] Nowak, M.A., Sigmund, K., Evolutionary dynamics of biological games, Sci-ence 303 (2004), 793–799.[13] P. Poletti, Human behaviour in epidemic modelling, Ph.D. thesis, Universityof Trento, 2010, http://eprints-phd.biblio.unitn.it/422/1/tesi.pdf .[14] P. Poletti, B. Caprile, M. Ajelli, A. Pugliese, S. Merler, Spontaneous be-havioural changes in response to epidemics, J. Theoret. Biol. 260 (2009), 31–40.[15] P. Poletti, M. Ajelli, S. Merler, Risk perception and effectiveness of uncoordi-nated behavioral responses in an emerging epidemic, Math. Biosci. 238 (2012),8089.[16] Verelst F, Willem L, Beutels P. 2016 Behavioural change models for infectiousdisease transmission: a systematic review (20102015). J. R. Soc. Interface 13:20160820. http://dx.doi.org/10.1098/rsif.2016.0820.
Appendix A. Matlab routines
The file findsingorbit.m is used to find a singular orbit. Parametervalues are entered in the file epimconstants.m. The files entryexitint0.mand entryexitint1.m are used by findsingorbit.m to evaluate entry-exitintegrals in x = 0 and x = 1 respectively. epimconstants.m % Constants used by other functions.epsilon = 0.005;betan = 0.5;betaa = 0.1;gam = 1/6;k = 0.3;mn = 5;ma = 2; entryexitint0.m % Entry-exit function in the plane x=0. function y = entryexitint0(S0,S1,v0)% In the plane x=0, evaluates the integral from S0 to S1, with parameter v0,% that is used to define the entry-exit function. Given S0 we will want to% find S1 such that the integral is 0.syms S;y = vpaintegral(top(S,v0)/bottom(S,v0),S0,S1);returnendfunction y1 = main(S,v0)epimconstants;syms S;y1 = v0-S+(gam/betaa)*log(S);returnendfunction y2 = top(S,v0)epimconstants;syms S;y2 = k-(mn-ma)*main(S,v0);returnendfunction y3 = bottom(S,v0);epimconstants;syms S;y3 = -betaa*S*main(S,v0);returnend entryexitint1.m function y = entryexitint1(S0,S1,v0)% In the plane x=1, evaluates the integral from S0 to S1, with parameter v0,% that is used to define the entry-exit function. Given S0 we will want to% find S1 such that the integral is 0.syms S;y = vpaintegral(top(S,v0)/bottom(S,v0),S0,S1); NALYSIS OF AN EPIDEMIC MODEL 29 returnendfunction y1 = main(S,v0)epimconstants;syms S;y1 = v0-S+(gam/betan)*log(S);returnendfunction y2 = top(S,v0)epimconstants;syms S;y2 = k-(mn-ma)*main(S,v0);returnendfunction y3 = bottom(S,v0);epimconstants;syms S;y3 = betan*S*main(S,v0);returnend findsingorbit.m function singorbit = findsingorbit(S0,I0,x0)epimconstantssyms usingorbit = [S0 I0 x0];while 0 < S0 & 0 < I0 & S0+I0 < =1 & 0 < x0 & x0 < ~ =k/(mn-ma)if I0 > k/(mn-ma)disp('solution is attracted to x=0')% In this case the solution is attracted to the plane x=0 and will% follow a solution in that plane.v0=I0+S0-(gam/betaa)*log(S0);% The solution will arrive at (S1,I1) with I1=0 where u=S1% satisfies: eqn1 = u-(gam/betaa)*log(u)==v0;S1=vpasolve(eqn1,u,[0,S0]);% The solution in SI-space will arrive at (S2,I2) with I2=k/(mn-ma)% where u=S2 satisfies:eqn2 = u-(gam/betaa)*log(u)==v0-k/(mn-ma);% Find where the solution in SI-space crosses the line I2=k/(mn-ma).disp('solution leaves the plane x=0')S2=vpasolve(eqn2,u,[S1,S0]);syms uEnd;eqn3 = entryexitint0(S0,uEnd,v0)==0;u4=vpasolve(eqn3,uEnd,[S1,S2]);% Continue solution from following point.S0=u4;I0=v0-S0+(gam/betaa)*log(S0);x0=0.1;singorbit = [singorbit;[S0 I0 x0]];elsedisp('solution is attracted to the plane x=1')% In this case the solution is attracted to the plane x=1 and will% follow a solution in that plane.v0=I0+S0-(gam/betan)*log(S0);% The solution will arrive at (S1,I1) with I1=0 where u=S1% satisfies:eqn1 = u-(gam/betan)*log(u)==v0;S1=vpasolve(eqn1,u,[0,S0]);% The following is the value of v0 for which the integral curve in% SI-space has its max on the line I=k/(mn-ma).v0tan = (gam/betan)-(gam/betan)*log(gam/betan)+k/(mn-ma);if v0 < = v0tandisp('entire solution lies below I=k/(mn-ma), solution terminates')% In this case the solution will not leave the plane x=1.S0=S1;I0=0;x0=1;singorbit = [singorbit;[S0 I0 x0]];elsedisp('solution crosses the line I=k/(mn-ma)')% In this case the solution in SI-space meets the line NALYSIS OF AN EPIDEMIC MODEL 31 % I=k/(mn-ma) in two points with S-values S3 > S2.eqn2 = u-(gam/betan)*log(u)==v0-k/(mn-ma);S3=vpasolve(eqn2,u,[S1,gam/betan]);int2 = entryexitint1(S0,S3,v0);if int2 < =0disp('nevertheless solution terminates')% In this case the solution will not leave the plane x=1.S0=S1;I0=0;x0=1;singorbit = [singorbit;[S0 I0 x0]];elsedisp('solution leaves the plane x=1')S2=vpasolve(eqn2,u,[gam/betan,S0]);syms uEndeqn3 = entryexitint1(S0,uEnd,v0)==0;u4=vpasolve(eqn3,uEnd,[S3,S2]);% Continue solution from following point.S0=u4;I0=v0-S0+(gam/betan)*log(S0);x0=0.9;singorbit = [singorbit;[S0 I0 x0]];endendendendreturnend Department of Mathematics, North Carolina State University, Box8205, Raleigh, NC 27695 USA
E-mail address ::