© The authors. Published by EDP Sciences, 2020

1 Introduction

COVID-19 is a rapidly spreading disease caused by a novel coronavirus, SARS-CoV-2, transmitted between people via respiratory droplets. Typical symptoms are fever, dry cough, and shortness of breath. In severe cases, infection leads to pneumonia and acute respiratory distress syndrome. There is currently no vaccine or specific antiviral treatment. The virus was first reported in Wuhan, Hubei, China, in December 2019. The disease had spread widely over the province of Hubei, but was contained after imposing strict quarantine. Reported peak of the daily confirmed cases of about 4000 was reached on February 3, ten days after the quarantine had been imposed. The number of daily new cases decreased below 100 in early March. Currently (as of March 20, 2020), despite some efforts, nearly exponential growth with more than 20 000 cumulative cases is observed in Italy (nearly 50 000 cumulative cases), Spain, Germany, and the US The World Health Organization (WHO) declared COVID-19 a pandemic.

To characterize the dynamics of COVID-19 spread and be able to consider hypothetical long-term scenarios, we developed a simple Susceptible–Exposed–Infectious–Removed model. The model has only 3 parameters: average incubation period, τ, estimated to be approximately 5 days [11, 14], and two parameters that are affected by nonpharmaceutical protective measures: the infectious–susceptible contact rate, β, and the exclusion rate, γ. Isolation (quarantine) reduces β, while testing increases the rate at which infectious individuals are hospitalised or isolated, γ. These two rates determine the basic reproduction number, ${\mathcal{R}}_0=\beta /\gamma$. The model is applied to estimate the fold change of β required to contain the epidemic at its early stage. Since initial government-mandated restrictions appear to be insufficient to contain the spread of the disease at its early stage in European countries, we use the model to put forward long-term scenarios, in which we assume time-dependent β that reflects anticipated changes in social behaviour rather than the impact of administratively-imposed control measures.

2 Results

2.1 Model

We consider a variant of the Susceptible–Exposed–Infectious–Removed (SEIR) model, itself an extension of the commonly used SIR model [1, 8]. As intypical SEIR models, we assume that susceptible individuals get infected by infectious individuals at the rate proportional to the contact rate, β, and in this way become exposed (and infected), and later, after the latent period, become infectious. Ultimately, at the exclusion rate γ, infectious individuals turn into removed individuals (which is an umbrella class comprising confirmed cases currently quarantined, handled by the healthcare, recovered with acquired resistance, or the deceased). Removed individuals do not infect.

In contrast to most SEIR models, instead of a single exposed (infected) state, we have introduced multiple exposed states (E₁, …, E₅ in Fig. 1). With k = 5 exposed states such representation is equivalent to the assumption that the latent period is distributed according to the Erlang distribution with the “shape” parameter k = 5 and the rate parameter λ = k∕τ, where τ is the mean latent period. Since there are no available estimates of the mean latent period and the disease spreads mainly through airborne droplets disseminated by the coughing or sneezing infected persons, we assume that during the incubation period the infected individuals are effectively noninfectious and thus the mean latent period can be set equal to the estimated average incubation period, τ ≈ 5 days [11, 14].

The serial interval, defined as the time between onset of symptoms in a primary case and onset of symptoms in its secondary cases, that can be seen as an emergent property of the model, is distributed according to the hypoexponential distribution resulting from a series of k = 5 exponential distributions, each originating from one step of the incubation process, having mean times τ∕k, and one exponential distribution, originating from the exclusion process, with mean time 1∕γ. At τ∕k = 1 day and 1∕γ = 3 days, at which the rates of the hypoexponential distribution are $\text{Λ}=\{1,1,1,1,1,\frac{1}{3}\}$, the model implicates a distribution of the serial interval that agrees perfectly with the distribution obtained from epidemic data (see Fig. 2B in Ref. [12]).

All model parameters—β, γ, and τ (together with k, as demonstrated further)—are thus constrained by available data. Deterministic model dynamics is described by a system of 8 ordinary differential equations (4.1a)–(4.1e) given and analysed in more detail in Section 4.

