Age-differentiated incentives for adaptive behavior during epidemics produce oscillatory and chaotic dynamics

Heterogeneity in contact patterns, mortality rates, and transmissibility among and between different age classes can have significant effects on epidemic outcomes. Adaptive behavior in response to the spread of an infectious pathogen may give rise to complex epidemiological dynamics. Here we model an infectious disease in which adaptive behavior incentives, and mortality rates, can vary between two and three age classes. The model indicates that age-dependent variability in infection aversion can produce more complex epidemic dynamics at lower levels of pathogen transmissibility and that those at less risk of infection can still drive complexity in the dynamics of those at higher risk of infection. Policymakers should consider the interdependence of such heterogeneous groups when making decisions.


Introduction
One of the principal failings of the attempts to model and predict future trends and dynamics of infectious disease epidemics has been the lack of incorporation of human behavior into these models [1].The social drivers of infectious disease dynamics have been relatively neglected in the academic literature compared with the vast resources and attention paid to the biomedical and, to a lesser extent, the ecological drivers of disease [2,3].COVID-19, however, has shone a spotlight on this discrepancy as socio-cultural factors have played an important role in the pandemic-from the political polarization of risk perception in the United States [4] to the social, cultural, and demographic factors associated with vaccine hesitancy [5].Models should incorporate relevant social phenomena, as well as their interactions, as the study of each phenomenon in isolation may not be informative or useful.For example, adaptive age-specific preventative behaviors motivated by differentiated risk of COVID-19 mortality should be included.The resulting contact patterns and associated infection dynamics are expected to interact in their influence on epidemic trajectories.
The mortality rate of COVID-19 is highly age-specific [6] with a log-linear increase of infection fatality ratio by age among individuals over 30 [7].This contributes to skewed risk perceptions across the population's age structure and economic activity associated with COVID-19 risk [8,9].Heterogeneity in contact patterns, mortality rates, and transmissibility among and between different age classes is known to have significant effects on epidemic outcomes [10,11].In Germany, for example, younger adults and teenagers were likely the main drivers of COVID-19 transmission dynamics during the first three pandemic waves [12].Goldstein et al. [13] found that vaccinating the elderly first against COVID-19 would save the most lives, although the model neglected key features of transmission dynamics [14].Further, Acemoglu and colleagues showed that optimal lockdown policies during COVID-19 are age-specific, with strict lockdown policies on the oldest group and reduction of interactions between age classes being most effective [15].Their model, however, assumed exogenous targeted policies and did not incorporate adaptive behavior, a cornerstone of COVID-19 dynamics.
As risk of transmission of a dangerous infection increases during an epidemic, individuals and governments tend to react by mitigating that risk with adaptive behavior.Adaptive behavior has enjoyed growing attention in the disease modeling literature [16,17], with techniques using game theory [18], fear-infection parallel contagions [19,20], and network-and agentbased approaches [21].Early work by Capasso et al. [22] experimented with the introduction of a negative feedback mechanism in the traditional susceptible-infected-recovered (SIR) model [23] and represented reduced contacts as a function of the number of infected.Philipson formalized the economic theory of adaptive behavior showing that rational agents following dynamic incentives may lead to oscillations around an indifference point, or equilibrium [24].Adaptive behavior can cause oscillatory dynamics because a system with a negative feedback can continuously overshoot adjustments [25,26].
Adaptive behavior may also lead to complex dynamics that are characteristic of deterministic, chaotic systems, especially with delays in adaptive response [27].Empirical evidence from the COVID-19 pandemic suggests that the behavior of the epidemic has been chaotic in a majority of countries [28].An investigation of the second derivative of infections over time during COVID-19 found that the pandemic qualitatively met Henrı ´Poincare ´'s criteria for chaos in deterministic dynamical systems: a large number of solutions, dynamic sensitivity, and numerical unpredictability [29].Measles models show chaotic dynamics that may also characterize the observed dynamics [30,31].In the SEIR framework with a periodically varying contact rate that represented seasonal changes, sustained oscillations [32] and period-doubling bifurcations [33] were found, and one case with a particularly high degree of contact led to chaotic dynamics [34].While age structure heterogeneity in susceptibility and seasonal variance in contact rates due to school attendance characterize measles modeling [35,36], age structure heterogeneity in dynamic, adaptive contact rates have not been modeled for COVID- 19.
Here, we model an infectious disease with adaptive behavior incentives of the form used in a previous model [27] and include different mortality rates for three age classes: the "young" (ages 20-49), the "middle-aged" (50-64), and the "old" (65+) (à la Acemoglu et al. [15]).

