It is well-known that prophylactic behavior affects disease transmissions by reducing rates of effective contacts between infectious and susceptible individuals in the target population (direct physical or indirect through shared media in which pathogens can survive) [23, 42]. In our modeling framework, change in disease dynamic affects back prophylactic behavior through observable disease prevalence (which depend, for instance, on the existence and effectiveness of a disease surveillance mechanism and testing capacity to detect infectious individuals), but disease related information is processed by a system of health beliefs and related attitudes. In other words, the feedback from disease dynamic to prophylactic behavior is mediated by health beliefs and normative attitudes in the individual's interpretation domain which determines risk tolerance. However, engagement in prophylactic behavior is also determined by how many individuals are already behaving prophylactically, and a social group can pay more attention to and mimick those behaving prophylactically than those behaving non-prophylactically, or vice-versa.

To describe disease dynamics, we consider an extended SEIR compartmental model framework [43–45] with differentiated infective states. The population size N is given at time t by

where S, E, I_a, I_s, I_d, and R are susceptible, exposed, asymptomatic infectious, symptomatic infectious, detected (tested positive and reported), and removed individuals, respectively. The disease-dependent compartments E, I_a, I_s, I_d, and R are considered homogeneous whereas the class of susceptible individuals (S) is further differentiated into two groups based on risk tolerance:

where S_i denotes susceptibles with standard of risk i∈A, A={-1,1}. For instance, in a population which only responds to current risk, S₋₁ may represent individuals with low standard of evidence (the most responsive to disease prevalence), while S₁ corresponds to individuals with high standard of evidence (least responsive to disease prevalence). However, difference in standard of risk may be qualitative rather than quantitative: S₋₁ may represent individuals responsive to only historical risk while S₁ corresponds to individuals sensible to both historical and predicted risk.

The primary source of disease-related information in an epidemic context is the timely number of new positive cases, denoted C. The number of individuals in the class I_d of detected infectious individuals (see Equation 1a) is related to C by

where the dot notation indicates the first derivative with respect to time (i.e., İ_d(t) = dI_d(t)/dt), ρ_d is the removal rate of individuals from I_d, with non-negative initial condition I_d(0) = I_d0. Solving the differential Equation (2) for I_d gives

We consider the perceived disease prevalence, denoted P, and the relative rate of change of new positive cases, denoted Q, as the basic pieces of information on which susceptible individuals will decide to be prophylactic. From Equation (3), the perceived disease prevalence is given at time t by

where the constant τ > 0 represents a time delay in the acquisition of information on detected and reported infectious. The relative rate of change Q is the quotient of the rate of change (increase or decrease) of the number of new detected cases C to the number of new detected cases, delayed by τ time units:

Note that Q(t) = 1 means that the timely number of new positive cases is doubling per unit time, whereas Q(t) = −1/2 means that the timely number of new positive cases is halving per unit time.

We assume that based on the information pieces P and Q, each group of susceptible S_i makes up an information aggregate, denoted η_i, and which satisfies η_i = 0 before disease outbreak (t < τ, P = Q = 0). This information aggregate is defined as a quadratic function of P and Q:

where a_i, b_i, c_i, d_i, and e_i are non-negative real coefficients expressing the weights of linear, quadratic, and interaction components of the prevalence P and the rate of change Q. It is worthwhile noticing that Equation (4c) is intended as an approximate summary of the available information that drives the decision to wear a mask. The quadratic form is indeed used as an approximation to the actual, likely non-linear mechanism by which the susceptible group S_i processes disease-related information. The non-negative signs imposed on the coefficients in η_i (a_i, b_i, c_i, d_i and e_i) ensure that when the timely number of new positive cases is non-decreasing (Q ≥ 0), η_i is not only non-negative for any perceived disease prevalence P but also non-decreasing in P.