Figure 1 Scheme of the model with notation guide and default parameter values. The model is fully characterized by parameters β, γ, and τ. The latent period, assumed to be substitutable with the incubation period, results from the inclusion of k = 5 intermediate exposed states and follows the Erlang distribution with average τ = 5 days. The serial interval follows the hypoexponential distribution with rates Λ = {k∕τ, k∕τ, k∕τ, k∕τ, k∕τ, γ}.

2.2 Dynamics of infection spread with respect to all model parameters

The exponential phase of the epidemic growth, during which the number of removed (confirmed) individuals grows as R(t) = R(t = 0) × exp(α t), can be analysed based on an implicit formula (see Sect. 4 for derivation): $$\frac{\beta }{\gamma +\alpha }={\left(1+\frac{\alpha \tau }{k}\right)}^k\text{,}$$(2.1)

which, when the number of intermediate exposed steps is very large (formally, when k → ∞), can be replaced by $$\frac{\beta }{\gamma +\alpha }=e^{\alpha \tau }\text{.}$$(2.2)

In the subsequent analysis we will use the daily multiplication coefficient: $$\theta :=e^{\alpha }\text{,}$$(2.3)

that describes the day-to-day increase of the new confirmed cases (for example, when the number of new cases is 25% higher than the number of new cases on the previous day, then θ = 1.25).

In Figure 2 we show the dependence of θ on parameters β, γ, and τ, according to (2.1). Either a decrease of the contact rate, β, or increase of the exclusion rate, γ, result in reduction of θ below 1 and containment of the epidemic. Practically, β is reduced by quarantine or isolation of individuals, which both lower the potential number of contacts per individual per day, while γ is increased by performing more diagnostic tests, which reduces the expected time during which an infected individual may infect susceptible ones. Coefficient θ decreases also with increasing τ (Fig. 2C), but this parameter is not controllable by protective measures.

