Parameter Estimation from ICC curves

This article introduces the concept of ICC (Incidence vs Cumulative Cases) curves, which are nonlinear transformations of the traditional ‘EPI curves’ (plots of disease incidence versus time) commonly used by epidemiologists. The main message of the present work is that describing an outbreak in terms of its associated ICC curve is not only natural from a dynamical point of view, but also advantageous for parameter identification and forecasting purposes. Specifically, we provide a method for the practical identification of epidemiological parameters of the SIR (Susceptible - Infected - Removed) compartmental model, directly from its associated ICC curve. Because the result is robust to under-reporting, this approach gives a simple way to estimate the basic reproduction ratio of a disease from reported incidence data. The present analysis also provides a theoretical justification for the forecasting methodology put forward in [1], which was shown to lead to useful estimates of the characteristics (expected duration, final number of cases, and peak intensity) of ongoing outbreaks.

We develop the concept of ICC curves in the context of the classical SIR model [2, 3], which is the simplest compartmental model that captures the basic phenomenon of disease transmission: before recovering (or being removed through disease-induced death), infected individuals transmit the disease to susceptible individuals, who in turn become infected. The process repeats itself until the epidemic has followed its course and no additional infections occur. Since many outbreaks have a time scale that is much shorter than the scale at which the total population $N$ changes, a common simplification is to assume that $N$ is constant. It is then natural to define the time-dependent quantities $s=S/N,$ $i=I/N,$ and $r=R/N$ as the expected sizes of the susceptible ($S$), infected ($I$), and recovered ($R$) compartments relative to $N,$ so that $s+i+r=1.$ The resulting SIR model is then given by

where $\beta$ is the contact rate of the disease and $\gamma$ is its recovery rate. Scaling time by $\tau=1/\gamma$ turns the above equations into a one-parameter ($R_{0}$) model, where $R_{0}=\beta/\gamma$ is the basic reproduction ratio [4, 5, 6] of the modeled epidemic. The disease spreads if $R_{0}>1$ and dies out when $R_{0}<1.$ It is also clear from the right-hand side of the second SIR equation that for $i\neq 0,$ the proportion of infected individuals grows as a function of time as long as $s>\gamma/\beta=1/R_{0}.$

In a compartmental model, an outbreak is a trajectory that starts in the vicinity of the disease-free equilibrium and ends at the endemic equilibrium, when the number of infected individuals approaches zero. We define an ICC curve, where ICC stands for Incidence vs. Cumulative Cases, as a representation of the dynamics associated with an outbreak, in the plane of coordinates $C$ (cumulative cases) and $\mathcal{I}$ (incidence). The following theorem, established in Section 2, gives an exact formulation of the ICC curves of the SIR model.

Details of the derivation of this theorem are given in Section 2.1. Its significance is illustrated in Figure 1, which shows a numerical simulation of the SIR model and a plot of the corresponding ICC curve for $R_{0}=2$. Instead of considering the time dependence of $S$, $I$, and $R$ as shown in the left panel of Figure 1, the course of the outbreak is captured by a single curve: a plot of incidence versus cumulative cases. As explained in Section 2.1, knowledge of the ICC curve is sufficient to reconstruct the temporal behavior of $S$, $I$, and $R$. In other words, the ICC curve contains the entire information associated with any outbreak described by the SIR model.

An important motivation for studying ICC curves is that they can easily be fitted to case data. Specifically, we show in Section 3.1 that, for fixed $N$, there exists a single set of parameters $(\beta,\gamma,\kappa)$ that minimizes the root mean square error between (noisy) simulated outbreak data and the ICC curve ${\mathcal{I}}=G_{\kappa,N}(C)$ given by Equation (1). Additionally, a numerical investigation indicates that such an estimation is robust to noise. As a consequence, the present approach provides a method to estimate the parameters of the SIR model, thereby making this model practically identifiable.

The rest of this manuscript is organized as follows. Section 2 presents the derivation of Equations (1) and 2, and discusses the robustness of ICC curves to systemic under- (or over-) reporting. Section 3 explains how model parameters may be retrieved from ICC curves and assesses the robustness of the proposed procedure to reporting noise. Section 4 briefly discusses parameter estimation as an outbreak unfolds. Section 5 shows a few examples of application to real data, for outbreaks of gastroenteritis and COVID-19. Future extensions and applications of the methodology proposed in this article are reviewed in Section 6.

Structural identifiability (see e.g. [8] for a review) relates to the ability of uniquely identifying model parameters from knowledge of the model output. In this context, $N$ is assumed to be known since it represents the population of the system to which the model is applied. Given $C(t),$ one can calculate ${\mathcal{I}}(t)=dC/dt=G_{\kappa,N}(C(t)).$ From Equation (1), it is clear that knowledge of $C(t)$ and ${\mathcal{I}}(t)$ uniquely determines $\beta,$ $\gamma,$ and $\kappa.$ Indeed, if two sets of parameters, $(\beta^{(1)},\gamma^{(1)},\kappa^{(1)})$ and $(\beta^{(2)},\gamma^{(2)},\kappa^{(2)})$, led to the same output $C(t)$ (and thus to the same values of its derivative ${\mathcal{I}}(t)$), we would have

