Issue 
Math. Model. Nat. Phenom.
Volume 15, 2020
Coronavirus: Scientific insights and societal aspects



Article Number  28  
Number of page(s)  12  
DOI  https://doi.org/10.1051/mmnp/2020011  
Published online  21 April 2020 
Dynamics of COVID19 pandemic at constant and timedependent contact rates
^{1}
Department of Biosystems and Soft Matter, Institute of Fundamental Technological Research, Polish Academy of Sciences,
02106
Warsaw, Poland.
^{2}
Faculty of Mathematics, Informatics and Mechanics, University of Warsaw,
02097
Warsaw, Poland.
^{*} Corresponding author: tlipnia@ippt.pan.pl
Received:
25
March
2020
Accepted:
5
April
2020
We constructed a simple Susceptible−Exposed–Infectious–Removed model of the spread of COVID19. The model is parametrised only by the average incubation period, τ, and two rate parameters: contact rate, β, and exclusion rate, γ. The rates depend on nontherapeutic interventions and determine the basic reproduction number, R_{0} = β/γ, and, together with τ, the daily multiplication coefficient in the early exponential phase, θ. Initial R_{0} determines the reduction of β required to contain the spread of the epidemic. We demonstrate that introduction of a cascade of multiple exposed states enables the model to reproduce the distributions of the incubation period and the serial interval reported by epidemiologists. Using the model, we consider a hypothetical scenario in which β is modulated solely by anticipated changes of social behaviours: first, β decreases in response to a surge of daily new cases, pressuring people to selfisolate, and then, over longer time scale, β increases as people gradually accept the risk. In this scenario, initial abrupt epidemic spread is followed by a plateau and slow regression, which, although economically and socially devastating, grants time to develop and deploy vaccine or at least limit daily cases to a manageable number.
Mathematics Subject Classification: 92D30
Key words: basic reproduction number / novel coronavirus
© The authors. Published by EDP Sciences, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
COVID19 is a rapidly spreading disease caused by a novel coronavirus, SARSCoV2, 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 COVID19 a pandemic.
To characterize the dynamics of COVID19 spread and be able to consider hypothetical longterm 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, . The model is applied to estimate the fold change of β required to contain the epidemic at its early stage. Since initial governmentmandated 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 longterm scenarios, in which we assume timedependent β that reflects anticipated changes in social behaviour rather than the impact of administrativelyimposed 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_{1}, …, E_{5} 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 , 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): (2.1)
which, when the number of intermediate exposed steps is very large (formally, when k → ∞), can be replaced by (2.2)
In the subsequent analysis we will use the daily multiplication coefficient: (2.3)
that describes the daytoday 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
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: (2.4)
Although the daily multiplication coefficient, θ, is not a function of solely (Fig. 2), the value of coefficient is critical for assessment of propagation of the epidemic: implies exponential progression, 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,
In Figure 3A we show the dependence 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 . Since to terminate the exponential growth phase, has to be reduced to 1, Figure 3A shows to which extent 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 (Fig. 3A). Thus, to “stabilise” the epidemic, China had to reduce 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 and effective reduction of about 25 times. We should notice that several authors [9, 19], also based on epidemiological data from China, obtained much lower estimates of ; we will come back to this important issue in the Section 2.4.
The ~10day 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 from about 12 to 2.6, which is less than 5 times. To just contain the epidemic, that is, reduce to 1 and stabilise the number of daily cases, Italy should further reduce the contact rate more than twofold. 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 lastweek θ ≥ 1.19 (as of March 20, 2020), to contain the epidemic should reduce the contact rate more than threefold (in the case of the US, eightfold). In Figure 3B we give current values of θ and corresponding values of 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 . Dotted lines guide the eye to relate and β coming from data for China before and after quarantine (‘Q’) and for several European countries and the US based on lastweek data. (B) Epidemic dynamics in selected countries with high number of cumulative confirmed cases of COVID19 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 lastweek 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 is computed assuming that 1∕γ is in between 3 and 5 days. 
2.4 Dependence of the basic reproduction number, , on the number of exposed intermediates, k
Value of for COVID19 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 (95% CI: 1.15–4.77), whereas Tang et al. [19] estimated that (95% CI: 5.71–7.23). Both estimates were based on SEIRtype 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 Erlangdistributed 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 estimated from θ. In Figure 4, we show 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 predicted using a model with one intermediate state is more than 1.5 times smaller than 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 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, , 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 longterm 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 longterm 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 “nonquarantined” 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 longterm quarantine.
Figure 5 Longterm 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 timedependent 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 administrativelyimposed 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 handwashing 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 (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 culturedependent. In Figure 6A we show longterm 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 β_{0}, 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 β_{0} × 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 (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 longterm 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 Longterm epidemic dynamics with timedependent 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 lastweek θ 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 COVID19 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, COVID19 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: and , 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 about 25fold (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 from about 12 to 2.6; thus, further ~3fold reduction is necessary to contain the exponential epidemic spread.
Higher reduction of 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 would allow only to limit growth of daily new cases, but further fourfold 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 weektomonth 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_{1} (t) + ⋯ + E_{k}(t) + I(t). The epidemic growth that is associatedwith large 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_{1}(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 .
Of concern is COVID19 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 θ^{n} 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^{10} and 1.28^{14}, 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 longerterm 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 50million 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 noncontainment 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:
where k = 5 is the assumed number of intermediate exposed states and f_{S} = S(t)∕N with constant N = S(t) + E_{1}(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:
with R_{0} being the initial cumulative number of confirmed cases and
where α is given implicitly: (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_{0}, R(t = 0) = R_{0}, and , where R_{0} is the given initial cumulative number of confirmed cases and E_{i0}(R_{0}) and I_{0}(R_{0}) are computed according to 4.3, assuring that the initial condition lies on the exponential trajectory. In the model variant with timedependent 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 threecompartment “SER” model with delay τ (equal to the average incubation period):
with f_{S}(t) = S(t)∕N where N = S(t) + E(t) + R(t) = const. The delay differential equationsbased formulation (4.5a–4.5c) may admit negativevalued 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: (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, lastweek 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 preepidemic 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
 R.M. Anderson and R.M. May, Infectious Diseases of Humans: Dynamics and Control. Oxford University Press, New York (1991). [Google Scholar]
 V. Capasso and G. Serio, A generalization of the Kermack–McKendrick deterministic epidemic model. Math. Biosci. 42 (1978) 43–61. [Google Scholar]
 COVID19 Community Mobility Reports. Available at: https://www.google.com/covid19/mobility/ (2020). [Google Scholar]
 S. Funk, M. Salathe and V.A.A. Jansen, Modelling the influence of human behaviour on the spread of infectious diseases: a review. J. R. Soc. Interface 7 (2010) 1247–1256. [CrossRef] [PubMed] [Google Scholar]
 G. Goh, Epidemic Calculator. Available at: https://gabgoh.github.io/COVID/index.html (2020). [Google Scholar]
 R.J. Hatchett, C.E. Mecher and M. Lipsitch, Public health interventions and epidemic intensity duringthe 1918 influenza pandemic. Proc. Natl Acad. Sci. USA 104 (2007) 7582–7587. [CrossRef] [Google Scholar]
 How many tests for COVID19 are being performed around the world? Available at: https://ourworldindata.org/covidtesting (2020). [Google Scholar]
 W.O. Kermack and A.G. McKendrick, Contributions to the mathematical theory of epidemics—I. Proc. R. Soc. 115A (1927) 700–721. [Google Scholar]
 A.J. Kucharski, T.W. Russell, C. Diamond, Y. Liu, J. Edmunds et al., Early dynamics of transmission and control of COVID19: A mathematical modelling study. To appear in: Lancet Infect. Dis. (2020). doi: 10.1016/S14733099(20)301444 [Google Scholar]
 C.C. Lai, T.P. Shih, W.C. Ko, H.J. Tang and P.R. Hsueh, Severe acute respiratory syndrome coronavirus 2 (SARSCoV2) and coronavirus disease2019 (COVID19): the epidemic and the challenges. Int. J. Antimicrob. Agents 55 (2020) 105924. [CrossRef] [PubMed] [Google Scholar]
 S.A. Lauer, K.H. Grantz, Q. Bi, F.K. Jones, Q. Zheng et al., The incubation period of 2019nCoV from publicly reported confirmed cases: estimation and application. Preprint medRxiv doi: 10.1101/2020.02.02.20020016 (2020). [Google Scholar]
 Q. Li, X. Guan, P. Wu, X. Wang, L. Zhou et al., Early transmission dynamics in Wuhan, China, of novel Coronavirus infected pneumonia. N. Engl. J. Med. 382 (2020) 1199–1207. [Google Scholar]
 R. Li, S. Pei, B. Chen, Y. Song, T. Zhang et al., Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARSCoV2). Science (2020) eabb3221. [Google Scholar]
 N.M. Linton, T. Kobayashi, Y. Yang, K. Hayashi, A.R. Akhmetzhanov et al., Incubation period and other epidemiological characteristicsof 2019 novel coronavirus infections with right truncation: a statistical analysis of publicly available case data. J. Clin. Med. 9 (2020) 538. [Google Scholar]
 Mathematica notebook. Available at: http://pmbm.ippt.pan.pl/model/covid19 (2020). [Google Scholar]
 Novel coronavirus (COVID19) cases, provided by JHU CSSE. Available at: https://systems.jhu.edu/research/publichealth/ncov (2020). [Google Scholar]
 N. Perra, D. Balcan, B. Goncalves and A. Vespignani, Towards a characterization of behaviordisease models. PLoS One 6 (2011) e23084. [CrossRef] [PubMed] [Google Scholar]
 J. Riou, A. Hauser, M.J. Counotte and C.L. Althaus, Adjusted agespecific case fatality ratio during the COVID19 epidemic in Hubei, China, January and February 2020. Preprint medRxiv doi: 10.1101/2020.03.04.20031104 (2020). [Google Scholar]
 B. Tang, X. Wang, Q. Li, N.L. Bragazzi, S. Tang et al., Estimation of the transmission risk of the 2019nCoV and its implication for public health interventions. J. Clin. Med. 9 (2020) 462. [Google Scholar]
 X. Wang, W. Tian, X. Lv, Y. Shi, X. Zhou et al., Effects of Chinese strategies for controlling the diffusion and deterioration of novel coronavirusinfected pneumonia in China. Preprint medRxiv doi: 10.1101/2020.03.10.20032755 (2020). [Google Scholar]
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 . Dotted lines guide the eye to relate and β coming from data for China before and after quarantine (‘Q’) and for several European countries and the US based on lastweek data. (B) Epidemic dynamics in selected countries with high number of cumulative confirmed cases of COVID19 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 lastweek 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 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, , 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 Longterm 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 Longterm epidemic dynamics with timedependent 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 lastweek θ for Italy. In all panels, population size N = 50 million and R(t = 0) = 500. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.