Equations (2.1) and (2.2) imply that $$\{\begin{array}{ll}\alpha >0\text{ (equivalently, }\theta >1\text{)}\hfill & \text{for }\beta >\gamma \text{,}\hfill \\ \alpha =0\text{ (equivalently, }\theta =1\text{)}\hfill & \text{for }\beta =\gamma \text{,}\hfill \\ \alpha <0\text{ (equivalently, }\theta <1\text{)}\hfill & \text{for }\beta <\gamma \text{.}\hfill \end{array}$$

Further, parameters β and γ determine the expected number of susceptible individuals that will be infected by each infectious individual, which in the initial phase of theepidemic is expressed by the basic reproduction number: $${\mathcal{R}}_0:=\frac{\beta }{\gamma }\text{.}$$(2.4)

Although the daily multiplication coefficient, θ, is not a function of solely ${\mathcal{R}}_0$ (Fig. 2), the value of coefficient ${\mathcal{R}}_0$ is critical for assessment of propagation of the epidemic: ${\mathcal{R}}_0>1$ implies exponential progression, ${\mathcal{R}}_0<1$ implies exponential regression.

Figure 2 Dependence of the daily multiplication coefficient, θ, on all model parameters: (A) contact rate, β, (B) exclusion rate, γ, and (C) average incubation period, τ. Default parameter values are marked with dotted rectangles.

2.3 Relation of the daily multiplication coefficient, θ, and the basic reproduction number, ${\mathcal{R}}_0$

In Figure 3A we show the dependence $\theta ({\mathcal{R}}_0)$ for two values of γ, 1/(5 days) and 1/(3 days), and in the limit of γ →∞. The figure should be read as follows. An initial θ that can be determined based on data for a considered country or region corresponds to some value of ${\mathcal{R}}_0$. Since to terminate the exponential growth phase, ${\mathcal{R}}_0$ has to be reduced to 1, Figure 3A shows to which extent ${\mathcal{R}}_0$ should be reduced to contain the epidemic. The reproduction number can be reduced by either reducing the contact rate or increasing the exclusion rate, but the latter is problematic when the number of cases grows rapidly.

Assuming γ = 1/(3 days), which is the default value in the model, the data showing epidemic dynamics in China (Fig. 3B) can be interpreted as follows. Based on the initial exponential phase (trends in one or two weeks marked in Fig. 3B) one can estimate that θ = 1.4 which corresponds to ${\mathcal{R}}_0\simeq 8.6$ (Fig. 3A). Thus, to “stabilise” the epidemic, China had to reduce ${\mathcal{R}}_0$ by a factor of 8.6. However, as indicated by the fall of the number of daily new cases, China reduced θ to approximately 0.89 = θ′ (see Sect. 4), which implies $\mathcal{R}=0.34$ and effective reduction of ${\mathcal{R}}_0$ about 25 times. We should notice that several authors [9, 19], also based on epidemiological data from China, obtained much lower estimates of ${\mathcal{R}}_0$; we will come back to this important issue in the Section 2.4.

The ~10-day time lag between introduction of the quarantine and the time of the peak number of daily new cases can be attributed to: (1) the incubation period that causes that infected individuals are registered about τ = 5 days after infection and (2) the fact that the quarantine does not prevent mixing of infectious and susceptible individuals in households, adding next τ = 5 days to onset of the disease among family members. Unfortunately, this implies that even if strict quarantine is imposed immediatelyin European countries it will have a delayed effect.

In the case of Italy, “soft quarantine” imposed on February 21 in the most affected Italian province, Lombardy, allowed to reduce θ from the initial value of 1.5 (calculated based on the week of February 23–29, 2020) to the current value of θ_{last week} = 1.14 (calculated based on the last week, March 14–20, 2020). This implies reduction of ${\mathcal{R}}_0$ from about 12 to 2.6, which is less than 5 times. To just contain the epidemic, that is, reduce ${\mathcal{R}}_0$ to 1 and stabilise the number of daily cases, Italy should further reduce the contact rate more than two-fold. Figure 3A indicates that several other large European countries (Germany, France, Spain, and United Kingdom) as well as the US, each having thousands of cases and last-week θ ≥ 1.19 (as of March 20, 2020), to contain the epidemic should reduce the contact rate more than three-fold (in the case of the US, eight-fold). In Figure 3B we give current values of θ and corresponding values of ${\mathcal{R}}_0$ for the mentioned countries (estimated for 1∕γ = 3 days and 1∕γ = 5 days).

A much different dynamics can be observed in South Korea, which is exiting the exponential phase (Fig. 3B). Among many preventive countermeasures, until March 20, 2020, South Korea performed 316 000 tests [7], likely substantially increasing the exclusion rate. The example of South Korea is promising as it suggests that, in addition to extensive contact tracing and surveillance, development of easily accessible, fast, and accurate tests may help to contain the spread of the disease.

Figure 3 Daily multiplication coefficient predicted by the model and estimated from data. (A) Dependence of the daily multiplication coefficient, θ, on the basic reproduction number ${\mathcal{R}}_0=\beta /\gamma$. Dotted lines guide the eye to relate ${\mathcal{R}}_0$ and β coming from data for China before and after quarantine (‘Q’) and for several European countries and the US based on last-week data. (B) Epidemic dynamics in selected countries with high number of cumulative confirmed cases of COVID-19 based on data from Johns Hopkins University (JHU) Center for Systems Science and Engineering (CSSE) [16]. Value of initial θ for China is estimated based on 7 or 14 days of the early exponential phase; current θ is estimated for countries being in the exponential phase based on last-week data (7 days since March 14 to March 20, 2020). Introduction of quarantine in China (January 23) is marked with ‘Q’. The range of values of ${\mathcal{R}}_0$ is computed assuming that 1∕γ is in between 3 and 5 days.

2.4 Dependence of the basic reproduction number, ${\mathcal{R}}_0$, on the number of exposed intermediates, k

Value of ${\mathcal{R}}_0$ for COVID-19 has been already estimated in several important studies (parameters of which are conveniently summarized on the webpage of an online SEIR model simulator [5]). In particular, Kucharski et al. [9] estimated that ${\mathcal{R}}_0=2.35$ (95% CI: 1.15–4.77), whereas Tang et al. [19] estimated that ${\mathcal{R}}_0=6.47$ (95% CI: 5.71–7.23). Both estimates were based on SEIR-type models. The main difference between our model and these models is the inclusion of a cascade of five intermediate exposed states. This is equivalent to assuming that (1) the incubation period is Erlang-distributed with parameters k = 5 and λ = k∕τ, and (2) the serial interval follows a hypoexponential distribution (in agreement with Li et al. [12]). Importantly, the simplifying assumption that there is only one intermediate state is equivalent to the assumption that the incubation period follows an exponential distribution. A bit surprisingly, this simplification substantially modifies ${\mathcal{R}}_0$ estimated from θ. In Figure 4, we show $\theta ({\mathcal{R}}_0)$ for two values of γ: 1∕3 (default) and 1∕5 day^{− 1}, calculated assuming that the number of intermediate states is 1, 2, 3, 5, 10, or ∞ (the last case has been determined based on an equivalent model represented by delay differential equations, see Sect. 4). One can observe that for θ = 1.4, the value of ${\mathcal{R}}_0$ predicted using a model with one intermediate state is more than 1.5 times smaller than ${\mathcal{R}}_0$ predicted for the model with k =5 intermediate states (5.4 vs. 8.6 at γ = 1/3 day^{− 1}; of note, the discrepancy would be larger for larger θ; the ratio of ${\mathcal{R}}_0$ computed for two different fixed values of k does not depend on γ). As shown in Figure 4C, the exponential distribution of the incubation period with rate 1∕τ, or even Erlang(k = 2, λ = k∕τ), are not supported by data that indicate much narrower distribution with its probability mass separated from zero [11, 14]. In conclusion, we think that our estimate of the basic reproduction number is closer to its true value.

Figure 4 (A, B)Dependence of the daily multiplication coefficient, $\theta ({\mathcal{R}}_0)$, and (C) dependence of the distribution of the incubation period, τ, on the number of intermediate exposed states. In the model it is assumed that k = 5. In panel C), Erlang distribution with k = 1 is equivalent to the exponential distribution with rate parameter 1∕τ.