Equation (4c) describes how disease risk information can be differently interpreted by different tolerance groups within the same population. For instance, a group of individuals may focus on the disease incidence rate while ignoring predictions of future disease risk. For such individuals, one or both parameters a_i and b_i in the information aggregate (Equation 4c) will be positive but parameters c_i, d_i, and e_i will be zero. Another group may pay attention to trend in disease incidence, opening the possibility of adopting prophylactic behavior earlier when the trend moves upwards, and also giving up such behavior when incidence decreases. Such a group will have at least one of the parameters c_i, d_i, and e_i being greater than zero. Between these two extreme situations, there are many possibilities with various combinations of small vs. large values of parameters a_i and b_i to reflect how much attention is paid to current disease incidence by a tolerance group, and parameters c_i, d_i, and e_i related to the relative importance of trends to the group.

For a susceptible class with standard of evidence i, we consider for simplicity two levels of prophylactic behavior: high prophylactic behavior (i.e., individuals properly wearing mask where and when this is recommended) vs. low prophylactic behavior. The overall prophylactic behavior in the class S_i can thus be summarized by the prophylactic proportion m_i ∈ [0, 1], i.e., the proportion of mask-wearers. We assume that change in the prophylactic proportion m_i is proportional to change in the information aggregate η_i and to the proportion of S_i individuals already wearing masks:

where α_i is a positive real which determines the nature and strength of in-group behavior, and we have taken the proportionality constant to be one to ensure that the coefficients of the linear components of η_i (a_i and e_i) in the model are identifiable from observed data. When α_i = 1, ∂m_i/∂t is proportional to both m_i and 1 − m_i, and S_i individuals give the same relative importance to both mask-wearers and non-mask-wearers: for a unit increase in the information aggregate, the highest increase in m_i occurs when m_i = 0.5, i.e., when half of the S_i individuals have engaged in prophylactic behavior. For general α_i values, the highest increase in m_i for a unit increase in η_i occurs when mi=(1+αi)-1/αi. It appears that for α_i ∈ (0, 1), the highest increase in m_i occurs when m_i < 0.5 (weak influence of in-group non-prophylactic behavior), and for α_i > 1, the highest increase in m_i occurs when m_i > 0.5 (strong influence of in-group non-prophylactic behavior). Hence, the larger α_i, the larger impact in-group non-prophylactic behavior will have, slowing down engagement in prophylactic behavior.

Since m_i depends on time only through the information aggregate η_i, we can interpret Equation (5) as a differentiation in chain and write ∂mi/∂ηi=mi(1-miαi)/αi, which appears to be Richards growth equation [38] with intrinsic growth rate equal to one. Solving for m_i yields the generalized logistic curve:

where δ_i is a constant related to the proportion m_i0 ∈ (0, 1) of S_i susceptibles who would hold a high prophylactic attitude even in the absence of any evidence of disease (i.e., when η_i = 0, which happens for t ≤ τ) by

Figure 1 shows the prophylactic proportion m_i as a function of time for a few selected parameter values, with a perceived disease prevalence P varying from zero to 33%. It appears that the coefficients b_i, c_i, d_i, and e_i in Equation (4c) have distinct effects on m_i and capture different reactions of susceptibles to disease risk. For instance, b_i and d_i can be described as response acceleration parameters for high and low prevalence values, respectively. Indeed, a susceptible group more responsive to large prevalence values than to low prevalence values corresponds to b_i > 0 (Figure 1A), whereas d_i > 0 corresponds to groups more responsive to low prevalence values than to high prevalence values (Figure 1C). Similarly, c_i > 0 corresponds to groups where engagement in prophylactic behavior is stronger (Figure 1B), and e_i > 0 to groups where engagement in prophylactic behavior is earlier (Figure 1D), but in both cases, disengagement also happens early, once the number of new detected case starts dropping (Q < 0).

