The paradox of neglecting changes in behavior: How standard epidemic models misestimate both transmissibility and final epidemic size

This is an uncorrected proof. Figures Abstract During epidemic outbreaks, populations adapt their behavior in response to disease burden, fundamentally altering transmission dynamics. Despite this, most compartmental models assume constant contact rates throughout outbreaks. To quantify biases from this assumption, we fitted a baseline SEIRD model with constant transmission and three behavioral variants—incorporating mortality-driven transmission reduction via exponential, rational, and mixed functional forms—to COVID-19 mortality data from 20 selected US locations during the first pandemic wave (March–July 2020). All three behavioral models achieved a lower median normalized sum of squared error in at least 18 of 20 locations, and Bayesian model selection favored them in at least 18 of 20 locations. More importantly, we identified systematic biases when behavioral responses are ignored: the baseline model consistently underestimated the basic reproduction number () while paradoxically overestimating the final epidemic size. Median estimates from the behavioral models exceeded the baseline estimates across all 20 locations, yet baseline models predicted larger cumulative infection burdens. Controlled synthetic experiments—where mortality trajectories were generated from behavioral models with known parameters—confirmed these biases result from model misspecification rather than data quality or stochastic variation. We prove analytically that for any fixed , the baseline model overestimates cumulative infections compared to behavioral models where mortality reduces transmission, regardless of functional form. This dual bias has potential implications for pandemic response: standard models may simultaneously underestimate pathogen contagiousness, which could contribute to delayed or insufficient early interventions while overestimating infection burden, which could bias planning for later epidemic phases. Our findings across 20 geographically diverse locations demonstrate that incorporating behavioral change substantially improves both model fit and estimation of epidemiological parameters relevant for public health policy. Author summary When diseases spread, people change their behavior—avoiding crowds, wearing masks, washing hands more frequently. Yet many mathematical models used to predict epidemics assume people behave the same way throughout an outbreak. We asked: What happens when models ignore these behavioral changes? Using COVID-19 data from 20 selected US locations, we compared a traditional model, which assumes constant transmission, against models where transmission decreases as deaths increase. We found a troubling paradox: models ignoring behavior consistently suggest diseases are less contagious than estimates from behavioral models, yet simultaneously predict that more people will ultimately become infected. This mismatch matters for decision-making: underestimating how easily a disease spreads could contribute to delayed or insufficient early interventions, while overestimating total infections could bias planning for later epidemic phases. We confirmed this paradox through computer simulations where we knew the true answer, and showed analytically that behavioral adaptation reduces final epidemic size when the basic reproduction number is held constant. Our work shows that incorporating human behavioral responses is not just a modeling refinement—it is important for epidemic predictions that inform public health decisions. Citation: Pant B, Lalovic M, Kiss IZ, Santillana M (2026) The paradox of neglecting changes in behavior: How standard epidemic models misestimate both transmissibility and final epidemic size. PLoS Comput Biol 22(9): e1014746. https://doi.org/10.1371/journal.pcbi.1014746 Editor: Nicola Perra, Queen Mary University of London, UNITED KINGDOM OF GREAT BRITAIN AND NORTHERN IRELAND Received: December 17, 2025; Accepted: August 20, 2026; Published: September 3, 2026 Copyright: © 2026 Pant et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The complete codebase used for simulation, inference, and figure generation is available at https://github.com/markolalovic/epi-behavior-models. No new datasets were generated for this study. The COVID-19 mortality data analyzed in this study are publicly available from the COVID-19 Data Repository maintained by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University: https://github.com/CSSEGISandData/COVID-19. Funding: Support for this research was provided to BP by an AMS-Simons Travel Grant. BP and MS were supported by cooperative agreement CDC-RFA-FT-23-0069 from the Centers for Disease Control and Prevention, Center for Forecasting and Outbreak Analytics. BP and MS received salary support from cooperative agreement CDC-RFA-FT-23-0069. ML was supported by a PhD studentship from Northeastern University London. The views expressed in this publication do not necessarily reflect the official policies of the Department of Health and Human Services/Centers for Disease Control and Prevention. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist. 1. Introduction Behavioral responses to epidemic threats fundamentally shape disease transmission. During the 1918 influenza pandemic, communities closed schools, banned public gatherings, and adopted mask-wearing to reduce transmission [1,2]. The HIV/AIDS epidemic saw widespread adoption of condom use, averting an estimated 100 million infections between 1990 and 2019 [3]. More recently, the COVID-19 pandemic triggered rapid behavioral shifts, including masking and social distancing, with mask mandates alone associated with significant reductions in cases and deaths [4,5]. Recent survey data from the United States suggests that human behavior, which varied spatially and temporally during COVID-19, was strongly associated with changes in COVID-19 mortality and hospitalization [6], indicating that observed disease severity indicators drive behavioral adaptation. These examples demonstrate that human populations do not passively experience epidemics—they adapt their behavior in response to perceived risk, fundamentally altering transmission dynamics. Compartmental models based on Susceptible-Exposed-Infected-Recovered (SEIR) frameworks have been widely used to understand epidemic trajectories and guide public health interventions [7]. Despite historical and contemporary evidence showing that people actively adapt their behavior in response to disease threats, most models assume populations passively experience outbreaks by maintaining constant contact patterns. Although the assumption of constant transmission rates is mathematically convenient, it ignores the dynamic feedback between disease burden and behavioral change. Recent modeling efforts have begun incorporating behavioral responses through various mechanisms, including prevalence-dependent [8–11], hospitalization-driven [12], and mortality-driven feedback [8,13–15]. LeJeune et al. [16] classify these approaches as either endogenous (behavior changes as a function of disease location variables within the model, creating a feedback loop) or exogenous (behavior changes through fixed or time-varying parameters informed by external data sources such as government policies and mobility patterns, or through modelers’ assumptions, rather than a location variable within a model itself). The critical question is whether systematic biases arise when models fail to account for behavioral change. If populations adapt their behavior in response to disease burden—as historical and contemporary evidence suggests—then models assuming constant transmission may systematically mischaracterize epidemic dynamics. Such mischaracterization could have serious consequences for public health decision-making: underestimating pathogen transmissibility might lead to insufficient interventions, while overestimating infection burden could result in misallocated resources. Despite these stakes, few studies have systematically investigated—either numerically or analytically—whether and to what extent omitting behavioral feedback introduces biases in epidemic models. Recent work has demonstrated that incorporating behavioral responses improves model fit to observed data [8,12,14,17,18], with behavioral models capturing epidemic trajectories better than their behavior-free counterparts. Behavioral models have also shown superior performance in out-of-sample validation [8,12]. Theoretical studies have examined how behavioral mechanisms—including vaccination behavior [19], awareness-driven protection [20], and incidence-dependent contact reduction [9]—alter epidemic outcomes such as cumulative infection burden; notably, models ignoring behavior change have systematically overshot final size projections [21]. Furthermore, using a 10-dimensional compartmental model fitted to a single US location, Pant et al. [8] show that models ignoring behavioral change underestimate the basic reproduction number relative to models that explicitly incorporate it; whether this finding generalizes across locations and simpler model structures motivates the present study. However, systematic comparison of key epidemiological quantities—such as basic reproduction number and final size—estimated through models with versus without behavioral feedback, using empirical data across many locations, remains lacking. We address this gap by systematically comparing compartmental models with and without behavioral feedback, fitted to COVID-19 mortality data from 20 US locations (19 U.S. states and the District of Columbia) during the first pandemic wave (March–July 2020). We implement three behavioral models where the transmission rate decreases as mortality increases, alongside a baseline model with constant transmission. Using Bayesian methods, we fit each model to observed mortality trajectories, estimate posterior distributions for the basic reproduction number and final epidemic size, and compare model performance both qualitatively and quantitatively. To test whether biases arise from model misspecification rather than data limitations, we conduct controlled synthetic experiments where true parameter values are known by design. We also analytically derive the relationship between final epidemic size in baseline versus behavioral models when the basic reproduction number is identical. Our analysis reveals three key findings with important public health implications. First, behavioral models consistently outperform baseline models across nearly all locations examined, as evidenced by both qualitative and quantitative comparisons — a finding consistent with existing works in the literature [8,12,14,18]. Second, baseline models systematically underestimate the basic reproduction number while counterintuitively overestimating the total proportion of the population infected compared to all three behavioral model formulations. We confirm this dual bias through synthetic experiments, demonstrating that these biases emerge when behavioral feedback is present in the data but models do not explicitly account for behavioral change — a finding that, to our knowledge, has not been previously demonstrated by fitting to empirical data across multiple locations or confirmed through controlled synthetic experiments. Third, for the class of models considered in this study, we prove analytically that for a fixed basic reproduction number, models that ignore behavior change overestimate the final epidemic size in comparison to model that explicitly accounts for behavior change — extending prior theoretical results for specific SIR formulations [21]. Overall, these findings demonstrate that incorporating behavioral dynamics is not merely a modeling refinement but essential for accurate epidemic parameter estimation and evidence-based pandemic response. This paper is organized as follows. In Section 2, we formulate the baseline model that does not explicitly account for behavioral change and three behavioral model variants, describe the COVID-19 mortality data used for calibration, and outline our Bayesian parameter inference approach. Section 3 presents our three key findings: behavioral models’ superior fit to observed data, systematic underestimation of the basic reproduction number by baseline models, and overestimation of infection burden when behavioral responses are ignored. We support these findings through both real data analysis across 20 US locations and controlled synthetic experiments. Section 4 discusses the public health implications and limitations of our work. 2. Methods 2.1 Model formulation We developed a compartmental SEIRD model incorporating behavioral responses to disease-induced mortality. The population is stratified into susceptible (S), exposed (E), infectious (I), recovered (R), and deceased (D) compartments, with total living population N(t) = S(t) + E(t) + I(t) + R(t) at time t. The model dynamics are governed by the following system of differential equations: (2.1)where is the incubation rate, is the recovery rate, and is the disease-induced mortality rate. Under homogeneous mixing, the transmission rate , where c(t) is the per-capita contact rate and p (assumed to be constant) is the probability of transmission per contact. Behavioral adaptation that reduces contacts thus enters the model directly as a reduction in . The transmission rate varies between model variants to capture different behavioral assumptions, whereby the severity of the outbreak is assumed to lead to behavioral adaptation aimed at reducing infection risk. Baseline model. The baseline model assumes constant transmission , representing traditional SEIRD dynamics without explicit behavioral feedback. Behavioral models. Behavioral models incorporate dynamic feedback through mortality-dependent transmission, with declining as deaths increase to reflect behavioral adaptation to perceived risk. The transmission rate depends on three quantities: , the baseline transmission rate in the absence of behavioral response; , the instantaneous death rate due to disease; and , which quantifies the population’s sensitivity to mortality. Larger values of correspond to stronger behavioral responses. We consider three functional forms of : - (a) Exponential form: where transmission rate decays exponentially with increasing mortality. - (b) Rational form: which produces similar behavior to the exponential form but saturates more slowly as deaths rise. - (c) Mixed form: which incorporates both exponential and rational forms. This hybrid form saturates faster than either individual form as mortality increases. All behavioral models reduce to the baseline model when . Hereafter, we refer to these as the Behavioral (Exponential), Behavioral (Rational), and Behavioral (Mixed) models. Fig 1 illustrates how the three functional forms respond to varying values of . Across all forms, increasing produces lower peak mortality and a more prolonged post-peak decline, while undergoes a sharp reduction during the mortality surge before gradually recovering as deaths subside. The forms differ in the degree and persistence of transmission suppression: as mortality rises and falls over the course of the outbreak, the mixed form attains the greatest reduction in most rapidly and recovers most slowly toward , while the rational form produces the weakest suppression and recovers fastest. The exponential form is intermediate in both respects. (a) Exponential form: , where transmission rate decays exponentially with increasing infections. (b) Rational form: , which saturates more slowly as infections rise. (c) Mixed form: , combining both exponential and rational responses. Panels (d)–(f) show the time evolution of corresponding to panels (a)–(c), respectively. The baseline model (, black line) represents traditional SEIRD dynamics without behavioral adaptation. Stronger behavioral responses (higher values) lead to lower peak mortality, more prolonged epidemic tails, and greater sustained suppression of . Model parameters: day-1, day-1, day-1, day-1, N = 106. Daily mortality. While individuals change behavior in response to instantaneous mortality , we compare model predictions to observed daily mortality data. Model-simulated daily mortality on day t is defined as: (2.5)representing the daily increment in cumulative deaths. The time series is compared to observed mortality data to estimate model parameters (Section 2.3). Basic reproduction number. The basic reproduction number, defined as the average number of secondary infections produced by a typical infectious individual in a wholly susceptible population (i.e., ), is given by: (2.6)Effective reproduction number. The effective reproduction number, denoted by , is the average number of new infections generated by a typical infected individual at time t [22]. The effective reproduction numbers for the baseline model, , and behavioral model, , are given respectively by: (2.7)Final epidemic size. The final epidemic size, denoted A, represents the total proportion of the population that becomes infected over the course of an outbreak: (2.8)where and denote the cumulative recovered and deceased individuals as . Numerically, we compute A from simulated SEIRD trajectories by integrating until the infectious compartment falls below a threshold (I(t) < 1) and calculating , where tend is the time at which the epidemic effectively terminates. 2.2 Data and location selection We calibrated the models to COVID-19 mortality data from the COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University [23]. Our analysis focused on the first wave of the pandemic in the United States, covering the 123-day period from March 1, 2020 to July 1, 2020. This period is chosen because it preceded the emergence of major SARS-CoV-2 variants of concern [24] and the availability of vaccines and standardized antiviral treatments, with transmission dynamics driven primarily by the ancestral strain in a largely immunologically naive population, and because non-pharmaceutical interventions (NPIs) — including stay-at-home orders, school closures, and mask mandates — represented the primary public health response, making behavioral adaptation the dominant mechanism modulating transmission [6,7]. Location selection criteria. From the 51 locations, we retained 20 for analysis by applying two exclusion criteria to the 7-day moving average of daily reported deaths, , for each location . The criteria were designed to ensure that each included trajectory carried a sufficient mortality signal for reliable model fitting and captured a sufficiently complete first mortality wave within the observation window. Here, a mortality wave refers to a mortality trajectory with an identifiable growth phase, a well-defined mortality peak, and a subsequent temporal decline within the observation window. - Insufficient mortality signal. Because reliable fitting of a deterministic model requires a sufficient mortality signal, we excluded locations whose total smoothed mortality signal over the analysis window, , was below 50 deaths, or whose peak smoothed daily mortality, , did not reach 10 deaths per day. - Incomplete mortality wave. Because good parameter estimation requires the mortality trajectory to carry information about both epidemic growth and decline, we excluded locations for which the first mortality wave was not sufficiently resolved within the observation window (i.e., 123 days). Specifically, locations were excluded if the mortality peak occurred fewer than 21 days before the end of the observation window (this allows us to exclude locations where mortality is still growing during the observation window), or if the mean smoothed daily mortality over the final 14 days exceeded 40% of the peak (this allows us to exclude locations where mortality plateaus after peak). Colorado (CO) was additionally excluded because of a reporting discontinuity in the cumulative mortality data that produced an artificial mortality spike. The 20 retained locations are: CT, DC, DE, IA, IL, IN, LA, MA, MD, MI, MN, NJ, NM, NY, OH, PA, RI, VA, WA, and WI. Details on the 31 excluded locations and their mortality trajectories are provided in Fig A1 and Table A2 in S1 Appendix. Full location names and standard two-letter abbreviations for all 51 locations are given in Table A1 in S1 Appendix. As a sensitivity analysis, we also report model-comparison results across all 51 locations, including those excluded from the main analysis, for the mixed behavioral model in Table A3 in S1 Appendix; analogous sensitivity tables for the exponential and rational behavioral models are provided in Tables A4 and A5 in S1 Appendix. These sensitivity results are included to show how the main qualitative conclusions vary under relaxed inclusion criteria, but estimates from excluded locations should be interpreted cautiously because those trajectories fail one or more data screening criteria. 2.3 Parameter inference We infer parameters using Approximate Bayesian Computation with Sequential Monte Carlo (ABC–SMC) [25–27], as implemented in the pyABC library [28], with the deterministic SEIRD model (2.1). The complete codebase for simulation, inference, and figure generation is available at https://github.com/markolalovic/epi-behavior-models. Fixed and inferred parameters are summarized in Table 1; additional implementation details, including ABC–SMC stopping criteria, perturbation kernels, and tolerance schedules, are provided in Section A3 in S1 Appendix. For each location and model, we simulate daily deaths over the 123-day window and apply the same 7-day moving average used for the observed data. We reparameterize in epidemiological terms, inferring , the initial prevalence , the disease-induced mortality rate , and the behavioral response parameter . (We note that is not directly inferred but can be calculated from the inferred based on (2.6).) The prior on is chosen to be weakly informative: The lower bound excludes subcritical or nearly subcritical regimes that are inconsistent with the selected mortality trajectories, all of which exhibit a clear first mortality wave. The upper bound is deliberately conservative, allowing for uncertainty in early SARS-CoV-2 transmissibility and for the fact that, in behavioral models, represents intrinsic transmissibility before behavioral adaptation reduces contacts. The priors on the initial prevalence and mortality rate are log-uniform, because both parameters are positive, expected to be small, and uncertain over several orders of magnitude. For the behavioral response strength we use Parameter is phenomenological rather than directly measured. Its scale is determined by the dimensionless product , where is the instantaneous disease-induced mortality rate. Thus controls the daily-mortality scale over which transmission is reduced. The prior includes , corresponding to no behavioral feedback, and extends to values producing substantial reductions in transmission at mortality levels observed in the retained locations. For example, under the exponential form, gives transmission multipliers , , and at mortality rates of 10, 50, and 100 deaths per day, respectively. Therefore, the prior spans negligible, moderate, and strong behavioral responses. To calibrate model-simulated to observed mortality trajectories, we use the Normalized Sum of Squared Errors (NSSE): (2.9)where is the model-simulated 7-day daily moving average mortality and is the observed 7-day moving average mortality over the T = 123 day observation window. This normalization ensures the distance is scale-invariant across locations with different mortality burdens. ABC–SMC uses 1000 particles per generation, a multivariate normal transition kernel adapted from the previous population, and a quantile-based tolerance schedule over up to 10 generations. Additional details are provided in Section A3 in S1 Appendix. 2.4 Synthetic data experiments To rigorously test whether behavioral effects can be recovered from mortality data and to quantify systematic biases in parameter estimation, we conducted controlled synthetic data experiments. We generated mortality trajectories from behavioral models using parameters calibrated to Massachusetts COVID-19 data, with varying values of the behavioral parameter . To approximate observational variability, we added small independent Gaussian noise to each simulated series: Negative values were truncated to zero before applying the 7-day moving average. We then fitted both the data-generating behavioral model and the baseline model to these synthetic mortality datasets. This approach allows us to assess: (i) whether the behavioral model correctly recovers the known underlying dynamics; (ii) whether the baseline model produces systematic biases when behavioral feedback is present but ignored; and (iii) how the magnitude of behavioral response (controlled by varying ) affects the ability to distinguish between model structures. Note that it is expected by construction that the behavioral model recovers the synthetic data better than the baseline; the purpose of these experiments is not to demonstrate this, but rather to confirm that systematic biases in and final epidemic size arise specifically from model misspecification. By knowing the true data-generating mechanism, we can definitively attribute any systematic biases in and final epidemic size observed in the empirical data to model misspecification, rather than to data noise (e.g., stochastic fluctuations in daily death counts), reporting irregularities (e.g., batch corrections, delayed notifications, or inconsistent cause-of-death classification across states), or other surveillance heterogeneities inherent to COVID-19 mortality data. 2.5 Model comparison We evaluate the relative performance of baseline and behavioral models using Bayesian model selection via ABC–SMC. We describe below the distance metric and Bayes factor calculation used for model comparison. Bayes factor. The Bayes factor is a measure of the relative evidence provided by the data in favor of one model over another, defined as the ratio of their marginal likelihoods [33]. Under equal model priors, it reduces to the posterior odds ratio. We performed ABC–SMC model selection [25–27] comparing the baseline model against each of the three behavioral model formulations (exponential, rational, and mixed) separately, assigning equal prior probabilities to each pair. To avoid extinction artifacts, we report results at the first stabilized generation , defined as the earliest g for which both model probabilities change by less than 0.02 over 3 consecutive generations and . At the stabilized generation , we summarize the evidence for a behavioral model via the Bayes factor under equal model priors: We interpret BFeq following the scale from [33]. Values of BFeq > 1 indicate evidence favoring the behavioral model over the baseline model, with strength categorized as: 1–3 (very weak), 3–20 (positive), 20–150 (strong), and > 150 (very strong). Values < 1 indicate evidence favoring the baseline model. 3. Results We fitted COVID-19 mortality data from 20 US locations (19 U.S. states and the District of Columbia) during the first pandemic wave (March–July 2020) to a baseline SEIRD model with constant transmission rate and three behavioral SEIRD models where transmission rate declines in response to disease-induced mortality. The 20 retained locations were selected using pre-specified data-quality criteria designed to isolate clear first-wave mortality trajectories suitable for single-wave inference (see Section 2 for detailed selection criteria). Our analysis reveals three key findings: (i) behavioral models provide significantly improved fits to observed mortality data across almost all locations examined, with both median normalized squared error and Bayes factors favoring behavioral models in 18 of 20 locations; (ii) baseline models that ignore behavioral change systematically underestimate the basic reproduction number (), with median from behavioral models consistently higher than from baseline models across all locations; and (iii) baseline models overestimate total infection burden despite underestimating transmissibility, creating a paradoxical dual bias with serious policy implications. 3.1 Behavioral models show superior fit to COVID-19 mortality data All three behavioral models (exponential, rational, and mixed) consistently outperformed the baseline model when fitted to COVID-19 mortality data from 20 locations during the first pandemic wave. Visually, the behavioral models capture the characteristic asymmetric mortality curves observed in many locations—sharp initial rises followed by gradual, extended declines—while the baseline model produces curves that decay too rapidly (Fig 2 shows the mixed behavioral model; Figs A4 and A5 in S1 Appendix show exponential and rational formulations, respectively). Posterior predictive fits comparing baseline and behavioral (mixed) models to 7-day averaged daily COVID-19 deaths during the first pandemic wave (March–July 2020). Black points show observed mortality data. Solid lines show posterior predictive medians and shaded regions show 90% credible intervals (5th–95th percentiles), computed from 1,000 weighted posterior samples drawn from the final ABC–SMC generation. The behavioral model accurately captures the characteristic asymmetric mortality curves with sharp initial rises and gradual extended declines observed across most locations. Locations are ordered from left to right and top to bottom by decreasing improvement in fit: those where the behavioral model outperforms the baseline model most dramatically (largest difference in median NSSE) appear first, progressing to locations where the two models perform most similarly. This ordering visually emphasizes that behavioral feedback provides substantial improvements in model fit across the vast majority of locations examined. This superior performance is confirmed quantitatively across all three formulations. For the mixed behavioral model, we observed lower median normalized sum of squared errors (NSSE) in 18 of 20 locations, with Bayesian model selection via ABC–SMC providing evidence for the mixed behavioral model in these same 18 locations (Table 2). Only Iowa and New York favored the baseline model. The exponential and rational formulations of behavioral models demonstrated comparable performance (see Tables A6 and A7 in S1 Appendix). Overall, across all three behavioral formulations, at least 18 of 20 locations favored behavioral models over baseline models based on NSSE and Bayes factor. Among the three behavioral formulations, the mixed model most consistently outperformed the baseline, particularly in locations with sharp, high-intensity mortality peaks such as DC. This likely reflects the mixed form’s faster saturation at high mortality levels (as depicted in Fig 1) - a property absent in the exponential and rational forms individually. In locations with more gradual mortality trajectories, the three forms produced comparable performance, suggesting that mortality signal strength mediates which functional form is most identifiable from data. The superior performance of behavioral models varies systematically across locations, with the largest improvements in locations exhibiting extended mortality decline phases. Fig 2 orders locations by decreasing difference in median NSSE between baseline and behavioral models, with locations showing the largest improvement appearing first. Locations like Washington and Indiana (top left) exhibit the most striking contrast—the behavioral models capture gradual, extended mortality declines that baseline models fail to reproduce. Moving toward Rhode Island and Delaware (bottom right), the performance gap narrows, with both model types achieving similar fit quality. Sensitivity analysis reveals a key feature of behavioral models: increasing the behavioral response parameter () produces lower peak mortality and slower decline in daily deaths following the peak (Fig 1), a mechanism absent in baseline models. This pattern is reminiscent of the “flattening the curve” phenomenon, where public health interventions reduce peak mortality while extending the duration of the outbreak. In all three behavioral formulations, exhibits a characteristic pattern: sharp decline during the mortality surge followed by gradual recovery as deaths decline (see Section A6 in S1 Appendix). This time-varying transmission rate better reflects the reality of behavioral adaptation during epidemics, consequently producing superior fits to observed mortality data compared to the constant transmission rate assumption of the baseline model. We validated this finding through controlled synthetic data experiments where mortality trajectories were generated from behavioral models with varying levels of behavioral sensitivity (), then fitted to both the data-generating behavioral model and the baseline model. Using the mixed behavioral formulation, behavioral models accurately recovered the true dynamics across all values, while the baseline model increasingly failed to capture mortality trajectories as behavioral feedback strengthened (Fig 4A). Results were consistent across all three functional forms examined (Figs A2A and A3A show exponential and rational forms in S1 Appendix, respectively), demonstrating that baseline models systematically mischaracterize epidemic trajectories when behavioral responses are present. 3.2 Baseline models systematically underestimate the basic reproduction number Baseline models that ignore behavioral change systematically underestimate the basic reproduction number () compared to behavioral models when fitted to the same COVID-19 mortality data across the 20 selected US locations (Fig 3A shows the comparison with the mixed formulation). For the behavioral mixed model, the baseline model inferred a lower median in all 20 retained locations. The cross-location difference, defined as , had median 1.26, interquartile range 0.91–1.97, and range 0.15–2.67. Furthermore, estimates from all three behavioral model variants were similar to each other, and the systematic underestimation by the baseline model was consistent across the exponential, rational, and mixed formulations (Figs A6 and A7 in S1 Appendix show comparisons for the exponential and rational forms, respectively). Posterior distributions of key epidemiological quantities estimated by fitting the baseline and behavioral (mixed) models to COVID-19 mortality data from 20 US locations during the first pandemic wave (March–July 2020). (A) Basic reproduction number . The behavioral model consistently estimates higher values than the baseline model across all 20 locations, indicating that models ignoring behavioral feedback systematically underestimate intrinsic transmissibility. (B) Final epidemic size (proportion of population infected). Despite estimating lower , the baseline model predicts larger cumulative infection burdens than the behavioral model in all 20 locations, demonstrating the paradoxical consequence of ignoring behavioral adaptation. Boxplots show median (center line), interquartile range (box), and whiskers extending to within the sample support computed from weighted posterior samples. Locations are ordered from left to right by decreasing difference in median NSSE between the baseline and behavioral (Mixed) models. The effective reproduction number for all four models is provided in Section A7 in S1 Appendix. In all cases, converges to values below 1 by late April–May across all locations (coinciding with peak infection), though behavioral models exhibit steeper initial declines corresponding to rapid behavioral responses during the early stages of the outbreak. Synthetic experiments where true was known by design confirmed this systematic underestimation. Baseline models consistently underestimated true across all behavioral sensitivity values, with bias increasing as behavioral feedback strengthened (increasing ), while behavioral models accurately recovered true across all values (Fig 4B for mixed; Figs A2B and A3B for exponential and rational forms in S1 Appendix, respectively). This confirms that underestimation is a systematic consequence of model misspecification. Synthetic data experiments where mortality trajectories were generated using the behavioral (mixed) model with known parameters and varying behavioral sensitivity , then fitted with both baseline and behavioral models. (A) Posterior predictive fits to synthetic mortality data (black points). Solid lines show posterior medians and shaded regions show 90% credible intervals. The behavioral (mixed) model accurately recovers the true dynamics across all values, while the baseline model increasingly fails to capture the extended post-peak decline as behavioral feedback strengthens with increasing values. (B) Posterior distributions of the basic reproduction number . The behavioral model correctly recovers the true (red dashed line) across all scenarios, while the baseline model systematically underestimates , with bias increasing as increases. (C) Posterior distributions of final epidemic size. The behavioral model accurately estimates the true final size (red dashed line), while the baseline model systematically overestimates infection burden despite underestimating transmissibility. Intuitively, this systematic underestimation of arises because baseline and behavioral models, although fitted to the same observed mortality data, assume fundamentally different transmission mechanisms. The baseline model uses a constant transmission rate throughout the epidemic, fitting the data with a moderate value that remains fixed. In contrast, behavioral models assume a time-varying transmission rate that declines as people adapt their behavior in response to rising mortality. To reproduce the same initial epidemic growth observed in the data, behavioral models require a substantially higher initial transmission rate —which then declines as behavioral responses take effect. As evident from Eq 2.6, is an increasing function of ; hence models that estimate higher values of necessarily yield higher estimates of . Thus, behavioral models, which require higher to account for subsequent behavioral adaptation, naturally produce higher estimates than baseline models. This higher initial transmission rate captured by behavioral models represents the intrinsic transmissibility of the pathogen before any behavioral adaptation occurs. The baseline model, lacking this behavioral feedback mechanism, systematically underestimates this initial transmissibility because it interprets the entire epidemic trajectory—including periods when behavior has already adapted—through the lens of a single, constant transmission rate. 3.3 Ignoring behavioral responses overestimates infection burden Despite underestimating , baseline models overestimate the total proportion of the population that becomes infected (Fig 3B shows the comparison with the mixed formulation; Figs A6B and A7B in S1 Appendix show comparisons with exponential and rational forms, respectively). For the behavioral mixed model, the baseline model produced a larger median final-size estimate in all 20 retained locations. The cross-location overestimation, defined as , had median 2.96 percentage points, interquartile range 1.58–6.22 percentage points, and range 0.64–7.87 percentage points. This overestimation was largest in WI, where the baseline model predicted a median final epidemic size that was 7.87 percentage points higher than the behavioral mixed model. Synthetic experiments confirmed this systematic overestimation of final epidemic size. Baseline models consistently overestimated true final epidemic sizes while behavioral models accurately recovered them across all values considered (Fig 4C shows comparison with mixed; Figs A2C and A3C show comparisons with exponential and rational forms in S1 Appendix, respectively). The bias in overestimating final epidemic size increases with stronger behavioral responses, demonstrating that ignoring behavioral dynamics leads to inflated infection burden estimates regardless of functional form. To gain theoretical insight into the relationship between behavioral feedback and epidemic outcomes, we analytically examine final epidemic size under a controlled scenario where is identical for both the baseline and a behavioral model. For analytical tractability, we assume a constant total population N = S + E + I + R + D, which differs from the time-varying living population used in our numerical simulations. We prove that for any fixed value of , behavioral models incorporating mortality-driven feedback produce final epidemic sizes less than or equal to those from baseline models. We first demonstrate this result for the exponential behavioral model. Theorem 3.1 (Final Size Reduction: Exponential Form). Consider an epidemic model with behavioral response function where . Furthermore, we assume that the total population size is fixed with N = S + E + I + R + D. Let A0 denote the final size of the baseline model (i.e., ). The final size of the behavioral model, , can be expressed as: where a1 < 0 is the first-order correction to the final epidemic size due to behavioral feedback. Hence, for small , the final size of the behavioral model is smaller than that of the baseline model. The proof of Theorem 3.1, based on asymptotic expansion, is given in Appendix A. This result extends beyond the exponential form. For any behavioral response function satisfying —including the exponential, rational, and mixed forms considered in this study—we show that behavioral models always produce smaller or equal final epidemic sizes compared to baseline models when is the same for both models. Theorem 3.2 (Final Size Reduction: General Form). Consider the SEIRD model (2.1) where N = S + E + I + R + D is the fixed total population size and the transmission rate decreases with increasing mortality through a behavioral response function satisfying , with being the maximum transmission rate and . If a behavioral model (with ) and the baseline model (with and constant transmission ) have equal basic reproduction numbers, given by , then the final epidemic size of the baseline model (A0) is always larger than or equal to that of the behavioral model (). The proof of Theorem 3.2 is given in Appendix B. We obtain final size relationships for the baseline model with constant transmission rate and for behavioral models with . By defining an auxiliary function based on these final size equations and exploiting its properties, including its unique maximum and concavity, we show that the final size of the behavioral model is always less than or equal to that of the baseline model when both share the same basic reproduction number . 4. Discussion Our analysis of COVID-19 mortality data from 20 US locations during the first pandemic wave reveals systematic biases that arise when epidemic models fail to account for behavioral responses to disease burden. Across all three behavioral models, at least 18 of 20 locations favored behavioral models—where transmission rates decline as mortality increases—over the traditional model with a constant transmission rate, based on both NSSE and Bayes factor. More critically, we identified a problematic paradox: models that ignore behavioral change systematically underestimate pathogen transmissibility while simultaneously overestimating total infection burden. We confirmed these biases through controlled synthetic experiments where the true dynamics were known by design, demonstrating that the discrepancies arise from model misspecification rather than data quality or stochastic effects. These findings held across three distinct mathematical formulations of behavioral models, strengthening confidence in our conclusions. The dual bias we identify has potential implications for public health decision-making. Models that ignore behavioral responses underestimate transmissibility as measured by the basic reproduction number, which may contribute to delayed or insufficiently aggressive early interventions. Conversely, overestimation of total infection burden could, in principle, contribute to excessive resource allocation for later epidemic phases that behavioral adaptation may prevent. While we do not formally evaluate allocation strategies or costs, these considerations suggest that the structural biases identified here are worth accounting for when translating model projections into policy — a direction that warrants dedicated future investigation. The mechanism underlying these biases is intuitive once recognized. Baseline models with constant transmission rates fit observed mortality data using relatively modest transmission parameters because the transmission rate remains fixed throughout the inferred mortality wave. In contrast, behavioral models require substantially higher initial transmission rates to fit the same mortality data precisely because transmission declines as mortality increases and people adapt their behavior. This difference directly manifests in the basic reproduction number: the baseline model settles on a lower value that suffices under constant transmission, while the behavioral model correctly identifies the higher intrinsic transmissibility that characterized the pathogen before behavioral responses took effect. The bias in final epidemic size follows naturally from this same mechanism. When populations reduce contact rates in response to rising deaths—through actions like avoiding crowds, wearing masks, or increasing hygiene—fewer people ultimately become infected compared to scenarios where behavior remains unchanged. This study has several important limitations. First, our models assume that behavioral changes are driven exclusively by location-specific reported mortality, whereas behavior in reality may also respond to mortality signals from nearby regions, national trends, policy announcements, or media coverage [6]. In addition, reported COVID-19 mortality may vary across locations because of differences in notification timing, testing availability, and cause-of-death classification, including the possibility of unrecognized COVID-19 deaths recorded under other causes [34]. Second, using mortality as a proxy for public risk perception does not capture the complex array of influences on behavior, including policy actions, misinformation, informal information channels, and individual risk tolerance [15]. This simplification may also overlook delays between changes in actual risk and behavioral adaptation, and neglects the possibility that current behavior may be influenced by past epidemic conditions (such as cumulative mortality over a preceding time window), neither of which our instantaneous mortality-driven feedback formulation captures [35,36]. Compartmental frameworks that explicitly encode richer behavioral determinants — including perceived susceptibility, perceived severity, fatigue — offer more mechanistic representations of behavioral adaptation [8,37], and comparing their performance against simpler phenomenological feedback formulations such as ours represents a valuable direction for future work. Third, we assume homogeneous populations within each location, ignoring substantial variation in risk tolerance, contact patterns, and protective behavior adoption across age groups, socioeconomic strata, and communities. Real populations exhibit heterogeneous responses to epidemic threats, with some groups adopting protective measures rapidly while others respond more slowly or not at all [6,35,36,38]. Fourth, our analysis focuses exclusively on the first pandemic wave (March–July 2020), when behavioral responses were likely strongest due to novelty and uncertainty. Behavioral patterns may differ substantially in later waves as pandemic fatigue sets in or as populations develop more nuanced risk assessments [6]. We note further that behavioral suppression of an initial wave may leave a larger susceptible population, potentially contributing to epidemic resurgence in subsequent waves — a dynamic not captured by our single-wave framework. Fifth, our models imply that extended post-peak mortality declines reflect stronger behavioral responses, but this relationship does not hold universally. New York experienced one of the earliest and most severe outbreaks [39–41] with stringent public health measures and well-documented behavioral changes [6,42,43], yet the mortality curve rose and fell relatively symmetrically without the prolonged tail our models associate with strong behavioral feedback. This disconnect suggests that the temporal signature of behavioral adaptation depends on factors we do not explicitly capture, such as the timing and stringency of policy interventions, underlying contact network structures, population density, or the rapidity with which protective behaviors were adopted. More broadly, our mortality-driven behavioral feedback implicitly aggregates multiple drivers of behavior change, including voluntary risk-aversion and policy-induced compliance, and cannot distinguish between these mechanisms at the aggregate level. Locations with sharp mortality declines might reflect either weak behavioral responses or, conversely, such rapid behavioral adaptation that transmission was quickly suppressed. Our models cannot distinguish between these possibilities, which limits our ability to infer behavioral mechanisms from mortality data alone. Finally, our behavioral models assume a specific functional relationship between mortality and transmission reduction. While we tested three different functional forms, real behavioral responses may follow more complex patterns involving thresholds, saturation effects, or time-varying sensitivity to disease burden. Relatedly, this study does not consider comparisons between baseline and behavioral models where the two differ in complexity by more than a single additional parameter. This is because having even a single additional parameter allows the model an additional degree of freedom to fit the data; hence, it is not unexpected that the behavioral model typically fits the data better than the baseline model, which can be obtained as a special case by turning off the behavioral response (i.e., setting ). Furthermore, even within the class of models considered here, the sensitivity of our qualitative and quantitative findings to the fixed epidemiological parameters remains unexplored. Because this study focuses on COVID-19, for which the incubation and infectious periods are relatively well characterized, we fixed both parameters at literature-informed values. Testing the sensitivity of our qualitative and quantitative findings to alternative values of the incubation and infectious periods is therefore a natural direction for future work, alongside testing sensitivity to different functional forms of behavioral feedback and to models of varying structural complexity. Future work should address these limitations through several avenues. Incorporating spatial coupling between locations could capture how mortality in nearby regions influences local behavior [44,45]. Integrating data on policy interventions, mobility patterns, and direct behavioral surveys could provide more nuanced representations of behavioral change beyond mortality-driven responses. In particular, incorporating data on policy stringency — for example, from the Oxford COVID-19 Government Response Tracker [46] — alongside mortality signals could help disentangle the contributions of voluntary risk-aversion and policy-induced compliance to the behavioral signal recovered from empirical data. While age stratification is absent from our framework — with all variability in disease severity absorbed into the fitted parameter , this is a deliberate modeling choice: the homogeneous mixing assumption affords the mathematical tractability required to derive the analytical results of Theorems 3.1 and 3.2, and ensures that observed biases can be cleanly attributed to behavioral misspecification rather than confounded by heterogeneous age-dependent mortality, contact patterns, or behavioral responses. Age-structured models could account for heterogeneous contact patterns, risk perceptions, recovery rates, and disease severity across demographic groups [22]. As such, future studies should explore the effects of including an age-dependent infection fatality ratio (IFR) inference strategy as it has been well-documented that COVID-19 infection fatality ratios are known to depend heavily on age. Multi-wave analyses could reveal how behavioral responses evolve over the course of a pandemic as novelty effects diminish and risk perceptions stabilize. Comparing fits to multiple data streams—cases, hospitalizations, and deaths—could help disentangle behavioral responses to different epidemic indicators and their respective time lags. More fundamentally, developing frameworks that explicitly model the cognitive and social processes underlying behavioral change—such as risk perception dynamics, information diffusion, and social norm formation—could provide deeper mechanistic understanding beyond phenomenological feedback functions [15,47–50]. Such frameworks could better predict behavioral responses to novel pathogens and inform more effective risk communication strategies. Finally, the infectious period (10 days) and incubation period (3 days) are held fixed and identical across all models; while the latter differs from some published estimates for ancestral SARS-CoV-2, this choice ensures that differences in model performance are attributable solely to the presence or absence of behavioral feedback rather than to differences in assumed epidemiological parameters, and a comprehensive sensitivity analysis with respect to these quantities — including the choice of distributional assumptions for infectious periods — remains an interesting direction for future work. Severing the coupling between pathogen transmission and human behavior produces models that systematically misrepresent reality. Infectious disease outbreaks unfold through feedback: transmission shapes behavior, and behavior alters transmission. Models that assume constant contact rates sacrifice this coupling for mathematical convenience, producing contradictory estimates of key epidemiological quantities—underestimating pathogen transmissibility while overestimating infection burden. These structural failures cascade through every decision that relies on model projections, from intervention timing to resource allocation. Our findings demonstrate that incorporating behavioral dynamics is not merely a modeling refinement but essential for accurate epidemic parameter estimation and evidence-based pandemic response. Across 20 geographically and demographically diverse US locations, behavioral models consistently provided superior fits to observed data and more plausible estimates of both transmissibility and final epidemic size. The analytical results further establish that this dual bias is an inevitable mathematical consequence of ignoring behavioral feedback when it is present in the data. Future pandemic preparedness must prioritize the development and deployment of models that explicitly account for the dynamic interplay between disease transmission and human behavioral responses. A Proof of Theorem 3.1 Proof. Let , where combines the behavioral sensitivity parameter with the mortality rate , and denotes the infectious population trajectory corresponding to . Furthermore, we assume that the total population size is fixed with N = S + E + I + R + D. The final epidemic size corresponding to is defined as (A.1)From the S–equation of model (2.1), we have Dividing both sides by S(t) and integrating in time from 0 to (noting S(t) > 0 for all finite time until extinction): (A.2)Since the removed individuals at infinity equal the initial population minus susceptibles at infinity, we have and hence, using , (A.2) can be rewritten as (A.3)This holds for all . When (i.e., ), we have with A0 and I0(t) denoting the final size and infectious population corresponding to the baseline model. Hence, the classical implicit final-size equation: (A.4)Small- expansion. We now expand for small . Let (A.5)where, a1 and J(t) are the first-order corrections to the final epidemic size and infectious trajectory, respectively. We expand both sides of (A.3) using the Taylor series. Right-hand side expansion. Using and from (A.5), we expand: (A.7)We expand the internal term in (A.3) using (A.7) as: Plugging in the expansion of and from (A.5) into (A.1) and comparing terms gives: Hence, the right-hand side of (A.3) can be finally expanded as: (A.8)Matching terms. Now we compare the matching terms of left-hand side expansion, given in (A.6), with right-hand side expansion, given in (A.8). At order O(1), we recover the classical relation (A.4). At order we have: which can be rewritten as: Therefore, solving for a1 gives: (A.9)Since , and A0 corresponds to the final size of the baseline model, to show that the final size of the behavioral model is smaller than that of the baseline model, it is sufficient to show that a1 < 0. Since we consider small , the terms are negligible compared to the first-order term , so the sign of a1 determines whether or . The numerator of (A.9) is always positive. The denominator of (A.9) is positive provided: which can be rewritten as: which is true since (see (B.9) for details). This concludes the proof. □ B Proof of Theorem 3.2 Proof. Let , where combines the behavioral sensitivity parameter with the mortality rate , and denotes the infectious population trajectory corresponding to . For example, . We require that , and that . We assume that the total population size is fixed with N = S + E + I + R + D. Furthermore, initial conditions are: S(0) = S0 > 0, , , R(0) = D(0) = 0 and . In this proof we work with the cumulative number infected, so the final epidemic size is defined as (B.1)The final-size proportion used in the main text is , so the inequality is unchanged after normalization. From the S–equation of model (2.1), we have Dividing both sides by S(t) and integrating in time from 0 to (noting S(t) > 0 for all finite time until extinction): (B.2)Since removed individuals at infinity equal the initial population minus susceptibles at infinity, we have and hence (B.2) can be rewritten as (B.3)Equation (B.3) is exact for any time-dependent and any solution of the SEIRD model (2.1) satisfying the stated assumptions. Let A0 denote the final size when . Recall the pointwise inequality assumption The above implies Substituting inequality of in (B.3) gives: (B.4)For the baseline case with (no behavioral response), we have and the infectious trajectory reduces to I0(t). Substituting in (B.3) yields: (B.5)Equations (B.4) and (B.5) cannot be directly compared because they involve both different final epidemic sizes on the left hand side ( versus A0) and different infectious population trajectories on the right hand side ( versus I0(t)). To enable comparison, we apply the relation from equation (B.1), which transforms inequality (B.4) into: (B.6)and for baseline case, equation (B.5) can rewritten as: (B.7)We define the function (B.8)Notice that for the baseline model f(A0) = 0 while for behavioral models . f(A) has two zeros, the trivial solution A = 0 (no epidemic) and the nontrivial solution A = A0 (baseline epidemic). The proof proceeds in three steps. First, we show f(A) is strictly concave and compute its unique maximum at . Second, since A0 is defined only implicitly through a transcendental equation (B.7), we cannot directly compare A0 with . We therefore normalize the susceptible population as s = S/N and express the baseline final size through the transcendental equation , where . By analyzing the auxiliary function , we show that it attains a unique maximum at . Furthermore, we prove that , which translates to . Third, combining with the concavity of f(A) and the inequality , we conclude . Details follow. Differentiating f(A) with respect to A gives (where is a constant): Thus f is strictly concave on [0,N) and can have at most one critical point (which will be a global maximum). Solving yields , which gives . Thus the unique critical point is where . Because , is the global maximizer of f on [0,N). Relation to A0. Introduce the normalized susceptible fraction s: = S/N. We use the standard approximation , consistent with the assumption of a nearly fully susceptible population at epidemic onset. It follows from (B.7) that the classical final-size equation for can be written in normalized form as Define the auxiliary function Then is a zero of , i.e., . Compute the derivative: Hence for , , and for . Thus attains its unique maximum at . Evaluate at the endpoints and at the maximum: and by the standard inequality (with equality only at x = 1), with strict inequality when . Thus the maximum value is nonnegative (positive for ) while . Because is continuous on (0,1], has a unique maximum at with , and , the equation has exactly two solutions in (0,1] when : one at s = 1 and the other strictly to the left of . The epidemiologically relevant solution is the nontrivial root . Because and , this nontrivial root must satisfy (B.9)Conclude . Translate back to A. Since , the inequality is equivalent to hence Therefore the unique maximizer of f(A), given by (B.8), indeed lies strictly to the left of the classical final size A0 (for , i.e., in the nontrivial epidemic regime). Since satisfies inequality (B.6), we conclude that must lie to the left of A0. Hence, . □ Supporting information S1 Appendix. Supplementary methods and results. https://doi.org/10.1371/journal.pcbi.1014746.s001 (PDF) References - 1. Luckingham B. To mask or not to mask: a note on the 1918 Spanish influenza epidemic in Tucson. J Ariz Hist. 1984;25(2):191–204. pmid:11635217 - 2. Dolan B. Unmasking history: Who was behind the anti-mask league protests during the 1918 influenza epidemic in San Francisco? Perspect Med Humanit. 2020;5(19). - 3. Stover J, Teng Y. The impact of condom use on the HIV epidemic. Gates Open Res. 2022;5:91. pmid:35355689 - 4. Adjodah D, Dinakar K, Chinazzi M, Fraiberger SP, Pentland A, Bates S, et al. Association between COVID-19 outcomes and mask mandates, adherence, and attitudes. PLoS One. 2021;16(6):e0252315. pmid:34161332 - 5. Kwon S, Joshi AD, Lo C-H, Drew DA, Nguyen LH, Guo C-G, et al. Association of social distancing and face mask use with risk of COVID-19. Nat Commun. 2021;12(1):3737. pmid:34145289 - 6. Urmi T, Pant B, Dewey G, Quintana-Mathe A, Lang I, Druckman J, et al. Characterizing population-level changes in human behavior during the COVID-19 pandemic in the United States. Proc Natl Acad Sci U S A. 2025;122(37):e2500655122. pmid:40932771 - 7. Ferguson NM, Laydon D, Nedjati-Gilani G, Imai N, Ainslie K, Baguelin M, et al. Report 9: Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand, vol. 16. London: Imperial College London; 2020. - 8. Pant B, Safdar S, Santillana M, Gumel AB. Mathematical assessment of the role of human behavior changes on SARS-CoV-2 transmission dynamics in the United States. Bull Math Biol. 2024;86(8):92. pmid:38888744 - 9. Del Valle S, Hethcote H, Hyman JM, Castillo-Chavez C. Effects of behavioral changes in a smallpox attack model. Math Biosci. 2005;195(2):228–51. pmid:15913667 - 10. Perra N, Balcan D, Gonçalves B, Vespignani A. Towards a characterization of behavior-disease models. PLoS One. 2011;6(8):e23084. pmid:21826228 - 11. Mahmud MS, Eshun S, Espinoza B, Kadelka C. Adaptive human behavior and delays in information availability autonomously modulate epidemic waves. PNAS Nexus. 2025;4(5):pgaf145. pmid:40432904 - 12. Oveson A, Girvan M, Gumel AB. Modeling the impact of hospitalization-induced behavioral changes on the spread of COVID-19 in New York City. Infect Dis Model. 2025;10(4):1055–92. pmid:40547859 - 13. Weitz JS, Park SW, Eksin C, Dushoff J. Awareness-driven behavior changes can shift the shape of epidemics away from peaks and toward plateaus, shoulders, and oscillations. Proc Natl Acad Sci U S A. 2020;117(51):32764–71. pmid:33262277 - 14. Gozzi N, Perra N, Vespignani A. Comparative evaluation of behavioral epidemic models using COVID-19 data. Proc Natl Acad Sci U S A. 2025;122(24):e2421993122. pmid:40504146 - 15. Jo Y, Sharbayta SS, Buonomo B. Behavioral change models for infectious disease transmission: a systematic review (2020-2025). arXiv:260216633 [Preprint]. 2026. - 16. LeJeune L, Ghaffarzadegan N, Childs LM, Saucedo O. Formulating human risk response in epidemic models: Exogenous vs endogenous approaches. Eur J Oper Res. 2025;324(1):246–58. - 17. He D, Dushoff J, Day T, Ma J, Earn DJD. Inferring the causes of the three waves of the 1918 influenza pandemic in England and Wales. Proc R Soc B: Biol Sci. 2013;280(1766):20131345. pmid:23843396 - 18. LeJeune L, Ghaffarzadegan N, Childs LM, Saucedo O. Mathematical analysis of simple behavioral epidemic models. Math Biosci. 2024;375:109250. pmid:39009074 - 19. Bauch CT, Earn DJD. Vaccination and the theory of games. Proc Natl Acad Sci U S A. 2004;101(36):13391–4. pmid:15329411 - 20. Funk S, Gilad E, Watkins C, Jansen VAA. The spread of awareness and its impact on epidemic outbreaks. Proc Natl Acad Sci U S A. 2009;106(16):6872–7. pmid:19332788 - 21. Eksin C, Paarporn K, Weitz JS. Systematic biases in disease forecasting - the role of behavior change. Epidemics. 2019;27:96–105. pmid:30922858 - 22. Pant B, Gumel AB. Mathematical assessment of the roles of age heterogeneity and vaccination on the dynamics and control of SARS-CoV-2. Infect Dis Model. 2024;9(3):828–74. pmid:38725431 - 23. Dong E, Du H, Gardner L. An interactive web-based dashboard to track COVID-19 in real time. Lancet Infect Dis. 2020;20(5):533–4. pmid:32087114 - 24. Korber B, Fischer WM, Gnanakaran S, Yoon H, Theiler J, Abfalterer W, et al. Tracking changes in SARS-CoV-2 spike: evidence that D614G increases infectivity of the COVID-19 virus. Cell. 2020;182(4):812-827.e19. pmid:32697968 - 25. Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MPH. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. J R Soc Interface. 2009;6(31):187–202. pmid:19205079 - 26. Sunnåker M, Busetto AG, Numminen E, Corander J, Foll M, Dessimoz C. Approximate Bayesian computation. PLoS Comput Biol. 2013;9(1):e1002803. pmid:23341757 - 27. Minter A, Retkute R. Approximate Bayesian Computation for infectious disease modelling. Epidemics. 2019;29:100368. pmid:31563466 - 28. Schälte Y, Klinger E, Alamoudi E, Hasenauer J. pyABC: Efficient and robust easy-to-use approximate Bayesian computation. JOSS. 2022;7(74):4304. - 29. Lauer SA, Grantz KH, Bi Q, Jones FK, Zheng Q, Meredith HR, et al. The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application. Ann Intern Med. 2020;172(9):577–82. pmid:32150748 - 30. He X, Lau EHY, Wu P, Deng X, Wang J, Hao X, et al. Temporal dynamics in viral shedding and transmissibility of COVID-19. Nat Med. 2020;26(5):672–5. pmid:32296168 - 31. Wölfel R, Corman VM, Guggemos W, Seilmaier M, Zange S, Müller MA, et al. Virological assessment of hospitalized patients with COVID-2019. Nature. 2020;581(7809):465–9. pmid:32235945 - 32. Puhach O, Meyer B, Eckerle I. SARS-CoV-2 viral load and shedding kinetics. Nat Rev Microbiol. 2023;21(3):147–61. pmid:36460930 - 33. Kass RE, Raftery AE. Bayes factors. J Am Stat Assoc. 1995;90(430):773–95. - 34. Kiang MV, Li ZR, Wrigley-Field E, Raquib RV, Lundberg DJ, Paglino E, et al. Applying machine learning to identify unrecognized COVID-19 deaths recorded as other causes of death in the United States. Sci Adv. 2026;12(12):eaef5697. pmid:41849604 - 35. Gadarian SK, Goodman SW, Pepinsky TB. Partisanship, health behavior, and policy attitudes in the early stages of the COVID-19 pandemic. PLoS One. 2021;16(4):e0249596. pmid:33826646 - 36. Kaashoek J, Testa C, Chen JT, Stolerman LM, Krieger N, Hanage WP, et al. The evolving roles of US political partisanship and social vulnerability in the COVID-19 pandemic from February 2020-February 2021. PLOS Glob Public Health. 2022;2(12):e0000557. pmid:36962752 - 37. Ryan M, Brindal E, Roberts M, Hickson RI. A behaviour and disease transmission model: incorporating the Health Belief Model for human behaviour into a simple transmission model. J R Soc Interface. 2024;21(215):20240038. pmid:38835247 - 38. Perlis RH, Lunz Trujillo K, Safarpour A, Quintana A, Simonson MD, Perlis J, et al. Community mobility and depressive symptoms during the COVID-19 pandemic in the United States. JAMA Netw Open. 2023;6(9):e2334945. pmid:37755830 - 39. Stadlbauer D, Tan J, Jiang K, Hernandez MM, Fabre S, Amanat F. Repeated cross-sectional sero-monitoring of SARS-CoV-2 in New York City. Nature. 2021;590(7844):146–50. - 40. Davis JT, Chinazzi M, Perra N, Mu K, Pastore Y Piontti A, Ajelli M, et al. Cryptic transmission of SARS-CoV-2 and the first COVID-19 wave. Nature. 2021;600(7887):127–32. pmid:34695837 - 41. Alroy KA, Crossa A, Dominianni C, Sell J, Bartley K, Sanderson M, et al. Population-based estimates of coronavirus disease 2019 (COVID-19)-like illness, COVID-19 illness, and rates of case ascertainment, hospitalizations, and deaths-Noninstitutionalized New York City residents, March-April 2020. Clin Infect Dis. 2021;73(9):1707–10. pmid:33458740 - 42. Kariya T. Rapid spread of COVID-19 in New York and the response of the community. Glob Health Med. 2020;2(2):123–6. pmid:33330789 - 43. Rǎdulescu A, Ballard S, Gonzalez K, Linton J. Dynamic coupling between the COVID-19 epidemic timeline and the behavioral response to PAUSE in New York State counties. PLoS One. 2021;16(8):e0255236. pmid:34347810 - 44. Colizza V, Vespignani A. Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: theory and simulations. J Theor Biol. 2008;251(3):450–67. pmid:18222487 - 45. Balcan D, Colizza V, Gonçalves B, Hu H, Ramasco JJ, Vespignani A. Multiscale mobility networks and the spatial spreading of infectious diseases. Proc Natl Acad Sci U S A. 2009;106(51):21484–9. pmid:20018697 - 46. Hale T, Angrist N, Goldszmidt R, Kira B, Petherick A, Phillips T, et al. A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker). Nat Hum Behav. 2021;5(4):529–38. pmid:33686204 - 47. Wang Z, Andrews MA, Wu Z-X, Wang L, Bauch CT. Coupled disease-behavior dynamics on complex networks: a review. Phys Life Rev. 2015;15:1–29. pmid:26211717 - 48. De Gaetano A, Barrat A, Paolotti D. Modeling the interplay between disease spread, behaviors, and disease perception with a data-driven approach. Math Biosci. 2024;378:109337. pmid:39510244 - 49. d’Andrea V, Gallotti R, Castaldo N, De Domenico M. Individual risk perception and empirical social structures shape the dynamics of infectious disease outbreaks. PLoS Comput Biol. 2022;18(2):e1009760. pmid:35171901 - 50. Ye Y, Zhang Q, Ruan Z, Cao Z, Xuan Q, Zeng DD. Effect of heterogeneous risk perception on information diffusion, behavior change, and disease transmission. Phys Rev E. 2020;102(4–1):042314. pmid:33212602