Mapping Incidence and Prevalence Peak Data for SIR Forecasting Applications

Compartmental models have seen broad usage at the onset of several disease outbreaks in the last century as a means to project expected numbers of infected and to drive healthcare response. Broadly speaking, a compartmental model describes the dynamics of a disease spread by breaking a population into set categories and modeling the process by which individuals transition through these categories. Perhaps the most basic is the Kermack–McKendrick model, often called the SIR or Susceptible-Infectious-Recovered model, which models the movement of subjects from being Susceptible, to Infected (and contagious), and then to Removed (Kermack, McKendrick and Walker, 1927). Since the development of the SIR model, further research has extended the idea to include additional compartments, such as the Exposed category in the SEIR model – where a subject has been exposed to the disease but is not yet infectious – and the ability to move back to the Susceptible category in the SEIS model (see Walter and Contreras (1999); Hethcote (2000) for an overview of each). Compartmental models have been applied to projection tasks such as the 2014-15 Ebola epidemic (Chretien, Riley and George, 2015), the 2009 A/H1N1 influenza pandemic (Nsoesie et al., 2014), several HIV outbreaks (Anderson, 1988; Golub, Gorr and Gould, 1993; Nyabadza, Mukandavire and Hove-Musekwa, 2011), the recent COVID-19 pandemic (Zhao and Chen, 2020; Cooper, Mondal and Antonopoulos, 2020; Zhang et al., 2022), and many other epidemiological projection tasks in the last century (see, for instance, Guanghong et al. (2004); LaDeau et al. (2011); Zhan et al. (2019) and references therein).

The elegance of compartmental models is in their succinct ability to describe the state of a population in terms of how subjects transfer to and from the different categories, and thus fitting these models requires estimation of interpretable quantities such as rates of infection and recovery. For instance, denote the proportion of the population in the susceptible, infectious, and removed compartments by $S_{t},I_{t},$ and $R_{t},$ respectively, such that $S_{t}+I_{t}+R_{t}=1$ for all $t$. Then the SIR model is determined by the equations

where $\beta>0$ is the disease transmission rate and $\gamma>0$ is the rate of recovery. If one knew these two rates, and the initial number of individuals in each category – $S_{0},I_{0}$ and $R_{0}$ – the numbers $S_{t},I_{t},\text{and}\;R_{t}$ could be numerically simulated for any time-point $t$. While modern computational tools do make simulating these quantities feasible, the need to simulate the entire system numerically to get the exact counts in each compartment is challenging. That is, fitting these models to data and can be computationally expensive for more intricate compartmental structures. Researchers have determined analytic solutions for the entire SIR model (explicit forms/approximations to the number of susceptible, infected, and removed at a specific time) that do not require numeric simulation (Harko, Lobo and Mak, 2014; Schlickeiser and Kröger, 2021; Carvalho and Gonçalves, 2021). However, these solutions involve reparameterizing the time axis, and explicit calculations back onto the original time axis require numeric integration or approximation methods. Similar approaches are also used to obtain exact and approximate solutions to the more complicated compartmental models, such as the SEIR model (Wang, Wei and Zhang, 2014; Piovella, 2020) and the SIRS model (Acedo, González-Parra and Arenas, 2010).

Analytic maps from the initial starting parameters of an SIR curve to quantities of interest (QoIs) – such as the value of the peak of the infected curve and the limiting number of susceptible individuals – were developed as early as the mid-1900s (Kendall, 1956; Bailey, 1957; Hethcote, 1976). The focus on these quantities has generally been motivated by their usefulness to making public health decisions. For instance, in papers such as Hethcote (2000) and Weiss (2013), the maximum value of the infected compartment (not the time of maximum infection) was studied for the purpose of informing public health officials of the maximum number of infections after the initial estimation of the disease transmission and recovery rates, since knowing this maximum quantity informs how many hospital beds might be needed. The authors in Castro et al. (2020) point out that the time of peak infection is also informative for healthcare officials, and they develop an analytic form for the peak time of infection using an approximation of the SIR curve. An exact form for the peak infection time is available in Kendall (1956) and Deakin (1975), albeit in terms of an integral without a closed form. Several modern papers study fast approximations to the integral form for peak infection time (Cadoni and Gaeta, 2020; Turkyilmazoglu, 2021), using the previously-mentioned analytic solutions to the SIR model. These different approximations for peak infection time are evaluated and compared in Kröger, Turkyilmazoglu and Schlickeiser (2021).