The in-group behavior parameter α_i allows additional flexibility in m_i by controlling how the aggregated information is jointly used with how many individuals are already behaving prophylactically in a susceptible group. Whereas b_i and d_i can be viewed as parameters inducing behavioral response acceleration with respect to change in risk information (P and Q, respectively), α_i is an intrinsic acceleration/deceleration parameter, i.e., the acceleration of m_i happens not because of change in risk information, but rather in response to the current (low or high) value of m_i itself. Figure 2 shows m_i curves for a few selected parameter values. It can be observed that ceteris paribus, a larger α_i value, implies an overall weaker prophylactic behavioral response. Indeed, the derivative of m_i with respect to α_i, given by

is negative (m_i decreases with α_i) for η_i > 0 (i.e., when prophylactic proportion is above the disease-free level m_i0). However, if information aggregate reaches zero (η_i = 0), ∂m_i/∂α_i (Equation 7) vanishes, and if η_i becomes negative, ∂m_i/∂α_i > 0. This happens around t = 96 days in Figure 2D, where e_i > 0 (i.e., in a population where Q is given much attention) allows η_i < 0 after Q becomes negative, and as the epidemic dies out (P → 0), Q dominates the information aggregate: the ordering of prophylactic proportions switches such that a lower α_i value corresponds to a weaker prophylactic behavioral response for η_i < 0. This can be interpreted as a return of the early engagement of a group with e_i > 0 (Figure 2D, m1(0.1)). In this respect, α_i appears as a parameter which exaggerates behavioral response regardless of the sign of information aggregate. But as the disease-free prophylactic proportion m_i0 is typically low (m_i0 = 0.05 in Figure 2), there will generally be less room for this exaggeration when the epidemic dies out (η_i < 0) than at disease outbreak (η_i > 0).

We assume for simplicity that the detected infectious individuals are isolated (e.g., hospital and home) and do not mix actively with other classes. Using the “quarantine-adjusted” incidence mechanism [46] yields the force of infection λ_i (the average number of adequate contacts of one S_i susceptible person with infectives per unit time):

where β_ia and β_is are rates of effective contacts with asymptomatic and symptomatic infectious individuals, respectively. The effective contact rate β_ij depends on a baseline contact rate β₀ (possibly restricted by public health policies), the prophylactic attitude of S_i susceptibles, the average efficiency κ ∈ (0, 1) of prophylactic behaviors in reducing transmissions, and the probability ϕ_j ∈ (0, 1) of disease transmission on contact with I_j infections:

For a target population, we consider a period of study short enough for both disease-related death and natural demographic rates (births, net immigration, and deaths) to be negligible relative to the total population size N (assumed large but finite). As a result, the population size in Equation (1a) remains constant and equal to an initial size N₀ = N(0). Joining the behavior and disease dynamics mechanisms described in Sections 2.3, 2.4 gives the Behavior-SEIR model depicted on the flow diagram in Figure 3 with parameters described in Table 1. After sufficient contacts with infectious individual(s), a susceptible individual enters an incubation period (class E) lasting 1/θ time units on average, in a non-infectious state, and without any disease symptom. Some of these exposed individuals are early detected with probability π thanks to contact tracing or systematic tests on target groups, and enter the class I_d. In the non-early detected exposed group, 100σ% develop symptoms and enter the class I_s, and 100(1−σ)% remain asymptomatic and enter the class I_a. Individuals in the class I_s are then identified at a high rate γ_s. Thanks to contact tracing or systematic tests on target groups again, some asymptomatic individuals in the class I_a are identified at a lower rate γ_a. All exposed individuals eventually recover from the disease, entering the class R (removals).

The Behavior-SEIR dynamics model is described at time t by the following system of non-linear differential equations:

with the non-negative initial conditions S_i(0) = S_i0, E(0) = E₀, I_a(0) = I_a0, I_s(0) = I_s0, I_d(0) = I_d0, and R(0) = R₀ such that N₀ = S₋₁₀+S₁₀+E₀+I_a0+I_s0+I_d0+R₀. The proposed model is a system of neutral delay differential equations [37] where the force of infection λ_i (Equation 8a) depends through m_i (Equation 8b) on the information available at time t, precisely P (dependence on state variables) and Q (dependence on first derivatives of state variables), which are delayed by τ time units.

