Complex dynamics induced by delayed adaptive behavior during epidemics

The behavioral response of individuals to an epidemic, and the timing of their behavior change, can have important consequences for epidemic dynamics. While the importance of behavior-change for epidemic dynamics is broadly recognized, the approaches to modeling the coupled dynamics of epidemics and behavior change have been limited. An important mode of behavior change in epidemics is a reduction in potentially infectious contacts. We develop a model for endogenous change of the effective contact rate of Susceptible-Infectious-Susceptible (SIS) epidemic model by positing a dynamic utility function associated with the expected number of contacts people have in the population. This utility function trades off the ideal number of social contacts with the expected cost of becoming infected. Our analysis of this simple, deterministic model reveals the existence of a non-zero endemic equilibrium, oscillatory dynamics under some parametric conditions, and complex dynamic regimes that shift under small parameter perturbations. Time-lag between current epidemic conditions and behavioral response is essential for the more exotic dynamics. These results suggest that reducing this lag may reduce both the final-size of the epidemic and the uncertainty associated with epidemic projections.


Introduction
Human behavior-change in response to an epidemic, whether autonomously adopted by individuals or externally directed by institutions, affects the dynamics of infectious diseases [1,2]. Prominent examples of important behavior change in response to infectious disease dynamics include measles-mumps-rubella (MMR) vaccination choices [3], social distancing in influenza outbreaks [4], condom purchases in HIV-affected communities [5], and non-pharmaceutical interventions in the SARS epidemic [6]. Social choices, driven by incentives, can affect rates and structure of contact, how people interact with the environment, and patterns of space, mobility, and migration, with important consequences for infectious disease dynamics [7]. Changes in personal health decision-making can be psychologically motivated by risk perception, from personal susceptibility to disease severity (i.e., morbidity and mortality), and a perception of self-efficacy, namely the awareness of potential behavior changes and the ability to make them [8].
Behavior should be considered endogenous to an infectious disease system because it is, in part, a consequence of the prevalence and severity of the disease, which also respond to changes in behavior [9]. Individuals have greater incentive to change behavior as prevalence increases; conversely they have reduced incentive as prevalence decreases [10,11]. Endogenous behavioral response may then theoretically produce a non-zero endemic equilibrium of infection. This happens because, at low levels of prevalence, the cost of avoidance of a disease may be higher than the private benefit to the individual, even though the collective, public benefit in the long-term may be far greater. However, in epidemic response we typically think of behavior change as an exogenously-induced intervention. While guiding positive change is an important intervention, neglecting to recognize the endogeneity of behavior can lead to the reversal of those changes when they no longer appear necessary from an individual perspective. This can frustrate efforts to eradicate a disease.
Although there is growing interest in the role of dynamic human behavior in infectious disease dynamics, there is still a lack of general understanding of the most important properties of such systems [1,9,12]. For example, behavior change was a notable, though intentional, omission in a prominent CDC model that predicted outcomes of the recent West African Ebola epidemic in 2014. This model predicted a final size of the epidemic two orders of magnitude greater than what actually occurred [13] and this excess was due, in large measure, to behavior change. Behavior is difficult to measure, quantify, or predict [9], in part, due to the complexity and diversity of human beings who make different choices under different circumstances. It is also highly variable depending on the relative time scales of the disease and behavioral responses.
Despite these challenges, modelers have adopted a variety of strategies to investigate the importance of behavior (see Funk et al., 2010 [2] for a general review). An early expansion of the Kermack-McKendrick Susceptible-Infectious-Removed (SIR) model simply allowed the transmission parameter (β) to be a negative function of the number infected, effectively introducing an intrinsic negative feedback to the infected class that regulates the disease [14]. Modelers have used a variety of tools, including agent-based modeling [15], network structures for the replacement of central nodes when sick [16] or for behavior change as a social contagion process [17], game theoretic descriptions of rational choice under changing incentives as with vaccination [4,11,18], branching process for heterogeneous agents and the effect of behavior during the West Africa Ebola epidemic in 2014 [19], and adaptive network-structured sexual behavior that can produce oscillatory dynamics [20].
A common approach to incorporating behavior into epidemic models is to track coevolving dynamics of behavior and infection [17,21,22], where behavior represents an i-state of the model [23]. In a compartment model, this could mean separate compartments (and transitions therefrom) for susceptible individuals in a state of fear and those not in a state of fear (e.g., Epstein et al. [17]). In contrast, we develop an approach in which overall behavior is modulated dynamically and endogenously simple as a function of the current fraction of the population infected.
Interactions between behavior-mediated transmission and the prevalence of disease depend on the nature and timing of an outbreak of infectious disease. Given enough time for response, behavior is prevalence-elastic. If the behavior in question is slow-changing, it may effectively be regarded as a parameter rather than a variable in a fast-acting outbreak [9], although it may still be important to the long-term interaction between economics, the environment, and disease [24]. Continuous-time modeling of epidemics and adaptive behavior generally assumes perfect instantaneous information and immediate response and does not allow for delay of information or reactions. However, in reality, information about the state of an epidemic can take time to collect and distribute, and the reaction to this information may also be delayed. For example, if behavior responds to mortality rates, there will inevitably be a lag with an average duration of of the incubation period plus the life expectancy upon becoming infected. In a tightly interdependent system, reacting to outdated information can result in an irrational response.
A delay in a behavioral response to an epidemic can take a variety of forms. For example, in the 2014-15 Ebola epidemic, there were delays in the identification of the virus, delays in the acquisition of reliable information on suspected and confirmed cases, delays in the mobilization of international relief efforts by the World Health Organization and affected governments [25], delays in the establishment of trust between local affected communities and health authorities [26], delays in the launching of community information and communication campaigns, and delays in laying the physical infrastructure for the treatment and isolation of an increasing number of patients. Recognizing the effects of such delays as a part of the endogenous relationship between disease and behavior is important to improved understanding of the complexity of disease dynamics.
To address the gap in the literature, we develop and analyze a simple model based on the principle of endogenous behavior change during epidemics of infectious diseases. We introduce a dynamic utility function that determines the population's effective contact rate at a particular time period. This utility function is based on information about the epidemic size that may not be current. Results from the model show that the system approaches an equilibrium in most cases, although small parameter perturbations can lead the dynamics to enter qualitatively distinct regimes. This dynamical behavior is similar to models of ecological population dynamics, and a useful mathematical parallel is drawn between these systems.