Some recent papers study SIR curve QoIs for inferential tasks related to modeling an epidemic (Miller, 2009; Lang et al., 2018). In Amaro (2023), the Gumbel distribution is suggested as a good approximation of the infection curve in an SIR model, and maps between the SIR curve parameters and peak value/time are used to develop Method of Moments estimators of the Gumbel parameters. The Gumbel distribution is then used to approximate the exact solution of the SIR model. In Osthus et al. (2017), the relationship between peak infection value and the SIR parameters is used to incorporate historical data on epidemic peaks into the inferential task of fitting an SIR curve during the early stages of an epidemic. The authors point out that an SIR curve is sensitive to small perturbations in the transmission and recovery rates, and that incorporating these data discourages models that well-represent early data but drastically over-predict the peak infection value. In McAndrew et al. (2024), both historic and surveyed QoIs are used to constrain an SEIRH model (where the added “H” category refers to people hospitalized at a given time) in much the same way as Osthus et al. (2017).

Both Amaro (2023) and Osthus et al. (2017) use information on peak infection value and time to inform the modeling of a pandemic and epidemic, respectively. Since most modern monitoring systems approximate the daily number of new cases of a disease, called incidence, this is not the most practically useful development. The number of infected individuals at any given time, called prevalence, is typically an unobserved quantity (Noordzij et al., 2010). For examples of all quantities in an SIR curve on a fixed time axis, see Figure 1. Note that peak prevalence and peak incidence need not occur at the same time, and that it is mathematically possible for the incidence to be greater than prevalence (as visualized on the right side of Figure 1). This being said, in most disease outbreaks of note, prevalence is typically greater than incidence.

The model in McAndrew et al. (2024) incorporates incidence QoI data by extending the parameter space to include peak incidence value (PIV) and peak incidence time (PIT) and by defining a prior on these values using either historic or survey data. While this model achieves a similar goal to the ultimate goal of this paper, its main distinction is that that it requires an importance sampling scheme to fit, which we outline here for completeness: priors on the rate parameters and initial values for the compartmental model are used to sample proposal values, the entire system is numerically simulated to determine the (PIV, PIT), then the likelihood of this system is determined via the prior probability of (PIV, PIT). Note that while this need not be true for the SEIRH model considered in McAndrew et al. (2024) (since it is not an SIR model), defining both a prior on the rate parameters and on PIV and PIT for an SIR model is redundant given the initial values $S_{0}$, $I_{0}$, and $R_{0}$, since a given set of starting parameters should immediately determine these QoIs.

Incorporating historical peak values and times of an outbreak into a Bayesian compartmental forecasting model via prior specifications is highly useful to constrain forecasts¹¹1For a review of the differences between disease forecasts and disease projections, see Massad et al. (2005)., especially early in the outbreak (see Osthus et al. (2017), Osthus et al. (2019), and McAndrew et al. (2024)). Relating these historical QoIs requires determining a relationship between historical peak values and times as well as the parameters and initial conditions of the compartmental model (e.g., the SIR model). The challenge lies in this relationship. Analytically, this relationship is known for observed peak prevalence value (PPV) and peak prevalence time (PPT). However, almost all infectious disease data is of incidence. Despite being principally unsound, the data application in Osthus et al. (2017) treated incidence data as if it were prevalence in order to use the known analytic relationships. In this paper, we make three contributions. First, we develop the methods to map peak values and times of incidence to the parameters of the SIR model, given the initial conditions. Second, we demonstrate how to incorporate these new mappings into the model of Osthus et al. (2017) and provide the code used to do so. Third, we compare the forecasts of the modeling framework of Osthus et al. (2017) using both incidence data as incidence, and using the misspecified prevalence data. In addition to comparing these forecasts, we show that if SIR parameter inference is desired, then mistaking incidence data for prevalence data will result in biased estimation.