2.5 Hypothetical long-term scenarios

2.5.1 A scenario with constant contact rate, β

The above analysis of the epidemiological situation in several European countries and in the US suggests that the protective measures instituted in these countries failed to contain the epidemic. It is thus worth to consider long-term scenarios in which the cumulative number of removed cases, R, becomes comparable to the population size, N, substantially decreasing the susceptible fraction, f_S. In such case, initial exponential growth ceases, the number of daily new cases culminates in a maximum and then drops due to depletion of susceptible individuals.

When the reproduction number remains constant, the maximum number of daily new cases is fully determined by initial θ and exceeds 0.3 million for θ = 1.05 and 1 million for θ = 1.1 in a population of N = 50 million people (Fig. 5A). It is clear that even the lower of the two numbers exceeds the healthcare capacity of a country of tens of millions inhabitants such as Italy (~60 million), South Korea (~50 million), or Poland (~40 million). This, in turn, implies higher case fatality rate, which may approach the rate of severe cases (10–15% [10]). Correspondingly, with θ changing from 1.05 to 1.1, time to reach the maximum of daily new cases shortens nearly twice: from 7.5 to 4 months, assuming the initial number of removed cases, R(t = 0), equal 500. We should notice that achieving θ = 1.05 as well as θ = 1.1 require quarantine; in “non-quarantined” countries, such as currently France, Germany, Spain, United Kingdom, and the US, θ ≳ 1.2 (Fig. 3B).

In Figure 5B, we simultaneously show the number of the maximum weekly and monthly new cases and the time to reach the maximum. Parameter θ > 1.1 implies thatwithin one (worst) month more than half of population will be infected, whereas θ < 1.1 require long-term quarantine.

Figure 5 Long-term epidemic dynamics at a constant contact rate. (A) Daily new cases for two values of the initial daily multiplication coefficient, θ, and the initial number of removed (confirmed) cases, R(t = 0), equal 500. (B) Time to the worst day for three initial numbers of removed (confirmed) cases and new cases in the worst week and month. In both panels, population size N = 50 million,1∕γ = 3 days.

2.5.2 Scenarios with time-dependent contact rate, β

We consider the scenario with the constant reproduction parameter not realistic. Initially, when the number of cases is relatively low, the only way to reduce the contact rate is an administratively-imposed quarantine. However, one could expect that when the number of daily new cases exceeds a threshold, people will begin reducing their contacts voluntarily [2] or taking other protective measures such as careful hand-washing or physical distancing to reduce their number of effective contacts (reviewed in Ref. [4]). To account for this behavioural aspect, we assume that the contact rate is a decreasing function of daily new cases. We propose to assume that $$\beta ={\beta }_0\times \frac{H}{H+(\text{daily new cases})}\text{.}$$(2.5)