Model specifications
To represent endogenous behavior change, we start with a discrete-time Susceptible-Infected-Susceptible (SIS) model [27], which, when incidence is relatively small compared to the total population [28,29], can be written in terms of the recursions where at time t, S t represents the number of susceptible individuals, I t the infected individuals, and N t the number of individuals that make up the population, which is assumed fixed in a closed population. We can therefore remove the time subscript and write N for the population size. Here γ is the rate of removal from I to S due to recovery. This model in its simplest form assumes random mixing, where the parameter b represents a composite of the average contact rate and the disease-specific transmissibility given a contact event.
In order to introduce human behavior, we substitute a dynamic b t for b, which is a function of both b 0 , the probability that disease transmission takes place on contact, and a dynamic social rate of contact c t whose optimal value, c * t , is determined at each time period t as in economic epidemiological models [30], namely where c * t represents the optimal contact rate, defined as the number of contacts per unit time that maximize utility for the individual. Here, c * t is represented as a function of the number of infected in the population according to the perceived risks and benefits of social contacts, which we model as a utility function, assuming there is a constant utility independent of contact, a utility lost associated with infection, and a utility derived from the choice of number of daily contacts with a penalty for deviating from the choice of contacts which would yield the most utility.
This utility function is thus assumed to take the form Here U represents utility for an individual at time t given a particular number of contacts per unit time c t , α 0 is a constant that represents maximum potential utility achieved at a target contact rateĉ. The second term, α 1 (c t −ĉ) 2 , is a scaled squared deviation to represent the penalty for deviating fromĉ. The third term, is the cost of infection, α 2 , multiplied by the probability of infection over the course of the time unit. The time-delay ∆ represents the delay in information acquisition and the speed of response to that information. We note that (1 − I N b 0 ) ct can be approximated by We thus assume I N is small, as in the beginning stages of an epidemic, and approximate U (c t ) in Eqn. 2.5 using Eqn. 2.6. Eqn. 2.5 assumes a strictly negative relationship between number of infecteds and contact.
We assume an individual will balance the cost of infection, the probability of infection, and the cost of deviating from the target contact rateĉ to select an optimal contact rate c * t , namely the number of contacts, which takes into account the risk of infection and the penalty for deviating from the target contact rate. We assume all agents select c t regardless of infection status, such as in scenarios where individuals are not aware of their own status. This optimal contact rate can be calculated by finding the maximum of U with respect to c t from Eqn. 2.5 with substitution from Eqn. 2.6, namely which vanishes at the optimal contact rate, c * t , given by Therefore, total utility will decrease as I t increases and c * t also decreases. Utility is maximized at each time step, rather than over the course of lifetime expectations. In addition, Eqn. 2.9 assumes a strictly negative relationship between number of infecteds at time t − ∆ and c * . While behavior at high degrees of prevalence or expectations of same have been shown to be non-linear and fatalistic [31,32], in this model, prevalence (i.e., I N ) is assumed small consistent with Eqn. 2.6.
We introduce the new parameter α = α 2 2α 1 b 0 , so that We can now rewrite the recursion from Eqn. 2.2, using Eqn. 2.4 and replacing c t with c * t as defined by Eqn. 2.10, as When ∆ = 0 and there is no time delay, f (I) is a cubic function, given by (2.12) 3 Analytical Results

