key: cord-0494538-hqi3w2a1
authors: Palomino, L. A. Rodr'iguez; Budini, Adri'an A.
title: Minimal epidemic model considering external infected injection and governmental quarantine policies: Application to COVID-19 pandemic
date: 2020-05-05
journal: nan
DOI: nan
sha: a5ba3391b36ba7e55331d6d08b0b06a25ac11910
doc_id: 494538
cord_uid: hqi3w2a1

Due to modern transportation networks (airplanes, cruise ships, etc.) an epidemic in a given country or city may be triggered by the arrival of external infected agents. Posterior government quarantine policies are usually taken in order to control the epidemic growth. We formulate a minimal epidemic evolution model that takes into account these components. The previous and posterior evolutions to the quarantine policy are modeled in a separate way by considering different complexities parameters in each stage. Application of this model to COVID-19 data in different countries is implemented. Estimations of the infected peak time-occurrence and epidemic saturation values as well as possible post-quarantine scenarios are analyzed over the basis of the model, reported data, and the fraction of the population that adopts the quarantine policy.

The COVID-19 pandemic started in the Chinese state Hubei at the end of 2019, and since then it propagated to almost all countries. Contrarily to an epidemic developing in a local closed region, the external propagation was triggered by the arrival of infected travelers that return from the other countries. In fact, at a continental level, after Asia, the epidemic propagated to Europe (with focus on Italy and Spain) and later to America (USA and Latino American countries) and Africa.

After the initial propagation, different administrations adopted quarantine policies. The main goal is to slow down the velocity of infected people increase, which in turn should avoid the collapse of the national health systems. This issue is of central importance because the disease mortality is decreased under appropriate medical attention.

At the present time, for most of the countries, the epidemic did not reach its maximal value. From a mathematical modeling perspective, the main challenge is to estimate the posterior epidemic evolution taking into account the quarantine policies and the available data. Magnitudes of special interest are the maximal (accumulative) number of infected agents as well as the epidemic peak and its time of occurrence. A qualitative description of possible post-quarantine scenarios is also of interest. In fact, while strict quarantine policies are effective for controlling the epidemic propagation, its sustainability and economical cost are aspects under discussion.

The mathematical modeling of epidemic propagation and evolution is a well established topic [1, 2] . One of the simplest descriptions is given by the so called Susceptible-Infectives-Recovered (SIR) models. The three variables * Electronic address: rodrigl@cab.cnea.gov.ar † Electronic address: adrianbudini@gmail.com of interest (S, I and R) count the number of persons in each possible state. A non-linear contribution (the product of S and I) governs the transition from the healthy to the infected group. There exists a vast and growing literature dealing with the COVID-19 pandemic description, which is based on SIR like models and related ones (see for example (non-exhaustive) Refs. [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] ). In contrast to standard SIR models, where saturation of an epidemic occurs due to a dominant number of recovered agents, here saturation and epidemic decay are mainly produced (in a first stage) by the quarantine policies.

The goal of this contribution is to formulate a simple and minimal epidemic mathematical model that takes into account both the external injection of infected agents as well as the quarantine government policies. In contrast to other approaches, we consider that the initial free epidemical growth and its posterior constrained evolution (due to quarantine) must be described in intrinsic separate ways. In fact, both stages are determined by different mechanisms whose underlying complexities may be quite dissimilar. The proposed evolution is defined with a minimal number of parameters whose estimation follows from data fitting and also public information. In the pre-quarantine period the dynamic is linear while in the post-quarantine stage the model can be related to a linearized SIR model. The approach is applied to available data of different European countries and Argentina. The maximal number of infected agents as well as the size of the epidemic peak are estimated in each case. A qualitative description of post-quarantine scenarios in terms of the fraction of isolated population is also proposed.

The manuscript is outlined as follows. In Sec. II we present the model and its derivation. In Sec. III it is applied in the description of COVID-19 data for different countries: Italy, Spain, France, United Kingdom, and Argentina. In Sec. IV post-quarantine scenarios are discussed. The Conclusions are provided in Sec. V. A mapping with SIR models is presented in the Appendix.

Due to the connectivity between diverse parts of the world (airplanes, cruise ships, etc.), in most countries the beginning of the epidemic was produced by travelers returning from an infected region. In this first period, where no local governmental policy has been taken yet, the growth of the (total) infected population n(t) results from external injection plus local contagions. These features depend on a complex way on flight connectivity and local passenger displacements, which in turn may lead to very different average behaviors. For example, the infected population may grow in an exponential (Latino American countries) or sub-exponential way (some European countries). Given the increase of infected people, in a second period local governments take different policies such as total or partial quarantine (avoiding local spread), jointly with the cancellation (partial or full) of infected external injection. In our approach, due to their intrinsic different origins, these two stages are modeled in a separate way.