This paper is outlined as follows. In Section 2, we review existing methods for mapping SIR parameters to peak/time of prevalence. In Section 3, we develop identical maps for incidence, then provide methods for inverting these maps for both incidence and prevalence. In Section 4, we develop and improve upon the modeling framework in Osthus et al. (2017) to incorporate historic incidence peak/time data when fitting an SIR curve. We apply this model to influenza data in Section 5.

The system (1a)-(1c) can be solved analytically for all time by reparameterizing the time axis to be instead in terms of the number of individuals removed from the system beyond the initial amount in the Removed category $R_{0}$,

where $\tau_{t}$ is the inverse of the map

where $[x]^{+}:=\max\{x,0\}.$ The above form for the SIR dynamics is particularly useful since it allows one to calculate the number of individuals in each category without numerically simulating the entire system. The major drawback of this form is that mapping values back to an interpretable time requires one to approximate the integral in (3), since this integral is nonelementary. There exist other methods that derive an analytic solution to the SIR system (Harko, Lobo and Mak, 2014). We use this formulation since it is well-known (having been originally developed in Kendall (1956)), and because it simplifies the calculation of SIR curve QoIs, such as PPV and PPT. For instance, the maximum of (2b) with respect to $\tau$ occurs at

which can be derived via straightforward calculus. PPT can then be approximated from (3), and PPV calculated from (2b).

The process above of calculating PPV and PPT, denoted $(I_{t_{\tau^{*}}},t_{\tau^{*}}),$ from the initial starting parameters can be reversed for fixed $S_{0},I_{0},R_{0}$. That is, given $(I_{t_{\tau^{*}}},t_{\tau^{*}})$, we propose a method to calculate $(\beta,\gamma)$. Define $\rho$ as

The quantity $S_{0}\rho$ is more widely known as the basic reproduction number, which measures the average number of additional infections generated by a single new infection (Cadoni and Gaeta, 2020). Plugging (4) into (2b) gives

From here, setting $I_{t}:=I_{t_{\tau^{*}}}$ and solving for $\rho$ is a (piecewise) convex optimization problem, which is typically fast and accurate computationally. Note that a change of variables for (3) gives

and that $\beta\tau^{*}=\log(\rho S_{0})$. Thus, once the value of $\rho$ is approximated, (7) can be evaluated to get $\beta t_{\tau^{*}}$. Since $t_{\tau^{*}}$ is assumed known a priori, the values of $\beta$ and $\gamma$ can be obtained via arithmetic.

We introduce the above derivations not only because they will be used in the next Section, but also to motivate the type of calculation we aim to develop in this paper. While the above approach still requires computational methods to map from PPV and PPT to the SIR model parameters, this approach is much more direct than the brute force method of simulating several $(\beta,\gamma)$ combinations until a prevalence curve with a peak sufficiently near $(I_{t_{\tau^{*}}},t_{\tau^{*}})$ is discovered (Prangle, 2016).

As mentioned above, a primary focus of this work will be to develop maps between PIV and PIT and the initial SIR curve parameters. In this direction, one can derive an equivalent formulation of the SIR dynamics in (1a) - (1c) by replacing (1b) with

Putting (8) into words, prevalence is equal to the prevalence at the last time step, plus those in the infected category that infect those in the susceptible category, minus the number in the infected category that are removed from the system. The form for prevalence in (8) is particularly useful for this application since it contains a term that explicitly models incidence at time $t$: $\beta S_{t-1}I_{t-1}$.

In this section, we introduce the Dirichlet-Beta state-space model (DBSSM) from Osthus et al. (2017) and update it to incorporate historic PIV and PIT data. The original formulation of the DBSSM was to answer an observed issue associated with using the SIR model for early-pandemic forecasting tasks. Namely, that two SIR curves that reasonably fit early count data may lead to drastically different PPV predictions. In a simulation example, Osthus et al. (2017) show two such SIR curves that have peaks that differ by 30% of the entire population, even though they have a nearly-identical fit to the early-pandemic data observations (see Figure 3 in the cited paper). To address this stability issue, the DBSSM incorporates historic PPV and PPT data into the prior specifications to discourage SIR curve fits with peak values that are greatly above reasonable expectations. As mentioned previously, this incorporation of prevalence data is not the most practical approach, since incidence data is generally the observed quantity. A further shortcoming in the original DBSSM formulation is that it learns a map between PPT and the SIR curve parameters, rather than using an analytic map. After introducing this model, we will identify ways that the methodology developed in this paper will improve these issues for the DBSSM.