Equilibria
To determine the dynamic trajectories of this model without time delay, we solve for the fixed point(s)Î of I t (i.e. value or values of I such that f (I t ) = I t ) from Eqn. 2.12; that is, we solve Solving this quadratic yields two solutions. The first iŝ and the second isÎ (3.4) The discriminant of both of these equilibria is given by ((α −ĉ) 2 + 4 α b 0 N γ) and is therefore always positive. These non-zero equilibria are denoted byÎ 1 andÎ 2 , respectively.
These equlibria are considered legitimate when 0 ≤Î k ≤ N for k = [0, 1, 2] as the number of infected cannot be negative and cannot exceed the total population.Î 0 is therefore a legitimate equilibrium.Î 1 is always positive because its constituent components are all positive. The condition under whichÎ 1 ≤ N is ThusÎ 1 cannot be less than N, and it can equal N only if γ = 0, a degenerate case as it would mean there is no recovery and therefore any number of contacts would lead to the whole population becoming infected.
The condition under whichÎ 2 ≤ N is given by

Stability
Assessing global asymptotic stability in epidemic models is an important task of mathematical epidemiology [33,34]. In this model, when ∆ = 0, there are three equilibria for the recursion equation Eqn. 2.11: the disease-free equilibrium,Î 0 , which is stable whenÎ 2 > 0, Î 1 , which is unstable and greater than N, andÎ 2 , which is stable and positive whenÎ 0 is unstable. Recall the condition for the existence ofÎ 2 when γ > 0 isĉb 0 > γ N (Eqn. 3.7). These two scenarios for the three equilibria are illustrated in Figs. 1-2. Calculations for stability can be found in the appendix.
The basic reproduction number, the number of secondary cases produced by the introduction of one infected person into an entirely susceptible population, R 0 , is a fundamental tool of epidemiology because it delineates a threshold between an epidemic and eradication [35,36]. R t , the reproduction number at time t, is the eigen-value of this system or I t+1 It [37,38]. When R t is greater than 1, the epidemic will continue to grow. When it is less than 1, it will shrink. Reducing R t to less than 1 is therefore a central goal of epidemic intervention. The effect of endogenous behavior change on R t in this model is illustrated in Fig. 3. Figure 1: Return map for system with ∆ = 0 whereÎ 2 is a stable equilibrium and between 0 and N .Î 0 is unstable andÎ 1 > N . This is valid whenĉb 0 > γ N and γ > 0. R t > 1 if I 0 <Î 2 unless I t becomes greater thanÎ 2 at which point R t < 1. The system thereby regulates itself to approach theÎ 2 equilibrium when ∆ = 0.

