The Complex Interplay Between Risk Tolerance and the Spread of Infectious Diseases

Abstract

Risk-driven behavior provides a feedback mechanism through which individuals both shape and are collectively affected by an epidemic. We introduce a general and flexible compartmental model to study the effect of heterogeneity in the population with regards to risk tolerance. The interplay between behavior and epidemiology leads to a rich set of possible epidemic dynamics. Depending on the behavioral composition of the population, we find that increasing heterogeneity in risk tolerance can either increase or decrease the epidemic size. We find that multiple waves of infection can arise due to the interplay between transmission and behavior, even without the replenishment of susceptibles. We find that increasing protective mechanisms such as the effectiveness of interventions, the number of risk-averse people in the population, and the duration of intervention usage reduces the epidemic overshoot. When the protection is pushed past a critical threshold, the epidemic dynamics enter an underdamped regime where the epidemic size exactly equals the herd immunity threshold. Lastly, we can find regimes where epidemic size does not monotonically decrease with a population that becomes increasingly risk-averse.

Introduction

Recent outbreaks such as the COVID-19 pandemic, the 2014 Ebola outbreak, the 2009 influenza A (H1N1) pandemic, and the 2002 SARS epidemic brought to light many of the challenges of mounting an effective and unified epidemic response in a country as large and as diverse as the United States. Particularly during the COVID-19 pandemic, people were split in opinion on questions such as the origin of the virus [maxmen_covid_2021], whether they would social distance or wear a mask [betsch_social_2020, fischer_mask_2021, yang_sociocultural_2022], or whether the country should even have a pandemic response at all [morens_concept_2022]. As time progressed, the situation became more dire and the death toll accumulated. People then had a new battery of questions to address, such as whether or not they would adhere to mandatory lock-downs [wong_paradox_2020, brzezinski_science_2021, kleitman_comply_2021] or whether they felt comfortable using the new mRNA vaccines [machingaidze_understanding_2021, fedele_covid-19_2021]. Compounding the issue were the multiple streams of information and potential misinformation spread through social media and other channels [gallotti_assessing_2020, loomba_measuring_2021, offer-westort_battling_2024, towers_mass_2015]. People’s stances on the questions and issues were diverse, arising from the milieu of differences in culture, geography, scientific education, sources of information, political leanings, and individual identity [murray_chapter_2016, kramer_infection_2021, lu_collectivism_2021].

Taken altogether, these differences within the population reflect a spectrum of people’s risk tolerances to a circulating infectious disease. For any given intervention, such as social distancing, wearing a mask, or taking a vaccine, each person in the population falls somewhere on a spectrum of willingness to adopt the intervention. Given a threat level of an infectious disease in the population, some people will readily wear masks, whereas other people will refuse to.

In this study, we aim to analyze the impact of heterogeneity in risk tolerance and the resulting behavioral response on the dynamics of epidemics. We seek to add to a burgeoning literature on the impact of human behavior in epidemic response [bansal_when_2007, funk_modelling_2010, fenichel_adaptive_2011, edmunds_evaluating_1999, weitz_awareness-driven_2020, wagner_economic_2020, tyson_timing_2020, espinoza_asymptomatic_2021, wagner_modelling_2022, espinoza_heterogeneous_2022, qiu_understanding_2022, traulsen_individual_2023, saad-roy_dynamics_2023, smith_covid-19_2023], which the recent pandemic highlighted as an area for further exploration in preparation for the next large scale global health crisis [morse_prediction_2012, osterholm_preparing_2020, bergstrom_human_2024]. To study the impact of heterogeneity in risk tolerance on epidemic dynamics, we introduce a simple and flexible modeling framework based on ordinary differential equations that can be used for different interventions and an arbitrary partitioning of the population with regard to risk tolerance and behavioral responses. We will examine and discuss potential interesting outcomes that can arise from coupling individual-level preferences and population-level epidemiology.

Results

Model of Adaptive Intervention Usage under Heterogeneous Risk Tolerance

Here we assume people’s risk aversion manifests as the rate at which they adopt individual interventions in response to an infectious disease outbreak. The intuition underlying this paradigm is that more risk-averse individuals are more sensitive to becoming sick and thus will adopt interventions at a faster rate than more risk-tolerant people. We consider the following SPIR compartmental model of a population with $n$ differing levels of risk tolerance (1-4). This model features four types of classes: unprotected susceptible ($S$), protected susceptible ($P$), infectious ($I$), and recovered with permanent immunity ($R$). Since there are $n$ differing levels of risk tolerance, we subdivide the susceptible population into $n$ discrete groups indexed by $i$, where $i\in\{1,2,...,n\}$. Each tolerance level is characterized by an intervention adoption rate parameter ($\lambda_{i}$) and an intervention relaxation rate parameter ($\delta_{i}$). Transitions of susceptibles between their unprotected class ($S_{i}$) and their corresponding protected class ($P_{i}$) are governed by the corresponding parameters of the same index ($\lambda_{i},\delta_{i}$). Overall, the system is governed by $3+2n$ parameters: a transmission rate parameter ($\beta$), a recovery rate parameter ($\gamma$), an intervention effectiveness parameter ($\epsilon$), and an intervention adoption rate ($\lambda_{i}$) and intervention relaxation rate ($\delta_{i}$) for each tolerance level.

$$\frac{dS_{i}}{dt} =-\beta S_{i}I-\lambda_{i}S_{i}I+\delta_{i}P_{i}$$ (1)

$$\frac{dP_{i}}{dt} =-(1-\epsilon)\beta P_{i}I-\delta_{i}P_{i}+\lambda_{i}S_{i}I$$ (2)

$$\frac{dI}{dt} =-\gamma I+\sum_{i=1}^{n}(\beta S_{i}I+(1-\epsilon)\beta P_{i}I)$$ (3)

$$\frac{dR}{dt} =\gamma I$$ (4)

The transition from the unprotected susceptible state to the protected susceptible state represents individuals implementing an intervention that confers them protection against disease transmission from an infected individual. The rate at which intervention adoption occurs may be driven by individuals considering information such as the epidemic incidence rate (e.g. cases per day), the total number of infected individuals in the population (e.g. total number of active cases), and mortality rate (e.g. deaths per day) [weitz_awareness-driven_2020]. Here we assume that individuals have knowledge about the total number of infected individuals ($I$) and respond accordingly. Parameterizing each person’s individual risk tolerance by $\lambda_{i}$, we assume each individual person adopts an intervention at a rate $\lambda_{i}I$. Then, if there are $S_{i}$ number of people that behave exactly the same (i.e. have the same level of risk-aversion), then at the population scale there is a collective adoption rate of $\lambda_{i}S_{i}I$. The same reasoning holds for each of the $n$ tolerance levels. We also consider a model where the adoption rate is driven by individuals reacting to the incidence rate (Supplemental Materials); while this produces a more complex mathematical model, the results are qualitatively similar.

The effectiveness of the intervention being used is captured by the parameter $\epsilon$, which linearly scales down the transmission rate between infected and protected susceptibles. In the limit of $\epsilon=1$, the intervention is perfectly effective and protected individuals cannot become infected. In the limit of $\epsilon=0$ then the intervention is completely ineffective, which reduces the model to an SIR model without interventions. For simplicity, we assume each epidemic features only a single type of intervention (whether that be masking, social distancing, vaccines, etc.) and that the effectiveness of an intervention is identical across the population. In reality, multiple interventions may be available concurrently, which would drive additional variation in behavior due to differences in risk sensitivity across the population.

This model allows for protected individuals to relax their usage of interventions, becoming unprotected in the process. Here, individuals in the protected class can relax back to the unprotected class through two means, either through a rate that is dependent on the quantity of infections present or through a rate that is independent of the number of infections. The infection-dependent rate is implicitly captured through the $\lambda SI$ term, which can be thought of a net rate that can be further decomposed into adoption and relaxation rate terms (i.e. $\lambda SI=\lambda_{adoption}SI-\lambda_{relaxation}SI$). Here we have assumed the adoption term to always be greater than the relaxation term, otherwise individuals would never adopt an intervention. The infection-independent rate is governed by the intervention relaxation rate parameter ($\delta_{i}$) for each tolerance level. In the limit of $\delta_{i}=0$, an intervention is irreversible, which would represent an intervention such as vaccines with permanent immunity. When $\delta_{i}$ is non-zero, individuals are using interventions such as masking or social distancing. The infection-independent rate is motivated by factors such as psychological fatigue of social distancing [franzen_fatigue_2021, jorgensen_pandemic_2022] and physical discomfort with wearing masks [cheok_appropriate_2021]. In general, we will consider the regime where the relaxation rate $\delta$ is of comparable scale or smaller than the transmission scale (i.e. $\delta\leq\beta$). This reflects intuition that people are likely to continue to protect themselves with interventions even beyond an initial outbreak [barak_experience_2022].