for all values of $\frac{C}{N}\in\left[0,\frac{C_{\infty}}{N}\right]$, which leads to $\beta^{(1)}=\beta^{(2)}$, $R_{0}^{(1)}=R_{0}^{(2)}$, and $\kappa^{(1)}=\kappa^{(2)}.$ Therefore, all of the parameters of the SIR model, including initial conditions, are structurally identifiable from knowledge of $C(t),$ the cumulative number of cases, and $N$, the size of the susceptible population. Similar results were obtained in [9, 10, 11]. If $N$ is unknown, setting

for all values of $\frac{C}{N}\in\left[0,\frac{C_{\infty}}{N}\right]$, additionally leads to $N^{(1)}=N^{(2)}$, making the size of the population $N$ a structurally identifiable parameter as well.

Whereas structural identifiability depends on the nature of the model itself, practical identifiability relates to parameter identifiability in the presence of measurement noise, as well as (hopefully small) deviations between model and reality (see e.g. [8]). This question may be addressed numerically by simulating an exact solution of the model, adding noise to its output, and estimating the model parameters that best fit the perturbed output. Two recent articles [11, 12] discuss the practical identifiability of a variety of compartmental models and reach different conclusions regarding parameter identification of the SEIR model from cumulative data. In [11], the cumulative data is multiplied by $1+\epsilon,$ where $\epsilon$ is normally distributed with zero mean, whereas in [12], Poisson-distributed noise of mean equal to the current incidence value [13] is used as a substitute for the incidence data. Whether or not parameters are identifiable in practice therefore depends on the type of observational noise that affects surveillance reports and thus the EPI curve.

Adding independently-distributed noise to the incidence rather than the cumulative data is more realistic since errors on $C$ are considered to accumulate and not be independent [14]. On the one hand, assuming Poisson-type of uncertainty is natural for data that represent counts, but the effect of replacing incidence data by a Poisson random variable of same mean will lead to significant relative deviations only if $N$ is small (a few hundreds). For larger values of $N$ and thus of ${\mathcal{I}},$ the relative change $\delta{\mathcal{I}}=|{\mathcal{I}}-\hbox{Poisson}({\mathcal{I}})|/{\mathcal{I}}$ scales like $1/\sqrt{\mathcal{I}}.$ Indeed,

With Stirling’s formula, we obtain

If it is believed that the diminishing relative influence of the Poisson-distributed noise for large values of $N$ is not realistic for a specific epidemic, then multiplying $\mathcal{I}$ by $(1+\epsilon)$ where $\epsilon$ is normally distributed with mean zero, circumvents this issue. We consider both types of noise below.

To test practical identifiability, we simulate the SIR model and generate a discrete time series for the cumulative number of cases of the simulated outbreak. This time series is used as a reference and represents the unperturbed (or true) data. We then add noise to the associated discrete incidence data (defined as the difference of cumulative values between consecutive time points), following the procedure described in [13], with the additional realistic constraint that the final number of cases is the same for all noisy realizations of the same outbreak. Two examples of simulated outbreaks generated in this fashion are shown in Figure 2 for different values of $R_{0}$ (associated with different orders of magnitude of the incidence $\mathcal{I}$) and different discretization time steps $\delta_{t}$. For $R_{0}=2$, we use $\delta_{t}=1$. For $R_{0}=10.5$, we use $\delta_{t}=2.$ As shown in Figure 2, this latter choice leads to a discretized ICC curve with very few points in regions of large incidence, and simulates a situation where incidence data is of poorer frequency. Each time series, consisting of noisy incidence data, is used to calculate an associated time series for the corresponding cumulative data. The resulting set of points is plotted in the ($C$, $\mathcal{I}$) plane, and defines an empirical ICC curve. In both panels of Figure 2, $C$ and $\mathcal{I}$ are estimated from discrete incidence data as described in Equation (5) below.

In principle, the formulas for $\beta_{c}$ and $\gamma_{c}$ given in Appendix A may be used with any number of points $M$, as long as the data points represent the entire ICC curve. Is is natural to ask how the method performs when applied to a developing outbreak. To test this, we simulate 50,000 possible outbreaks using the same reference simulations as before, and with the two different types of noise introduced above to represent reporting error. We then follow the evolution of the distributions of the estimated parameters as time increases. For $R_{0}=2$, the outbreak has run its course after 60 units of time, and we include up to $M_{max}=40$ data points, with consecutive values separated by 1 unit of time (which could correspond to a day or a week depending on the disease). For $R_{0}=10.5$, we pick $M_{max}=80$, with data points separated by 2 units of time. Figure 6 shows estimated values of $\beta$, $\gamma$, and $R_{0}=2$ for the two types of reporting noise. In both cases, estimates of $\beta$ and $\gamma$ are close to their actual values near $t=20$, which is before the peak of the outbreak (see left panel of Figure 1). The value of $R_{0}$ converges more slowly (circles in the bottom row of Figure 6) but the ratio of the estimated values of $\beta$ and $\gamma$ (stars) is close to the actual value of $R_{0}$ near $t=12$, fairly early in the outbreak.