This simple formula implies, in particular, that when the number of daily new cases equals H, people reduce their contacts twice. The value of parameter H can be culture-dependent. In Figure 6A we show long-term epidemic dynamics assuming H = 10 000. In the considered scenario, the number of daily new cases relatively quickly reaches 15 000 and then decreases slowly (over years) with a decreasing susceptible fraction, f_S. The scenario shown in Figure 6A results in a manageable number of new cases and grants time to develop, produce, and distribute a vaccine. For a given β₀, the maximum number of daily cases is an increasing function of H, whereas the time in which the number of daily new cases decreases to the half of the maximum value, t_↘max∕2, decreases with H; a bit surprisingly, t_↘max∕2 depends mainly on β₀ × H (Fig. 6B).

Finally, we consider a scenario in which a secondary behavioural effect modifies the primary effect discussed above. As before, we assume that people reduce their contacts in response to an increasing risk, but in the course of their (self)quarantine they gradually accept an increasing risk [17]. (As an aside, when protective measures are lifted rapidly, such behaviour can lead to epidemic waves observed during the 1918–1919 A/H1N1 influenza pandemic [6].) Here we propose to model this effect by assuming that $$H(t)=H_0\times \frac{t+T_{\text{negl}}}{T_{\text{negl}}}\text{,}$$(2.6)

where t is (somewhat arbitrary) time from the beginning of the epidemic and T_negl is the time scale on which H (and risk acceptance, or “negligence”) doubles. In Figure 6C we show long-term epidemic dynamics assuming H = 10 000 and three values of T_negl: 3, 6, and 9 months. Interestingly, the dynamics is characterized by a relatively broad plateau, which is desirable due to limited capacity of the healthcare system.

Figure 6 Long-term epidemic dynamics with time-dependent contact rate. (A, B) Contact rate reduced in response to the surge in the daily number of cases. (C) Contact rate reduced in response to the surge in the daily number of cases and modulated by the duration of the quarantine. In panels A and C, initial θ = 1.14, which is equal to last-week θ for Italy. In all panels, population size N = 50 million and R(t = 0) = 500.

3 Discussion

In this study, we proposed a Susceptible–Exposed–Infectious–Removed model to analyse how the dynamics of the COVID-19 spread depends on the contact rate, β, and the exclusion rate, γ—two model parameters that can be modulated by nontherapeutic interventions. The remaining model parameter is the average incubation period, assumed equal 5 days [11, 14]. A decrease of the contact rate, that can be achieved by imposing quarantine, results in a decrease of the daily multiplication coefficient θ that describes growth of the cumulative as well as the daily new cases in the exponential phase. Coefficient θ can be reduced also by an increase of the exclusion rate that in turn depends on efficiency of the healthcare. Large exclusion rate is achieved by fast identification and isolation of infectious individuals. However, when the number of cases is very high, achieving large exclusion rate is nearly impossible; in this case, reduction of the contact rate remains the only option.

In the early phase, before preventive measures were introduced, COVID-19 had been spreading in China and then in Italy with θ = 1.4 and θ = 1.5, respectively. This allowed us to estimate the initial basic reproduction number in these countries: ${\mathcal{R}}_0\simeq 9$ and ${\mathcal{R}}_0\simeq 12$, respectively. Of note, even without prevention the reproduction number depends on numerous different factors such as sociological context, hygiene, and weather.

Data describing epidemic dynamics in China suggest that very strict quarantine and other preventive efforts imposed in the most affected China province, Hubei, [9, 20] allowed to reduce ${\mathcal{R}}_0$ about 25-fold (which is on par with the estimated reduction of transmissibility of 97–100% [18]). This allowed the Chinese to reduce coefficient θ from about 1.4 in the exponential growth phase to about θ′ = 0.89 in the regression phase. In Italy, however, in response to the mild restrictions, θ was reduced from about 1.5 in the early phase to 1.14 observed currently (March 14–20, 2020), which corresponds to the reduction of ${\mathcal{R}}_0$ from about 12 to 2.6; thus, further ~3-fold reduction is necessary to contain the exponential epidemic spread.