For simplicity, we consider the model for the case when $n=1$ and $n=2$. A schematic for these two cases is shown in Figure 1. However, the framework is general and can be extended to any discrete number of groups.

Figure 1: Flow diagram for an SIR model with adaptive interventions for either (a) a population with homogeneous risk tolerance or (b) a heterogeneous population with two different levels of risk tolerance. Susceptible individuals can access a more protected susceptible state through usage of interventions. The transition rate to the protected state depends on the incidence level. The protected state offers a $1-\epsilon$ reduction in transmission rate over the normal susceptible state.

For convention, when there are two susceptible classes, we assume the first susceptible class ($S_{1}$) has a lower risk tolerance for becoming infected (i.e. more risk-averse). As a result, these individuals more readily adopt the intervention (i.e. $\lambda_{1}>\lambda_{2}$), making individuals in this class transition more rapidly to the protected susceptible state ($P_{1}$). The second susceptible class ($S_{2}$) is more risk tolerant (i.e. more risk-taking), and thus is less eager to use the intervention, making individuals in this class transition more slowly to their protected susceptible state ($P_{2}$).

Adaptive Adoption of Interventions Can Produce Damped Oscillations

The coupling of intervention usage to the incidence rate and the resulting adaptive changes enables the epidemic dynamics to display a much richer set of behavior over the simple SIR model. From Figure 2, we see this particular set of conditions can deterministically produce multiple waves of infection, even when vital dynamics (i.e. birth and death processes) are not considered.

Figure 2: Time series for population with homogeneous risk tolerance and adaptive intervention usage and the corresponding phase space trajectory indicate the presence of damped oscillations.

Evidence for cycling of individuals between using interventions and not using interventions during the COVID pandemic can be seen in longitudinal usage [mathieu_coronavirus_2020, hale_global_2021, rader_mask-wearing_2021, salomon_us_2021] data. The possibility for these oscillations highlight the intimate connection between individual human behavior and intervention usage in shaping the dynamics of epidemics, while also be affected by the collective decision of everyone in the population. The coupling of behavior and epidemiology here provides a feedback mechanism where an increasing incidence rate prompts more individuals to adopt an intervention, which lowers the overall incidence rate; however, as the epidemic wanes and factors such as fatigue or discomfort set in, people begin dropping their usage of interventions, which may eventually lead to another wave of outbreaks if enough people become unprotected while infected individuals still remain, and then the cycle can be repeated.

Protective Mechanisms Saturate in Underdamped Regime that Eliminates Epidemic Overshoot

One might have the intuition that having more people that will more readily adopt an intervention (i.e. mask, social distance, or vaccinate) or increasing the effectiveness of the intervention in reducing transmission will further decrease the size of the epidemic. While we find this intuition to be mostly correct, we unexpectedly find that the protection conferred by either of these mechanisms can saturate once a critical parameter threshold has been passed.

In Figure 3$,left$, we see that increasing the effectiveness of the intervention or increasing the fraction of the population that are risk-averse monotonically decreases the epidemic size. However, in the dark blue region (which we will refer to as the underdamped regime) where both protection mechanisms are at their highest, we see no further reduction in the epidemic size. This regime corresponds to an epidemic where the epidemic size exactly equals the herd immunity threshold.

The orbits of the dynamics from different areas of this parameter space are shown in Figure 3$,right$. We see that even though the final epidemic size is the same throughout the underdamped region, the trajectories to reach the same final epidemic size can look qualitatively different.

Figures S4-S5 show a larger sampling of trajectories if one fixes either the fraction of the population with low risk tolerance or the intervention effectiveness respectively. It becomes clear that at the border of the underdamped region, we can see a clear change in the qualitative behavior of the trajectories as the threshold is crossed. Under some assumptions, one can prove that that the epidemic overshoot is eliminated in the underdamped regime (Supplemental Information).

Figure 3: Left. Epidemic size as a function of varying the fraction of the population that are low-risk tolerance (i.e. those with higher $\lambda$). Right. Corresponding orbits in the I versus S+P plane for the sampled points in parameter space.

However, we should make the point that this is not evidence that highly effective interventions are a waste or that the overall population should tolerate risky behavior. As this is a model with a large parameter space, we cannot visualize all of it. If we could, we would find many parameter regimes where the protection mechanisms never reach a critical threshold, which implies the conventional intuition of increasing intervention effectiveness and having more risk-averse people always being beneficial applies.

The Threshold to the Underdamped Regime is Reduced when Intervention Usage is Prolonged

The transition to the underdamped regime is more easily accessed when the usage of interventions is prolonged (or equivalently when the rate at which protected individuals relax back into the regular susceptibility classes decreases). Consider the following scenario which is identical to the previous setup, except now the intervention reversion rate ($\delta P_{i}$) has been reduced by an order of magnitude (Figure 4). This corresponds to a scenario where people continue to use the intervention (i.e. such as wearing masks or social distancing) on a timescale significantly longer than the transmission timescale ($\delta>>\beta$).

(a) $\epsilon=0.5,\delta_{1}=\delta_{2}=1$

(b) $\epsilon=0.5,\delta_{1}=\delta_{2}=0.1$

(c) $\epsilon=0.9,\delta_{1}=\delta_{2}=1$

(d) $\epsilon=0.9,\delta_{1}=\delta_{2}=0.1$

Figure 4: Final epidemic size versus fraction of population that are risk-averse ($S_{1}$). Simulations in the left column have a higher $\delta$ than simulations in the right column. The other simulation parameters and initial conditions are $\lambda_{1}=100,\lambda_{2}=1,I=10^{-7},P_{1}=P_{2}=0,R=0$.

This suggests that increasing the timescale at which individuals continue to use interventions decreases the number of risk-averse individuals needed to achieve the same epidemic size. This is reflected in the horizontal shift of the transition region to the left when comparing figures in the left column and the right column (Figure 4).

Increasing Heterogeneity in Risk Tolerance can Either Increase or Decrease the Epidemic Size

The literature generally suggests that increasing heterogeneity in the population through increasing the variation in their contact patterns, age, etc. results in a reduction in the epidemic size [miller_epidemic_2007, britton_mathematical_2020, gomes_individual_2022, allard_role_2023].

We find in this model of heterogeneous behavior that it is possible to switch from a regime where increasing the heterogeneity in risk tolerance results in a decrease in epidemic size to a regime where increasing the heterogeneity in risk tolerance results in an increase in epidemic size. This result also does not have to necessarily be due to a dramatic shift in parameters. From Figure 5, we see this shift can arise from solely varying the fraction of the population with low risk tolerance by a small amount.

Figure 5: Epidemic size under differing levels of heterogeneity in the adoption rate for interventions. The mean adoption rate of the two groups (i.e. geometric average of $\lambda_{1},\lambda_{2}$) is compared to the difference in the two adoption rates as parameterized by a homogeneity index (see Methods for definition). $Left$ is when the fraction of the population with low risk tolerance ($x_{1}$) is 0.2, $center$ is when $x_{1}=0.35$, $right$ is when $x_{1}=0.5$.

The intuition underlying this result is that when the average adoption rate ($\lambda_{average}$), which is expressed as a (geometric or arithmetic) weighted mean of the adoption rates of the two groups, is fixed at a constant level, then the epidemic size can be suppressed through either varying the fraction of the population in each group or through varying each group’s adoption rate. When risk-averse people make up a smaller fraction of the population than risk-taking people, then it would be more beneficial in reducing epidemic size to have the adoption rates of the two groups be more similar (i.e. more homogeneous) as that would imply risk-taking people (which are then the majority of the population) would have a similar adoption rate to risk-averse people. As an increasingly larger fraction of the population is composed of risk-averse people, then it becomes increasingly beneficial in reducing epidemic size to have the adoption rates of the two groups be more different (i.e. more heterogeneous) as the deleterious effects of highly risk-taking people (which are then the minority of the population) can be mitigated by the large presence of risk-averse people.

Epidemic Size Does Not Necessarily Decrease with an Increasing Number of Risk-Averse People