Similar results for $R_{0}=10.5$, together with plots showing the evolution of the distributions of $\beta_{c}$, $\gamma_{c}$, and $R_{0}$, are provided as supplementary material. Parameter estimates for $\beta$ and $\gamma$ are close to their true values near $t=22$, which is near the peak of the outbreak, with the ratio $<\beta>/<\gamma>$ giving a reasonable estimate of $R_{0}$ starting at $t=18$, even in the presence of large uncertainty.

So far, we have assumed the size of the susceptible population, $N$, to be known. In practice, the value of $N$ that should be used to fit the ICC curve to outbreak data is often unknown since it may reflect under-reporting (as discussed in Section 2.2) or spatial effects (existence of localized clusters to which the spread is restricted) associated with the outbreak. Moreover, it could be argued that the SIR model is too simplistic to represent real-life outbreaks. The examples of [1] and those discussed below however show there is merit to the present approach.

In what follows, we fit ICC curves to surveillance reports by picking a range of values of $N$ for which the root mean square error (RMSE) between the ICC curve and the data is within 2% of the minimum RMSE value. The data points used in the evaluation of the RMSE are found from the reported cumulative data after smoothing and interpolation. This latter procedure estimates missing incidence values and removes negative incidence reports, if any. The resulting time series is expected to be a better approximation of the ‘true’ (i.e. as close as possible to the output of a compartmental model) evolution of the cumulative data, and is used for the parameter estimation procedure discussed in this article. Our first example is a gastroenteritis outbreak in a nursing home in Mallorca, Spain, and is therefore spatially localized. The other examples are estimations of the basic reproduction ratio for COVID-19 outbreaks in Hubei Province (China), the Republic of Korea, and France, based on data available until 4/15/2020.

This article introduces a parameter estimation method for disease outbreaks that bypasses the numerical integration of epidemiological models. The approach relies on the use of ICC curves, also introduced here. ICC curves relate incidence to the cumulative number of cases $C$ and may be viewed as nonlinear transformations of the traditional epidemiological curves, in which the time variable is replaced by $C$, a monotonically increasing but nonlinear function of time. For each single-wave outbreak, the ICC curve has a simple concave-down shape that crosses the horizontal axis at the origin and at the final value of $C$.

The formulas presented in Section 2, which extend the parabolic approximation suggested in [1], are exact for the SIR model. The numerical experiments of Sections 3 and 4 show the method is robust to noise and may be used for parameter estimation as an outbreak unfolds, with the understanding that reasonably accurate estimates can only be reached shortly before, or after, incidence peaks. This is not a shortcoming of the present approach and was also observed in [11] when fitting synthetic prevalence time series.

Because of the simplicity of Equation 1, and the existence of a unique set of parameters that best fits any collection of data points, parameter distributions may be generated directly from epidemiological data. In the case of the SIR model, this methodology presents a fast and novel alternative to more traditional parameter estimation strategies, such as particle Makov Chain Monte Carlo methods (PMCMC; see [21] for a review). This may be beneficial in pandemic situations where epidemiological estimates need to be updated daily and in many locations simultaneously. For instance, recent work on COVID-19 data from Wuhan suggests that an SIR model can better capture the information contained in case reports than an SEIR model [22]. For more complex compartmental models, distributions generated by the present approach may be used as priors for PMCMC estimations of the contact and recovery rates of a disease.

When applying the proposed method to surveillance data, an estimate of $N$ may not always be available. The examples of Section 5 suggest that this parameter may be identified by optimizing the fit between the ICC curve and a smoothed and interpolated version of the reported data. However, any variability in the selected value of $N$ will be associated with variability in estimates of the basic reproduction number $R_{0}$. Indeed, for an outbreak that has completed its course, any good fit of the reported data by an ICC curve will produce consistent values of $C_{\infty}$, the final number of cases for the outbreak (see for instance the plots of Figure 8). Because of Equation 3, if $C_{\infty}$ is known, the estimated value of $R_{0}$ is only a function of $N$. This uncertainty is inherent to the SIR model and cannot be avoided. The existence of a value $N_{m}$ of $N$ that best fits the data is encouraging, but more robust estimates are likely to be obtained if additional information is available. In the absence thereof, a range of values of $N$ close to the optimal value $N_{m}$ should be used to produce credible intervals for the basic reproduction ratio $R_{0}$.

As previously mentioned, empirical (see [1]) and numerical observations by the author suggest the concave-down shape of the ICC curve is ubiquitous in outbreak data and in compartmental epidemiological models. It would therefore be interesting to explore methods which, like the discussion of Section 2, lead to exact formulations of ICC curves. In particular, knowledge of how model parameters affect the shape of ICC curves could provide simple means to visualize the consequences of mitigation effects. Combined with data assimilation, ICC curves may thus present a convenient paradigm for forecasting the course of an outbreak.