Modelling fatality curves of COVID-19 and the effectiveness of intervention strategies

Data source

Data used in this study were downloaded from the database made publicly available by the Johns Hopkins University (JHU, 2020), which lists in an automated fashion the number of confirmed cases and deaths per country. We have also compared the JHU data with the corresponding data from Worldometer (2020) for eventual data inconsistency and data redundancy. In all cases considered here we have used data up to May 8, 2020. In the present study we considered the mortality data of COVID-19 from the following countries: China, Brazil, France, Germany, Iran, Italy, South Korea, and Spain.

Confirmed cases vs. mortality data

Because a large proportion of COVID-19 infections go undetected (Li et al., 2020), it is difficult to estimate the actual number of infected people within a given population. As many carriers of the virus are either asymptomatic or develop only mild symptoms, they will not be detected unless they are tested. In other words, the number of confirmed cases for COVID-19 is a poor proxy for the total number of infections. Furthermore, the fraction of confirmed cases relative to the total number of infections depends heavily on the testing policy of each country, which makes it problematic to compare the evolution of confirmed cases among different countries. In contrast, the number of deaths attributed to COVID-19 is a somewhat more reliable measure of the advance of the epidemic and its severity. The official numbers of deaths attributed to COVID-19 have, of course, uncertainties of their own, as countries may differ as to the criteria and protocols for recording deaths related to the disease. For instance, some countries’ death figures do not include (or only later in the epidemic started to include) deaths outside hospitals, which naturally leads to under-reporting. There may also be delays in the reporting of deaths, also leading to uncertainties as to the number of deaths at a given time. Furthermore, other factors, such as the age structure of a population and quality of care, may affect the fraction of deaths relative to the number of confirmed cases. Nevertheless, taking all these factors into consideration, it is still reasonable to assume that the evolution of the number of confirmed deaths bears a relation to the dynamics of the number of infections (Famulare, 2020). Under these circumstances and in the absence of more reliable estimates for the number of infected cases for COVID-19, we decided here to seek an alternative approach and model the fatality curves, defined as the cumulative number of deaths as a function of time, rather than the number of confirmed cases, as is more commonly done.

Mathematical models

We model the time evolution of the number of cases in the epidemic by means of the Richards growth model (RGM), defined by the following ordinary differential equation (ODE) (Wang, Wu & Yang, 2012; Hsieh, 2009): (1) $$\frac{dC}{dt}=rC(t)[1-{\left(\frac{C(t)}{K}\right)}^{\alpha }]$$where C(t) is the cumulative number of cases at time t, r is the growth rate at the early stage, K is the final epidemic size, and the parameter α measures the asymmetry with respect to the S-shaped dynamics of the standard logistic model, which is recovered for α = 1.

It is worth to point out that the Richards model has an intrinsic connection to the SIR epidemic model, see, for example, Wang, Wu & Yang (2012). As a matter of fact, by identifying the variable C(t) of the Richards model with the cumulative number of deaths of a modified SIRD model, akin to the SIR model of Wang, Wu & Yang (2012), it is possible to establish a sort of ‘map’ between the parameters (α, r) of the Richards model and the parameters (β, γ₁, γ₂) of the SIRD model, where β is the transmission rate, γ₁ is the recovery rate and γ₂ is the death rate (Macêdo et al., 2020). The advantage, however, of phenomenological models, such as the Richards model, is that they allow for exact solutions (Chowell et al., 2016), which makes the data analysis much simpler. Furthermore, by avoiding ‘the description of biological mechanisms that may be difficult to identify,’ especially in an ongoing epidemic, they ‘can be utilised for efficient and rapid forecasts with quantified uncertainty’ (Bürger, Chowell & Lara-Daz, 2019). It is also worth pointing out that phenomenological growth models have been successfully applied to other epidemics, such as Zika (Chowell et al., 2016) and influenza (Bürger, Chowell & Lara-Daz, 2019), which makes these models good candidates for describing the ongoing COVID-19 epidemic, where there is still substantial uncertainty in the epidemiological parameters.

