key: cord-104133-d01joq23
authors: Arthur, Ronan F.; Jones, James H.; Bonds, Matthew H.; Ram, Yoav; Feldman, Marcus W.
title: Adaptive social contact rates induce complex dynamics during epidemics
date: 2020-07-14
journal: bioRxiv
DOI: 10.1101/2020.04.14.028407
sha: 
doc_id: 104133
cord_uid: d01joq23

The COVID-19 pandemic has posed a significant dilemma for governments across the globe. The public health consequences of inaction are catastrophic; but the economic consequences of drastic action are likewise catastrophic. Governments must therefore strike a balance in the face of these trade-offs. But with critical uncertainty about how to find such a balance, they are forced to experiment with their interventions and await the results of their experimentation. Models have proved inaccurate because behavioral response patterns are either not factored in or are hard to predict. One crucial behavioral response in a pandemic is adaptive social contact: potentially infectious contact between people is deliberately reduced either individually or by fiat; and this must be balanced against the economic cost of having fewer people in contact and therefore active in the labor force. We develop a model for adaptive optimal control of the effective social contact rate within a Susceptible-Infectious-Susceptible (SIS) epidemic model using a dynamic utility function with delayed information. This utility function trades off the population-wide contact rate with the expected cost and risk of increasing infections. Our analytical and computational analysis of this simple discrete-time deterministic model reveals the existence of a non-zero equilibrium, oscillatory dynamics around this equilibrium under some parametric conditions, and complex dynamic regimes that shift under small parameter perturbations. These results support the supposition that infectious disease dynamics under adaptive behavior-change may have an indifference point, may produce oscillatory dynamics without other forcing, and constitute complex adaptive systems with associated dynamics. Implications for COVID-19 include an expectation of fluctuations, for a considerable time, around a quasi-equilibrium that balances public health and economic priorities, that shows multiple peaks and surges in some scenarios, and that implies a high degree of uncertainty in mathematical projections. Author summary Epidemic response in the form of social contact reduction, such as has been utilized during the ongoing COVID-19 pandemic, presents inherent tradeoffs between the economic costs of reducing social contacts and the public health costs of neglecting to do so. Such tradeoffs introduce an interactive, iterative mechanism which adds complexity to an infectious disease system. Consequently, infectious disease modeling typically has not included dynamic behavior change that must address such a tradeoff. Here, we develop a theoretical model that introduces lost or gained economic and public health utility through the adjustment of social contact rates with delayed information. We find this model produces an equilibrium, a point of indifference where the tradeoff is neutral, and at which a disease will be endemic for a long period of time. Under small perturbations, this model exhibits complex dynamic regimes, including oscillatory behavior, runaway exponential growth, and eradication. These dynamics suggest that for epidemic response that relies on social contact reduction, secondary waves and surges with accompanied business re-closures and shutdowns may be expected, and that accurate projection under such circumstances is unlikely.