General intuition would suggest that as one increases the number of risk-averse people in the population that the overall epidemic size would go down. However, the introduction of the adaptive behavior mechanism allows for regimes where this is no longer strictly the case. Thus, it is no longer a guarantee that decreasing the population’s overall risk tolerance will always improve epidemic outcomes.

Consider Figure 6(a), where we find a small region after the transition to the underdamped regime where there is an increase in the epidemic size when increasing the fraction of the population that are risk-averse.

(a)

(b)

Figure 6: (a) Final epidemic size as a function of the proportion of the population that are risk averse ($S_{1}$). Model parameters and initial conditions for the simulation are $\epsilon=1,\lambda_{1}=10,\lambda_{2}=0.5,\delta_{1}=0.1,\delta_{2}=0.1,I(0)=10^{-7},x_{S_{2}}(0)=1-x_{S_{1}}-I(0),P_{1}(0)=P_{2}(0)=0,R(0)=0$. (b) Same as in (a) except now the effectiveness of the intervention ($\epsilon$) is allowed to vary.

If we also vary the effectiveness of the intervention as an additional axis (Figure 6b), we observe that there is a small trench in the threshold region surrounding the plateau area. This double descent suggests that the landscape can potentially be quite complicated when risk tolerance in the population is partitioned into even more groups.

Discussion and Conclusions

In this paper, we have proposed a simple model to model heterogeneity in risk tolerance levels in the population. We find that including a behavioral mechanism for adopting interventions that adapts with the level of infections greatly expands the variety in epidemic dynamics and outcomes that can occur.

The general picture from the findings suggest that epidemic dynamics under adaptive intervention adoption fall into either an underdamped regime or an overdamped regime. The underdamped regime has a special property in which the epidemic size equals the herd immunity threshold exactly, which means no epidemic overshoot occurs. The system can be driven into this regime when protection mechanisms (such as numbers of risk-averse people, intervention effectiveness, and duration of intervention usage) are increased to a sufficiently high level. This regime is also marked by damped oscillations in the phase space of infecteds and susceptibles. In direct contrast, the overdamped regime closely resembles the dynamics of a simple $SIR$ model without behavior, in which there are no oscillations and a non-zero overshoot, which makes the epidemic size greater than the herd immunity threshold.

We have looked for some evidence in the historical data on outbreaks for these damped oscillations due to cycles in adoption and relaxation of interventions. While such data is in very limited supply, previous analysis suggested that relaxation of social distancing measures may have led to multiple waves of infection in the Spanish flu of the early 20th century [hatchett_public_2007, caley_quantifying_2007]. Dating back to the time of the bubonic plague, there is data from an outbreak in 1636 in the parish of St. Martin in the Fields, which showed how relaxation of quarantining and isolation measures lead to a smaller secondary wave of infections [newman_shutt_2012]. We also see some evidence from time-series data from early in the COVID pandemic on masking policy [hale_global_2021]. Taking mask policy as a proxy for the general level of mask usage [hale_global_2021], we can observe the policy level in the United States as it relates to the overall incidence in COVID infections during that same period [mathieu_coronavirus_2020], which we take as a rough proxy for the number of infected people. We see that a relaxation in the mask policy from its strictest level coincided with the rise of the Omicron variant soon afterward (Figure 7). While there are many confounding variables at play here including viral evolution that make it difficult to disentangle the contribution of individual factors, the timing suggests that behavior plays a key role in shaping the epidemic landscape as well.

Figure 7: Time series data for masking policy and confirmed new cases in the USA for the first two years of the COVID pandemic. Data extracted from [mathieu_coronavirus_2020, hale_global_2021].

The results here are reminiscent of feedback control systems commonly studied in control theory. Here the set point is the herd immunity threshold, which is determined by the basic reproduction number ($R_{0}$). The ability for the population to reach this set point for epidemic size without additional overshoot depends on the effectiveness of the feedback mechanism from coupling intervention usage to the number of infected people. In the model presented, the adoption of interventions is a continuous process, in which the different groups are constantly reacting to the level of infections without requiring any notion of time or thresholds. In contrast, existing research on mitigation have considered more active control where activation and intervention timing play a key role [bussell_applying_2019, lauro_optimal_2021, morris_optimal_2021]. Future work may explore how to synergistically utilize both active and continuous mechanisms for control.

The inclusion of heterogeneity in risk tolerance and adaptive adoption of interventions leads to several unexpected conclusions. We find that increasing heterogeneity in risk tolerance levels in the population can lead to either an increase or decrease in the epidemic size. The direction of the trend depends nonlinearly on the composition of the population in terms of the ratio of risk-averse to risk-taking individuals and their respective intervention adoption rates. This adds to a small literature that demonstrates how heterogeneity can actually lead to a larger epidemic [volz_effects_2011, espinoza_heterogeneous_2022].

Interestingly, these results on heterogeneity also can be used to address the question of whether distributed or centralized control of mitigation results in a smaller outbreak. Control of factors such as mobility may be more practically achieved in a more centralized and unified fashion [bonaccorsi_economic_2020, espinoza_mobility_2020], whereas a distributed approach may be more appropriate when considering factors such as speed and individual agency. In a centralized scenario, a single entity controls the dynamics. That situation has an exact correspondence to the homogeneous population considered here, where all individuals respond in unison. Distributed control allows for more localized control, such as individuals or small groups deciding if they want to mask or social distance. This corresponds to the heterogeneous scenarios considered here, where there are multiple groups each with differing risk tolerance level. The results suggest that in some scenarios a single, coordinated response would be better for mitigation, whereas in other parameter regimes, a more decentralized strategy would be more optimal.

We also find that increasing overall protection mechanisms does not always result in a monotonic decrease in epidemic size. In scenarios when the adoption rate begins to approach the transmission rate, near the critically damped boundary a nonmonotonicity can arise. This suggests that when intervention usage and effectiveness are tenuous, the dynamics become more complex and predicting what epidemic outcomes will result becomes significantly more difficult. Understanding how these nonlinear effects combine with other biological and behavioral heterogeneities will be important to explore in future work.

Methods

Defining the $\lambda$ Homogeneity Index

The $\lambda$ homogeneity index is defined as follows. We will assume the initial condition that at the beginning of the dynamics, the total population is composed of the fraction of the population in the low risk tolerance group $x_{1}$, the high risk tolerance group $x_{2}$, or the infected class. We will assume the fraction of the population initially infected is sufficiently small so that size of the two susceptible compartments is given by $x_{1}$ and $1-x_{1}$ respectively.

We can define the average adoption rate as being either a geometric or arithmetic mean of the two adoption rates. The choice one makes is arbitrary, so we present prescriptions for both routes. In both cases, we will map the level of homogeneity to the unit interval.

Geometric Average

Let the geometric average of the two adoption rates be given by $\lambda_{\text{Geometric Average}}$.

$$\lambda_{\text{Geometric Average}}=\lambda_{1}^{x_{1}}\lambda_{2}^{1-x_{1}}$$ (5)

Define $\lambda_{2}$ as a fraction between 0 and 1 of the average $\lambda$.

$$\lambda_{2}=c\lambda_{\text{Geometric Average}},c\in(0,1]$$ (6)

These two equations combine to define $\lambda_{1}$.

$$\lambda_{1}=\frac{\lambda_{\text{Geometric Average}}}{c^{\frac{x_{1}}{1-x_{1}}}},c\in(0,1]$$ (7)

An index value of 1 indicates $\lambda_{1}=\lambda_{2}$, while decreasing the index value towards 0 increases the difference between $\lambda_{1}$ and $\lambda_{2}$.

Arithmetic Average

Let the arithmetic average of the two adoption rates be given by $\lambda_{\text{Arithmetic Average}}$.

$$\lambda_{\text{Arithmetic Average}}=x_{1}\lambda_{1}+(1-x_{1})\lambda_{2}$$ (8)

Define $\lambda_{2}$ as a fraction between 0 and 1 of the average $\lambda$.

$$\lambda_{2}=c\lambda_{\text{Arithmetic Average}},c\in(0,1]$$ (9)

These two equations combine to define $\lambda_{1}$.

$$\lambda_{1}=\frac{\lambda_{\text{Arithmetic Average}}(1-x_{1}c)}{1-x_{1}},c\in(0,1]$$ (10)

Again, an index value of 1 indicates $\lambda_{1}=\lambda_{2}$, while decreasing the index value towards 0 increases the difference between $\lambda_{1}$ and $\lambda_{2}$.