As already mentioned, in the present article we shall apply the RGM to the fatality curves of COVID-19, so that C(t) will represent the cumulative numbers of deaths in a given country at time t, where t will be counted in days from the first death. Nevertheless, in the interest of generality, in this and in the following section we shall refer to C(t) simply as the number of cases.

Equation (1) must be supplemented with a boundary condition. Here it is convenient to choose (2) $$\ddot{C}(t_c)=0$$for some given t_c, where dots denote time derivatives, that is $\ddot{C}(t)=d^{2}C(t)/d{t}^2$. A direct integration of Eq. (1) subjected to condition Eq. (2) yields the following explicit formula: (3) $$C(t;r,\alpha ,K,t_c)=\frac{K}{{\{1+\alpha \exp [-\alpha r(t-t_c)]\}}^{1/\alpha }}$$where in the left-hand side we have explicitly denoted the dependence of the solution of the RGM on the four parameters, namely r, α, K, and t_c. In fitting Eq. (3) to empirical data it is convenient to set C(0) = C₀, where C₀ is the number of deaths recorded at the first day that a death was reported. Using that $C(0)=K/{\left[1+\alpha \exp \left(\alpha r{t}_c\right)\right]}^{1/\alpha }$, we can eliminate t_c in favour of the other parameters, so that we are left with only three free-parameters, namely (r,α, K), to be numerically determined.

We note, however, that the RGM is not reliable in situations where the epidemic is in such an early stage that the available data is well below the estimated inflection point t_c, that is, when the epidemic is still in the so-called exponential growth regime (Wu et al., 2020). In this case, it is preferable to use the so-called generalised growth model (Wu et al., 2020; Chowell, 2017), which is defined by the following ODE: (4) $$\frac{dC}{dt}=r{[C(t)]}^q$$where the parameter q provides an interpolation between the sub-exponential regime (0 < q < 1) and the exponential one (q = 1). The solution of Eq. (4) is (5) $$C(t;r,q,A)={[A+(1-q)rt]}^{1/(1-q)}=A^{1/(1-q)}{e}_q(rt/A)$$where the function e_q(x) = [1 + (1 − q)x]^{1/(1 − q)} is known in the physics literature as the q-exponential function (Tsallis, 1988; Picoli et al., 2009). Here the parameter A is related to the initial condition, that is A = C(0)^{(1 − q)}, but we shall treat A as a free parameter to be determined from the fit of Eq. (5) to a given dataset.

Intervention strategy and efficiency