The parameter τ accounts for two potential sources of information delay: (i) reporting delay, i.e., the delay between the moment exposed or infectious individuals are detected and isolated from the mixing population, and the moment the number of detected cases is publicized and can be considered by susceptible individuals to assess their risk, and (ii) reaction time, i.e., the delay between the moment the number of detected cases is made public and the moment susceptible individuals actually consider the information to adjust their adherence to preventive prophylactic behavior. We focus on reporting delay which can be included in public health policy design [47]. Indeed, official statistics are often reported with a time delay that may arise from a desire for thorough verification [48]. However, reporting delays can produce the dangerous illusion of an improving epidemic situation since the most recent days have the least cases accounted [48]. Including the time delay parameter τ in model (Equation 9) allows to investigate the extent to which reporting delay can affect the evolution of an epidemic through the behavioral response to delayed information. For our simulation experiments (see Section 2.8), we consider information delays ranging from 1 day to 1 week (τ = 1, 3, 5, 7 days).

We compute the effective reproductive number based on the Behavior-SEIR model (Equation 9) using the next-generation matrix approach [49]. Starting from any disease-free state Xc=(S-10,S10,0,0,0,0,R0)⊤, the basic reproductive number R(0) for system Equation (9) is given by

where N₀ = S₋₁₀+S₁₀+R₀ > 0 and Ro (Equation 10b) is the reproductive number when there is no differential evidentiary group, R₀ = 0, and all prophylactic proportions are zero (m_i0 = 0). As expected, the basic reproductive number R(0) (Equation 10a) depends on both the distribution of the population between evidentiary groups, and the prophylactic proportion of each group when no disease evidence is available (both perceived disease prevalence P and rate of change Q are zero). Along an epidemic, the effective reproductive number is then given by

where Ri(t) (Equation 10d) is the effective reproductive number for S_i susceptibles such that R(t) is the average of Ri(t) over all groups of susceptibles.

We would normally conduct here a stability analysis of the model, including the bifurcation diagram for interesting model parameters such as the time delay τ and the level of protection κ. However, these investigations are out of the scope of this conceptual analysis which focuses on the behavioral response to an outbreak and its impact on the dynamic of an epidemic. Such investigations can however be carried out after extending the model to represent more realistic scenarios, including, for instance, vital rates (births and immigration, natural deaths, and disease-related death), flow between susceptible compartments (S₋₁ and S₁), and immunity lost (flow from R back into S₋₁ and S₁).

To allow for comparison between various epidemic scenarios, we define some measures to quantify epidemic severity and overall disease dynamic.

We carried out a simulation experiment to explore combinations of model parameters related to disease transmission (β₀), early detection of exposed individuals (π), delayed acquisition or reaction to risk information (τ), behavioral reaction (a_i, b_i, c_i, d_i, and e_i), in-group behavior (α_i), and efficiency of prophylactic behaviors (κ) that lead to patterns where the course of an epidemic depends more or less on differences between risk perception and related feedback loops.

Our simulation results indicate important variations in epidemic severity measures between levels of in-group pressure. Indeed, an increase in in-group pressure for pro-prophylactic behavior (decreasing α_i) increases the average number of secondary waves, but decreases the epidemic peak size and final size (Table 2).