Numerical Solutions and Code

Author Contributions

designed research, performed research, and wrote and reviewed the manuscript.

Additional information

Supplemental Material: The Complex Interplay Between Risk Tolerance and the Spread of Infectious Diseases

Analysis of the Model When Individuals Respond to the Incidence Rate

The rate at which intervention adoption occurs may be driven by individuals considering information such as the epidemic incidence rate (e.g. cases per day), the total number of infected individuals in the population (e.g. total number of active cases), and mortality rate (e.g. deaths per day) [weitz_awareness-driven_2020]. Here we will consider the first case, where individuals adopt interventions based on the incidence rate for infections. Recall that the incidence rate is given by $\sum_{i=1}^{n}(\beta S_{i}I+(1-\epsilon)\beta P_{i}I)$. Parameterizing each person’s individual risk tolerance by $\lambda_{i}$, let us assume each individual adopts an intervention at rate $\lambda_{i}\sum_{i=1}^{n}(\beta S_{i}I+(1-\epsilon)\beta P_{i}I)$. Then, if there are $S_{i}$ number of people that behave exactly the same (i.e. have the same level of risk-aversion), then at the population scale there is a collective adoption rate of $\lambda_{i}S_{i}\sum_{i=1}^{n}(\beta S_{i}I+(1-\epsilon)\beta P_{i}I)$. The same reasoning holds for each of the $n$ tolerance levels. The corresponding equations for this model are given by (11-14).

$$\frac{dS_{i}}{dt} =-\beta S_{i}I-\lambda_{i}S_{i}\sum_{i=1}^{n}(\beta S_{i}I+(1-\epsilon)\beta P_{i}I)+\delta_{i}P_{i}$$ (11)

$$\frac{dP_{i}}{dt} =-(1-\epsilon)\beta P_{i}I-\delta_{i}P_{i}+\lambda_{i}S_{i}\sum_{i=1}^{n}(\beta S_{i}I+(1-\epsilon)\beta P_{i}I)$$ (12)

$$\frac{dI}{dt} =-\gamma I+\sum_{i=1}^{n}(\beta S_{i}I+(1-\epsilon)\beta P_{i}I)$$ (13)

$$\frac{dR}{dt} =\gamma I$$ (14)

When comparing the final epidemic size and epidemic trajectories between this model and the model in the main text where individuals adopt interventions at a rate that is based on the total number of infected people (Figures S4-S5), the results are indistinguishable (Figures S1-S2).

Figure S1: Left. Epidemic size as a function of varying the fraction of the population that are low-risk tolerance (i.e. those with higher $\lambda$) when individuals react to incidence rate. Right. Corresponding orbits in the I versus S+P plane for the sampled points in parameter space when the fraction of the population that are low-risk tolerance has been fixed.

Figure S2: Left. Epidemic size as a function of varying the fraction of the population that are low-risk tolerance (i.e. those with higher $\lambda$) when individuals react to incidence rate. Right. Corresponding orbits in the I versus S+P plane for the sampled points in parameter space when the intervention effectiveness has been fixed.

Thus, given the equivalence in results, we have gone with the mathematically simpler and cleaner model based on total number of infected people in the main text.

The Herd Immunity Threshold is Set by a Complex Interplay Between Transmission ($R_{0}$), Behavior, and Intervention Effectiveness ($\epsilon$).

While we could make some analytical calculations for the epidemic size in the homogeneous model, the heterogeneous two group case requires a numerical approach to find the plateau region. In general, it is set by a highly nonlinear interaction between the transmission ($R_{0}$), behavior as determined by the fraction of the population that are risk-averse and risk-taking, and the effectiveness of the intervention ($\epsilon$).

Consider the following progression of figures (Figure S3), where in each subsequent figure, the effectiveness of the intervention in blocking transmission ($\epsilon$) is increasing.

(a) $\epsilon=0$

(b) $\epsilon=0.3$

(c) $\epsilon=0.5$

(d) $\epsilon=0.7$

(e) $\epsilon=0.9$

(f) $\epsilon=1$

Figure S3: Final epidemic size versus fraction of population that are risk-averse ($S_{1}$) with a progressive increase in intervention effectiveness ($\epsilon$). The other simulation parameters and initial conditions are $\lambda_{1}=100,\delta_{1}=1,\lambda_{2}=1,\delta_{2}=1,I(0)=10^{-7},P_{1}(0)=P_{2}(0)=0,R(0)=0$.

Comparing Orbits Inside and Outside of the Plateau Region of Herd Immunity

The following figures sample more orbits for the Figure considered in the main text.

Figure S4: Left. Epidemic size as a function of varying the fraction of the population that are low-risk tolerance (i.e. those with higher $\lambda$). Right. Corresponding orbits in the I versus S+P plane for the sampled points in parameter space when the fraction of the population that are low-risk tolerance has been fixed.

Figure S5: Left. Epidemic size as a function of varying the fraction of the population that are low-risk tolerance (i.e. those with higher $\lambda$). Right. Corresponding orbits in the I versus S+P plane for the sampled points in parameter space when the intervention effectiveness has been fixed.

Proving Underdamped Regime Eliminates Epidemic Overshoot

In general, the nonlinear feedback between the protected classes and infected individuals make it difficult to make analytical calculations in the full model. Under some simplifications however, we can make some progress. In this section, we derive the epidemic size at which the protection saturates in a restricted case of the homogeneous model ($n=1$).

Let us consider the homogeneous model in the limit of an intervention with perfect effectiveness (i.e. $\epsilon=1$). We will consider the case where the recovery rate from infection and relaxation rate for interventions are comparable (i.e. $\gamma=\delta$). Since the equation for recovered individuals can be ignored since the population is closed ($S+P+I+R=1$), this reduces the dynamics to the following system:

$$\frac{dS}{dt} =-\beta SI-\lambda SI+\gamma P$$ (15)

$$\frac{dP}{dt} =\lambda SI-\gamma P$$ (16)

$$\frac{dI}{dt} =\beta SI-\gamma I$$ (17)

$$S(0) =f,P(0)=0,I(0)=\alpha,R(0)=1-f-\alpha$$ (18)

The initial conditions sets the number of individuals initially susceptible to be $f$, the number of initially infected is assumed to be small $I(0)<<1$, and the remainder of the population is already immune to infection (i.e. recovered). To ensure that the epidemic initially grows in size, we assume that $S(0)>\frac{1}{R_{0}}$. This follows from comparing the incidence term ($\beta SI$) to the recovery term ($\gamma I$) in 17.

To find the asymptotic behavior for $S$, we will attempt to eliminate $P$ from (15). We start first by seeking an equation that relates the $P$ and $I$ compartments. Consider the following ansatz that considers the difference between the two compartments:

$$x=\frac{\beta}{\lambda}P-I$$ (19)

Differentiating this equation with respect to time and using (16)-(17) yields:

$$\frac{dx}{dt}=\gamma\left(I-\frac{\beta}{\lambda}P\right)$$ (20)

Using (19), this simplifies to $\frac{dx}{dt}=-\gamma x$, which has a critical point at $x=0$. At $x=0$, we obtain that $P=\frac{\lambda}{\beta}I$, indicating a regime where the behavior of $P$ and $I$ scale linearly with each other. Combining this with (15) yields:

$$\frac{dS}{dt} =\left(-(\beta+\lambda)S+\frac{\gamma\lambda}{\beta}\right)I$$ (21)

Now that we have an equation that is linear in $I$, we are in a good position to find a final size relationship for the number of susceptibles. To start we take the ratio of (17) and (21).

$$\frac{dI}{dS} =\frac{\beta SI-\gamma I}{\left(-(\beta+\lambda)S+\frac{\gamma\lambda}{\beta}\right)I}$$ (22)

Using the partial fractions $\frac{-\beta}{\beta+\lambda}$ and $\frac{\frac{\gamma\beta}{\beta+\lambda}}{((\beta+\lambda)S-\frac{\gamma\lambda}{\beta})}$, we get upon indefinite integration of (22) that
$k=\frac{\beta}{\beta+\lambda}\left(\frac{\gamma}{\beta+\lambda}\text{ln}\left((\beta+\lambda)S-\frac{\gamma\lambda}{\beta}\right)-S\right)-I$, where $k$ is a constant that holds throughout the trajectory of the dynamics. Thus, considering the values of $S$ and $I$ at the beginning of the epidemic ($t=0$) and the end of the epidemic ($t=\infty$) and using the conditions that $I(\infty)=0$ and $I(0)\approx 0$ yields the following transcendental equation for the final epidemic size.