Let $y_{t}$ be the observed proportion in a population that tested positive for some disease at some timepoint $t$, and let $\theta_{t}=(S_{t},I_{t},R_{t})^{\prime}$. Then the DBSSM is defined as,

where $\phi=\{S_{0},I_{0},R_{0},\beta,\gamma,\lambda,\iota\}$, $\lambda,\iota$ are variance control parameters, and $f$ is a map that propagates the SIR system determined by $(\theta_{t-1},\beta,\gamma)$ forward one step according to (1a)-(1c). Note that, by this set up, the set of parameters $\theta_{0:t^{\prime}}=\{\theta_{0},\theta_{1},\dots,\theta_{t^{\prime}}\}$ is a first-order Markov chain, and that for all $s\neq t$, the data observations $y_{s}$ and $y_{t}$ are independent given $\theta_{t}$. The variable $\text{In}_{t}$ denotes the incidence of the system at time $t$; in the original formulation of this model, the prevalence – $I_{t}$ – was used here. The incidence at time $t$ is directly calculated using $\theta_{t}$ and (8).

The conditional expectations of the model described by (17) and (18) are unbiased,

while their respective variances reduce to zero as $\lambda,\iota\to\infty$. Of course, the conditional mean in (20) is dependant upon the accuracy of $f$ in propogating the latent space $\theta_{t-1}$ forward one time step. The authors in Osthus et al. (2017) used a fourth-order Runge-Kutta approximation and observed reasonable accuracy.

We review the full Bayesian framework of the DBSSM and provide the prior specification in Appendix A. The main innovation of this model is that the parameter space is expanded to include the latent variable $z=(PPT,PPV)$, and this latent variable is given a prior that incorporates historic PPT and PPV data. In the following, we review this prior, $\pi(z|\theta_{0})$, and the mechanism by which this prior informs the SIR parameters $\beta$ and $\gamma$. These priors are then each updated according to the theory developed in this paper.

We recreate the data application from Osthus et al. (2017), using the updated model from Section 4. The aim of this application is not to improve the forecasting in the original formulation of the DBSSM (see Remark 1). Rather, we will demonstrate that the forecasting capabilities of this model remain the same, while we also observe different estimates for the infection rate $\beta$, the recovery rate $\gamma$ and the basic reproduction number $\rho$.

The source data modified and then used for this application are counts of patients seen in the US with an influenza-like illness (ILI), where ILI is defined as having a temperature of at least 100 degrees Fahrenheit, a cough and/or a sore throat, and no known cause for those symptoms other than influenza (CDC, 2024). These data are collected weekly, where more than 3400 outpatient healthcare providers report to the CDC the number of patients with ILI they treated (CDC, 2023).

The number of patients reported as having ILI will naturally also include cases of respiratory illnesses other than influenza. Following the approach of Shaman et al. (2013), we use virologic surveillance data (where patients are actually tested for influenza) to estimate the proportion of ILI patients with influenza, then multiply ILI data by this proportion. This corrected data is referred to as ILI+. For more details on this adjustment, see Shaman et al. (2013). Note that the ILI+ data estimates the weekly incidence of influenza cases – not prevalence.

To fit both versions of the DBSSM, we use ten influenza seasons: the seasons that started in the years 2002-2007, and the seasons the started during 2010-2013. The years 2008 and 2009 were omitted to be consistent with Osthus et al. (2017); these two years correspond to a pandemic and the focus of that work was to forecast seasonal influenza. Each season is defined as 35 consecutive weeks starting on roughly the first week of October (epidemiology week 40, treated as $t=1$).

For a estimated proportion of individuals in a population infected with influenza at timepoints $\mathbb{T}=\{1,\dots,T\}$, suppose only the ILI+ data up through $t^{\prime}\in\mathbb{T}$ are observed. Given this, we simulate 62500 from the posterior $\pi(\boldsymbol{\theta_{{1:t^{\prime}}}},\boldsymbol{\phi}|y_{1:t^{\prime}})$ for four separate chains, discarding the first 12500 as burn-in and thinning out all but every tenth observation in the remaining samples. Given these draws from the posterior distribution, the posterior predictive density, $\pi(y_{(t^{\prime}+1):T}|y_{1:t^{\prime}})$, is used to estimate “future” observations of ILI+ data.