Table 3 shows summaries of epidemic severity measures comparing populations with the reference response profile 0 and profiles A–B, averaged over all levels of in-group pressure (α_i ∈ {0.1, 1, 2, 3}). All investigated parameter settings result into an epidemic, the disease quickly dying out in only 1% of settings with the homogeneous profile B (i.e., populations with a reaction to rate of change Q proportional to perceived prevalence P). When there is an epidemic, one secondary epidemic wave occurs in 9% of settings with the reference profile 0. Secondary epidemic waves are the most likely (19%) under profile A [homogeneous populations with a strong (quadratic) reaction to P only], and the least likely (2%) under profile B. The epidemic peaks after about 64 days on average in populations responsive to P only (profile 0) to ~3,209 new reported cases. Slightly lower average peak time and size result from a stronger reaction to P (3,050 case after 62 days for profile A) or an additional response to Q (3,126 cases after 60 days for profile D). Ignoring the sign of Q increases the average peak time to about 72 days (profile C) while decreasing the peak size to ~32,748 cases, with intermediate result for profile B (2,840 cases after 66 days). The final epidemic size is typically large, 91% on average, under profile 0 (and C). A stronger reaction to P (profile A) slightly reduces the final size to 90% whereas an additional response to Q (profiles B and D) increases the final size to 94%. In summary, a stronger reaction to P essentially leads to more secondary waves whereas an additional response to Q hastens the epidemic peak and increases the final size, with lower peak sizes in both scenarios. A response to Q proportional to P also increases average final size, but delays the peak (with a lower size).

Heterogeneity in the behavioral reaction to risk information generally leads to intermediate results half way between the two combined homogeneous profiles (Table 3). Exceptions include the time to curb the first epidemic wave (larger than expected) under profiles 0 × C and 0 × D, and the peak time (larger than expected) and size under profile 0 × D (smaller than expected). Figure 4 shows the timely number of new positive cases under profiles 0, D, and 0 × D, exhibiting an interactive effect between in-group pressure and heterogeneity. Indeed, it appears that in an heterogeneous pro-prophylactic behavior population (Figure 4, α = 0.1) where half of individuals are responsive to P and the other half is additionally responsive to Q, the epidemic peak size (1,380 new detected cases) is smaller than in homogeneous populations with profile 0 (1,819 new detected cases) or D (1,653 new detected cases).

When we restrict attention to populations where in-group pressure is neutral (α_i = 1), we observe that a stronger reaction to P leads to one secondary epidemic wave in 24% of settings, as opposed to 19% across α_i values (Supplementary Table S3). Similarly, an additional response to Q increases the peak time and does not increase the final epidemic size (90% under profile D), unless the reaction to Q is proportional to P (93% under profile B). These discrepancies point out to important variations between levels of in-group pressure.

The results of fitted models shown in Table 4 corroborate our observations for the number of secondary waves (model coefficient estimates < 0 for α₋₁ and Δα), the peak size, and the final epidemic size (coefficient estimates > 0 for α₋₁ and Δα). It also appears that among the profiles A–D, only profile B (population with a response to the rate of change proportional to prevalence) leads to a decrease in the number of secondary epidemic waves and an increase in the final epidemic size, as compared to the reference profile 0. In addition, the model results indicate that for an heterogeneous population, the expected number of secondary epidemic waves or the expected final epidemic size is intermediate between the outcomes for the two corresponding homogeneous populations, except when half the population has profile D. For the latter, after controlling for disease dynamic and in-group pressure parameters, the expected number of secondary epidemic waves is 15.9% (100 × 0.159) higher for the homogeneous profile D, but 31.8% higher for the heterogeneous profile 0 × D, as compared to an homogeneous profile 0.

The fitted models also indicate that a 1-day increase in risk information delay (τ) leads to a 2.6% average increase (100 × 0.026) in the average odd ratio for a random individual to get infected over the course of the epidemic, a 3.3% increase in the overall peak size, and a 6.5% decrease in the number of secondary waves (Table 4). A discussion of the biological interpretation of these statistics is provided in Section 4.1. As for other varied model parameters (β₀, π, and κ), apart from in-group pressure, the variations of the number of secondary waves are mostly driven by the level of protection by prophylactic behavior, to the extent that a 1%-point increase in κ results into a 15.5% increase in the number of secondary waves. In other words, if 80% level of protection yields one epidemic wave on average, then increasing the level of protection by 10% (i.e., from 80 to 90%) results in an expected 2.55 waves (1+1 × 0.155 × 10). Both the peak size and the final epidemic size are mainly determined by the baseline transmission rate (β₀) and the probability of early detection of exposed individuals (π).

