Age classes stratified by risk and adaptive behavior during epidemics

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₀ 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, α₁ represents the utility gained by achieving ĉ contacts and α₂ 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₁, …, 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:

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,

3. Analytical Results

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₁ = 5000, N₂ = 5000, I₀ = 1, γ₁ = 0.1, γ₂ = 0.1,, ĉ₁ = 0.02, ĉ₂ = 0.01, c₁₂ = 0.005, α₁ = 1, α₂₁ = 20, α₂₂ = 40, b₀ = 0.009.

By increasing the transmissibility parameter b₀ 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₁₂, 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).

• Download figure
• Open in new tab

Figure 1:

Disease dynamics for 2-population age-structured model using default parameters and varying transmissibility b₀. A) b₀ = 0.009, Convergence to stable equilibria; B) b₀ = 0.013 =, 2-point cycle; C) b₀ = 0.02, 4-point cycle; D) b₀ = 0.021, chaos; E) b₀ = 0.0235, 6-pt cycle; F) b₀ = 0.025, Collapse.

• Download figure
• Open in new tab

Figure 2:

Bifurcation diagram of equilibria, oscillatory dynamics, and chaotic behavior as a function of transmissibility b₀.

• Download figure
• Open in new tab

Figure 3:

Disease dynamics for 2-population age-structured model using default parameters and varying between-group contact rate c₁₂. A) c₁₂ = 0.005, Convergence to stable equilibria; B) c₁₂ = 0.025, 2-point cycle; C)c₁₂ = 0.04, 6-point cycle; D) c₁₂ = 0.05, chaos; E) c₁₂ = 0.08, 6-pt cycle; F) c₁₂ = 0.15, 8-pt cycle for Old and 2-pt cycle for Young.

• Download figure
• Open in new tab

Figure 4:

Bifurcation diagram of equilibria, oscillatory dynamics, and chaotic behavior as a function of inter-group contact rate c₁₂.

• Download figure
• Open in new tab

Figure 5:

Heat maps of dynamic regimes as a function of the relationship between aversion to infection in the youth (α₂₁) and in the elderly (α₂₂) with varying transmissibility (b₀): A) b₀ = 0.009; B) b₀ = 0.013; C) b₀ = 0.02; D) b₀ = 0.021; E) b₀ = 0.0235; F) b₀ = 0.025.

4.2 The 3-population model

For our 3-age-class model, A₁ represents the youth, A₂ represents the middle-aged, and A₃ 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₁ = 5000, N₂ = 5000, N₃ = 5000, I₀ = 1, γ₁ = 0.1, γ₂ = 0.1, γ₃ = 0.1, ĉ₁ = 0.02, ĉ₂ = 0.015, ĉ₃ = 0.01, c₁₂ = 0.005, c₁₃ = 0.003, c₂₃ = 0.007, α₁ = 1, α₂₁ = 20, α₂₂ = 30, α₂₃ = 40, b₀ = 0.01.

By increasing the transmissibility b₀, 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₁₂, c₁₃, and c₂₃, 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.

• Download figure
• Open in new tab

Figure 6:

Disease dynamics for 3-population age-structured model using default parameters and varying transmissibility b₀. A) b₀ = 0.01, Convergence to stable equilibria; B) b₀ = 0.013, 2-point cycle; C)b₀ = 0.022, Chaos into a 2-pt cycle cycle; D) b₀ = 0.05, Collapse.

• Download figure
• Open in new tab

Figure 7:

Disease dynamics for 3-population age-structured model using default parameters and varying between-group contact rates c₁₂, c₁₃, and c₂₃.

c₁₂ = 0.005, c₁₃ = 0.003, c₂₃ = 0.007, Convergence to stable equilibria; B) c₁₂ = 0.008, c₁₃ = 0.01, c₂₃ = 0.01, 2-point cycle; C) c₁₂ = 0.022, c₁₃ = 0.02, c₂₃ = 0.024, 4-point cycle; D) c₁₂ = 0.032, c₁₃ = 0.03, c₂₃ = 0.034, Chaos; E) c₁₂ = 0.05, c₁₃ = 0.048, c₂₃ = 0.052, 5-pt cycle; F) c₁₂ = 0.06, c₁₃ = 0.058, c₂₃ = 0.062, 6-pt cycle.

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