Higher reduction of ${\mathcal{R}}_0$ is required in four other considered European countries, France, Germany, Spain, and United Kingdom, with the number of cases exceeding 10 000 and current θ ≳ 1.2 (March 14–20, 2020). Such reduction of ${\mathcal{R}}_0$ would allow only to limit growth of daily new cases, but further four-fold reduction is necessary to bring the epidemic to the regression phase with θ′ = 0.89, as was reported in Hubei. Higher θ′∈ (0.89, 1.0) may also allow to eliminate the disease but the required quarantine time would be longer. It appears that most of European countries follow the relatively mild regulations adopted in Italy, which (as of March 20, 2020) proved to be insufficient to contain the epidemic, and, correspondingly, with a week-to-month delay, follow the Italian trajectory of the cumulative daily new cases.

It is important to notice that the current number of removed (confirmed) individuals may be much lower than the total number of current exposed and infectious individuals: E₁ (t) + ⋯ + E_k(t) + I(t). The epidemic growth that is associatedwith large ${\mathcal{R}}_0$ implies a proportional large number of not yet confirmed cases that in the next few days will contribute to the number of cumulative confirmed cases. The confirmed fraction f_confirmed(t) = R(t)∕[E₁(t) + ⋯ + E_k(t) + I(t) + R(t)], that remains constant in time in the exponential phase, can be calculated based on (4.3) and (4.4) as $f_{\text{confirmed}}=1/{\mathcal{R}}_0$.

Of concern is COVID-19 lethality, which in quantitative terms is usually expressed as the case fatality rate, assessed as the ratio of deaths to confirmed cases. This definition, used by WHO, works well in the steady statewhen the number of active cases remains constant or when historical data of previous epidemics are analysed. However, in the exponential growth phase, the case fatality rate underestimates the real rate of deaths due to the disease. When the average time between the onset of symptoms or positive test and death is n days, the current cumulative number of deaths should be divided by the cumulative number of cases recorded m days before, which is θⁿ times lower than the current number. This observation may to some extent explain puzzling large differences between fatality rates across European countries. In Germany, having the lowest fatality rate (among countries with a high number of cases), identification of people by testing at early stage of the disease and high standard of public healthcare may result in large n, which together with high θ = 1.28 (March 14–20, 2020) may result in a significant underestimation of the case fatality rate. Assuming n between 10 and 14 days, one can find that the case fatality rate may be underestimated by a factor between 1.28¹⁰ and 1.28¹⁴, that is, 12–32 times. Finally, we should notice that in noncontainment scenarios, the expected death toll may not be predicted based on the case fatality rate, because a significant number of cases may be asymptomatic (and some individuals could possibly turn out resistant) [13].

Since current data suggest that in Europe, US, and Iran the epidemic will not be contained in its early phase, we have analysed longer-term scenarios. The model with the constant contact and exclusion rates predicts that even for a relatively small exponential growth phase coefficient, θ = 1.05, the number of daily new cases (in a 50-million population) will exceed 0.3 million, well above the capacity of the healthcare system. Additionally, keeping the coefficient θ at such low level requires stringent limitation of personal contacts, while the time to the peak of the epidemic would be about 7 months. One of potential consequences of exceeding the healthcare system capacity is the increase of case fatality rate due to the lack of necessary medical staff and equipment.

We think that the more realistic scenarios are these in which the contact rate varies over time. In early phase of the epidemic, the contact rate may be reduced only by forcing people to stay home; in the latter phase, when the number of daily cases exceeds a threshold, people isolate themselves to reduce the risk. In such scenario, the number of daily new cases reaches peak proportional to an assumed “fear” threshold and then slowly decreases due to the decreasing fraction of susceptible individuals. Such scenario seems more realistic and, although devastating for both the economy and social life, grants time to develop and administer vaccine. Historical data on the 1918–1919 A/H1N1 influenza pandemic suggest also that this “fear” threshold may not be constant in time, because people suffering from prolonged quarantine may tend to accept higher risk. When this “negligence” effect is included, one may obtain trajectories for which fast growth is followed by a plateau and then relatively fast decrease of daily cases. A bit surprisingly, this scenario, although not resulting from centrally imposed preventive policies, may be the most plausible non-containment scenario balancing health and economic costs.

4 Methods

4.1 Model equations