Figures 5, 6 show the cumulative number of cases detected up to 1,000 days after the outbreak for various values of κ, π, and β₀. For κ ∈ {0.5, 0.6}, no secondary epidemic wave was observed (Figure 5). It appears that if the disease surveillance mechanism for early detection and removal is loose (π = 0.25), disease dynamic is barely sensitive to behaviors (Figure 5A). When disease surveillance is more effective (π≥ 0.5), disease dynamic (peak and observed final size) becomes more sensitive to behaviors (Figures 5B–D), especially when the baseline disease transmission rate is low (β₀ = 0.25).

When prophylactic behavior offers (almost) perfect protection (κ ∈ {0.95, 1}), secondary epidemic waves (one or two) were observed in 28% of these settings. Figure 6 shows the cumulative numbers of detected cases when there are secondary waves (see Supplementary Figure S1 for unique wave scenarios). It can be observed that at such high levels of protection by prophylactic behavior, disease dynamic is highly sensitive to behavioral changes, even when the baseline disease transmission rate is high (β₀ = 3).

When the efficiency of prophylactic behavior is between these two extremes (κ ∈ {0.75, 0.9}), secondary epidemic waves (one or two) were observed in 4% of these settings (Supplementary Figures S2, S3): with these levels of protection by prophylactic behavior, disease dynamic is already sensitive to behavioral changes, especially when the baseline disease transmission rate is β₀ < 2, and the higher the probability of early detection (π), the higher the number of secondary epidemic waves.

In this study, we introduce a new behavior-disease compartmental model where risk perception is a function of both perceived disease dynamic and an interpretation domain (in-group pressure). That is, the risk information derived from disease dynamics can include predicted incidence (as expressed by the rate of change Q of new positive cases) in addition to historical incidence (as expressed by the perceived prevalence P given the disease surveillance system), and the actual risk perceived by an individual can arise not only from disease dynamics but also from the pro/non-prophylactic behavior of others in the individual's social group.

Under the social influence-based model which uses true disease incidence as risk information, Tyson et al. [14] found that populations more responsive to risk information can experience more severe epidemics in terms of final size and undergo multiple epidemic waves, although the epidemic peak sizes will be smaller.

We observe similar trends when matching increase in social influence with increase in pro-prophylactic behavior in-group pressure. Considering populations responsive to historical risk information, our results regarding the number of epidemic waves and the peak size are consistent with Tyson et al.'s [14] findings, i.e., a stronger reaction to perceived disease prevalence produces more epidemic waves but smaller peak size. Indeed, a strong reaction to prevalence results in transient prophylaxis that slows down and stops disease progression (smaller peak) while many susceptible individuals are still in the population. Low disease prevalence then leads to relaxation of prophylactic behaviors and subsequently to a new outbreak (multiple waves).

A new aspect of behavior-disease dynamic captured by our model is the importance of response profiles. For instance, we considered two kinds of populations, ones that are only responsive to historical risk and ones that are additionally responsive to predicted risk. For risk predicting populations, the epidemic peak size does not monotonically decrease with increasing pro-prophylactic in-group pressure (the epidemic peak size is smallest when in-group behavior is neutral). In other words, in the context of risk information overload [36], in particular, when disease evolution curves are overly discussed on mass media and social media, stronger social influence or in-group pressure can lead to more severe epidemic outcome, at least as measured by the epidemic peak size.