We perform two separate fits of this posterior model, on the first thirteen days ($t^{\prime}=13$) and on the first twenty two days ($t^{\prime}=22$), for the considered ILI+ data for the 2010 influenza season in the United States. These fits are performed both using the original formulation of the DBSSM, which naively uses incidence data directly in place of prevalence data, and using the new formulation developed in this paper, which uses the new maps developed in Section 3 for a more principled treatment of incidence data. These fits and forecasts are outlined in Figure 2. The dark shaded grey regions prior to $t^{\prime}$ mark the 95 percentiles of the posterior density, while the lighter grey shaded regions after $t^{\prime}$ make the 95% prediction intervals. The forecast using incidence data up until $t^{\prime}=13$ has a narrower prediction interval than the one using prevalence data; each of the forecasts that use data up through $t^{\prime}=22$ are laregly comparable.

In addition to the slight improvements on forecasting we observe in Figure 2, this method also has strong implications for the interpretability of $\beta,\gamma$ for the fitted model. Indeed, only the updated version of the DBSSM developed in this paper leads to realizations of these parameters that can accurately be interpreted as the infection rate ($\beta$), recovery rate ($\gamma$), and the basic reproduction number ($S_{0}\rho$).

The main contribution of this paper is the development of methods to map the time and value of peak incidence to the SIR curve parameters, and vice versa, for the purpose of forecasting tasks and inference on disease rate parameters. We do this by computationally solving a system of equations ((9) and (10)). There are several impactful uses of these maps in the context of previous literature. First, much like how the peak prevalence value (PPV) and time (PPT) are useful for public health response to an epidemic (Weiss, 2013), the analogous quantities for incidence are also useful, since they describe the influx new patients entering the hospital system on a given day. Second, this work improves upon existing work that uses historical prevalence data to model epidemics by creating a map from PIT and PIV to the SIR parameters, since incidence is typically the data that is available for ongoing epidemics (Osthus et al., 2017; Amaro, 2023). In the case of the application in Osthus et al. (2017), where incidence data were used in place of prevalence data without justification, we have shown that this leads to biased SIR parameter estimates (see Figure 3). Furthermore, our results indicate that forecasts performed using the erroneous data specification leads to larger prediction intervals earlier on in the outbreak, although this forecast is largely comparable for the correct specification later on in the outbreak (see Figure 2 and Remark 1). Of course, it remains more correct in principle to use incidence data appropriately when fitting compartment models for forecasting with ongoing incidence data (Nsoesie, Mararthe and Brownstein, 2013; Chowell et al., 2016; Abolmaali and Shirzaei, 2021). We have provided a modeling framework that incorporates this data appropriately (see Section 4). Lastly, while the methods discussed in McAndrew et al. (2024) do incorporate incidence data to fit forecasting models, this paper uses a Bayesian framework and importance sampling. Since the maps developed here are deterministic, they can be used in both a Bayesian and a Frequentist framework.

As a direction for future work, it would be useful to investigate better approximations for the solution to (13). The Taylor Approximation in Section 3.2.1 was by far the fastest computationally, but it came with the highest error on PIT. Finding a fast and accurate approximation to this equation would greatly increase the runtime for applications where the map between the SIR parameters and PIT/PIV must be evaluated several hundreds of thousands of times. However, for most applications (including the one in this paper), the Compute Integral approximation is sufficiently fast.

As a second direction for future work, it would be interesting to investigate the analogous maps for more complicated compartmental models, such as the SEIR model, the SIRS model, and the SEIRH model.

All software used to perform the simulations and studies in this paper are publicly available at https://github.com/lanl/precog.

Research presented in this article was partially supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20240066DR. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

This research was partially funded by NIH/NIGMS under grant R01GM130668-01 awarded to Sara Y. Del Valle.

The authors have no non-financial or other financial competing interests to declare that are relevant to the content of this article other than the aforementioned declared funding sources.