$$\frac{\beta}{\beta+\lambda}\left(\frac{\gamma}{\beta+\lambda}\text{ln}\left((\beta+\lambda)S(0)-\frac{\gamma\lambda}{\beta}\right)-S(0)\right) =\frac{\beta}{\beta+\lambda}\left(\frac{\gamma}{\beta+\lambda}\text{ln}\left((\beta+\lambda)S(\infty)-\frac{\gamma\lambda}{\beta}\right)-S(\infty)\right)$$ (23)

Since $S(0)>\frac{1}{R_{0}}$, then the argument of the logarithm on the left hand side must be positive, and subsequently the left hand side evaluates to a real number. Due to the equality, the right hand side must also evaluate to a real number, implying the argument of the logarithm on the right hand side must also be positive. Positivity implies the following inequality for the lower bound for the final number of susceptibles:

$$S(\infty)>\frac{\gamma\lambda}{\beta(\beta+\lambda)}$$ (24)

An upper bound can be given by simply noting that in the long time limit, the recovery term ($\gamma I$) must be at least as large as the incidence term ($\beta S(\infty)I$) in (17), otherwise the epidemic would still be growing. This implies $S(\infty)\leq\frac{1}{R_{0}}$.

To summarize, when the interventions are perfectly effective, the rate of relaxation from the protected class is equal to the rate of recovery from infection, and the number in the protected class scales linearly with the number of infected, then the final fraction of susceptibles is bounded as follows:

$$\frac{\lambda}{R_{0}(\beta+\lambda)}<S(\infty)\leq\frac{1}{R_{0}}$$ (25)

We see in the parameter limit of when the adoption rate of interventions is very fast compared to the transmission rate (i.e. $\lambda>>\beta$), that the lower bound reduces to $\frac{1}{R_{0}}$. Since both bounds now coincide, then $S_{\infty}$ must equal that value. Interestingly this corresponds to the herd immunity threshold of the standard SIR model. As the overshoot is the excess number of cases beyond the herd immunity threshold, we see that in this parameter limit there is no overshoot.

This analysis for the homogeneous case also carries over to the heterogeneous case for two groups when the adoption rate between the two groups is significantly different (i.e. $\lambda_{1}>>\lambda_{2}$). This results in a separation of time scales in which the faster adopters quickly transition to the protected state and can essentially be treated as immune over the course of the remaining epidemic over the slow adopters. This amounts to effectively reducing the dynamics to the homogeneous model considered here where $f$ and $1-f$ fractions of the population in the susceptible and recovered respectively correspond to the fraction of the population in the slow adopter ($\lambda_{2}$) and fast adopter groups ($\lambda_{1}$).

Heterogeneity in Risk Tolerance through Arithmetic Averaging

There are two ways for calculating the difference (heterogeneity) in adoption rates, either through geometric or arithmetic averaging. Both can be justified, and we presented the geometric formulation in the main text. We find that the arithmetic formulation gives qualitative similar results under a suitable parameter shift (Figure S6).

Figure S6: Epidemic size under differing levels of heterogeneity in the adoption rate for interventions. The mean adoption rate of the two groups (i.e. arithmetic average of $\lambda_{1},\lambda_{2}$) is compared to the difference in the two adoption rates as parameterized by a homogeneity index (see Methods for definition). $Left$ is when the fraction of the population with low risk tolerance ($x_{1}$) is 0.2, $center$ is when $x_{1}=0.5$, $right$ is when $x_{1}=0.8$.

Code to Generate Figures

Code executed in MATLAB R2023a.

⬇

1%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2

3function totalRecovered = plotSIRDynamics(susceptibleFrac, R0, transmissionReduction, PlotOption, gamma, lambda1, delta1, lambda2, delta2, modelType)

4

5beta = R0*gamma;

6

7% Initial conditions

8initialInfected = 0.00000001;

9initialSusceptible1 = susceptibleFrac;%(1-initialInfected)/2;

10initialSusceptible2 = 1-initialInfected-initialSusceptible1;

11initialMasked1 = 0;

12initialMasked2 = 0;

13initialRecovered = 0;

14

15% Time vector

16tspan = [0 10000];

17

18% SIR ODE system

19

20%% Reversible w/ total level response

21sirODE_level = @(t, y) [

22 -beta*y(1)*y(5) - lambda1*y(1)*y(5) + delta1*y(2);

23 -beta*(1-transmissionReduction)*y(2)*y(5) - delta1*y(2) + lambda1*y(1)*y(5);

24 -beta*y(3)*y(5) - lambda2*y(3)*y(5) + delta2*y(4);

25 -beta*(1-transmissionReduction)*y(4)*y(5) - delta2*y(4) + lambda2*y(3)*y(5);

26 beta*(y(1)+y(3))*y(5) + beta*(1-transmissionReduction)*(y(2)+y(4))*y(5) - gamma*y(5);

27 gamma * y(5)];

28

29%% Reversible w/ incidence rate response

30

31sirODE_incidence = @(t, y) [

32 -beta*y(1)*y(5) - lambda1*y(1)*(beta*(y(1)+y(3))*y(5) + beta*(1-transmissionReduction)*(y(2)+y(4))*y(5)) + delta1*y(2);

33 -beta*(1-transmissionReduction)*y(2)*y(5) - delta1*y(2) + lambda1*y(1)*(beta*(y(1)+y(3))*y(5) + beta*(1-transmissionReduction)*(y(2)+y(4))*y(5));

34 -beta*y(3)*y(5) - lambda2*y(3)*(beta*(y(1)+y(3))*y(5) + beta*(1-transmissionReduction)*(y(2)+y(4))*y(5)) + delta2*y(4);

35 -beta*(1-transmissionReduction)*y(4)*y(5) - delta2*y(4) + lambda2*y(3)*(beta*(y(1)+y(3))*y(5) + beta*(1-transmissionReduction)*(y(2)+y(4))*y(5));

36 beta*(y(1)+y(3))*y(5) + beta*(1-transmissionReduction)*(y(2)+y(4))*y(5) - gamma*y(5);

37 gamma * y(5)];

38

39% Solve ODE

40if modelType == 1

41 [t, y] = ode89(sirODE_level, tspan, [initialSusceptible1; initialMasked1; initialSusceptible2; initialMasked2; initialInfected; initialRecovered]);

42else

43 [t, y] = ode89(sirODE_incidence, tspan, [initialSusceptible1; initialMasked1; initialSusceptible2; initialMasked2; initialInfected; initialRecovered]);

44end

45

46if PlotOption == 1

47% Plot results

48 figure;

49 plot(t, y(:, 1), ’k-’, ’LineWidth’, 2, ’DisplayName’, ’Susceptible - Type 1’);

50 hold on;

51 plot(t, y(:, 3), ’b-’, ’LineWidth’, 2, ’DisplayName’, ’Susceptible - Type 2’);

52 plot(t, y(:, 2), ’g-’, ’LineWidth’, 2, ’DisplayName’, ’Protected - Type 1’);

53 plot(t, y(:, 4), ’m-’, ’LineWidth’, 2, ’DisplayName’, ’Protected - Type 2’);

54 plot(t, y(:, 5), ’r-’, ’LineWidth’, 2, ’DisplayName’, ’Infected’);

55 plot(t, y(:, 6), ’b-’, ’LineWidth’, 2, ’DisplayName’, ’Recovered’);

56 plot(t, y(:, 1)+y(:, 2)+y(:, 3)+y(:, 4), ’k.’, ’LineWidth’, 2, ’DisplayName’, ’Total Susceptible + Protected’);

57 xlabel(’Time’);

58 ylabel(’Proportion of Population’);

59

60 ylim([0 1])

61 title([’\beta=’,num2str(beta),’, \gamma=’,num2str(gamma),’, \epsilon=’,num2str(transmissionReduction), ’, \lambda_1=’,num2str(lambda1),’, \delta_1=’,num2str(delta1),’, \lambda_2=’,num2str(lambda2),’, \delta_2=’,num2str(delta2),’, x_{S_1}=’,num2str(susceptibleFrac)]);

62 legend(’Location’, ’eastoutside’);

63 grid on;

64 hold off;

65

66end

67

68if PlotOption == 2

69 hold on

70 plot(y(:, 1)+y(:,2)+y(:, 3)+y(:,4), y(:,5), ’b-’)