At the beginning, the epidemic increase is produced by an external injection of infected agents, which in turn may produce a local spreading of the disease. At this stage, we propose the evolution

The function a(t) > 0 measures the velocity of external infected injection (i.e. airplanes). Furthermore, the function r(t) > 0 gives the rate of local (regional) contagions. This evolution is taken between the initial time t = 0, the day before the arrival of the first infected agent, up to the time t q where a government policy is taken. Time is measured in day. We analyze different possibilities for the functions a(t) and r(t).

One assumption is to consider that the local spreading of the disease gives a small or null contribution. Thus, we consider the functions

Here, the rate constants a and c have units 1/(day) ν , while ν > 0 is a dimensionless positive parameter. This functional form of a(t) is justified from the solution it gives rise. Eq. (1) can be integrated as n(t) = n(0) + t 0 dt ′ a(t ′ ), with n(0) = 0. Therefore,

Alternatively, we may take a nonvanishing contribution due to local propagation, r(t) = 0. We assume

where as before the rate constants a and c have units 1/(day) ν , while ν > 0 is again a dimensionless positive parameter. With these definitions, the solution of Eq. (1) reads [n(0) = 0]

We notice that the models (2) and (4) lead to exactly the same solutions, Eqs. (3) and (5). Hence, their underlying different mechanisms cannot be discerned from n(t). On the other hand, the solutions present the desired time behavior. In fact, in a short time regime n(t) ≃ (a/ν)t ν , in an intermediate regime,

, while in the posterior long time regime n(t) ≃ (a/c) exp(ct ν /ν). The parameter ν allows us to cover very different growing behaviors. A standard exponential growth corresponds to ν = 1.

One may argue that the equality of the previous solutions only applies to the special chosen functions a(t) and r(t) [Eqs. (2) and (4)]. This is certainly true. Nevertheless, in real experimental data it is hard to separate between imported and local contributions. Thus, while it is possible to chose other functional forms of a(t) and r(t), from the data is very hard to distinguish between them. Thus, we model the epidemic growth as

This function corresponds to Eq. (3) or Eq. (5). As a solution of compromise, the value of the parametersã,c, and ν are determined from the experimental data. These parameters take in an effective way both the external injection as well as the local propagation of the disease. In a roughly way, the parameter ν can be read as a measure of the flight and local spread networks complexities. The time t q is taken from public information corresponding to each country. Given that any change in social behavior does not occur instantaneously, this time has an intrinsic and unavoidable uncertainty.

Given the global network information structure, in diverse countries (or cities) the local governments dictated quarantine policies in order to change the free epidemic growth [Eq. (6) ] and try to get saturation in the total number of reported cases at the minimum number of cases possible. In order to take into account this underlying change in the epidemic dynamics, we split the total number of infected as

The term n R (t) gives the number of infected agents that recovered their health or died. Thus, n I (t) is the number of infected agents able to propagate the disease. For times t ≤ t q , for the data we analyze in the next section, it is possible to approximate

where n(t) is given by Eq. (6) .

After the quarantine, we assume the evolution

Here, the timeτ is a parameter to determine. It can be read as the average number of days that an infected agent is able to propagate the disease in the quarantine period. The rate of local contagions r q (t) > 0 is proportional to the number of susceptible population able to be infected in the quarantine period. We notice that the structure of Eq. (9) [with an exponential r q (t)] can be derived from a linearization of a non-linear SIR like model [see Appendix]. For n R (t) we take the evolution

which guarantees that in a long time limit all the population falls in this group. From Eqs. (9) and (10), and the relation (7), the total number of infected people evolves as

Given that quarantine leads to saturation of n(t), the function r q (t) must fulfill lim t→∞ r q (t) = 0. (12) In addition, at time t = t q we demand the continuity of the derivative of n(t), which leads to

where the left and right terms follow from Eqs. (6) and (11) respectively.

The post-quarantine dynamics is completely defined after knowing the rate r q (t). In order to fulfill condition (12) we propose the function

whereγ has units 1/(day) µ , while µ > 0 is a dimensionless positive parameter. As before, this parameter can be associated with different complexities of the postquarantine dynamics. Under the condition (13), taking