We define an intervention strategy in the context of the RGM by considering that the corresponding measures, as applied to the actual population, induce at some time t₀ a change in the parameters of the model. In this way, the solution for the whole epidemic curve acquires the following piecewise form: (6) $$C(t)=\{\begin{array}{ll}C(t;r,K,\alpha ,t_c),& t\le t_0\\ C(t;r\mathrm{\prime },K\mathrm{\prime },\alpha \mathrm{\prime },t_c\mathrm{\prime }),& t>t_0\end{array}$$where we obviously require that K′/K < 1. Furthermore we impose the following boundary conditions at t₀: (7) $$C(t_0;r,K,\alpha ,t_c)=C(t_0;r\mathrm{\prime },K\mathrm{\prime },\alpha \mathrm{\prime },t_c\mathrm{\prime })$$ (8) $$\dot{C}(t_0;r,K,\alpha ,t_c)=\dot{C}(t_0;r\mathrm{\prime },K\mathrm{\prime },\alpha \mathrm{\prime },t_c\mathrm{\prime })$$

Note that condition Eq. (8) takes into account, albeit indirectly, the fact that intervention measures take some time to affect the epidemic dynamics. In other words, the trend (i.e. the derivative) one sees at a given time t reflects in part the measures taken at some earlier time (or lack thereof). Thus, imposing continuity of the derivative of the epidemic curve at time t₀ in our ‘intervention strategy’ seeks to capture this delay effect.

In our intervention model above, we assume that the net result of the intervention is to alter the parameters r and α of the RGM after the time t₀; see Eq. (6). In other words, we assume that the parameters of the RGM can capture (albeit in an effective and simplified manner) some basic aspects of the underlying epidemic dynamics, so that changes in the mechanisms of the disease propagation—owing, say, to the introduction of intervention measures—could be described in terms of variations in the model parameters. It is in this sense that we associate actual interventions with possible changes in the model parameters, assuming of course that the growth model is still valid after the interventions. Admittedly, the difficult part is to estimate how a particular set of intervention measures (e.g. social distancing, contact tracing and quarantine, school closures, etc.) would influence these parameters. This link between actual interventions and the RGM parameters cannot, of course, be obtained within the context of the Richards model alone. It requires, for example the use more complex models, such as agent-based or compartmental models, to implement specific interventions, after which one can use the RGM to fit the resulting epidemic curves; see, for example Chowell (2017) where a similar approach was adopted albeit in the context of quantifying the uncertainty in epidemiological parameter estimates. Comparison of the RGM parameters after the intervention with those for the reference curve (i.e. without intervention) could thus shed light on how actual interventions can be mimic within the RGM approach. We are currently pursuing this strategy—namely, using the RGM formula to fit simulations from both agent-based and SIR-type models—to gain a better understanding of how interventions can be reflected in the parameters of the RGM. Such an analysis, however, is still ongoing (Macêdo et al., 2020) and is beyond the scope of the present article.

Note that, as discussed above, the time t₀ is not the actual time of adoption of the intervention but rather the time after which the corresponding measures start to affect the epidemic dynamics, as reflected in a change in the evolution of the number of cases. Nevertheless we shall for simplicity refer to t₀ as the intervention ‘adoption time.’ We remark, furthermore, that as far as interventions go, there are basically two aims: (i) to reduce the speed of an epidemic, which is essentially to ‘flatten the curve’ of daily deaths and (ii) to reduce the epidemic size, that is the total number of deaths due to the epidemic. In this article we address the latter case, since the interventions modelled by the strategy Eq. (6) are designed to yield K′ < K. In other words, the aim here is to ‘bend’ the cumulative curve of deaths as soon as possible, so as to reach a lower plateau at the end of the epidemic (for more details, see “Discussion”).

As defined in Eq. (6), an intervention strategy adopted at time t₀ can be viewed as a map (r, K, α, t_c) → (r′, K′, α′, t_c′) in the parameter space of the RGM, which results in the condition K′/K < 1. Let us therefore define the intervention efficiency as the relative reduction (expressed in percentage) of the total number of cases: η(t₀) = (K − K′)/K, where it is assumed that η(t₀) > 0. Using conditions Eqs. (7) and (8) in Eq. (3), one obtains that (9) $$\eta (t_0)=1-\frac{y}{{[1-x(1-y^{\alpha })]}^{1/\alpha \mathrm{\prime }}}$$where y = C(t₀)/K and x = r/r′. Later we will exemplify the above measure of intervention efficiency, using as input the parameters r and α obtained from the fatality curve of the COVID-19 from Italy. This will allow us to investigate how the efficiency of different strategies (i.e. for different choices of r′ and α′) varies as a function of the adoption time t₀.

Statistical fits

All statistical fits in the article were performed using the Levenberg–Marquardt algorithm (Moré, 1978) to solve the corresponding non-linear least square optimisation problem. In the case of the Richards model, we set C(0) = C₀, where C₀ is the number of deaths recorded at the first day that a death was reported, so that there remain three parameters, namely (r, K, α), to be determined. For the q-exponential growth model we also need to determine three parameters (r, q, A). The fitting procedures were implemented in the opensource software QtiPlot, which was also used to produce the corresponding plots in Figs. 1–3. The plots in Fig. 4 were produced with the data visualisation library Matplotlib for Python.

Fatality curves

In Fig. 1 we show the official cumulative number of deaths (blue symbols) attributed to COVID-19 for China, where it is clearly visible the jump on day 84 from the first death, when the death toll was revised upward by almost 50%. Since officials in Wuhan, China, informed only that this increase ‘reflected updated reporting and deaths outside hospitals’ (BBC News, 2020), it does not seem possible to reconstruct the actual fatality curve for China; nor does it make much sense to fit any model to the data prior to this correction, owing to their unreliability. It is possible nonetheless to render the Chinese data amenable to a statistical fit, if only as a test of the model, if one somehow smooths the revised data. In this spirit, we adopted the following ad hoc procedure to obtain a smooth ‘empirical’ curve for China: we multiplied all data points prior to the revision date by a factor corresponding to the ratio between the numbers of deaths after and before the revision date, meaning that the excess deaths due to the data correction was uniformly redistributed by the same proportion for all dates before the revision. This is admittedly an arbitrary procedure (perhaps a worst case scenario in terms of change of shape of the unknown true curve), which is intended only as a numerical way of ‘welding’ the two sides of the curve at the jump. This ‘renormalised’ fatality curve (red circles) for China is shown in Fig. 1, superimposed with the corresponding statistical fit of the RGM curve (black solid line), where the fit parameters are shown in the inset of the figure. One sees from the figure that the fatality curve in this case, where the epidemic has apparently levelled off, is well adjusted by the RGM formula.

The good performance of the RGM for the renormalised data from China, exhibited in Fig. 1, encouraged us to apply the model to other countries at earlier stages of the epidemic. Here, however, care must be taken when estimating the model parameters from small time series, since it is well known that the Richards model (Wang, Wu & Yang, 2012; Hsieh, 2009) and its variants (Wu et al., 2020) are susceptible to the problem of over fitting, owing to the redundancy of the parameters. This may lead, for example, to estimation of certain parameters that are outside of biologically or otherwise reasonable ranges. For example, when applied to the number of infected cases in an epidemic, the parameter α should be constrained within the interval (0,1) (Wang, Wu & Yang, 2012). Here we apply the RGM instead to the number of deaths, but we assume that the same constraint should be observed. In other words, fits that return α outside this interval are disregarded as not reliable. Similarly, we restrict the acceptable values of r to the range (0,1), as we observed that values of r outside this interval tends to be an indication that the RGM is not quite suitable to fit the data. In other words, for our purposes here we assume that the restrictions 0 < r < 1 and 0 < α < 1 are useful empirical criteria for the validity of the RGM, which can also be verified from the map between the Richards model and the modified SIRD model (Macêdo et al., 2020). The unsuitability of the RGM is particularly evident when the available data does not encompass the inflection point t_c (Chowell, 2017; Wu et al., 2020), although as more data is accumulated the model is expected to become more accurate. As an empirical criteria, we thus consider here that the RGM is only acceptable if t_c is smaller than the time of the last data point; if this criteria is not fulfilled by a particular dataset, we then apply the generalised growth model as given in Eq. (5).

With these considerations in mind, we applied the RGM to the mortality data of COVID-19 from several other countries. In Fig. 2, we show the cumulative deaths for Italy, Spain, France, Germany, Iran, and South Korea, together with the respective RGM fits. Here again the RGM seems to be able to provide a reasonably good fit to the data for all the above countries. In the case of Brazil, where the epidemic curve has not yet reached its inflection point, the RGM is not justifiable. In such cases, it is more advisable to use a simpler growth model, such as the q-exponential. In Fig. 3 we show the fit of the q-exponential curve Eq. (5) to the Brazilian data, where one sees that the agreement is very good. From the fit we get q = 0.72, indicating that the fatality curve is in a sub-exponential growth regime.