## Abstract

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 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.

## 1 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 agent-based 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í Poincaré’s criteria for chaos in deterministic dynamical systems: a large number of solutions, dynamic sensitivity, and numerical unpredictability [29]. Measles models show characteristically 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,” the “middle-aged,” and the “old” (‘a la Acemoglu et al. [15]).

## 2 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 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 gained by achieving *ĉ* 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 2.4 with respect to c yields 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 The utility function (as in Eq 2.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. 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]),
where represents the within-group contact of *A*_{i}, optimized at time step *t* to maximize utility in *A*_{i} (as in Eq 2.4), and *c*_{ij} represents the contact between *A*_{i} and *A*_{j} for *i* ≠ *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 (2.6), the transition between susceptible and infected disease states for age class *A*_{i}, as in Eqs. 2.1-2.3, can be expressed as
We assume the following constraints:

The aversion to infection is greatest in the elderly and least in the young, i.e.,

The target contact rate is greatest in the young and least in the elderly, i.e.,

The recovery rate from the infected category to the susceptible is the same for all age classes, namely, for all

*i, j*,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,

## 3. Analytical Results

### 3.1 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 2.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. 2.6-2.13, with *k* = 2, become
with *S*_{1t} + *I*_{1t} = *N*_{1}, *S*_{2t} + *I*_{2t} = *N*_{2}, and *N* = *N*_{1} + *N*_{2}.

For simplicity, we assume the population sizes of all groups are equal

Then, substituting from Eq 3.1 and *α*_{1} from Eq 3.2, 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
and

### 3.2 Stability

To analyze local stability, we calculate the Jacobian of the differential equations corresponding to Eqs 3.8-3.9, 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},

## 4 Computational Results

### 4.1 The 2-population model

We first examine the numerical iteration of the discrete time SIS recursions (Eqs 3.8-3.9) without time-delay, and set default values for all parameters, namely:

*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,3,5). 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 2,4).

### 4.2 The 3-population model

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 2.6-2.13, for *k* = 3, are
with
Equilibria for the system defined by Eqs. 4.2-4.4 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, *N*_{3} = 5000, *I*_{0} = 1, *γ*_{1} = 0.1, *γ*_{2} = 0.1, *γ*_{3} = 0.1, *ĉ*_{1} = 0.02, *ĉ*_{2} = 0.015, *ĉ*_{3} = 0.01, *c*_{12} = 0.005, *c*_{13} = 0.003, *c*_{23} = 0.007, *α*_{1} = 1, *α*_{21} = 20, *α*_{22} = 30, *α*_{23} = 40, *b*_{0} = 0.01.

By increasing the transmissibility *b*_{0}, the model goes from simple convergence to a 2-point cycle to chaos into a 2-pt cycle to collapse (Fig 6). By uniformly increasing the between-group contact rates *c*_{12}, *c*_{13}, and *c*_{23}, the model exhibits convergence, a 2-pt cycle, a 4-pt cycle, chaos, a 5-py cycle, and a 6-pt cycle. The amplitude of oscillations and the variance of chaotic dynamics are greatest for the elderly population and least for the youth population.

## 5 Discussion

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 5). 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. 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 within-group 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. Further, we constrained the model with 4 mathematical assumptions (Eqs 2.10-2.13), some of which may not always hold in real-world scenarios. The relaxation of the above assumptions may yield different model outcomes. However, while the thresholds between dynamical regimes may shift, the regimes themselves are likely robust. We note these assumptions were made in order to construct the simplest possible mathematical model that includes adaptive behavior and age-structured risk perceptions in an epidemic.

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. In Germany, youths were found to drive COVID-19 transmission dynamics in the first three pandemic waves [12]. 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 sub-populations with differentiated characteristics or reactions to epidemic dynamics. It is worth noting that the justification for the bifurcated structure need not be restricted to age-based 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 het-erogeneous 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 5).

In our previous work [27], we compared the complex dynamics associated with endogenous behavior change during epidemics to similar theoretical behavior in ecological systems. Often, when a single population is modeled in ecology (e.g. in a fishery [42]), 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.

## Data Availability

All data produced in the present work are contained in the manuscript.