The solution of Eq. (9) can be written as

where the auxiliary functionn(t) readŝ

For shortening the expression we defined the ratesc ν ≡ ct ν−1 q , andγ µ ≡γt µ−1 q . In Eq. (15) the fitting parameters areγ,τ , and µ. The initial condition n I (t q ) = n(t q ) follows from the fitting function Eq. (6) .

The functional form of n I (t) always assumes an extremum maximal value, whose time of occurrence is denoted as t M , (d/dt)n I (t)| t=tM = 0. This time can be obtained in an analytical way in the case µ = 1. From Eq. (15) we get

The solution of Eq. (10) can be written as

where the initial condition is n R (t q ) = 0. In general, the integral can be performed in a numerical way [16] . Complementarily, a series solution can be obtained for µ = 1.

The total number of infected n(t) can be determined from Eq. (7) and the solutions defined by Eqs. (15) and (17) .

Its saturation value is lim t→∞ n(t) = lim t→∞ n R (t). For µ = 1, from Eq. (18) it follows lim t→∞ n(t) = n(t q )ec ν /γ 1−c ν /γ 1 +γτ +· · · , µ = 1. (19) Consistently, the saturation value depends on both the pre-and post-quarantine parameters.

In this section, on the basis of the proposed evolution, we analyze the (2020) COVID-19 pandemic in different countries. For times before and after quarantine the predicted behaviors are given by Eqs. (6) and (7) respectively. The number of infected and dead in each day are taken from "Our World in Data" [17] , while the number of recovered are taken from [18] , which are based on official data provided by each country. Furthermore, the (approximate) times of the quarantine onset t q are taken from public information. As a result of standard fitting techniques, in a first step, the data previous to the quarantine implementation allows us to obtain the parameters c,ã, and ν. In a second step, the parametersγ,τ , and µ are chosen for fitting the data after the quarantine imposition. The starting date of the epidemic is considered to be 12/31/19 in Wuhan-China, after which it becomes a pandemic spreading around the entire world. The data of disease reported in this work is until 04/11/2020.

We consider the European countries Italy, Spain, France, United Kingdom (UK), and Argentina. In these cases the quarantine has been applied strictly as a unified national policy. Among these countries, Italy and Spain reported collapse of the national health systems, which may be related with a delay in the implementation of the quarantine. In Table I (month/day) in each country is provided. It was taken as the date where an effective increase is observed in the reported data. Furthermore, the quarantine times t q were 6) and (7). The red and blue colors indicates times before and after the quarantine respectively. obtained from public information. We also include the results of the fitting parameters {c,ã, ν} associated to Eq. (6). Table II summarizes the fitting parameters {γ, τ , µ} associated to Eq. (7), with solutions (15) and (17).

For the European countries, in Fig. 1 it is shown the (available) data and the performed fits for the total number of infected n(t). Even in a logarithmic scale, the model provides a very good fitting of the reported data. From Table I , it follows that the dynamics of virus propagation in the first pre-quarantine stage [Eq. (6)] is not exponential, ν = 1. For Italy, Spain and France, the growth is proportional to t 1/2 , while for UK is proportional to t 0.78 . This last value is consistent with a very pronounced increase of the number of infections at the end of the data. On the other hand, in the post quarantine period we find that the best fitting of the available data is with µ ≈ 1 [ Table II ], corresponding to the solutions (15) and (18). This feature is valid for all countries, which shows that the epidemic after quarantine assumes the same dynamics.

The previous results provide a strong support to our modeling, which in a first stage is able to capture departures with respect to a pure exponential growth, as well as a posterior development of a universal dynamics.

In a linear time-scale, in Fig. 2 we show the predicted time evolution of n(t) and n I (t) jointly with the available data. It can clearly be seen that the quarantine turns out to be an effective policy that reduces infections considerably. In fact, with quarantine (full line) n I (t) develops a maximal value and departs considerably from the unconstrained (free) evolution (dotted line) Eq. (6) . For Italy, Spain and France, the peak of the epidemic is estimated to develop in 55-60 days after the epidemic beginning. Instead, for the UK case, the peak emerges after 90 days. These estimations follow from Eq. (16) . On the other hand, in all cases, the saturation value of n(t) can be very well approximated by the first term of Eq. (19) .