Computational Results
Computer iterations were used to interrogate the parameter space of this model and the associated dynamics. Parameters include: • N , the total population • I 0 , the initial number of infected individuals • ∆, the time delay in behavioral response • b 0 , the transmissibility parameter • γ, the probability that an infected individual recovers to the susceptible class •ĉ, the optimal contact rate when I t = 0 • α 1 , the utility associated with choosing a contact rate equal toĉ • α 2 , the negative utility associated with becoming infected Under different parameter sets, the system shows five qualitatively distinct kinds of dynamical behaviors: 1) monotonic convergence to the disease-free equilibriumÎ 0 (number of infected goes to 0); 2) monotonic convergence to theÎ 2 equilibrium; 3) damped oscillation that converges toÎ 2 ; 4) perpetual oscillation; and 5) a collapse, or sudden extinction of the epidemic when the number of infecteds crosses 0. In Fig. 3 and Fig. 4 each of these dynamical classes is demonstrated with small differences in b 0 and ∆, respectively.
In order to identify the critical thresholds between the dynamical classes of results, default parameters were chosen to remain fixed while each of the other parameters was varied. The default parameters were N = 10000, I 0 = 1, b 0 =0.05, γ = 0.08,ĉ = 0.0015, α 1 = 0.02, α 2 = 0.3, ∆ = 3 We note that γ, b 0 ,ĉ, and ∆ are important to the qualitative landscape of the dynamics. While I 0 , α 1 , and α 2 may change the location of the non-zero fixed point and may influence the time it takes to reach equilibrium, they do not affect which fixed point the model is attracted to or the nature of the dynamics exhibited. The system is dynamically sensitive to γ, b 0 , and ∆ and includes critical transitions into distinct classes of dynamical behavior of increasing complexity as illustrated in Fig. 3.
Tables 1-4 list sets of runs for varying values of b 0 , γ, ∆, andĉ respectively. Critical thresholds for each of these parameters are calculated, beyond which a regime shift to a distinct qualitative dynamic classification occurs. These thresholds are not fixed for each of the parameters -they are mathematically interdependent. An example of values for each parameter that produces each system class is given.  Table 1): a) when b 0 = 0.005, there is monotonic convergence to the disease-free equilibrium atÎ 0 = 0 where R 0 < 1; b) when b 0 = 0.01, there is monotonic convergence to theÎ 2 equilibrium when R 0 > 1, but R t approaches 1 as the system goes to equilibrium; c) when b 0 = 0.033 there is damped oscillation that overshoots the equilibrium with progressively smaller amplitude until it reaches the equilibrium where R = 1; d) when b 0 = 0.04, there is perpetual oscillation without consistent period or amplitude with R t fluctuating frequently below and above 1 results; and e) when b 0 = 0.05, a collapse of the system occurs.  (Table 3), which are set such thatÎ 2 is stable at ∆ = 0: a) when ∆ = 0, there is monotonic convergence to theÎ 2 equilibrium; b) when ∆ = 1, there is damped oscillation that converges toÎ 2 (see Fig. 5 for a return map of this system); c) when ∆ = 2, there is perpetual oscillation with a complex pattern as shown in Fig. 6; and d) when ∆ = 3, there is collapse, or extinction of the epidemic.  Thê I 2 equilibrium is indicated. This regime is perpetual oscillation and shows a complex pattern of oscillations with a bounded range. More than one I t+1 value exists for each I t value. The importance of both values makes the system difficult to predict when ∆ is unknown.