The COVID-19 pandemic had infected almost 9 million people and caused over 450,000 17 deaths worldwide as of June 23, 2020 [1] . In the absence of effective therapies and 18 vaccines [2] , many governments responded with lock-down policies and social distancing 19 laws to reduce the rate of social contacts and curb transmission of the virus. Prevalence 20 of COVID-19 in the wake of these policies in the United States indicates they may have 21 been successful at decreasing the reproduction number (R t ) of the epidemic [1] . 22 However, they have also led to economic recession with an unemployment rate at an 23 80-year peak, the stock market in decline, and the federal government forced to borrow 24 heavily to financially support businesses and households. Solutions to these economic 25 crises may conflict with public health recommendations. Thus, governments worldwide 26 must decide how to balance the economic and public health consequences of their 27 epidemic response interventions. 28 Behavior-change in response to an epidemic, whether autonomously adopted by 29 individuals or externally directed by governments, affects the dynamics of infectious 30 diseases [3, 4] . Prominent examples of behavior-change in response to infectious disease 31 prevalence include measles-mumps-rubella (MMR) vaccination choices [5] , social 32 distancing in influenza outbreaks [6] , condom purchases in HIV-affected communities [7], 33 and social distancing during the ongoing COVID-19 pandemic [2] . Behavior is 34 endogenous to an infectious disease system because it is, in part, a consequence of the 35 prevalence of the disease, which in turn responds to changes in behavior [8, 9] . 36 Individuals and governments have greater incentive to change behavior as prevalence 37 increases; conversely they have reduced incentive as prevalence decreases [10, 11] . 38 Endogenous behavioral response may then theoretically produce a non-zero endemic 39 equilibrium of infection. This happens because, at low levels of prevalence, the cost of 40 avoidance of a disease may be higher than the private benefit to the individual, even 41 though the collective, public benefit in the long-term may be greater. However, in 42 epidemic response we typically think of behavior-change as an exogenously-induced 43 intervention without considering associated costs. While guiding positive change is an 44 important intervention, neglecting to recognize the endogeneity of behavior can lead to 45 a misunderstanding of incentives and a resurgence of the epidemic when behavior 46 change is reversed prematurely. 47 Although there is growing interest in the role of adaptive human behavior in 48 infectious disease dynamics, there is still a lack of general understanding of the most 49 important properties of such systems [3, 8, 12] . Behavior is difficult to measure, quantify, 50 or predict [8] , in part, due to the complexity and diversity of human beings who make simply allowed the transmission parameter (β) to be a negative function of the number 56 infected, effectively introducing an intrinsic negative feedback to the infected class that 57 regulated the disease [13] . 58 Modelers have used a variety of tools, including agent-based modeling [14] , network 59 structures for the replacement of central nodes when sick [15] or for behavior-change as 60 a social contagion process [16] , game theoretic descriptions of rational choice under 61 changing incentives as with vaccination [6, 11, 17] , and a branching process for 62 heterogeneous agents and the effect of behavior during the West Africa Ebola epidemic 63 in 2014 [18] . A common approach to incorporating behavior into epidemic models is to 64 track co-evolving dynamics of behavior and infection [16, 19, 20] , where behavior 65 represents an i-state of the model [21] . In a compartmental model, this could mean 66 separate compartments (and transitions therefrom) for susceptible individuals in a state 67 of fear and those not in a state of fear [16] .

68 Periodicity (i.e. multi-peak dynamics) has long been documented empirically in 69 epidemiology [22, 23] . Periodicity can be driven by seasonal contact rate changes (e.g. 70 when children are in school) [24] , seasonality in the climate or ecology [25] , sexual 71 behavior change [26] , and host immunity cycling through new births of susceptibles or a 72 decay of immunity over time. Some papers in nonlinear dynamics have studied delay 73 differential equations in the context of epidemic dynamics and found periodic 74 solutions [27] . Although it is atypical to include delay in modeling, delay is an 75 important feature of epidemics. For example, if behavior responds to mortality rates, 76 there will inevitably be a lag with an average duration of the incubation period plus the 77 life expectancy upon becoming infected. In a tightly interdependent system, reacting to 78 outdated information can result in an irrational response and periodic cycling. The original epidemic model of Kermack and McKendrick [28] was first expressed in 92 discrete time. Then by allowing "the subdivisions of time to increase in number so that 93 each interval becomes very small" the famous differential equations of the SIR epidemic 94 model were derived. Here we begin with a discrete-time Susceptible-Infected-Susceptible 95 model that is adjusted on the principle of endogenous behavior-change through an 96 adaptive social-contact rate that can be thought of as either individually motivated or 97 institutionally imposed. We introduce a dynamic utility function that motivates the 98 population's effective contact rate at a particular time period. This utility function is 99 based on information about the epidemic size that may not be current. This leads to a 100 time delay in the contact function that increases the complexity of the population 101 dynamics of the infection. Results from the discrete-time model show that the system 102 approaches an equilibrium in many cases, although small parameter perturbations can 103 lead the dynamics to enter qualitatively distinct regimes. The analogous continuous-time model retains periodicities for some sets of parameters, but numerical 105 investigation shows that the continuous time version is much better behaved than the 106 discrete-time model. This dynamical behavior is similar to models of ecological 107 population dynamics, and a useful mathematical parallel is drawn between these 108 systems. To represent endogenous behavior-change, we start with the classical discrete-time 112 susceptible-infected-susceptible (SIS) model [28] , which, when incidence is relatively 113 small compared to the total population [29, 30] , can be written in terms of the recursions 114