In Fig. 3 we show the Argentinean case, country that also follows a severe quarantine. Contrarily to the European cases, here the epidemic growth before the quarantine is exponential (ν = 1, Table I ). During the quarantine, a considerable decrease in the growing velocity is observed. Given that this country implemented the quarantine when the number of reported cases hadn't increase considerably, a small level of infection is predicted at the epidemic peak when compared with the previous countries. The maximal number of infected cases can also be very well approximated by the first term of Eq. (19) .

With the information of Tables I and II , it is possible to estimate an important parameter that characterizes the propagation dynamics of the disease, this is the so-called reproduction number R(t) ≡ n(t + α)/n(t) [Eq. (6) ]. The time interval α is the number of days for which R(t) acquires a certain value. In the case of Argentina, in the pre-quarantine period, taking for ex- ample R(t) = 2 (doubling the number of infections) it is possible to obtain α ≈c −1 ln(2) = 3.5 days, which is in accordance with the values of 3-4 days obtained from the infections reported at the beginning of the spread of the disease. Making a similar analysis in the postquarantine stage, the new value of α for which R(t) = 2 is given by α ≈γ −1 ln(2) = 12.6 days. This result is also in agreement with the values of 10-11 days obtained from reported data. The change in the value of R(t) (3.5 → 12.6 days), clearly reflects the beneficial effect of the quarantine in decreasing the spread of the epidemic. For European countries a similar analysis can be done.

The quarantine policy is one of the most appropriate mechanisms to control an epidemic in a social way. Nevertheless, this kind of solution cannot be maintained during large periods of time. Therefore, establishing simple criteria that allow softening tightness of the quarantine without implying a collapse of the health system, is a problem whose characterization is demanded.

In a roughly way, each quarantine scenario can be defined by the fraction f ∈ (0, 1) of the total population that is completely isolated. Thus, a quarantine softening can be related to a change in this parameter. In our minimal model, the parameters associated to the quarantine period are the time-dependent rate r q (t) and the timeτ [see Eqs. (9) and (10)]. Given thatτ defines the time scale for the transition n I (t) → n R (t), we consider that it has a weak dependence on f. Nevertheless, it has a strong regional dependence [see Table II ].

For the analyzed countries r q (t) is defined by Eq. (14) with µ = 1. Therefore, the unique parameter that can be related with f is the rate parameterγ, which defines the time scale of the decay of r q (t). As an ansatz, we relate both parameters as follows

Here, t I is the average number of days it takes a person to develop health signals of the disease. Consistently a higher f leads to a higherγ, which from Eq. (19) implies a lowering in lim t→∞ n(t). In Table II we summarize the values of f that follow from Eq. (20) taking t I = 14 days [19] . These values give a roughly estimation of the quarantine tightness implemented by each country. The relation (20) allows us to analyze in a qualitative way the changes in the epidemic dynamics under a governmental softening of the quarantine policy. Each scenario corresponds (or can be related) to a different value of f → f p < 1, which now measures departures (or relaxation) from a more strict quarantine policy, f p < f.

The inverse case f p > f corresponds to an increasing of the quarantine tightness. The problem is to estimate the new epidemic peak (if it develops) and new epidemic saturation value as a function of f p and the time t p at which the quarantine relaxation begins (t p > t q ).

We assume that for times t ≥ t p the evolution remains the same [Eqs. (9) and (10)], being defined with updated parameters. Thus, the expressions for n I (t) and n R (t) are given by Eqs. (15) (µ = 1) and (18) respectively, under the replacements

The replacementc ν →c ν e −γ(tp−tq) guarantees that the evolution remains the same when the translation t q → t p is introduced. In a second step, in all places the rate parameterγ is updated as follows

which introduces the change in the quarantine policy. It is simple to check that when f p = f, the epidemic dynamics remains unaltered or invariant. By introducing the replacements (21) and (22) into Eq. (16), after some algebra it is simple to obtain the time t Mp at which n I (t) may develop (or not) a second peak,

In fact, considering that the condition t p < t Mp must be fulfilled, a second peak emerges if the following inequality is satisfied

Thus, when this inequality is satisfied a second epidemic peak is observed in the interval t > t p . When the condition (24) is not satisfied the epidemic does not develop a second peak in the interval t > t p . After quarantine softening, the saturation value of the total number of infected follows again from the relation lim t→∞ n(t) = lim t→∞ n R (t). From Eq. (18), by using the update rules (21) and (22) we get

In order to enlighten the previous results, in Fig. (4) we analyze different possible post-quarantine scenarios for the Argentinean case [ Fig. 3 ]. The characteristic parameters before and after the quarantine are those of Table I and II respectively. In particular, f = 0.77.