71 xlabel(’S+P’)

72 ylabel(’I’)

73 ylim([0 .15])

74end

75

76totalRecovered = y(end,6);

77end

78

79%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

80

81%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

82

83%%%%% SIR model parameters

84R0 = 3; % Basic reproduction number

85gamma = 1; % Recovery rate

86beta = R0*gamma; % Transmission rate of unmasked susceptibles

87transmissionReduction = 0.8; % Effectiveness of intervention (\epsilon)

88lambda1 = 10; % Rate parameter for unprotected susceptibles to adopt intervention

89delta1 = .01; % Rate parameter for protected individuals to remove intervention (mask)

90

91PlotOption = 1; % Plotting parameter. If 1, then plot graph, otherwise none.

92

93% Initial conditions

94initialInfected = 0.000001; % Seed fraction of population that are initially infected

95initialSusceptible1 = 1-initialInfected;

96initialMasked1 = 0;

97initialRecovered = 0;

98

99% Time vector

100tspan = [0 1500];

101

102% SIR ODE system

103%% Reversible w/ total level response

104sirODE_level = @(t, y) [

105 -beta*y(1)*y(3) - lambda1*y(1)*y(3) + delta1*y(2);

106 -beta*(1-transmissionReduction)*y(2)*y(3) - delta1*y(2) + lambda1*y(1)*y(3);

107 beta*y(1)*y(3) + beta*(1-transmissionReduction)*y(2)*y(3) - gamma*y(3);

108 gamma * y(3)];

109

110% Solve ODE

111[t, y] = ode89(sirODE_level, tspan, [initialSusceptible1; initialMasked1; initialInfected; initialRecovered]);

112

113if PlotOption == 1

114 % Plot results

115 figure;

116 plot(t, y(:, 1), ’k-’, ’LineWidth’, 2, ’DisplayName’, ’Susceptible - Type 1’);

117 hold on;

118 plot(t, y(:, 2), ’g-’, ’LineWidth’, 2, ’DisplayName’, ’Protected - Type 1’);

119 plot(t, y(:, 3), ’r-’, ’LineWidth’, 2, ’DisplayName’, ’Infected’);

120 plot(t, y(:, 4), ’b-’, ’LineWidth’, 2, ’DisplayName’, ’Recovered’);

121 plot(t, y(:, 1)+y(:, 2), ’c-’, ’LineWidth’, 2, ’DisplayName’, ’Total Susceptible + Protected’);

122 xlabel(’Time’);

123 ylabel(’Proportion of Population’);

124 ylim([0 1])

125 title([’\beta=’,num2str(beta),’, \gamma=’,num2str(gamma),’, \epsilon=’,num2str(transmissionReduction), ’, \lambda_1=’,num2str(lambda1),’, \delta_1=’,num2str(delta1)]);

126 legend(’Location’, ’eastoutside’);

127 grid on;

128 hold off;

129%

130% saveas(gcf,[’NPImodel_’,’beta’,num2str(beta),’_gamma’,num2str(gamma),’_epsilon’,num2str(transmissionReduction), ’_lambda1-’,num2str(lambda1),’_delta1-’,num2str(delta1),’.png’])

131

132 figure

133 hold on

134 plot3(y(:, 1),y(:,2), y(:,3), ’b-’)

135 xlabel(’S’)

136 ylabel(’P’)

137 zlabel(’I’)

138end

139

140totalRecovered = y(end,4);

141

142

143%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

144susceptFracVec = 0:0.02:1;

145R0vec = 2;

146colorVec = jet(length(R0vec));

147transmissionReduction = 0:0.02:1;

148

149% SIR model parameters

150gamma = 1; % Recovery rate

151lambda1 = 100; % Rate of susceptibles with lower tolerance to switch to masked class

152delta1 = .1; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

153lambda2 = 1; % Rate of susceptibles with higher tolerance to switch to masked class

154delta2 = .1; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

155

156recoveredFracVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

157xVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

158effectivenessVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

159

160for z = 1:length(R0vec)

161 for y = 1:length(transmissionReduction)

162 for x = 1:length(susceptFracVec)

163 recoveredFracVec(x,y,z) = plotSIRDynamics(susceptFracVec(x), R0vec(z),transmissionReduction(y), 0, gamma, lambda1, delta1, lambda2, delta2, 1);

164 xVec(x,y,z) = susceptFracVec(x);

165 effectivenessVec(x,y,z) = transmissionReduction(y);

166 end

167 end

168end

169

170figure

171

172subplot(1,2,1)

173[x_new, e_new] = meshgrid(susceptFracVec, transmissionReduction);

174contourf(x_new,e_new,recoveredFracVec(:,:,z)’, ’LevelStep’, 0.01)

175xlabel(’Fraction of Population with Low Risk Tolerance’, ’FontSize’, 14)

176ylabel(’Effectiveness of Intervention (\epsilon)’, ’FontSize’, 14)

177h = colorbar;

178title(h, ’Epidemic Size’, ’FontSize’, 14)

179title([’R_0=’, num2str(R0vec(z)), ’, \lambda_1=’, num2str(lambda1), ’, \delta_1=’, num2str(delta1), ’, \lambda_2=’, num2str(lambda2), ’, \delta_2=’, num2str(delta2)]);

180

181hold on

182points_x = [0.05, 0.2,.9, 0.5, .95, .75];

183points_y = [0.8, 0.1, .2,0.9, .95, .5];

184

185scatter(points_x, points_y, ’r’, ’filled’)

186

187labels = {’A’, ’B’, ’C’, ’D’, ’E’, ’F’};

188text(points_x, points_y, labels, ’Color’, ’red’, ’VerticalAlignment’, ’bottom’, ’HorizontalAlignment’, ’right’)

189

190subplot(1,2,2)

191

192hold on

193for a = 1:length(points_x)

194 plotSIRDynamics(points_x(a), R0vec(z),points_y(a), 2, gamma, lambda1, delta1, lambda2, delta2, 1);

195end

196

197points_x = [0.45,0.4, .36, .85, .95, 0.72];

198points_y = [0.025,0.13, 0.02,.01, .002, 0.02];

199

200scatter(points_x, points_y, ’Marker’, ’none’)

201labels = {’A’, ’B’, ’C’, ’D’, ’E’, ’F’};

202text(points_x, points_y, labels, ’Color’, ’red’, ’VerticalAlignment’, ’bottom’, ’HorizontalAlignment’, ’right’)

203

204

205%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure 4 and S3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

206% Repeat for different parameter conditions to generate each panel

207

208susceptFracVec = 0:0.01:1;

209R0vec = 1:0.5:10;

210colorVec = jet(length(R0vec));

211transmissionReduction = 1;

212

213% SIR model parameters

214gamma = 1; % Recovery rate

215lambda1 = 100; % Rate of susceptibles with lower tolerance to switch to masked class

216delta1 = .1; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

217lambda2 = 1; % Rate of susceptibles with higher tolerance to switch to masked class

218delta2 = .1; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

219

220

221recoveredFracVec = zeros(length(R0vec),length(susceptFracVec));

222figure

223hold on

224for y = 1:length(R0vec)

225 for x = 1:length(susceptFracVec)

226 recoveredFracVec(y, x) = plotSIRDynamics(susceptFracVec(x), R0vec(y),transmissionReduction, 0, gamma, lambda1, delta1, lambda2, delta2, 1);

227 end

228 hold on

229 plot(susceptFracVec, recoveredFracVec(y,:), ’Color’, colorVec(y,:), ’LineWidth’,2.5)

230end

231

232xlabel(’Fraction of Population that are Risk Averse’, ’FontSize’, 14)

233ylabel(’Epidemic Size’, ’FontSize’, 14)

234title([’\epsilon=’,num2str(transmissionReduction), ’, \lambda_1=’,num2str(lambda1),’, \delta_1=’,num2str(delta1),’, \lambda_2=’,num2str(lambda2),’, \delta_2=’,num2str(delta2)]);

235hLegend = legend(string(R0vec), ’Location’, ’eastoutside’, ’FontSize’, 12);

236title(hLegend, ’R_0’)

237

238%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure 5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

239 % Repeat for different parameter conditions to generate each panel

240

241figure

242susceptFracVec = 0.2;

243R0vec = 2;

244colorVec = jet(length(R0vec));

245transmissionReduction = .7;

246

247% SIR model parameters

248cVec = [0.01:0.01:1];