Model specifications
We begin with a susceptible-infected-susceptible (SIS) compartmental disease model [37], which includes an adaptive contact behavior that maximizes a utility function, as in Arthur et al. [27].With susceptible and infectious individuals, denoted by S t and I t , respectively, at time t, and N t the population size, we have: where b 0 represents transmissibility given contact, γ represents the removal rate, and c * t represents a contact rate at time t chosen to maximize utility U(c) as a function of contacts: Here, c represents contacts per unit time, ĉ represents the optimal contact rate when the disease risk is non-existent, α 1 represents the utility lost for deviating from ĉ contacts and α 2 represents the utility lost if infected.Here, Δ represents delayed information, such that an individual may base their perception of infection risk on prevalence during past time periods, rather than the current one.Maximizing Eq 4 with respect to c yields c * t , the contact rate chosen at each time step to maximize utility.
To incorporate age structure, we stratify the population into k discrete age classes, represented by A 1 , . .., A k .The number in sub-population A i is given by N i .We assume the size of each age class does not change: N it = N i and let P k i¼1 N i ¼ N. The utility function (as in Eq 4) for each age class i is given by where α 1i represents the utility loss of reduced contacts for A i , ĉi represents the target contact rate for A i , and α 2i represents the utility loss (i.e., aversion) to infection based on a delayed perception of population-level prevalence for A i .Here it is assumed that each age class perceives its risk according to the disease prevalence of the whole population (i.e., I tÀ D ¼ P k i¼1 I i;tÀ D ), rather than just of their own group.
Interactions between and within age classes at time t are defined in terms of a dynamic contact matrix M t (censu Ram & Schaposnik, 2021 [38]), ; ð6Þ where c * it represents the within-group contact of A i , optimized at time step t to maximize utility in A i (as in Eq 4), and c ij represents the contact between A i and A j for i 6 ¼ j.For simplicity, institutional behavior change is assumed to only affect the within-group contact (e.g., via school and workplace closures and nursing home quarantines), and thus, between-group contact rates are assumed fixed and not adaptive to changing infection risks.It is assumed that c ij = c ji .
Using the contact matrix (Eq 6), the transition between susceptible and infected disease states for age class A i , as in Eqs 1-3, can be expressed as We assume the following constraints: 1.The aversion to infection is greatest in the elderly and least in the young, i.e., 2. The target contact rate is greatest in the young and least in the elderly, i.e., 3. The recovery rate from the infected category to the susceptible is the same for all age classes, namely, for all i, j, 4. The number of infecteds in any age class A i may never be greater than the population size of that class or less than zero, i.e., for all i and t, Analytical results