The dynamics of the system depicted in Fig. 1 is governed by the following system of 8 ordinary differential equations: $$\begin{array}{lll}\frac{\text{d}S}{\text{d}t}\hfill & =\hfill & -\beta f_{\text{S}}(t)I(t)\text{,}\hfill \\ \frac{\text{d}{E}_1}{\text{d}t}\hfill & =\hfill & \beta f_{\text{S}}(t)I(t)-\frac{k}{\tau }{E}_1(t)\text{,}\hfill \\ \frac{\text{d}{E}_j}{\text{d}t}\hfill & =\hfill & \frac{k}{\tau }{E}_{j-1}(t)-\frac{k}{\tau }{E}_j(t)\text{ for }2\le j\le k\text{,}\hfill \\ \frac{\text{d}I}{\text{d}t}\hfill & =\hfill & \frac{k}{\tau }{E}_k(t)-\gamma I(t)\text{,}\hfill \\ \frac{\text{d}R}{\text{d}t}\hfill & =\hfill & \gamma I(t)\text{,}\hfill \end{array}$$

where k = 5 is the assumed number of intermediate exposed states and f_S = S(t)∕N with constant N = S(t) + E₁(t) + ⋯ + E_k(t) + I(t) + R(t) being the population size (with included deaths).

In the early phase of the epidemic, f_S ~ 1. Under the assumption that f_S = 1, the systems has the exponential solution: $$\begin{array}{lll}{E}_j(t)\hfill & =\hfill & E_{j0}\exp (\alpha t)\text{for}\text{1}\le \text{j}\le \text{k},\hfill \\ I(t)\hfill & =\hfill & I_0\exp (\alpha t),\hfill \\ R(t)\hfill & =\hfill & R_0\exp (\alpha t),\hfill \end{array}$$

with R₀ being the initial cumulative number of confirmed cases and $$\begin{array}{lll}{I}_0(R_0)\hfill & =\hfill & R_0\frac{\alpha }{\gamma }\text{,}\hfill \\ E_{j0}(R_0)\hfill & =\hfill & R_0\frac{\beta }{\gamma }\frac{\alpha \tau }{k}{\left(1+\frac{\alpha \tau }{k}\right)}^{-j}\text{ for }1\le j\le k\text{,}\hfill \end{array}$$

where α is given implicitly: $$\beta =\left(\gamma +\alpha \right){\left(1+\frac{\alpha \tau }{k}\right)}^k.$$(4.4)

When the population cannot be considered as entirely susceptible, that is, when N − S(t) ≪ N is not satisfied, the system becomes essentially nonlinear and has no analytical solutions. This system is solved numerically using Wolfram MATHEMATICA (Wolfram Research, Inc., Champaign, IL, US; notebook is available online [15]). All numerical simulations start from an initial condition in which E_i(t = 0) = E_i0, I(t = 0) = I₀, R(t = 0) = R₀, and $S(t\text{}=\text{}0)=S_0=N-{\sum }_{j=1}^k{E}_{j0}-I_0-R_0$, where R₀ is the given initial cumulative number of confirmed cases and E_i0(R₀) and I₀(R₀) are computed according to 4.3, assuring that the initial condition lies on the exponential trajectory. In the model variant with time-dependent contact rates, daily new cases are computed as γ I(t). Simulations of this model variant start from the same initial condition that lies on the exponential trajectory.

It is noteworthy that in the limit of k →∞ the considered model converges to the following three-compartment “SER” model with delay τ (equal to the average incubation period): $$\begin{array}{lll}\frac{\text{d}S}{\text{d}t}\hfill & =\hfill & -\beta f_{\text{S}}(t)E(t-\tau ),\hfill \\ \frac{\text{d}E}{\text{d}t}\hfill & =\hfill & \beta f_{\text{S}}(t)E(t-\tau )-\gamma E(t-\tau ),\hfill \\ \frac{\text{d}R}{\text{d}t}\hfill & =\hfill & \gamma E(t-\tau ).\hfill \end{array}$$

with f_S(t) = S(t)∕N where N = S(t) + E(t) + R(t) = const. The delay differential equations-based formulation (4.5a–4.5c) may admit negative-valued solutions in the epidemic regression phase.

4.2 Estimation of θ