Discussion
This simple model of endogenous behavior change contains a disease-free equilibrium and an endemic non-zero equilibrium, both of which are stable under different conditions. Several parameters, including b 0 , the transmissiblity parameter, Γ, the removal rate, ∆, the time delay, andĉ, the target contact rate, are each dynamically influential and can cause a regime shift through five distinct classes of dynamics found in the computational results, including monotonic convergence to the disease-free equilibrium, monotonic convergence to the non-zero equilibrium, damped oscillation to equilibrium, perpetual oscillation, and collapse.
Our model makes a number of simplifying assumptions. We assume, for example, that all individuals in the population respond in the same fashion. We assume that the individuals choose their contact rates according to an optimized utility function, which is homogeneous across all individuals in the population. Finally, we assume that the utility function is symmetric around the optimal number of contacts so that increasing or decreasing contacts above or below the target contact rate respectfully yield the same reduction in utility. These assumptions allowed us to 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 parameterized and predictive, and strategic models, which aim to be as simple as possible to highlight phenomenelogical generalities, this is a strategic model [39].
We note that the five distinct kinds of dynamical trajectories (Fig. 3) seen in these computational experiments come from a simple set of purely deterministic equations. This means that oscillations and even erratic, near-chaotic dynamics and collapse in an epidemic may not necessarily be due to seasonality, complex agent-based interactions, changing or stochastic parameter values, demographic change, host immunity, or socio-cultural idiosyncracies. This dynamical behavior in number of infecteds can result from mathematical properties of a simple deterministic system with homogeneous endogenous behavior change, similar to complex population dynamics of biological organisms [40]. The mathematical consistency with population dynamics suggests a parallel in ecology, that the indifference point for human behavior functions in a similar way to a carrying capacity in ecology, below which a population will tend to grow and above which a population will tend to shrink. For example, the Ricker Equation [41], commonly used in population dynamics to describe the growth of fish populations, exhibits similar complex dynamics and qualitative state thresholds. The existence of a non-zero, stable equilibrium in our model is consistent with economic epidemiology theory: if individuals are incentivized to change their behavior to protect themselves, they will, and they will cease to do this when they are not [10]. If we couple this with delayed information, the results can lead to limit-cycle dynamics, consis-tent with other negative feedback mechanisms with time delays [42,43]. This is because the system is reacting to conditions that were true in the past, but not necessarily true in the present.
Observed epidemic curves of many transient disease outbreaks typically inflect and go extinct, as opposed to this model that oscillates perpetually or converges to an endemic disease equilibrium. Our model is meant to demonstrate what effect personal incentives can have on infectious disease dynamics. Including institutional and public efforts that are incentivized to eradicate, rather than to optimize personal utility trade-offs, would alter the dynamics to look more like real-world epidemic curves. This may have a useful implication for policy. For example, beyond infectious diseases that remain endemic to society, outbreaks may also flare up once or multiple times, such as the double-peaked outbreaks of SARS in 3 countries in 2003 [44]. There may be many causes for such double-peaked outbreaks, one of which may be a lapse in behavior change after the epidemic begins to die down due to decreasing incentives, as represented in our simple theoretical model [17]. This is consistent with findings that voluntary vaccination programs suffer from decreasing incentives to participate as prevalence decreases [45,46].
One of the responsibilities of infectious disease modelers is to predict and project forward what epidemics will do in the future in order to better assist in the proper and strategic allocation of preventative resources. An insight from this model, however, is that prediction becomes very difficult if we take into account endogenous behavior change because the system is qualitatively sensitive to small differences in values of key parameters. These parameters are very hard to measure precisely; they change depending on the disease system and context and their measurement is generally subject to large errors. To maximize the ability to predict and minimize loss of life or morbidity, outbreak response should not only seek to minimize the reproduction number, but also the length of time taken to gather and distribute information.
In our model, the existence of complex dynamic regimes require a time delay. If behavior change arises from fear and fear is triggered by high local mortality, such delays seem plausible since death is a lagging epidemiological indicator. Lags mean that people can respond sluggishly to an unfolding epidemic crisis, but they also mean that people can abandon protective behaviors prematurely. Developing approaches to incentivize protective behavior throughout the duration of any lag introduced by the natural history of the infection (or otherwise) should be a priority in applied research. This paper represents a first step in understanding endogenous behavior change and time-lagged protective, and we anticipate further developments along these lines. Important future questions include: How and in what context do these phenomena occur empirically? What is the relative influence of insti-tutional interventions? How can we respond to epidemics in the future to best leverage the personal incentives to change behavior already inherent in this complex health landscape?

Stability
To determine stability in this model in the case ∆ = 0, we compute the derivative of the recursion function given by Eqn. 2.12 at the equilibria d dI (f (I)) For the disease-free equilibrium,Î 0 , then When the right side of Eqn. 6.2 is greater than 1, this equilibrium is unstable. The condition for instability then isĉ b 0 > γ N .
This condition also ensures thatÎ 2 is legitimate from Eqn. 3.7. This suggests that the disease-free equilibrium is unstable whenÎ 2 exists and stable when it does not.