Parameter estimation from ICC curves

Abstract

Incidence vs. Cumulative Cases (ICC) curves are introduced and shown to provide a simple framework for parameter identification in the case of the most elementary epidemiological model, consisting of susceptible, infected, and removed compartments. This novel methodology is used to estimate the basic reproduction ratio of recent outbreaks, including those associated with the ongoing COVID-19 pandemic.

Keywords:

• Outbreak
• parameter identification
• compartmental models

1. Overview

This article introduces the concept of ICC (Incidence vs. Cumulative Cases) curves, which are nonlinear transformations of the traditional epidemiological (EPI) curves (plots of disease incidence versus time). 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 systematic 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 [11], 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 [8, 9], 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 $\frac{\mathrm{d}s}{\mathrm{d}t}=-\beta si;\frac{\mathrm{d}i}{\mathrm{d}t}=\beta si-\gamma i;\frac{\mathrm{d}r}{\mathrm{d}t}=\gamma i,$where β is the contact rate of the disease and γ 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 [3, 4, 17] of the modelled 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\ne 0,$ the proportion of infected individuals grows as a function of time as long as $s>\gamma /\beta =1/R_0.$

In a general 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 curve of the SIR model.

Theorem 1.1

For a total susceptible population of size N and each initial condition $\kappa =\frac{S(0)}{N}=1-\frac{C(0)}{N}$, the ICC curve of the SIR model is given by (1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(1) where $C=I+R\in [0,C_{\mathrm{\infty }}]$ and $C_{\mathrm{\infty }},$ the final number of cases, is the positive solution of the transcendental equation (2) $C_{\mathrm{\infty }}+\frac{N}{R_0}\ln (1-\frac{C_{\mathrm{\infty }}}{N})-\frac{N}{R_0}\ln (\kappa )=0.$(2)

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 behaviour of S, I, and R. In other words, the ICC curve contains the entire information associated with any outbreak described by the SIR model.

Figure 1. (Color online). Simulation of the SIR model showing the time-dependence of the S, I, and R compartments (left), for $R_0=2$. The corresponding ICC curve is shown on the right. The exact expression given by Equation (1(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(1) ) is plotted as a solid black line; the numerically estimated ICC curve is shown as a thick cyan line.

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(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(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(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(1) ) and (2(2) $C_{\mathrm{\infty }}+\frac{N}{R_0}\ln (1-\frac{C_{\mathrm{\infty }}}{N})-\frac{N}{R_0}\ln (\kappa )=0.$(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.

2. ICC curves of the SIR model

2.1. Proof of Theorem 1.1

Epidemiologists use the epidemiological (EPI) curve, a plot of incidence as a function of time, to visualize the course of an epidemic. In the SIR model, incidence is defined as $\mathcal{I}=\beta S\frac{I}{N}=\frac{\mathrm{d}C}{\mathrm{d}t},$where C = R + I is the cumulative number of cases at time t. Note that C includes those who have recovered (and thus were a case at some point in time), as well as those who are currently infected. To derive Equation (1(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(1) ), write $\mathrm{d}s/dr=-\beta s/\gamma =-R_{0}s,$ which can be integrated to give $s=\kappa \exp (-R_{0}r),$ where $\kappa \in (0,1)$ is a constant that depends on initial conditions. Typically, $\kappa \simeq 1$ (see e.g. [1]) since for an outbreak in a naive population we normally have $r(0)=0,$ $i(0)<<1,$ and thus $\kappa =s(0)\simeq 1.$ However here we do not make this assumption and keep the parameter κ in the following discussion. With $c=C/N$ and $s=1-(i+r)=1-c,$ we obtain $r=-\frac{1}{R_0}\ln (s/\kappa )=-\frac{1}{R_0}\ln ((1-c)/\kappa ),$and consequently $\frac{\mathrm{d}c}{\mathrm{d}t}=\beta is=\beta (c-r)(1-c)=\beta (c+\frac{1}{R_0}\ln ((1-c)/\kappa ))(1-c).$Moving back to dimensional variables gives the expression in Equation (1(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(1) ), i.e. $\mathcal{I}=\mathrm{d}C/\mathrm{d}t=G_{\kappa ,N}(C)$.

An approximate expression of the SIR ICC curve, assuming $\kappa =1$ and C/N sufficiently small was given in [11]. Since $G_{\kappa ,N}\ge 0$ on $[0,C_{\mathrm{\infty }}]$, for a given initial condition $\kappa =1-C(0)/N$, the part of the ICC curve traced out as the outbreak unfolds will be such that $C(t)\ge C(0)>0$. However, the ICC curve is well defined for C in the range $[0,C_{\mathrm{\infty }}]$, which is the definition we adopt here. Even though the time variable does not appear in the ICC curve, the value of κ is not arbitrary: $C(0)$ is the value of C when R = 0, where R represents the removed compartment.

ICC curves are particularly useful for forecasting the course of an outbreak because they contain information on quantities that are of interest to public health practitioners, namely the final number of cases ($C_{\mathrm{\infty }}$) and peak incidence (${\mathcal{I}}_{\max }$). The latter is simply the maximum of $G_{\kappa ,N}$; the former may be found as follows. Since C is the cumulative number of cases, it is a monotonic function of time and its derivative, $G_{\kappa ,N}(C),$ must remain non-negative. For C/N small, $G_{\kappa ,N}(C)=G_{\kappa ,N}(0)+C(\beta -\gamma +\gamma \ln (\kappa ))+\mathcal{O}(C^2)$where $G_{\kappa ,N}(0)=-N\gamma \ln (\kappa )>0.$ By continuity, $G_{\kappa ,N}(C)$ will remain non-negative on the interval $[0,C_{\mathrm{\infty }}],$ where $C_{\mathrm{\infty }}$ is the unique positive solution of Equation (2(2) $C_{\mathrm{\infty }}+\frac{N}{R_0}\ln (1-\frac{C_{\mathrm{\infty }}}{N})-\frac{N}{R_0}\ln (\kappa )=0.$(2) ). By writing Equation (2(2) $C_{\mathrm{\infty }}+\frac{N}{R_0}\ln (1-\frac{C_{\mathrm{\infty }}}{N})-\frac{N}{R_0}\ln (\kappa )=0.$(2) ) as $N=C_{\mathrm{\infty }}+\kappa N\exp (-R_0{C}_{\mathrm{\infty }}/N),$ where the right-hand-side is the sum of two non-negative monotonic functions of $C_{\mathrm{\infty }}$, one increasing and the other decreasing, and since $0<\kappa <1,$ one can see that such a solution exists and is unique. Asymptotically, when the disease has followed its course and all infections have occurred, I = 0 and $C=R.$ Then, $S=\kappa N\exp (-R_{0}R/N)=\kappa N\exp (-R_{0}C/N)$ and the conservation law N = S + I + R becomes $N=\kappa N\exp (-R_{0}C/N)+C.$ In other words, $C_{\mathrm{\infty }}$ defined in Equation (2(2) $C_{\mathrm{\infty }}+\frac{N}{R_0}\ln (1-\frac{C_{\mathrm{\infty }}}{N})-\frac{N}{R_0}\ln (\kappa )=0.$(2) ) is simply the final (asymptotic) number of cases associated with the outbreak. This concludes the proof of Theorem 1.1.

To see that the entire dynamics of an outbreak in the SIR model is captured by the associated ICC curve, note that knowledge of $C(t)$ leads to knowledge of $S(t)=N-C(t),R(t)=-\frac{N}{R_0}\ln \left(\frac{1}{\kappa }(1-\frac{C(t)}{N})\right),\mathrm{and}I(t)=C(t)-R(t).$In other words, the time course of S, I, and R may be obtained from the solution of the differential equation $\mathrm{d}C/\mathrm{d}t=G_{\kappa ,N}(C)$ prescribed by the ICC curve, with appropriate initial conditions (which in turn define the value of κ).

2.2. Robustness to under- or over-reporting

The quantity N appearing in Equation (1(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(1) ) refers to the total population taken into account in the model, which is conserved by the dynamics. In case of over- or under-reporting, the observed incidence expectation ${\mathcal{I}}_o$ is a fraction of the actual incidence $\mathcal{I}$, so that ${\mathcal{I}}_o=\alpha \mathcal{I},$ where we assume that α is constant or varying slowly enough for its derivative to be negligible. The ICC curve estimated from case report data will then be a graph of ${\mathcal{I}}_o$ as a function of ${\displaystyle C_o={\int }_0^t{\mathcal{I}}_o(s)\mathrm{d}s=\alpha C.}$ We thus have $\begin{array}{rl}\frac{\mathrm{d}{\mathcal{I}}_o}{\mathrm{d}t}& =\alpha \frac{\mathrm{d}\mathcal{I}}{\mathrm{d}t}=\alpha G_{\kappa ,N}(C)=\alpha \beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N})\\ & =\beta (C_o+\frac{\alpha N}{R_0}\ln (1-\frac{C}{N})-\frac{\alpha N}{R_0}\ln (\kappa ))(1-\frac{C}{N})\\ & =\beta (C_o+\frac{N_o}{R_0}\ln (1-\frac{C_o}{N_o})-\frac{N_o}{R_0}\ln (\kappa ))(1-\frac{C_o}{N_o})=G_{\kappa ,N_o}(C_o),\end{array}$where $N_o=\alpha N$. Therefore, the ICC curve in the presence of over-reporting ($\alpha >1$) or under-reporting ($\alpha <1$) has the same functional form, with the same parameter values, as the true ICC curve, provided $\mathcal{I}$, C, and N are rescaled by the factor α. For the same reason, $\kappa =S(0)/N$ is also independent of the value of α.

While systematic under- or over-reporting does not affect the ICC curve provided the size of each compartment is appropriately rescaled, the reproduction ratio depends on N through the ratio $C_{\mathrm{\infty }}/N$, as can be seen from Equation (2(2) $C_{\mathrm{\infty }}+\frac{N}{R_0}\ln (1-\frac{C_{\mathrm{\infty }}}{N})-\frac{N}{R_0}\ln (\kappa )=0.$(2) ), in which setting $\kappa =1$ leads to (3) $R_0\frac{C_{\mathrm{\infty }}}{N}+\ln (1-\frac{C_{\mathrm{\infty }}}{N})=0⟹R_0=-\frac{\ln (1-C_{\mathrm{\infty }}/N)}{C_{\mathrm{\infty }}/N}.$(3) Finally, rescaling time amounts to rescaling β and γ, and consequently the ICC curve.

3. Parameter identification

Structural identifiability (see e.g. [14] 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)=\mathrm{d}C/\mathrm{d}t=G_{\kappa ,N}(C(t)).$ From Equation (1(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(1) ), it is clear that knowledge of $C(t)$ and $\mathcal{I}(t)$ uniquely determines β, γ, and κ. 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 $0=({\beta }^{(1)}-{\beta }^{(2)})\frac{C}{N}+({\gamma }^{(1)}-{\gamma }^{(2)})\ln (1-\frac{C}{N})-({\gamma }^{(1)}\ln ({\kappa }^{(1)})-{\gamma }^{(2)}\ln ({\kappa }^{(2)})),$for all values of $\frac{C}{N}\in [0,\frac{C_{\mathrm{\infty }}}{N}]$, 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 [5, 6, 18]. If N is unknown, setting $\begin{array}{rl}0={\beta }^{(1)}\frac{C}{N^{(1)}}& -{\beta }^{(2)}\frac{C}{N^{(2)}}+{\gamma }^{(1)}\ln (1-\frac{C}{N^{(1)}})-{\gamma }^{(2)}\ln (1-\frac{C}{N^{(2)}})\\ & -({\gamma }^{(1)}\ln ({\kappa }^{(1)})-{\gamma }^{(2)}\ln ({\kappa }^{(2)})),\end{array}$for all values of $\frac{C}{N}\in [0,\frac{C_{\mathrm{\infty }}}{N}]$, 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. [14]). 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 [16, 18] 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 [18], the cumulative data is multiplied by $1+ϵ,$ where ϵ is normally distributed with zero mean, whereas in [16], Poisson-distributed noise of mean equal to the current incidence value [2] 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 [10]. 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}-\text{Poisson}(\mathcal{I})|/\mathcal{I}$ scales like $1/\sqrt{\mathcal{I}}.$ Indeed, $\begin{array}{rl}⟨\delta \mathcal{I}⟩& =\sum _{k=0}^{\mathcal{I}}\frac{{\mathcal{I}}^k{e}^{-\mathcal{I}}}{k!}\frac{\mathcal{I}-k}{\mathcal{I}}+\sum _{k=\mathcal{I}+1}^{\mathrm{\infty }}\frac{{\mathcal{I}}^k{e}^{-\mathcal{I}}}{k!}\frac{k-\mathcal{I}}{\mathcal{I}}\\ & =\sum _{k=0}^{\mathcal{I}}\frac{{\mathcal{I}}^k{e}^{-\mathcal{I}}}{k!}-\sum _{k=0}^{\mathcal{I}-1}\frac{{\mathcal{I}}^k{e}^{-\mathcal{I}}}{k!}+\sum _{k=\mathcal{I}}^{\mathrm{\infty }}\frac{{\mathcal{I}}^k{e}^{-\mathcal{I}}}{k!}-\sum _{k=\mathcal{I}+1}^{\mathrm{\infty }}\frac{{\mathcal{I}}^k{e}^{-\mathcal{I}}}{k!}=2\frac{{\mathcal{I}}^{\mathcal{I}}{e}^{-\mathcal{I}}}{\mathcal{I}!}.\end{array}$With Stirling's formula, we obtain $⟨\delta \mathcal{I}⟩\sim 2{\mathcal{I}}^{\mathcal{I}}{e}^{-\mathcal{I}}{\left(\sqrt{2\pi \mathcal{I}}{\mathcal{I}}^{\mathcal{I}}{e}^{-\mathcal{I}}\right)}^{-1}=\sqrt{\frac{2}{\pi \mathcal{I}}}\text{ as }\mathcal{I}\to \mathrm{\infty }.$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+ϵ)$ where ϵ 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 [2], 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(5) $\mathcal{I}(t_{i-\frac{1}{2}})=\frac{{\mathcal{C}}_i-{\mathcal{C}}_{i-1}}{{\delta }_t}=I_i,C(t_{i-\frac{1}{2}})=\frac{1}{2}({\mathcal{C}}_i+{\mathcal{C}}_{i-1})=C_i.$(5) ) below.

Figure 2. (Color online). Empirical ICC curves for two hypothetical outbreaks described by the SIR model with Poisson-distributed reporting noise. Incidence information is assumed to be available every ${\delta }_t$, with ${\delta }_t=1$ for $R_0=2$ and ${\delta }_t=2$ for $R_0=10.5$. In each panel, stars correspond to the simulated outbreak data plotted in the (C, $\mathcal{I}$) plane. Data points without added noise (circles) and the exact ICC curve (solid line) are also shown.

3.1. Fitting ICC curves to data

Given a simulated time series for the number of cases of an outbreak in a population of size N, we find the parameters of the ICC curve that best fits the data by least square minimization. This approach is different from parameter identification methods that have been proposed in the literature until now. Indeed, traditional methods aim to find the best combination of model parameters consistent with observed time series; they therefore require numerical integration of the compartmental model whose parameters need to be identified, in order to estimate and then minimize an associated cost function (see for instance [2, 16, 18]). Instead, the approach we put forward here consists in finding model parameters by directly fitting the relevant ICC curve to the data. We claim that this is a more robust approach because the functional that needs to be minimized has a unique extremizer in parameter space, which can be calculated exactly from the data points.

Given a time series of discrete incidences $\{{\mathcal{I}}_k,k=0,\dots ,M\},$ or of discrete cumulative cases $\{{\mathcal{C}}_k=\sum _{j=0}^k{\mathcal{I}}_j,k=0,\dots ,M\}$ recorded every ${\delta }_t$ units of time, we introduce a cost function that measures how much the data points depart from the ICC curve of the SIR model parametrized by β, γ, and $p=\ln (\kappa ).$ This function, ${\mathcal{E}}_e$, is defined by (4) ${\mathcal{E}}_e(\beta ,\gamma ,p)=\sum _{i=1}^M{(\frac{1}{N}\mathcal{I}(t_{i-\frac{1}{2}})-G(\beta ,\gamma ,p,C(t_{i-\frac{1}{2}})))}^2,$(4) where $G(\beta ,\gamma ,p,C)=(\beta \frac{C}{N}+\gamma \ln (1-\frac{C}{N})-\gamma p)(1-\frac{C}{N}),$and (5) $\mathcal{I}(t_{i-\frac{1}{2}})=\frac{{\mathcal{C}}_i-{\mathcal{C}}_{i-1}}{{\delta }_t}=I_i,C(t_{i-\frac{1}{2}})=\frac{1}{2}({\mathcal{C}}_i+{\mathcal{C}}_{i-1})=C_i.$(5) A Maple [13] calculation shows that the cost function ${\mathcal{E}}_e$ has a unique critical point, whose expression is given in Appendix. For outbreaks that start in a naive population, the formulas are hugely simplified by setting p = 0. The resulting cost functional then has a unique minimum, whose value may easily be calculated. With the following notation, ${\mathcal{J}}_i=\frac{I_i}{N},U_i=\frac{C_i}{N}(1-\frac{C_i}{N}),V_i=-(1-\frac{C_i}{N})\ln (1-\frac{C_i}{N}),$we have (after setting p = 0), $\frac{\mathrm{\partial }{\mathcal{E}}_e}{\mathrm{\partial }\beta }=-2\sum _{i=1}^M{U}_i({\mathcal{J}}_i-\beta U_i+\gamma V_i);\frac{\mathrm{\partial }{\mathcal{E}}_e}{\mathrm{\partial }\gamma }=2\sum _{i=1}^M{V}_i({\mathcal{J}}_i-\beta U_i+\gamma V_i).$Setting the right-hand sides to zero and solving the resulting linear system for β and γ leads to the following expression for the unique global minimizer of ${\mathcal{E}}_e$ when p = 0: $\begin{array}{rl}{\beta }_m& ={\displaystyle \frac{1}{h}(\left(\sum _{i=1}^M{V}_i^2\right)\left(\sum _{i=1}^M{U}_i{\mathcal{I}}_i\right)-\left(\sum _{i=1}^M{V}_i{\mathcal{I}}_i\right)\left(\sum _{i=1}^M{U}_i{V}_i\right)),}\\ {\gamma }_m& =\frac{1}{h}(\left(\sum _{i=1}^M{U}_i{V}_i\right)\left(\sum _{i=1}^M{U}_i{\mathcal{I}}_i\right)-\left(\sum _{i=1}^M{U}_i^2\right)\left(\sum _{i=1}^M{V}_i{\mathcal{I}}_i\right)),\end{array}$where $h=\left(\sum _{i=1}^M{U}_i^2\right)\left(\sum _{i=1}^M{V}_i^2\right)-{\left(\sum _{i=1}^M{U}_i{V}_i\right)}^2=\frac{1}{2}\sum _{i,j=1}^M(U_i{V}_j-V_i{U}_j{)}^2>0.$In what follows, we use the expressions given in Appendix, which can easily be coded up. We checked that for time series generated from an SIR model with $\kappa \simeq 1$, the simplified expressions above provide good estimates of the actual values of β and γ that were used to generate the reference outbreak. Importantly, we now have a means of associating a single set of values for the parameters β, γ, and $p=\ln (\kappa )$ to each outbreak time series.

3.2. Parameter identification in the presence of noise

We now use the expressions found in Section 3.1 to assess the robustness of the proposed parameter identification method to noise. As discussed above, the SIR model will be declared practically identifiable provided good estimates of the parameters may be obtained from noisy data. We therefore simulate 100,000 noisy time series from a ‘reference’ numerical integration of the SIR model, use the formulas of Appendix to calculate $p_c$, ${\beta }_c$, and ${\gamma }_c$ for each time series, and plot histograms of the estimated values of β and γ. We consider the two types of noise discussed above, both affecting the reported discrete incidence $I_i=\mathcal{I}(t_{i-\frac{1}{2}}).$ In Case A, we assume that at each time point $t_{i-\frac{1}{2}},$ the reported value of $\mathcal{I}$ has a Poisson distribution of mean ${\mathcal{I}}_{\mathrm{ref}}(t_{i-\frac{1}{2}}).$ In Case B, it is assumed that the reported value of $I_i$ at each time point is of the form $(1+ϵs_i){\mathcal{I}}_{\mathrm{ref}}(t_{i-\frac{1}{2}}),$ where ${\mathcal{I}}_{\mathrm{ref}}$ is the discrete incidence of the reference outbreak, ϵ is the noise amplitude, and $s_i$ is standard normal.

3.2.1. Case A: Poisson-distributed incidence

Figure 3 shows the distributions (in the form of properly scaled histograms) of ${\beta }_c,$ ${\gamma }_c,$ and $R_0=\frac{{\beta }_c}{{\gamma }_c}$ obtained from 100,000 noisy realization from a reference simulation of an SIR model with parameters $\beta =0.5$ and $\gamma =0.25$ ($R_0=2$). Figure 4 shows similar plots for $\beta =0.5$ and $\gamma \simeq 0.0476$, so that $R_0=10.5$. In each case, the size of the total population is N = 10, 000 individuals. The corresponding sample means and standard deviations are displayed in Table 1. These results indicate that the respective means of ${\beta }_c$ and ${\gamma }_c$ provide very good estimates of the parameters β and γ. Accuracy is not as good for the larger value of $R_0$, whose distribution shows a fairly long tail. This is to be expected because the actual value of γ in this case is close to zero. Moreover, the number of data points away from the end points of the outbreak is much smaller for $R_0=10.5$ than for $R_0=2$ (see Figure 2), leading to a potential loss of information.

Figure 3. (Color online). Numerically estimated distributions of the parameter values ${\beta }_c$, ${\gamma }_c$, and $R_0={\beta }_c/{\gamma }_c$ for 100,000 outbreaks corresponding to noisy (Poisson-distributed) realizations of a reference outbreak with parameters $\beta =0.5$ and $\gamma =0.25$ ($R_0=2$). Corresponding means and standard deviations are displayed in Table 1.

Figure 4. (Color online). Same as Figure 3, but with parameters $\beta =0.5$ and $R_0=10.5$ ($\gamma \simeq 0.0476$). Corresponding means and standard deviations are displayed in Table 1.

Table 1. Empirical mean (μ) and standard deviation (σ) of estimated parameter values in the presence of Poisson-distributed noise, for two values of $R_0$.

	$R_0=2$		$R_0=10.5$
	μ	σ	μ	σ
${\beta }_c$	0.4994	0.0060	0.4936	0.0120
${\gamma }_c$	0.2497	0.0030	0.0495	0.0070
$R_0$	2.0001	0.0017	10.15	1.29

3.2.2. Case B: normally distributed noise

At each time point, we now multiply the reference incidence by $1+ϵs$, where s is standard normal and the amplitude ϵ takes on one of three values: 0.05, 0.15, or 0.25. This type of noise is less realistic for incidence reports, whose variability is expected to be well described by a Poisson distribution, as in Case A above. We however perform the present tests to illustrate the robustness of the parameter estimation method introduced in this article.

For $R_0=2$, the parameter estimation procedure is robust at all values of ϵ considered. Means and standard deviations of ${\beta }_c$, ${\gamma }_c$, and $R_0={\beta }_c/{\gamma }_c$ are displayed in Table 2. The corresponding distributions are shown in the top row of Figure 5. As expected, estimates of β and γ are not as accurate for simulations with $R_0=10.5$, but return mean values for β and γ that are close to the exact values. The distributions of ${\beta }_c/{\gamma }_c$ however show very long tails (see bottom row of Figure 5), due to values of ${\gamma }_c$ close to zero. In such a case, estimating $R_0$ as the ratio $<{\beta }_c>/<{\gamma }_c>$ however gives acceptable values, equal to 9.96, 9.19, and 7.53, for ϵ equal to 0.05, 0.15, and 0.25 respectively.

Figure 5. (Color online). Similar to Figures 3 (top row) and 4 (bottom row) but for normally distributed noise of variable amplitude, set at 0.05, 0.15, or 0.25. Corresponding means and standard deviations are displayed in Table 2.

Table 2. Empirical mean (μ) and standard deviation (σ) of estimated parameter values in the presence of normally distributed noise of mean zero and amplitude ϵ, for two values of $R_0$.

	$R_0=2$
	$ϵ=0.05$		$ϵ=0.15$		$ϵ=0.25$
	μ	σ	μ	σ	μ	σ
${\beta }_c$	0.4994	0.0048	0.5001	0.0152	0.5013	0.0267
${\gamma }_c$	0.2497	0.0024	0.2500	0.0076	0.2506	0.0135
$R_0$	2.0001	0.0009	2.0002	0.0028	2.0002	0.0050
	$R_0=10.5$
	$ϵ=0.05$		$ϵ=0.15$		$ϵ=0.25$
	μ	σ	μ	σ	μ	σ
${\beta }_c$	0.4939	0.0201	0.5053	0.0577	0.5396	0.0938
${\gamma }_c$	0.04961	0.0104	0.0550	0.0293	0.0717	0.0465
$R_0$	10.3682	2.14	27.0950	769.77	29.9396	863.86

In summary, the above numerical simulations indicate that the ICC-based parameter estimation method presented here provides good results in the presence of reporting noise. Not surprisingly, the quality of the estimates depends on the rate at which incidence is sampled and, in the case of normally distributed noise, on the amplitude of the noise.

4. Parameter identification as the outbreak unfolds

In principle, the formulas for ${\beta }_c$ and ${\gamma }_c$ given in Appendix may be used with any number of points M, as long as the data points represent the entire ICC curve. It 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 β, γ, and $R_0=2$ for the two types of reporting noise. In both cases, estimates of β and γ 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 β and γ (stars) is close to the actual value of $R_0$ near t = 12, fairly early in the outbreak.

Figure 6. (Color online). Estimates of β, γ and $R_0$ as the outbreak unfolds for Poisson-distributed noise (left) and normally distributed noise of size 0.15 (right), in a situation where $R_0=2$. Estimates of $R_0$ show the evolution of $<\beta /\gamma >$ (circles) and $<\beta >/<\gamma >$ (stars), where $<\cdot >$ indicates sample mean. Error bars correspond to one standard deviation on each side of the mean.

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 information. Parameter estimates for β and γ 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.

5. Application to outbreak data

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 [11] 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 15 April 2020.

5.1. Gastroenteritis outbreak in mallorca, Spain

We use the same data set as in [11], estimated from the epidemiological curve provided in [12]. We fit the ICC curve to the interpolated data (i.e. find ${\beta }_c(N)$, ${\gamma }_c(N)$, and $p_c(N)$) for a range of values of N, identify the value $N_m$ of N that best fits the data (i.e. for which the RMSE is minimum), and select a range of values of N near $N_m$ that give a RMSE within 2% of its minimum. Then, for each value of N in the selected range, we apply the parameter estimation procedure of Section 3, with 10,000 oubreak realizations, obtained from adding Poisson-distributed noise to the smoothed and interpolated data. Figure 7 shows the resulting parameter distributions (top panel; $N\in [93,186]$) and a plot of the ICC curve for the optimal parameter values (bottom panel; N = 120, $\beta =1.34$, and $\gamma =0.88$). The value of $R_0$ for the ICC curve displayed is 1.52 and corresponds to an attack rate of 59.8%. This is near the upper boundary of the 95% confidence interval (38.5–61.3%) for the overall attack rate, but well within the 95% confidence interval of the attack rate among nursing home residents (42.1–72.2%) reported in [12].

Figure 7. (Color online). Parameter estimates for a gastroenteritis outbreak in a nursing home in Mallorca, Spain [12]. The top panel shows histograms (scaled to represent probability distributions) of ${\beta }_c$, ${\gamma }_c$, and $R_0$, for $N\in [93,186]$. The bottom panel shows the ICC curve for the optimal parameter values, corresponding to a value of $R_0=1.52$.

5.2. COVID-19

Figure 8 shows ICC curves for the first waves of COVID-19 outbreaks in a few countries. The data was obtained from the World Health Organization daily situation reports [19] as well as from the Wuhan Health Municipal Commission [21]. For Hubei Province, only laboratory confirmed cases were used. Consequently, daily increment values for 17–19 February 2020 (days 65–67 in the data set, when China reported combined laboratory and clinically confirmed cases for a brief period of time) were set to half of the reported number of confirmed cases; the number of cumulative cases was then obtained by summing daily reports of laboratory confirmed cases. For the Republic of Korea, only the first 58 points in the data set were used, since they correspond to the first wave of the outbreak. The third example is for a country (France), where the outbreak was not over as of 15 April 2020. For each of these examples, we show the ICC curve for the set of parameters that best fits the interpolated data. The optimal values of N are N = 51, 160 for Hubei Province, N = 10, 282 for the Republic of Korea, and N = 161, 142 for France. The corresponding optimal basic reproduction ratio values are $R_0=2.74$ for Hubei Province, $R_0=2.13$ for the Republic of Korea, and $R_0=1.96$ for France. A histogram of estimated values of $R_0$, scaled to represent a probability distribution function, is shown in the right panel of each row of Figure 8. These were obtained as described above in the case of the outbreak in Mallorca, except that no Poisson-distributed noise was added to the data. This is because incidence values are generally large and that, as discussed in Section 3, the relative effect of the added noise is minimal, leading to a level of variability of $R_0$ that is insignificant compared to the effect of changing N. Values of $R_0$ in the 2–3 range are consistent with early estimates of the basic reproduction ratio of COVID-19 published in the literature [20, 22].

Figure 8. (Color online). Optimal ICC curves and estimated distributions of $R_0$ values for early COVID-19 outbreaks in Hubei Province (top panel), the Republic of Korea (middle panel), and France (bottom panel).

The time course of the COVID-19 outbreak in Hubei Province, as described by the SIR model with optimal parameters identified by the present method ($\beta =0.401$, $\gamma =0.146$ ($R_0=2.74$), and N = 51, 160) is shown in Figure 9, together with the reported data. The solid curves are obtained by integration of the ICC curve, with initial conditions corresponding to day 10 of the outbreak, and alignment with the data at day 45. Similar plots for the other outbreaks discussed in this section are provided as supplementary information. The good visual agreement between simulation and data indicates that the SIR model, and therefore the ICC curve approach presented in this manuscript, are able to capture the overall dynamics of real-life outbreaks.

Figure 9. (Color online). Cumulative cases (left) and incidence (right) as functions of time for the COVID-19 outbreak in Hubei Province, China. The solid curves represent the predictions of the optimal ICC curve and the open circles are the data points. Time is measured in days from 14 December 2019 (day 0).

6. Conclusions

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 graph that crosses the horizontal axis near the origin and at the final value of C.

The formulas presented in Section 2, which extend the parabolic approximation suggested in [11], 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 [18] when fitting synthetic prevalence time series.

Because of the simplicity of Equation (1(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(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 [7] 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 [15]. 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 fits of the reported data by an ICC curve will produce similar values of $C_{\mathrm{\infty }}$, the final number of cases for the outbreak (see for instance the plots of Figure 8). Because of Equation (3(3) $R_0\frac{C_{\mathrm{\infty }}}{N}+\ln (1-\frac{C_{\mathrm{\infty }}}{N})=0⟹R_0=-\frac{\ln (1-C_{\mathrm{\infty }}/N)}{C_{\mathrm{\infty }}/N}.$(3) ), if $C_{\mathrm{\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 [11]) and numerical observations by the author suggest ICC curves of shape similar to Equation (1(1) $G_{\kappa ,N}(C)=\beta (C+\frac{N}{R_0}\ln (1-\frac{C}{N})-\frac{N}{R_0}\ln (\kappa ))(1-\frac{C}{N}),$(1) ) are 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.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

For codes and data sets used in this study, please see https://jocelinelega.github.io/EpiGro/.

Appendix. Critical parameter values

The function ${\mathcal{E}}_e$ defined in Equation (4(4) ${\mathcal{E}}_e(\beta ,\gamma ,p)=\sum _{i=1}^M{(\frac{1}{N}\mathcal{I}(t_{i-\frac{1}{2}})-G(\beta ,\gamma ,p,C(t_{i-\frac{1}{2}})))}^2,$(4) ) has a unique extremizer $({\beta }_c,{\gamma }_c,p_c)$ given by the following expressions: $\begin{array}{rl}{\beta }_c& =\frac{(-B_s{O}_s+F_s{A}_s)p_c^2+(2O_s{D}_s-F_s{N}_s-A_s{Q}_s)p_c+(N_s{Q}_s-O_s{P}_s)}{(A_s^2-B_s{L}_s)p_c^2+2(-A_s{N}_s+D_s{L}_s)p_c-L_s{P}_s+N_s^2},\\ {\gamma }_c& =\frac{(-A_s{O}_s+F_s{L}_s)p_c-L_s{Q}_s+N_s{O}_s}{(A_s^2-B_s{L}_s)p_c^2+2(-A_s{N}_s+D_s{L}_s)p_c-L_s{P}_s+N_s^2},\\ p_c& =\frac{(-N_s{Q}_s+O_s{P}_s)A_s+(L_s{Q}_s-N_s{O}_s)D_s+(-L_s{P}_s+N_s^2)F_s}{(A_s{O}_s-F_s{L}_s)D_s-(A_s^2-B_s{L}_s)Q_s-N_s(B_s{O}_s-F_s{A}_s)},\end{array}$where $\begin{array}{rl}{A}_s& =\sum _{i=1}^M\frac{C_i}{N}{(1-\frac{C_i}{N})}^2=\sum _{i=1}^M{U}_i{P}_i,B_s=\sum _{i=1}^M{(1-\frac{C_i}{N})}^2=\sum _{i=1}^M{P}_i^2,\\ D_s& =\sum _{i=1}^M\ln (1-\frac{C_i}{N}){(1-\frac{C_i}{N})}^2=-\sum _{i=1}^M{V}_i{P}_i,F_s=\sum _{i=1}^M\frac{I_i}{N}(1-\frac{C_i}{N})=\sum _{i=1}^M{\mathcal{J}}_i{P}_i,\\ L_s& =\sum _{i=1}^M{\left(\frac{C_i}{N}\right)}^2{(1-\frac{C_i}{N})}^2=\sum _{i=1}^M{U}_i^2,N_s=\sum _{i=1}^M{(1-\frac{C_i}{N})}^2\ln (1-\frac{C_i}{N})\frac{C_i}{N}=-\sum _{i=1}^M{U}_i{V}_i,\\ O_s& =\sum _{i=1}^M\frac{I_i}{N}(1-\frac{C_i}{N})\frac{C_i}{N}=\sum _{i=1}^M{\mathcal{J}}_i{U}_i,P_s={\displaystyle \sum _{i=1}^M{(\ln (1-\frac{C_i}{N}))}^2{(1-\frac{C_i}{N})}^2=\sum _{i=1}^M{V}_i^2,}\\ Q_s& =\sum _{i=1}^M\frac{I_i}{N}\ln (1-\frac{C_i}{N})(1-\frac{C_i}{N})=-\sum _{i=1}^M{\mathcal{J}}_i{V}_i,\end{array}$and we have used the notation ${\mathcal{J}}_i=\frac{I_i}{N},P_i=1-\frac{C_i}{N},U_i=\frac{C_i}{N}(1-\frac{C_i}{N}),V_i=-\ln (1-\frac{C_i}{N})(1-\frac{C_i}{N}).$