Parameter θ for the initial exponential growth case is estimated based on JHU CSSE data [16] of cumulative cases as: $$\theta ={\left(\frac{\text{cumulative number of cases in day }{d}_{\text{last}}}{\text{cumulative number of cases in day }{d}_{\text{first}}}\right)}^{\frac{1}{\text{Δ}t}},$$(4.6)

where for China day d_first is January 20 and day d_last is January 26 or February 2, 2020, with Δt = d_last − d_first − 1 = 6 or 13, respectively. For other countries, last-week data are used: d_first is March 14 and d_last is March 20, 2020; Δt = 6.

It should be noticed that as long as the epidemic progresses exponentially, coefficient θ estimated based on the increase of cumulative cases R(t) and coefficient θ′ estimated from the number of daily new cases, R(t) − R(t − (1 day)), should be equal (which is obscured by unavoidable fluctuations in daily data). However, when the epidemic ceases, θ calculated based on R(t) reaches 1, while θ′ becomes smaller than 1, reflecting an exponential fall of the number of the daily new (as well as the active) cases. Parameter θ′ in the epidemic regression phase in China has been estimated as the decay rate of the daily new cases in the time span from February 2 to March 2, 2020.

Note added in the proof.

As of April 4, 2020, the nonessential industry lockdown in Italy, Spain, and France, which led to an unprecedented decrease of workplace traffic to 1/3 of pre-epidemic level [3], allowed to terminate the exponential growth phase. These three countries appear to have reached (at least local) maxima of daily new cases.

Acknowledgements

This study was supported by the National Science Centre (Poland) grant 2018/29/B/NZ2/00668.

References

All Figures

Figure 1 Scheme of the model with notation guide and default parameter values. The model is fully characterized by parameters β, γ, and τ. The latent period, assumed to be substitutable with the incubation period, results from the inclusion of k = 5 intermediate exposed states and follows the Erlang distribution with average τ = 5 days. The serial interval follows the hypoexponential distribution with rates Λ = {k∕τ, k∕τ, k∕τ, k∕τ, k∕τ, γ}.

In the text

	Figure 2 Dependence of the daily multiplication coefficient, θ, on all model parameters: (A) contact rate, β, (B) exclusion rate, γ, and (C) average incubation period, τ. Default parameter values are marked with dotted rectangles.
In the text

Figure 3 Daily multiplication coefficient predicted by the model and estimated from data. (A) Dependence of the daily multiplication coefficient, θ, on the basic reproduction number ${\mathcal{R}}_0=\beta /\gamma$. Dotted lines guide the eye to relate ${\mathcal{R}}_0$ and β coming from data for China before and after quarantine (‘Q’) and for several European countries and the US based on last-week data. (B) Epidemic dynamics in selected countries with high number of cumulative confirmed cases of COVID-19 based on data from Johns Hopkins University (JHU) Center for Systems Science and Engineering (CSSE) [16]. Value of initial θ for China is estimated based on 7 or 14 days of the early exponential phase; current θ is estimated for countries being in the exponential phase based on last-week data (7 days since March 14 to March 20, 2020). Introduction of quarantine in China (January 23) is marked with ‘Q’. The range of values of ${\mathcal{R}}_0$ is computed assuming that 1∕γ is in between 3 and 5 days.

In the text

	Figure 4 (A, B)Dependence of the daily multiplication coefficient, $\theta ({\mathcal{R}}_0)$, and (C) dependence of the distribution of the incubation period, τ, on the number of intermediate exposed states. In the model it is assumed that k = 5. In panel C), Erlang distribution with k = 1 is equivalent to the exponential distribution with rate parameter 1∕τ.
In the text

Figure 5 Long-term epidemic dynamics at a constant contact rate. (A) Daily new cases for two values of the initial daily multiplication coefficient, θ, and the initial number of removed (confirmed) cases, R(t = 0), equal 500. (B) Time to the worst day for three initial numbers of removed (confirmed) cases and new cases in the worst week and month. In both panels, population size N = 50 million,1∕γ = 3 days.

In the text

Figure 6 Long-term epidemic dynamics with time-dependent contact rate. (A, B) Contact rate reduced in response to the surge in the daily number of cases. (C) Contact rate reduced in response to the surge in the daily number of cases and modulated by the duration of the quarantine. In panels A and C, initial θ = 1.14, which is equal to last-week θ for Italy. In all panels, population size N = 50 million and R(t = 0) = 500.

In the text