249lambdaAvg = [0.5, 1:1:100];

250gamma = 1; % Recovery rate

251delta1 = .5; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

252delta2 = .5; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

253

254recoveredFracVec = zeros(length(cVec),length(lambdaAvg));

255

256for x = 1:length(cVec)

257 for y = 1:length(lambdaAvg)

258 recoveredFracVec(x,y) = plotSIRDynamics(susceptFracVec, R0vec, transmissionReduction, 0, gamma, lambdaAvg(y)/(cVec(x)^(susceptFracVec/(1-susceptFracVec))), delta1, cVec(x)*lambdaAvg(y), delta2, 1);

259 end

260end

261

262[c_new, l_new] = meshgrid(cVec, lambdaAvg);

263contourf(c_new,l_new,recoveredFracVec’, ’LevelStep’, 0.01)

264xlabel(’\lambda Homogeneity Index’, ’FontSize’, 14)

265ylabel(’\lambda_{Geometric Average}’, ’FontSize’, 14)

266h = colorbar;

267title(h, ’Epidemic Size’, ’FontSize’, 14)

268title([’R_0=’, num2str(R0vec), ’, x_1=’, num2str(susceptFracVec), ’, \epsilon=’, num2str(transmissionReduction), ’, \delta_1=’, num2str(delta1),’, \delta_2=’, num2str(delta2)]);

269

270%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure 6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

271

272figure

273subplot(1,2,1)

274

275susceptFracVec = 0:0.01:1;

276R0vec = 1:0.5:7;

277colorVec = jet(length(R0vec));

278transmissionReduction = 1;

279

280% SIR model parameters

281gamma = 1; % Recovery rate

282lambda1 = 10; % Rate of susceptibles with lower tolerance to switch to masked class

283delta1 = .1; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

284lambda2 = .5; % Rate of susceptibles with higher tolerance to switch to masked class

285delta2 = .1; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

286

287

288recoveredFracVec = zeros(length(R0vec),length(susceptFracVec));

289figure

290hold on

291for y = 1:length(R0vec)

292 for x = 1:length(susceptFracVec)

293 recoveredFracVec(y, x) = plotSIRDynamics(susceptFracVec(x), R0vec(y),transmissionReduction, 0, gamma, lambda1, delta1, lambda2, delta2, 1);

294 end

295 hold on

296 plot(susceptFracVec, recoveredFracVec(y,:), ’Color’, colorVec(y,:), ’LineWidth’,2.5)

297end

298

299xlabel(’Fraction of Population that are Risk Averse’, ’FontSize’, 14)

300ylabel(’Epidemic Size’, ’FontSize’, 14)

301title([’\epsilon=’,num2str(transmissionReduction), ’, \lambda_1=’,num2str(lambda1),’, \delta_1=’,num2str(delta1),’, \lambda_2=’,num2str(lambda2),’, \delta_2=’,num2str(delta2)]);

302hLegend = legend(string(R0vec), ’Location’, ’eastoutside’, ’FontSize’, 12);

303title(hLegend, ’R_0’)

304

305susceptFracVec = 0:0.02:1;

306R0vec = 2.5; % Only pick one R0 to show the surface plot for. Use other plot to pick the value of interest

307colorVec = jet(length(R0vec));

308transmissionReduction = 0:0.02:1;

309

310recoveredFracVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

311xVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

312effectivenessVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

313

314for z = 1:length(R0vec)

315 for y = 1:length(transmissionReduction)

316 for x = 1:length(susceptFracVec)

317 recoveredFracVec(x,y,z) = plotSIRDynamics(susceptFracVec(x), R0vec(z),transmissionReduction(y), 0, gamma, lambda1, delta1, lambda2, delta2, 1);

318 xVec(x,y,z) = susceptFracVec(x);

319 effectivenessVec(x,y,z) = transmissionReduction(y);

320 end

321 end

322end

323

324[X,Y,Z] = meshgrid(susceptFracVec,transmissionReduction,R0vec);

325

326subplot(1,2,2)

327for z = 1:length(R0vec)

328 surf(X(:,:,z), Y(:,:,z), recoveredFracVec(:,:,z), recoveredFracVec(:,:,z), ’FaceAlpha’, 0.8, ’EdgeColor’, ’interp’);

329 hold on

330end

331

332xlabel(’Fraction of Population that are Risk-Averse’, ’FontSize’, 14)

333ylabel(’Effectiveness of Intervention (\epsilon)’, ’FontSize’, 14)

334zlabel(’Epidemic Size’, ’FontSize’, 14)

335title([’\lambda_1=’, num2str(lambda1), ’, \delta_1=’, num2str(delta1), ’, \lambda_2=’, num2str(lambda2), ’, \delta_2=’, num2str(delta2)]);

336

337

338%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure 7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

339

340maskPolicyDataVsCases = readtable(’maskPolicyDataVsCases.xlsx’);

341

342dates = datetime(maskPolicyDataVsCases.Date, ’InputFormat’, ’your_date_format’); % Replace ’your_date_format’ with the actual format

343policy = maskPolicyDataVsCases.MaskPolicyLevel;

344cases = maskPolicyDataVsCases.ConfirmedNewCases;

345

346figure;

347

348% Plot the first variable on the left y-axis

349yyaxis left;

350plot(dates, policy);

351ylabel(’Mask Policy Level’);

352grid off;

353ylim([0, 5]); % Uncomment and adjust if you need specific limits for the left y-axis

354

355% Plot the second variable on the right y-axis

356yyaxis right;

357plot(dates, cases);

358ylabel(’Confirmed New Cases in USA’);

359

360

361datetick(’x’, ’yyyy-mm-dd’);

362xlabel(’Date’);

363legend(’Mask Policy Level’, ’Confirmed New Cases’);

364

365

366%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure S1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

367

368susceptFracVec = 0:0.02:1;

369R0vec = 2;

370colorVec = jet(length(R0vec));

371transmissionReduction = 0:0.02:1;

372

373% SIR model parameters

374gamma = 1; % Recovery rate

375lambda1 = 100; % Rate of susceptibles with lower tolerance to switch to masked class

376delta1 = .1; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

377lambda2 = 1; % Rate of susceptibles with higher tolerance to switch to masked class

378delta2 = .1; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

379

380recoveredFracVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

381xVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

382effectivenessVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

383

384for z = 1:length(R0vec)

385 for y = 1:length(transmissionReduction)

386 for x = 1:length(susceptFracVec)

387 recoveredFracVec(x,y,z) = plotSIRDynamics(susceptFracVec(x), R0vec(z),transmissionReduction(y), 0, gamma, lambda1, delta1, lambda2, delta2, 0); % individuals adopt based on incidence

388 xVec(x,y,z) = susceptFracVec(x);

389 effectivenessVec(x,y,z) = transmissionReduction(y);

390 end

391 end

392end

393

394

395figure

396

397subplot(1,2,1)

398[x_new, e_new] = meshgrid(susceptFracVec, transmissionReduction);

399contourf(x_new,e_new,recoveredFracVec(:,:,z)’, ’LevelStep’, 0.01)

400xlabel(’Fraction of Population with Low Risk Tolerance’, ’FontSize’, 14)

401ylabel(’Effectiveness of Intervention (\epsilon)’, ’FontSize’, 14)

402h = colorbar;

403title(h, ’Epidemic Size’, ’FontSize’, 14)

404title([’R_0=’, num2str(R0vec(z)), ’, \lambda_1=’, num2str(lambda1), ’, \delta_1=’, num2str(delta1), ’, \lambda_2=’, num2str(lambda2), ’, \delta_2=’, num2str(delta2)]);

405

406hold on

407points_x = 0.6*ones(100,1);

408points_y = rand(100,1);

409

410scatter(points_x, points_y, ’r’, ’filled’)

411

412subplot(1,2,2)

413

414for a = 1:length(points_x)

415 plotSIRDynamics(points_x(a), R0vec(z),points_y(a), 1, gamma, lambda1, delta1, lambda2, delta2, 0); % individuals adopt based on incidence

416end

417

418

419%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure S2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

420

421susceptFracVec = 0:0.02:1;

422R0vec = 2;

423colorVec = jet(length(R0vec));

424transmissionReduction = 0:0.02:1;

425

426% SIR model parameters

427gamma = 1; % Recovery rate

428lambda1 = 100; % Rate of susceptibles with lower tolerance to switch to masked class

429delta1 = .1; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

430lambda2 = 1; % Rate of susceptibles with higher tolerance to switch to masked class