where at time t, S t represents the number of susceptible individuals, I t the infected 115 individuals, and N t the number of individuals that make up the population, which is 116 assumed fixed in a closed population. We can therefore write N for the constant 117 population size. Here γ, with 0 < γ < 1, is the rate of removal from I to S due to 118 recovery. This model in its simplest form assumes random mixing, where the parameter 119 b represents a composite of the average contact rate and the disease-specific 120 transmissibility given a contact event. In order to introduce human behavior, we 121 substitute for b a time-dependent b t , which is a function of both b 0 , the probability that 122 disease transmission takes place on contact, and a dynamic social rate of contact c t 123 whose optimal value, c * t , is determined at each time t as in economic epidemiological 124 models [31] , namely

where c * t represents the optimal contact rate, defined as the number of contacts per unit 126 time that maximize utility for the individual. Here, c * t is a function of the number of This utility function is assumed to take the form

Here U represents utility for an individual at time t given a particular number of 134 contacts per unit time c, α 0 is a constant that represents maximum potential utility 135 achieved at a target contact rateĉ. The second term, α 1 (c −ĉ) 2 , is a concave function 136 that represents the penalty for deviating fromĉ. The third term, the delay in information acquisition and the speed of response to that information. We 140 note that (1 − I N b 0 ) c can be approximated by

when I N b 0 is small and c I N b 0 << 1. We thus assume I N (b 0 ) is small, and approximate 142 U (c) in Eq. 5 using Eq. 6. Eq. 5 assumes a strictly negative relationship between 143 number of infecteds and contact. 144 We assume an individual or government will balance the cost of infection, the 145 probability of infection, and the cost of deviating from the target contact rateĉ to select 146 an optimal contact rate c * t , namely the number of contacts, which takes into account the 147 risk of infection and the penalty for deviating from the target contact rate. This 148 captures the idea that individuals trade off how many people they want to interact with 149 versus their risk of getting sick, or that authorities want to reopen the economy during a 150 pandemic and have to trade off morbidity and mortality from increasing infections with 151 the need to allow additional social contacts to help the economy restart. This optimal 152 contact rate can be calculated by finding the maximum of U with respect to c from Eq. 153 5 with substitution from Eq. 6, namely

Differentiating, we have

which vanishes at the optimal contact rate, c * , which we write as c * t to show its 156 dependence on time. Then

which we assume to be positive. Therefore, total utility will decrease as I t increases and 158 c * t also decreases. Utility is maximized at each time step, rather than over the course of 159 lifetime expectations. In addition, Eq 9 assumes a strictly negative relationship between 160 number of infecteds at time t − ∆ and c * . While behavior at high degrees of prevalence 161 has been shown to be non-linear and fatalistic [32, 33] , in this model, prevalence (i.e., 162 b0It N ) is assumed to be small, consistent with Eq. 6. 163 We introduce the new parameter α = α2

We can now rewrite the recursion from Eq. 2, using Eq. 4 and replacing c t with c * t as 165 defined by Eq. 10, as

When ∆ = 0 and there is no time delay, f (·) is a cubic polynomial, given by

July 10, 2020 5/17

For the susceptible-infected-removed (SIR) version of the model, we include the removed category and write the (discrete-time) recursion system as