In Figs. 4(a) to (d) we take a fixed t p = 60 days while f p = 0.6, 0.5, 0.4, and 0.3. As expected, in all cases, under the quarantine softening f → f p the maximal number [lim t→∞ n(t)] of total infected increases [Eq. (25)]. This change starts to be significant for smaller values of f p , where most of the population recovers its mobility. In addition, a diminishing in f p also leads to the appearance of extra epidemic peaks in the number of active infected n I (t), whose time of occurrence can be inferred from Eq. (23). Due to the condition (24), for the higher value of f p the extra peak does not develop.

Given the international network of public transportation, the epidemic growth in a given city or country may be triggered by external agents that arrive from an infected region. This first stage is not universal. After that, local governmental policies such as quarantine ones are taken in order to mitigate the epidemic growth. We presented a minimal model were all these components are taken into account, being defined by parameters that can be deduced from public information and fitting of the available reported data.

The first stage, that is the importation of external cases plus the beginning of local contagions, was modeled by a linear dynamics [Eq. (1) ]. The proposed solution [Eq. (6) ] is able to fit different sub, super, as well as standard exponential behaviors. The quarantine period was modelled by splitting the total number of infected in active and inactive (recovered and dead) [Eq. (7)], which dynamics can be read as a linearization of a SIR like model [Eqs. (9) and (10)].

The model was applied to the COVID-19 pandemic, analyzing reported data in different countries. For European countries such as Italy, Spain, France, and UK, the model provides a very good fitting of the available data. In addition, estimations for the corresponding epidemic peaks were obtained. Similarly to the case of Argentina, these countries established quarantine as a rigorous national policy. While in the pre-quarantine period a universal behavior is not observed, this is the case in the post-quarantine stage, where independently of the country the epidemic can be fitted with the same complexity parameter. These features support the splitting introduced in our model, as well as the estimations obtained from it.

The present approach also allows to analyze possible quarantine softening scenarios. The dynamics remains the same, being defined with a modified rate that depends on a factor that is proportional to the fraction of the total population under quarantine. The formalism furnishes roughly estimations for the changes in the epidemic peak and saturation values. Thus, we conclude that the present contribution may provide a valuable tool for taking relevant decisions about epidemic control. In particular, it allows to obtain qualitative estimations of the degree of quarantine softening in order to avoid a posterior collapse of a given health system.

Networks, and Introduction

Effective containment explains sub-exponential growth in confirmed cases of recent COVID-19 outbreak in Mainland China

Data analysis and modeling of the evolution of COVID-19 in Brazil

Modelling Recovered Cases and Death Probabilities for the COVID-19 Outbreak

Imperial College COVID-19 Response Team, Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand

Age profile of susceptibility, mixing, and social distancing shape the dynamics of the novel coronavirus disease 2019 outbreak in China

Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions

Evolving Epidemiology and Impact of Non-pharmaceutical Interventions on the Outbreak of Coronavirus Disease

Covid-19 -A simple statistical model for predicting ICU load in exponential phases of the disease

Mathematical modeling of the spread of the coronavirus disease 2019 (COVID-19) taking into account the undetected infections. The case of China

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (sars-cov2)

A modified SEIR model to predict the COVID-19 outbreak in Spain: Simulating control scenarios and multi-scale epidemics

Epidemic analysis of COVID-19 in China by dynamical modeling

Relaxing quarantine after an epidemic: A mathematical study of the spanish covid-19 case

For the analyzed data, we performed this calculus by using different integration algorithms

The authors thank Orlando Billoni for sending relevant references and information. Also to Mariano Bonifacio and Esteban Sorocinschi for a critical reading of the manuscript. L.A.R.P. and A.A.B. thanks support from CONICET, Argentina.

In the recent literature there are alternative generalizations of SIR models that take into account different governmental quarantine policies [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] . For example, from [3] , we writeS and I are the number of susceptible and infected persons respectively. The rate κ 0 measures the removed agents due to quarantine, while the rate κ measures the removed infected people. The constant α and β have the usual meaning in SIR models [1, 2] . The previous non-linear evolution can be linearized under the condition αSI ≪ κ 0 S. Thus, it follows that S ≃ s 0 exp(−κ 0 t). Under the mapping r q (t) = αs 0 e −κ0t andτ −1 = (β + κ 0 + κ), we recover the proposal Eq. (9), n I (t) ↔ I.