431delta2 = .1; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

432

433recoveredFracVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

434xVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

435effectivenessVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

436

437for z = 1:length(R0vec)

438 for y = 1:length(transmissionReduction)

439 for x = 1:length(susceptFracVec)

440 recoveredFracVec(x,y,z) = plotSIRDynamics(susceptFracVec(x), R0vec(z),transmissionReduction(y), 0, gamma, lambda1, delta1, lambda2, delta2, 0); % individuals adopt based on incidence

441 xVec(x,y,z) = susceptFracVec(x);

442 effectivenessVec(x,y,z) = transmissionReduction(y);

443 end

444 end

445end

446

447

448figure

449

450subplot(1,2,1)

451[x_new, e_new] = meshgrid(susceptFracVec, transmissionReduction);

452contourf(x_new,e_new,recoveredFracVec(:,:,z)’, ’LevelStep’, 0.01)

453xlabel(’Fraction of Population with Low Risk Tolerance’, ’FontSize’, 14)

454ylabel(’Effectiveness of Intervention (\epsilon)’, ’FontSize’, 14)

455h = colorbar;

456title(h, ’Epidemic Size’, ’FontSize’, 14)

457title([’R_0=’, num2str(R0vec(z)), ’, \lambda_1=’, num2str(lambda1), ’, \delta_1=’, num2str(delta1), ’, \lambda_2=’, num2str(lambda2), ’, \delta_2=’, num2str(delta2)]);

458

459hold on

460points_x = rand(100,1);

461points_y = 0.3*ones(100,1);

462

463scatter(points_x, points_y, ’r’, ’filled’)

464

465subplot(1,2,2)

466

467for a = 1:length(points_x)

468 plotSIRDynamics(points_x(a), R0vec(z),points_y(a), 1, gamma, lambda1, delta1, lambda2, delta2, 0); % individuals adopt based on incidence

469end

470

471%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure S4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

472

473susceptFracVec = 0:0.02:1;

474R0vec = 2;

475colorVec = jet(length(R0vec));

476transmissionReduction = 0:0.02:1;

477

478% SIR model parameters

479gamma = 1; % Recovery rate

480lambda1 = 100; % Rate of susceptibles with lower tolerance to switch to masked class

481delta1 = .1; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

482lambda2 = 1; % Rate of susceptibles with higher tolerance to switch to masked class

483delta2 = .1; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

484

485recoveredFracVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

486xVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

487effectivenessVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

488

489for z = 1:length(R0vec)

490 for y = 1:length(transmissionReduction)

491 for x = 1:length(susceptFracVec)

492 recoveredFracVec(x,y,z) = plotSIRDynamics(susceptFracVec(x), R0vec(z),transmissionReduction(y), 0, gamma, lambda1, delta1, lambda2, delta2, 1); % individuals adopt based on infection level

493 xVec(x,y,z) = susceptFracVec(x);

494 effectivenessVec(x,y,z) = transmissionReduction(y);

495 end

496 end

497end

498

499

500figure

501

502subplot(1,2,1)

503[x_new, e_new] = meshgrid(susceptFracVec, transmissionReduction);

504contourf(x_new,e_new,recoveredFracVec(:,:,z)’, ’LevelStep’, 0.01)

505xlabel(’Fraction of Population with Low Risk Tolerance’, ’FontSize’, 14)

506ylabel(’Effectiveness of Intervention (\epsilon)’, ’FontSize’, 14)

507h = colorbar;

508title(h, ’Epidemic Size’, ’FontSize’, 14)

509title([’R_0=’, num2str(R0vec(z)), ’, \lambda_1=’, num2str(lambda1), ’, \delta_1=’, num2str(delta1), ’, \lambda_2=’, num2str(lambda2), ’, \delta_2=’, num2str(delta2)]);

510

511hold on

512points_x = 0.6*ones(100,1);

513points_y = rand(100,1);

514

515scatter(points_x, points_y, ’r’, ’filled’)

516

517subplot(1,2,2)

518

519for a = 1:length(points_x)

520 plotSIRDynamics(points_x(a), R0vec(z),points_y(a), 1, gamma, lambda1, delta1, lambda2, delta2, 1); % individuals adopt based on infection level

521end

522

523

524%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure S5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

525

526susceptFracVec = 0:0.02:1;

527R0vec = 2;

528colorVec = jet(length(R0vec));

529transmissionReduction = 0:0.02:1;

530

531% SIR model parameters

532gamma = 1; % Recovery rate

533lambda1 = 100; % Rate of susceptibles with lower tolerance to switch to masked class

534delta1 = .1; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

535lambda2 = 1; % Rate of susceptibles with higher tolerance to switch to masked class

536delta2 = .1; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

537

538recoveredFracVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

539xVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

540effectivenessVec = zeros(length(susceptFracVec),length(transmissionReduction),length(R0vec));

541

542for z = 1:length(R0vec)

543 for y = 1:length(transmissionReduction)

544 for x = 1:length(susceptFracVec)

545 recoveredFracVec(x,y,z) = plotSIRDynamics(susceptFracVec(x), R0vec(z),transmissionReduction(y), 0, gamma, lambda1, delta1, lambda2, delta2, 1); % individuals adopt based on infection level

546 xVec(x,y,z) = susceptFracVec(x);

547 effectivenessVec(x,y,z) = transmissionReduction(y);

548 end

549 end

550end

551

552

553figure

554

555subplot(1,2,1)

556[x_new, e_new] = meshgrid(susceptFracVec, transmissionReduction);

557contourf(x_new,e_new,recoveredFracVec(:,:,z)’, ’LevelStep’, 0.01)

558xlabel(’Fraction of Population with Low Risk Tolerance’, ’FontSize’, 14)

559ylabel(’Effectiveness of Intervention (\epsilon)’, ’FontSize’, 14)

560h = colorbar;

561title(h, ’Epidemic Size’, ’FontSize’, 14)

562title([’R_0=’, num2str(R0vec(z)), ’, \lambda_1=’, num2str(lambda1), ’, \delta_1=’, num2str(delta1), ’, \lambda_2=’, num2str(lambda2), ’, \delta_2=’, num2str(delta2)]);

563

564hold on

565points_x = rand(100,1);

566points_y = 0.3*ones(100,1);

567

568scatter(points_x, points_y, ’r’, ’filled’)

569

570subplot(1,2,2)

571

572for a = 1:length(points_x)

573 plotSIRDynamics(points_x(a), R0vec(z),points_y(a), 1, gamma, lambda1, delta1, lambda2, delta2, 1); % individuals adopt based on infection level

574end

575

576

577

578%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure S6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

579% Repeat this script with corresponding parameters to generate each panel

580

581figure

582susceptFracVec = 0.2;

583R0vec = 2;

584colorVec = jet(length(R0vec));

585transmissionReduction = .7;

586

587% SIR model parameters

588cVec = [0.01:0.05:1];

589lambdaAvg = [0.5, 1:5:100];

590gamma = 1; % Recovery rate

591delta1 = .5; % Rate of masked class with lower tolerance for adoption to remove intervention (mask) to return to susceptible class

592delta2 = .5; % Rate of masked class with higher tolerance for adoption to remove intervention (mask) to return to susceptible class

593

594recoveredFracVec = zeros(length(cVec),length(lambdaAvg));

595

596for x = 1:length(cVec)

597 for y = 1:length(lambdaAvg)

598 recoveredFracVec(x,y) = plotSIRDynamics(susceptFracVec, R0vec, transmissionReduction, 0, gamma, lambdaAvg(y)*(1-susceptFracVec*cVec(x))/(1-susceptFracVec), delta1, cVec(x)*lambdaAvg(y), delta2, 1); %Arithmetic parameterization

599 end

600end

601

602[c_new, l_new] = meshgrid(cVec, lambdaAvg);

603contourf(c_new,l_new,recoveredFracVec’, ’LevelStep’, 0.01)

604xlabel(’\lambda Homogeneity Index’, ’FontSize’, 14)

605ylabel(’\lambda_{Arithmetic Average}’, ’FontSize’, 14)

606h = colorbar;

607title(h, ’Epidemic Size’, ’FontSize’, 14)

608title([’R_0=’, num2str(R0vec), ’, x_1=’, num2str(susceptFracVec), ’, \epsilon=’, num2str(transmissionReduction), ’, \delta_1=’, num2str(delta1),’, \delta_2=’, num2str(delta2)]);