the baseline contact rate and c * t specified by 169 Eq. 10. With b t = b, say, and not changing over time, Eqs. 13-15 form the discrete-time 170 version of the classical Kermack-McKendrick SIR model [28] . The inclusion of the 171 removed category entails thatĨ = 0 is the only equilibrium of the system Eqs. 13-15; 172 unlike the SIS model, there is no equilibrium with infecteds present. In general, since c * t 173 includes the delay ∆, the dynamic approach toĨ = 0 is expected to be quite complex. 174 Intuitively, since the infecteds are ultimately removed, we do expect that from any 175 initial frequency I 0 of infecteds all N individuals will eventually be in the R category.

Numerical analysis of this SIR model shows strong similarity between the SIS and SIR 177 models for several hundred time steps before the SIR model converges toĨ = 0 with 178 R = N . In the section "Numerical Iteration and Continuous-Time Analog" we compare 179 the numerical iteration of the SIS (Eq. 11) and SIR (Eqs. [13] [14] [15] and integration of the 180 continuous-time (differential equation) versions of the SIS and SIR models. To determine the dynamic trajectories of (11) without time delay, we first solve for the 184 fixed point(s) of the recursion (11) (i.e., value or values of I such that

From Eq. 16, it is clear that I = 0 is an equilibrium as no new infections can occur 187 in the next time-step if none exist in the current one. This is the disease-free 188 equilibrium denoted byĨ. Other equilibria are the solutions of

We label the solution with the + sign I * and the one with the − signÎ. 

It is important to note that under these conditionsÎ is an equilibrium of the for this to hold whenÎ is legitimate is

If inequalities (20) and Nĉb 0 > γ hold, thenÎ is locally stable. However, even if 210 both of these inequalities hold, the number of infecteds may not converge toÎ. It is well 211 known that iterations of discrete-time recursive relations, of which (12) is an example 212 (i.e., with ∆ = 0), may produce cycles or chaos depending on the parameters and the 213 starting frequency I 0 of infecteds.

214 Table 1 shows an array of possible asymptotic dynamics with ∆ = 0 found by 215 numerical iteration of (12) for a specific set of parameters and an initial frequency Table 1 are examples for which, beginning with a 217 single infected, the number of infecteds explodes, becoming unbounded; of course, this is 218 an illegitimate trajectory since I t cannot exceed N . However, in the case marked * ,Î is 219 locally stable and with a large enough initial number of infecteds, there is damped 220 oscillatory convergence toÎ. In the case marked * * , with I 0 = 1 the number of infecteds 221 becomes unbounded, but in this case,Î is locally unstable, and starting with I 0 close to 222 I a stable two-point cycle is approached; in this case df (I)/dI| I=Î < −1. Table 1 .

Stability analysis of the SIS model is more complicated when ∆ = 0, and in the 224 Appendix we outline the procedure for local analysis of the recursion (11) nearÎ. Local 225 stability is sensitive to the delay time ∆ as can be seen from the numerical iteration of 226 July 10, 2020 7/17 (11) for the specific set of parameters shown in Table 2 . Some analytical details related 227 to Table 2 are in the Appendix. Table 1 reports an array of dynamic trajectories for some choices of parameters and, 240 in two cases, an initial number of infecteds other than I 0 = 1. The first three rows show 241 three sets of parameters for which the equilibrium values ofÎ are very similar but the 242 trajectories of I t are different: a two-point cycle, a four-point cycle, and apparently 243 chaotic cycling above and belowÎ. In all of these cases, df (I)/dI| I=Î < −1. Clearly the 244 dynamics are sensitive to the target contact rateĉ in these cases. The fourth and eighth 245 rows show that I t becomes unbounded (tends to +∞) from I 0 = 1, but a two-point 246 cycle is approached if I 0 is close enough toÎ: df (I)/dI| I=Î < −1 in this case. For the 247 parameters in the ninth row, if I 0 is close enough toÎ there is damped oscillation intoÎ: 248 here −1 < df (I)/dI| I=Î < 0.

The fifth and sixth rows of Table 1 exemplify another interesting dynamic starting 250 from I 0 = 1. I t becomes larger thanÎ (overshoots) and then converges monotonically 251 down toÎ; in each case 0 < df (I)/dt| I=Î < 1. For the parameters in the seventh row, 252 there is oscillatory convergence toÎ from I 0 = 1 (−1 < df (I)/dI| I=Î < 0), while in the 253 last row there is straightforward monotone convergence toÎ.

A continuous-time analog of the discrete-time recursion (11), in the form of a 255 differential equation, substitutes dI/dt for I t+1 − I t in (11). We then solve the resulting 256 delay differential equation numerically using the VODE differential equation integrator 257 in SciPy [36, 37] (source code available at https://github.com/yoavram/SanJose). Using 258 the parameters in Table 2 Figure 1 with I 0 = 1.

In Figure 1 , with no delay (∆ = 0) and a one-unit delay (∆ = 1), the discrete and 266 continuous dynamics are very similar, both converging toÎ. However, with ∆ = 2 the 267 differential equation oscillates intoÎ while the discrete-time recursion enters a regime of 268 inexact cycling aroundÎ, which appears to be a state of chaos. For ∆ = 3 and ∆ = 4, 269 the discrete recursion "collapses". In other words, I t becomes negative and appears to 270 go off to −∞; in Figure 1 , this is cut off at I = 0. The continuous version, however, in 271 these cases enters a stable cycle aroundÎ. It is important to note that in Figure 1 for respectively. In Fig. S5S there appears to be convergence toÎ, but in Fig. S5L after 302 about 500 time units, in both discrete-and continuous-time SIR versions, the number of 303 infected begins to decline towards zero.

It is worth noting that if the total population size of N decreases over time, for 305 example, if we take N (t) = N exp(−zt), with z = 50b 0ĉ γ, then the short-term dynamics 306 of the SIS model in (11) begins to closely resemble the SIR version. This is illustrated 307 in Supplementary Fig. S5N , where b 0 ,ĉ, γ are, as in Figs. S5S and S5L, the same as in 308 Fig. 2, panel (a) . With N decreasing to zero, both S and I will approach zero in the Our model makes a number of simplifying assumptions. We assume, for example, 320 that all individuals in the population will respond in the same fashion to government 321 policy. We assume that governments choose a uniform contact rate according to an 322 optimized utility function, which is homogeneous across all individuals in the population. 323 Finally, we assume that the utility function is symmetric around the optimal number of 324 contacts so that increasing or decreasing contacts above or below the target contact 325 rate, respectively, yield the same reduction in utility. These assumptions allowed us to 326 create the simplest possible model that includes adaptive behavior and time delay. In

Holling's heuristic distinction in ecology between tactical models, models built to be 328 parameterized and predictive, and strategic models, which aim to be as simple as 329 possible to highlight phenomenological generalities, this is a strategic model [38] . 330 We note that the five distinct kinds of dynamical trajectories seen in these 331 computational experiments come from a purely deterministic recursion. This means 332 that oscillations and even erratic, near-chaotic dynamics and collapse in an epidemic 333 may not necessarily be due to seasonality, complex agent-based interactions, changing or 334 stochastic parameter values, demographic change, host immunity, or socio-cultural 335 idiosyncracies. This dynamical behavior in number of infecteds can result from 336 mathematical properties of a simple deterministic system with homogeneous endogenous 337 behavior-change, similar to complex population dynamics of biological organisms [39] .

The mathematical consistency with population dynamics suggests a parallel in ecology, 339 that the indifference point for human behavior functions in a similar way to a carrying 340 capacity in ecology, below which a population will tend to grow and above which a individuals are incentivized to change their behavior to protect themselves, they will, 346 and they will cease to do this when they are not [10] . Further, our results show certain 347 parameter sets can lead to limit-cycle dynamics, consistent with other negative feedback 348 mechanisms with time delays [41, 42] . This is because the system is reacting to 349 conditions that were true in the past, but not necessarily true in the present. In our 350 discrete-time model, there is the added complexity that the non-zero equilibrium may 351 be locally stable but not attained from a wide range of initial conditions, including the 352 most natural one, namely a single infected individual.

Observed epidemic curves of many transient disease outbreaks typically inflect and 354 go extinct, as opposed to this model that may oscillate perpetually or converge [43] , and surges in fluctuations in COVID-19 cases globally [1] . There may be 363 many causes for such double-peaked outbreaks, one of which may be a lapse in 364 behavior-change after the epidemic begins to die down due to decreasing incentives [16] , 365 as represented in our simple theoretical model. This is consistent with findings that 366 voluntary vaccination programs suffer from decreasing incentives to participate as 367 prevalence decreases [44, 45] . It should be noted that the continuous-time version of our 368 model can support a stable cyclic epidemic whose interpretation in empirical terms will 369 depend on the time scale, and hence on the meaning of the delay, ∆.

One of the responsibilities of infectious disease modelers (e.g. COVID-19 modelers) 371 is to predict and project forward what epidemics will do in the future in order to better 372 assist in the proper and strategic allocation of preventative resources. COVID-19 373 models have often proved wrong by orders of magnitude because they lack the means to 374 account for adaptive response. An insight from this model, however, is that prediction 375 becomes very difficult, perhaps impossible, if we allow for adaptive behavior-change 376 because the system is qualitatively sensitive to small differences in values of key 377 parameters. These parameters are very hard to measure precisely; they change 378 depending on the disease system and context and their inference is generally subject to 379 large errors. Further, we don't know how policy-makers weight the economic trade-offs 380 against the public health priorities (i.e., the ratio between α 1 and α 2 in our model) to 381 arrive at new policy recommendations.

To maximize the ability to predict and minimize loss of life or morbidity, outbreak 383 response should not only seek to minimize the reproduction number, but also the length 384 of time taken to gather and distribute information. Another approach would be to use a 385 predetermined strategy for the contact rate, as opposed to a contact rate that depends 386 on the number of infecteds.

In our model, complex dynamic regimes occur more often when there is a time delay. 388 If behavior-change arises from fear and fear is triggered by high local mortality and high 389 local prevalence, such delays seem plausible since death and incubation periods are 390 lagging epidemiological indicators. Lags mean that people can respond sluggishly to an 391 unfolding epidemic crisis, but they also mean that people can abandon protective 392 behaviors prematurely. Developing approaches to incentivize protective behavior 393 throughout the duration of any lag introduced by the natural history of the infection (or 394 otherwise) should be a priority in applied research. This paper represents a first step in 395 understanding endogenous behavior-change and time-lagged protective behavior, and we 396 anticipate further developments along these lines that could incorporate long incubation 397 periods and/or recognition of asymptomatic transmission. 

In the neighborhood of the equilibriumÎ, write I t =Î + ε t and I t−∆ =Î + ε t−∆ , where 437 ε t and ε t−∆ are small enough that quadratic terms in them can be neglected in the 438 expression for I t+1 =Î + ε t+1 . The linear approximation to (A1) is then

and in the case ∆ = 0, this reduces to

We focus first on ∆ = 0 and write (A3) as ε t+1 = ε t L(Î). Recall thatÎ satisfies Eq. 441 (17) , and substituting γ from (17) 

Now we turn to the general case ∆ = 0 and Eq. (A2), which we write as

where A and B are the corresponding terms on the right side of (A2 constants with respect to time. Local stability ofÎ is then determined by the properties 450 of recursion (A7), whose solution first involves solving its characteristic equation

In principle there are ∆ + 1 real or complex roots of (A8), which we represent as 452 λ 1 , λ 2 , . . . , λ ∆+1 , and the solution of (A7) can be written as

where c i are found from the initial conditions. Convergence to, and hence local stability 454 ofÎ, is determined by the magnitude of the absolute value (if real) or modulus (if complex) of the roots λ 1 , λ 2 , . . . , λ ∆+1 :Î is locally stable if the largest among the ∆ + 1 456 of these is less than unity.

In Table 2 , results of numerically iterating the complete recursion (11) are listed for 458 the delay ∆ varying from ∆ = 0 to ∆ = 4, all starting from I 0 = 1, with N = 10, 000 459 and the stated parameters. Figure 3 illustrates the discrete-and continuous-time 460 dynamics summarized in Table 2 

with complex roots 0.4999 ± 0.6461i whose modulus is 0.8169, which is less than 1. The 465 complexity implies cyclic behavior, and since the modulus is less than one, we see 466 locally damped oscillatory convergence toÎ.

For ∆ = 2, the characteristic equation is the cubic

which has one real root 0.6383 and complex roots 0.8190 ± 0.6122i. Here the modulus of 469 the complex roots is 1.0225, which is greater than unity so thatÎ is not locally stable. 470 In this case the dynamics depend on the initial value I 0 . If I 0 < 72, I t oscillates but not 471 in a stable cycle. If I 0 > 73, the oscillation becomes unbounded. 

World Health Organization. Coronavirus disease (COVID-19): situation report

Scientific and ethical basis for social-distancing 489

interventions against COVID-19. The Lancet Infectious Diseases

Social factors in epidemiology

Modelling the influence of human behaviour on 493 the spread of infectious diseases: a review

Evolving public perceptions and stability in 496 vaccine uptake

Game theory of social distancing in response to an epidemic

The responsiveness of the demand for condoms 500 to the local prevalence of AIDS

Nine 502 challenges in incorporating the dynamics of behaviour in infectious diseases 503 models

impact and behaviour: the importance of social forces to 506 infectious disease dynamics and disease ecology

Economic epidemiology and infectious diseases

Erratic 511 Flu Vaccination Emerges from Short-Sighted Behavior in Contact Networks

Capturing human behaviour

A generalization of the Kermack-McKendrick deterministic 516 epidemic model

A hybrid epidemic model: 518 combining the advantages of agent-based and equation-based approaches

Winter. IEEE

The effect of a prudent adaptive 521 behaviour on disease transmission

Coupled contagion 523 dynamics of fear and disease: mathematical and computational explorations

A general approach for population games with 526 application to vaccination

Ebola cases and health system demand in Liberia

The spread of awareness and its 532

Proceedings of the National Academy of Sciences

A review

The Dynamics of Physiologically Structured Populations

Periodicity in epidemiological models

Measles in England and Wales-I: an analysis of factors 544 underlying seasonal patterns

Seasonal and interannual cycles of endemic cholera in 547

Bengal 1891-1940 in relation to climate and geography

Etiology of Newly Emerging Marine Diseases

Epidemic cycles driven by host behaviour

Periodic solutions of delay differential equations 552 arising in some models of epidemics

A contribution to the mathematical theory of 555 epidemics

The Royal Society

Modeling infectious diseases in humans and animals

Princeton University Press

Time series modelling of childhood diseases

dynamical systems approach

Adaptive human behavior in epidemiological models

Choices, beliefs, and infectious disease dynamics

Higher disease prevalence can 568 induce greater sociality: a game theoretic coevolutionary model

Global stability of an SIR epidemic 571

Global stability for the SEIR model in epidemiology

Scipy-based delay differential equation (dde) solver

The strategy of building models of complex ecological systems

Simple mathematical models with very complicated dynamics

Journal of the Fisheries Board of Canada

Time-delay versus stability in population models with two and three 587 trophic levels

Time delays are not necessarily destabilizing

Different epidemic curves for severe acute respiratory 591

Rational epidemics and their public control

Group interest versus self-interest in smallpox 596 vaccination policy