Equilibria
To understand the dynamic behavior of the number of infecteds in each age class, we first examine conditions for the existence of equilibria.If I t b 0 is small relative to N, then on linearizing with respect to I in Eq 5, the optimal value of c i at time t for A i is found to be Define the parameter α i as We begin with a 2-age-class model, where A 1 represents the youth and A 2 represents the middle-aged and elderly, in which case, Eqs 6-13, with k = 2, become For simplicity, we assume the population sizes of all groups are equal (i.e., . Then, substituting c * 1t from Eq 14 and α 1 from Eq 15, with Δ = 0, we obtain By symmetry, A fixed point (i.e., equilibrium) exists for A i when I i(t+1) = I it for i = 1, 2. Thus, equilibria are the roots of the two simultaneous polynomial equations

Stability
To analyze local stability, we calculate the Jacobian of the differential equations corresponding to Eqs 21 and 22, where I 1(t+1) = f 1 (I 1t , I 2t ) and I 2(t+1) = f 2 (I 2t , I 2t ) and evaluate the Jacobian at the equilibrium.Differentiating each function with respect to each variable I 1 and I 2 , ð29Þ

Computational results
The The aversion to infection for the elderly and the youth given by α 22 and α 21 , respectively, may produce more complex dynamics at higher values, but the relationship between the two values is also important (Fig 3).Simpler dynamics (e.g., equilibrium and 2-point cycles) can be found closer to the line y = x (i.e., α 22 = α 21 ), while highly differentiated aversions to infection are more likely to produce chaos or oscillations of higher periodicity.; ð30Þ

The 3-population model
with Equilibria for the system defined by Eqs 31-33 solve I i(t+1) = I it for i = [1, 2, 3] and are the roots of three simultaneous polynomial equations.
We set default values for parameters and initial conditions, such that N 1 = 5000, N 2 = 5000,

Discussion
Our model has generalized previous work on adaptive behavior during epidemics to include age structure with heterogeneous aversion to infection and target contact rates.Results from the 2-and 3-population adaptive behavior models indicate that, for certain parameter values, stable equilibria can be reached for each sub-population i, including 0 for all age groups and an equilibrium between 0 and N i .In the 2-population case, the equilibrium for each population can be derived numerically.With increasing values of b 0 and c ij (the transmissibility and between-group contact rate, respectively), the system may converge to a non-zero equilibrium, oscillate perpetually with 2-, 4-, 6-, 8-, and 10-period intervals, become chaotic, and collapse.Complexity of dynamics can be found at lower levels of transmissibility with greater differentiation between aversions to infection (Fig 3).For both the 2-population and 3-population models, the younger population has a higher non-zero equilibrium size than the older population, and the older population has greater amplitude of oscillations and variance of chaotic dynamics than the younger population.
Our model is built on a number of simplifying assumptions.First, the model is of the SIS type and does not consider an Exposed class (E) or a Recovered class (R).The inclusion of E or R (in, e.g., SEIRS models) may decrease the number of susceptible individuals and thus slow some of the oscillatory or chaotic dynamics found here.The inclusion of R for SIR models depletes susceptibles and prevents sustained dynamic patterns.Risk tolerance and reactions to shifting prevalence are stratified by age class, but assumed homogeneous within each class.In reality, individuals would have heterogeneous risk tolerance according to their age, political affiliation, and other characteristics.We modeled adaptive behavior with dynamic withingroup contact rates, but assumed between-group contact rates were fixed.While this was justified by COVID-19 public health policies that controlled within-group contact more than between-group contact, between-group contact is also responsive to prevalence.Expanding this model to finer age categories than the ones used in this paper should include empiricallyderived contact rates and differentiated optimization for between-group contacts.Further, we constrained the model with four mathematical assumptions (Eqs 10-13), some of which may not always hold in real-world scenarios.The relaxation of the above assumptions may yield different model outcomes.We note the model was not trained on real-world data to produce parameters as empirical data on adaptive behavior and incentives is lacking.While the dynamical regimes found here exist in theory, we do not know where in practice real-world disease systems may be.However, COVID-19 exhibited oscillatory patterns and chaotic dynamics [28], suggesting that this theoretical work may be directly relevant to diseases with strong adaptive behavior components.
Without contact with other groups, the Elderly age class exhibits simpler dynamics (e.g., monotonic convergence) (Fig 6).While the older population has a higher risk associated with infection, the younger population's lower aversion and higher baseline contact rates affect the epidemic dynamics in the older population.Thus, transmission in the young may not only lead to transmission in the elderly, but also increase the variance of the elderly dynamics and destabilize them.Indeed, even at lower levels of infection aversion for the elderly, the more differentiated the aversions to infection between the two age classes, the more complex the dynamics (Fig 3).Children are known to be the primary drivers of Influenza transmission, although severe morbidity and mortality are mostly seen in older age groups [39].The importance of the youth in the dynamics of our model provides theoretical support for COVID-19 lockdown policies that reduce between-group interaction (i.e., c ij ), as found by Acemoglu et al. [15].Transmission within the household is a key point of risk for the elderly and middle-aged, and high levels of transmission in the youth are likely to significantly affect disease outcomes in other age categories.
Our model may be usefully compared to other adaptive behavior models that look at subpopulations with differentiated characteristics or reactions to epidemic dynamics.It is worth noting that the justification for the bifurcated structure need not be restricted to agebased differentiation, but can also include political party affiliation, income levels, or other demographic, social, behavioral, physical, or geographic differences.Some studies use structured 2-population models with varying homophily and outgroup aversion or varying awareness separation and mixing separation [40,41].Results from these models indicate that heterogeneous populations, even when simply structured compartmentally as two populations, can produce greater complexity in epidemic dynamics, including large second waves and interconnected dual epidemics.Our results broadly agree with this theme: a 2-population model with varying infection aversion can produce more complexity at lower levels of transmissibility, and the difference between aversions can also drive complexity (e.g., Fig 3).
In our previous work [27], we compared the complex dynamics associated with endogenous behavior change during epidemics to similar theoretical behavior in ecological systems.For example, the logistic map is an extensively-studied simple model of population growth with limited resources that exhibits similar dynamics to our adaptive behavior models: convergence, oscillations, chaos, and collapse [42].The driver of these dynamics in the logistic map is r, the growth rate.In our model, two of the principal drivers are b 0 , the transmissibility, and c ij , the between-group contact rate, as the epidemic's reproduction number R 0 is directly proportional to transmissiblity and contact rate.Often, when a single population is modeled in ecology (e.g., in a fishery [43]), the carrying capacity operates as an attractor, above which the population is attracted downwards and below which the population is attracted upwards.With endogenous behavior, the equilibrium, or indifference point, is also an attractor, below which the population is motivated to relax protective behaviors and above which the population is motivated to adopt protective behaviors.It may be productive to extend this parallel further in the case of multiple populations.For example, the predator-prey model considers two interdependent populations.When the prey population is high, the predator population grows, but when the prey are low in number, the predators decrease.While our epidemiological model does not include such direct competition or predation, the behavior is somewhat similar: when the young population has high prevalence of disease, the prevalence in the older population increases until the indifference point is crossed.When the young population has low prevalence, the prevalence in the older population follows suit.If a comparison between ecology and epidemiology is appropriate, it follows that careful study of the literature in theoretical ecology may provide insights relevant to epidemiology, a field with as yet comparatively little exploration of such complex system dynamics.
We recommend further theoretical work in both adaptive behavior modeling and complexity in epidemiology.A systems perspective may better represent the inherent complexities and heterogeneities of real-world epidemics.Social drivers of disease have been relatively neglected in the literature [2], but played an important role in COVID-19 outcomes.For example, the divergence of risk assessment by age class that characterizes our model was a social phenomenon, though biologically motivated.Other complex phenomena important to COVID-19 include simultaneous asynchronous epidemics, political bifurcation of attitudes and practices, and the co-evolution of the human immune system and the virus.As public health policies depend on our ability to forecast different scenarios under a high degree of uncertainty and complexity, such modeling will play an important role in improving policy and health outcomes in future epidemics.
2-population modelWe first examine the numerical iteration of the discrete time SIS recursions (Eqs 21 and 22) without time-delay, and set default values for all parameters.No empirical research has analyzed comparative utility lost or gained as a function of contacts.Thus, the parameter set has been chosen to demonstrate possible dynamic regimes, rather than to parameterize COVID-19 or other infectious diseases.These default parameters are:N 1 = 5000, N 2 = 5000, I 0 = 1, γ 1 = 0.1, γ 2 = 0.1, ĉ1 ¼ 0:02, ĉ2 ¼ 0:01, c 12 = 0.005, α 1 = 1, α 21 = 20, α 22 = 40, b 0 = 0.009.By increasing the transmissibility parameter b 0 from the default value, we see a progression of dynamical regimes across critical thresholds from simple convergence to cyclic behavior in 2, 4, and 6-point cycles, chaos, and collapse (Figs 1, 2 and 3).By increasing the contact rate between the 2 populations, c 12 , there is a progression from simple convergence to a 2-pt cycle, 6-pt cycle, chaos, back to a 6-pt cycle, and finally an asynchronous 8-pt and 2-pt cycle (Figs 4 and 5).Bifurcation diagrams in Figs 2 and 5 indicate sustained dynamic patterns for a continuous parameter change in either b 0 (Fig 2) or c 12 (Fig 5).When c 12 = 0 and c * t is only based on prevalence in the same age category, we see what would happen in a homogeneous population without age structure with the same parameters (Fig 6).Without this interdependent age-structure, Youth dynamics are even more complex, while Elderly dynamics have less complexity.

For our 3 -
age-class model, A 1 represents the youth, A 2 represents the middle-aged, and A 3 represents the elderly.Then the analogous equations to Eqs 6-13, for k = 3, are

Fig 2 .Fig 3 .
Fig 2. Bifurcation diagram of equilibria, oscillatory dynamics, and chaotic behavior as a function of transmissibility b 0 .Points in the figure represent sustained long-term results in the model; thus any b 0 value with, for example, two corresponding yaxis points represents a sustained oscillation between two values for infected individuals.https://doi.org/10.1371/journal.pcbi.1011217.g002

Fig 5 .
Fig 5. Bifurcation diagram of equilibria, oscillatory dynamics, and chaotic behavior as a function of inter-group contact rate c 12 .Points in the figure represent sustained long-term results in the model; thus any c 12 value with, for example, two corresponding y-axis points represents a sustained oscillation between two values for infected individuals.https://doi.org/10.1371/journal.pcbi.1011217.g005