Our results also indicate that erratic disease evolution curves can be explained by strong behavioral response to predicted disease curves shown on mass media or in social media. For a population paying attention to predicted risk information, a strong response reduces the time scale of the chain of reactions that leads to multiple epidemic waves, giving rise to very fast oscillations in the observed disease incidence curve (see Figure 4), both as the epidemic establishes (before peak) and as the epidemic is waning. This risk-prediction feature also makes our model quite different from others where each epidemic peak is necessarily followed by an almost disease-free interval before emergence of the following peak [14, 57]. Because of this, which is rather frequent is real epidemic data [50], our model requires a wave delimitation method to identify epidemic waves. Also note that, unlike in Aziz-Alaoui et al. [2], our model targets a short-term dynamic, and epidemic waves in our model framework are not related to immunity lost (no flow from R back to S_i in Figure 3) but fully generated by the behavioral response to the outbreak.

Our model allows us to explore the impact of heterogeneity in the behavioral response of a population to disease risk information. In general, two-group heterogeneity in response profile (pro- vs. contra-prophylactic behavior in-group pressure, or responsiveness to predicted risk or not) leads to an intermediate epidemic outcome as compared to the two sub-populations evolving separately. One interesting finding is that under strong pro-prophylactic in-group pressure, a population consisting of two same sized sub-populations where only one group is exposed/responsive to predicted risk information can experience a less severe epidemic as measured by peak size.

Finally, we investigated the effect of delayed risk information on the severity of an epidemic in our model framework. Our simulation experiments indicate that delayed risk information slows down the behavioral response to the progression of the epidemic, contributing to a more severe epidemic outcome, i.e., larger peak and final size. For an individual with a 50% average risk to get infected in a 1-day reporting delay context, the risk to get infected becomes 54% if the reporting delay is 1 week (6 days increase). Although a 4% increase may appear small at an individual level, it would represent 4,000 more infections, given the population size in our simulations (community of 100,000 individuals). Similarly, for an epidemic that peaks to 1,500 new detected cases (see, e.g., Figure 4) under a 1-day reporting delay scenario, a reporting delay of 1 week would result into 50 additional detected cases on the peak day. These simulated results are in accordance with the work of Gutierrez et al. [47], who found that the COVID-19 epidemic in Mexico progressed much faster when delays are larger, resulting in more severe epidemic outcomes (larger death peak size and cumulative death toll). A viable solution for policymakers to reduce information delay-related increase in epidemic severity is to use a nowcasting technique to adjust the daily number of confirmed new cases for delayed reporting [48], especially in populations reactive to predicted risk, since delayed risk information creates an illusion of a downward trend. We must remind the reader, however, that our primary interest is to do a conceptual analysis and none of these results should be taken literally without empirical validation.

Although the proposed model framework is quite general for coupling risk tolerance and disease dynamics, we have in our presentation limited attention to populations with relatively simple structures. In this section, we highlight and discuss the most important assumptions that may strongly affect our conclusions, and explore some potential routes to their relaxation. The prophylactic proportion defined in Equation (6) makes the strong assumption that risk information derived from in-group behavior is of no use if there is no change in risk information aggregate. However, in-group behavior can lead to change in prophylactic behavior even if risk information is constant, i.e., people can change behavior by inferring the need to engage (or the appropriateness of disengaging) based on how many individuals are already (or are still) behaving prophylactically, not because the risk has increased or decreased but because there is extant risk and in-group pressure has changed the risk perception.

Another strong assumption in Equation (6) is that the evolution of prophylactic behavior does not depend on how long people have already been engaging in prophylactic behavior, and how strong their engagement was during that time. However, prophylactic behavior is often subject to fatigue. Indeed, preventive behavior fatigue during secondary epidemic waves of a disease has been particularly well-documented since the COVID-19 pandemic event [58–75]. Failing to account for prophylactic behavior fatigue in behavior-disease model likely introduces bias in model outputs.

To account for the above two limitations of Equation (6), we propose to extend Equation (5) describing change in the proportion of prophylactic individuals in a susceptible group i as:

where ω_i ≥ 0 is the weight of “indirectly perceived change in risk” derived from in-group behavior (how informative in-group behavior is about perceived risk) and f_i ∈ [0, 1] represents prophylactic fatigue in group i. It follows from Equation (12) that Equation (6a) can be generalized as:

where we assumed that enthusiasm for prophylactic behavior decays exponentially as experience under disease-related restrictions increases, with ϵ ≥ 0 expressing the extent to which prophylactic behavior is exhausting (ϵ = 0 means that there is no prophylactic fatigue over time, f_i = 0), and the variable M_i quantifies how much effort the susceptible group i has invested in prophylactic behavior since the epidemic outbreak:

and M_i(0) = 0. Note from Equation (13a) that Ṁ_i can be rewritten as Ṁ_i = η_i(1 − f_i) + ζ_i. Equation (13a) is reduced to Equation (6a) when ω_i = 0 and f_i(t) = 0. To solve system (9) with the prophylactic proportion m_i given by Equation (13), we can extend the differential system Equation (9) with Equation (13d) for the pseudo-state variable M_i, and an additional pseudo-state variable ζ_i whose first derivative is given by the integrand in Equation (13c), i.e., ζ˙i=ηi[ḟi+ωi(1-fi)], which simplifies to

However, it appears that prophylactic fatigue, as introduced in Equations (12–14), affects equally engagement and disengagement in prophylactic behavior. A post-hoc but inelegant solution to that issue is to modify Equation (12) to have the form in Equation (15):

where g_i = 1 if ∂η_i/∂t ≥ 0, and 0 otherwise. We leave the implementation and exploration of such speculations for future work.

Another strong hypothesis of our SEIR framework is the well-mixture assumption, i.e., for each time point, the probability of interaction between two random individuals in the population is uniformly distributed [76]. While a low standard individual might tend to interact more with low standard individuals than high standard individuals, for instance, our model assumes an homogeneous mixture of low standard and high standard groups. More generally, heterogeneity in a population goes beyond risk tolerance groups, and may be related to other factors such as geographical location [77], age [78], and behavioral risk factors [79]. In particular, the mixture of individuals from various geographical regions is generally non-uniform, and explains the large disparities among different geographical locations in the COVID-19 pandemic context [80]. Thus, to be realistic, our model should be extended to include spatial components to reveal or account for the contribution of the spatial structure of individuals to an observed epidemic dynamic. This can be achieved using, for instance, non-autonomous coupling functions between adjacent areas [6], a diffusion process [81], an agent-based infection graph [82], or a location network [83].

Finally, we point out to the possibility of an evolving pathogen, as observed for the COVID-19 pandemic [84]. For a multi-strain pathogen, a realistic epidemic model should account for mutation process occurring during infection of individuals in the population [85, 86]. As such a future direction for our study is to consider many variants of a target disease and investigate behavioral feedback loops as a pathogen strain dominates the population and is then replaced by a new pathogen, appearing through evolutionary process or interactions with adjacent geographical regions.

In this study, we assessed the impact of differential behavioral response profile on epidemic outcomes. Our main contributions for understanding feedback loops between transient prophylaxis and disease dynamic include (i) the distinction of historical risk from predicted risk information overly discussed on mass media and social media, and (ii) the inclusion of an interpretation domain where the collected risk information is subjected to in-group pressure. It was known that the final size of an epidemic has a non-monotonic relation with the behavioral response of a population to risk information. Our results indicate that this non-monotonicity extends to epidemic peak size as a measure of epidemic severity, in populations under strong in-group pressure.

An obvious future direction is to assess the ability of this model to predict disease dynamics, using real epidemic data (from, e.g., the world health organization, https://covid19.who.int) along with behavioral change data from surveys (e.g., adherence to COVID-19 protective measures [87, 88]). The inclusion of more than two social groups, based, for instance, on age or spatial location may be integral parts for establishing some predictive power of the model. However, our aim here was not predictive, but rather to better understand how to represent assumptions about heterogeneous risk tolerance and in-group pressures, and then in turn study their potential effects on disease dynamics.