An adaptive social distancing SIR model for COVID-19 disease spreading and forecasting

The SIR model

The SIR model is one of the most epidemiological compartmental models used for the investigation of the spread of disease caused by a virus. This model was introduced by Kermack and McKendrick in 1927 to describe the evolution of the relative proportions of three disjoint classes which change with time t (Figure 2.1).

Figure 2.1:

Compartimental diagram for SIR model.

In the SIR model, the population is subdivided into the following three epidemiologically distinct types of individuals and these classes are be represented by

1. The susceptible class S(t) the individuals have not yet been infected by the disease of interest.

2. The infected class I(t) the individuals who are currently infected and capable to transmit the disease to others.

3. The removed class R(t) the individuals formally infected who are deceased, or have recovered and are permanently immune.

Many hypotheses are being assumed in the model, including a large and closed population, the outbreak is relatively short lived, we neglect the natural births-deaths, the absence of latent period and the healing after infection grant a lifetime immunity. It is assumed also that members of both susceptible and infectious populations are homogeneously distributed in space. In other words, every individual will be meeting every other individual per unit time with equal probability. Mathematically, the model can be formulated as a system of three equations (Kermack and McKendrick 1927):

Herein, the term β N S ( t ) I ( t ) represents the disease transmission rate by contact between susceptible and infected individuals. This means that the rate of come across between susceptible and infected individuals is proportional to the product of the population sizes. The parameter β is the infection rate (to be estimated) and γ is the recovery rate which can be computed with respect to the period of viral shedding. Adding all equations in the model, we observe that the total population N = S + I + R is a constant. Indeed, we see that

The system (2.1) are completed by the initial conditions representing the first infected individual appeared in the population, which correspond to S(t ₀) = N, I(t ₀) = 0 and R(t ₀) = 0. Note that, the case where the number of infected and recovered are neglected with respect to the total population the number of susceptible S, we obtain the following approximation: S(t) ∼ S(t ₀) = N. Observe that, the first Eq. (2.1) can be written as

Therefore, we obtain an ordinary differential equation with constant coefficients and the solution can be expressed as

The doubling time, which represents the number of infected people can be doubled over each time interval of length T _d = ln(2)/(β − γ). Based on available data, this exponential growth is observed in the first days of starting COVID-19 outbreak data from various countries, see the following figure (Figure 2.2):

Figure 2.2:

Reported infective data during first 40 days since the appears of the first case.

For more analysis about Classical SIR model, see for instance (Hethcote 2000; Murray 2002; WEISS 2009) and the references therein and the model has been proposed recently to forecasting COVID-19 pandemic in (Bertozzi et al. 2020).

The SIRD model

Now regarding system (2.1), we can extend this SIR model including an additional class D(t) of deceased individuals due to the pandemic and σ the mortality rate of the infected (Figure 2.3).

Figure 2.3:

Compartmental diagram for SIRD model.

The dynamic of SIRD model is governed by the following system of ODEs

We observe that the total population S(t) + I(t) + R(t) + D(t) = N is a constant for each time t. Now, thanks to (Murray 2002), the basic reproduction number of (2.4) can be computed by

Similarly to subsection (2.1), we arrive to I(t) ∼ I(t ₀)  exp ((β − γ − σ)t) for S(t ₀) ∼ N. In this case, it is easy to see that the number of infected population will grow for β > γ + σ (i.e. R 0 > 1 ) and when β < γ + σ (i.e. R 0 < 1 ) it will decay which plays the key role behind disease propagation. On the other hand, thanks to (2.4), the cumulative infections can be computed by I _c (t) = I(t) + R(t) + D(t) and the evolution equation of I _c can be written as

Integrating this equation, we observe that I c ( t ) ∼ β β − γ − σ I ( t 0 ) exp ( ( β − γ − σ ) t ) for S(t) ∼ N.

The non-linear SIRD model

In order to reduce the impact of COVID-19, government measures intended to modify human behavior, such as the social distancing strategy that can help to slow the spread of COVID-19. Various restrictive measures adopted for the coronavirus infection suggest to introduce some new models of the infection spreading. Inspired by (Eksina, Paarpornb, and Weitz 2019), our goal in this work is to introduce a new model taking account of human behavior to the propagation of diseases and which individuals reduce their contacts with the infected individuals. Now, we propose the following modified SIRD model where the social distancing is taken into account

Herein, L is a non-linear continuous function that describe the social distancing and the sensitivity to the disease. This function is given in (3.2) below. Note that, in the case where the function L ( t , r ) ∼ r for any r , t ∈ R + , we get the SIRD model (2.4). We would like to mention that after the pandemic starts to propagate, there is a large reduction in the population S(t) by social distancing measures. For that, we assume that L ( t , r ) ∼ r for t ≤ t _c and L ( t , r ) ∼ r h ( r ) for t ≥ t _c , where t _c is the time of beginning lockdown and h(r) is a non decreasing function. A typical examples of the function L are

where p ε δ ( ⋅ ) > 1 and q ^δ (⋅) ≥ 0 are continuous functions such that

and

Herein, ɛ ∼ 0 is a regularized parameter and δ represents a period for which the lockdown is taking effect. Regarding the social distancing function p ε δ ( t ) (res. q ^δ (t)), it is clear that when the social distancing starts the proportion of susceptible individuals is reduced. Note that, before of the beginning of the outbreak at time t _c (the time at which social distancing is started), we get easily from (3.1)–(3.4) p ε δ ( t ) ∼ 1 (res. q ^δ (t) ∼ 0) and we deduce the classical SIRD model (2.4). From time t _c to t _c + δ, is the interval when social distancing measures are started. After time t _c + δ, the social distancing measures are completely effective. Without loss of generality, we use the expression of L in the rest of this paper.

The reproduction number and phase plane

It is clear that our model has a non-trivial equilibrium of the system (3.1), given by (S ^⋆, 0, R ^⋆, D ^⋆) where S ^⋆ + R ^⋆ + D ^⋆ = N. Moreover, the pandemic will occur if the infected proportion increases, namely

Recall that, the total population S(t) + I(t) + R(t) + D(t) = N for each time t ≥ 0. In the next step we focus only on the equations in S and I. Following the works in (Driessche and Watmough 2002), the right side of the system (3.1) can be decomposed as F + V + + V − where

Now, we let F = β L t 0 , S ( t 0 ) N and V = −(γ + σ) to deduce that

where ρ is a spectral radius. Moreover, the number of susceptible individuals that can get infected by an individual in the population at time t, can be represented by the following effective reproduction number

Next, we study the phase plane of our model. We divide the first equation by the second in (3.1), to formulate the phase portrait to our system in the following equation

For L t , S N = S N p ε δ ( t ) , we integrate Eq. (3.9) with respect to S, to obtain

for any time t > t _c . Notice that

Therefore, the maximum I _max reach after t _c + δ (recall that I _max is obtained when d I d t = 0 ) for S ⋆ = N R 0 1 p ⋆ . Consequently, the maximum number of ineffective population is given by

Moreover, the first infected individual appears in the population corresponding to S(t ₀) + I(t ₀) = N and it follows

Now, dividing the first equation by the last one of (3.1), the phase portrait for the pandemic system, can be formulated as

Now, integrating this equation with respect to S, we get

Letting t → ∞, assuming R(t ₀) = D(t ₀) = 0 and S(t ₀) = N and thanks to the fact that I(∞) = 0, we obtain R(∞) + D(∞) = N − S(∞). Additionally, the proportion of susceptible S(∞) can be computed by solving the following non-linear equation

Numerical approches

Generally, explicit treatment of (2.1) means stricter conditions on time step and the parameters to achieve stability of time discretization schemes. To overcome this, we will use the implicit Euler method, discretizing the time variable as t = kτ where k is a positive number and τ is the time step. The discrete iterative scheme of the proposed model can be written as

where S _k+1, I _k+1, R _k+1 and D _k+1 are numbers of the susceptible, the infected and the recovered and death individuals in the specific day k + 1 respectively. First, assume that L t k + 1 , S k + 1 N = S k + 1 N for any k ≥ 0, the discrete system (3.11) implies

Now, we add the first and the second equations in (3.12) to obtain

Using this and the first equation in (3.12), we deduce

Replacing (3.14) in (3.13) and using the fact that S _k + I _k + R _k + D _k = N, we get

Consequently, we have

which can be computed directly. Now we take Θ k + 1 ≔ N S k + 1 L t k + 1 , S k + 1 N to obtain from (3.11)

Therefore, the solution of (3.17) can be obtained by using the following fixed point method

where Θ k + 1 n ≔ N S k + 1 n L t k + 1 , S k + 1 n N for k , n ∈ N . We shall see further that all results obtained for scheme (3.13)–(3.15) are also true, with convenient adaptations, for the fully implicit scheme. Now, we present the algorithm used to solve our proposed model:

Parameter estimation

In this subsection, we present the technic to fit the proposed model to COVID-19 data. We estimate the transmission rate β, the recovery rate γ, the death rate σ and the distancing parameter p ^⋆. The model parameters ϑ ≔ (β, γ, σ, p ^⋆) can be estimated via least-square fitting of the model solution to the observed data (Liu and Stechlinski 2017). This can be expressed by finding for the set of parameters ϑ ⋆ ≔ ( β ⋆ , γ ⋆ , σ ⋆ , p ⋆ ⋆ ) that minimizes the sum of squared differences between the real data (ζ ₁, ζ ₂, …, ζ _k ) and the corresponding model solution named by (ζ _ϑ (t ₁), ζ _ϑ (t ₂), …, ζ _ϑ (t _k )). Consequently, the objective function is written as

where ζ _ϑ (t) ≔ (I _C,ϑ (t), I _ϑ (t), R _ϑ (t), D _ϑ (t)). Herein, (I _ϑ (t), R _ϑ (t), D _ϑ (t)) is given by the solution of the system (2.1) and the cumulative infection is computed by I _C,ϑ (t) = I _ϑ (t) + R _ϑ (t) + D _ϑ (t). The time of observed data is denoted by t _i , where i is the number of data points. We accord more weight to smaller data points and the corresponding weight can expressed as w 1 , t j = 1 I C , j , w 2 , t j = 1 I j , w 3 , t j = 1 R j , w 4 , t j = 1 D j for j = 1, …, k − 1 and some alternative choices of the weight can be found in (Chowell 2017).

Figure 4.1:

Numerical simulation of our model with the parameters β, γ and σ which are computed by fitting SIRD model during the first 30 days. The first and the third rows represent the infected cases, the second and the last rows are the cumulative cases. We investigate the effect of changing the parameter p ^⋆ in the simulation of COVID-19 spreading for various countries (recall that p ^⋆ = 1 corresponds to the classical SIRD model).

Figure 4.2:

Fit of data from Germany, Spain, Italy, France, Algeria and Morocco with classical SIR and the proposed model. The black dots are the susceptible, infective, recovered and death reported cases.

Figure 4.3:

Simulation of the reproduction number R t p from (4.3) for the previous countries.

Numerical tests

Various countries has implemented social distancing measures to reduce the transmission of the outbreak between individuals. These types of measures have an impact on the susceptible numbers on the proposed non-linear proposed model (3.1): reducing the contact L . , S N I . To estimate the impacts on the non-linearity by these measures, we develop several numerical simulations of the effects of these measures by adjusting the estimates of the parameter p ^⋆.

Now we illustrate various numerical tests based on some real data. We consider the problem (3.1) with

In the following, we assume that ɛ = 10⁻⁶, δ = 8 days and t _c is the beginning time of lockdown announced by the authorities in each country. In the first test, we use our algorithm and we fit the SIRD model (2.4) with reported data in the period from the beginning of the outbreak until time t _c at which social distancing is started, and we compute the the parameters β, γ and σ. Then by using these parameters we simulate the infected and cumulative cases in our model (3.1) and we illustrate the effect social distancing by varying the parameter p ^⋆.

The previous tests illustrate the effect of social distancing represented by p ^⋆ on the number of infected population in various countries. In the absence of social distancing i.e. p ^⋆ = 1 the number of infected cases grows exponentially to reach a maximum infected number and also decreases exponentially. The case where p ^⋆ > 1 corresponds to various social distancing measures for e.g. closing schools and universities, isolation and home quarantine. We observe that when p ^⋆ increases, the impact of the pandemic is reduced by flattening the curve, reducing the peak incidence and overall the infected numbers (Figure 4.1).

In the next test, we fit the observed data and we compute the appropriate parameters of the proposed model. Additionally, we compare the fitting results with the classical SIRD model (see Table 2). This permits to evaluate the social distancing measures in each country. In Figure 4.2, we can see which of SIRD models is the best by fitting them to data. We observe that our model-fitting SIRD system with p ^⋆ > 1 is better compared to the classical case p ^⋆ = 1. Moreover, the adaptive nature of our model gives a considerable results. We plan to update these results regularly on our platform (http://vps730710.ovh.net:3000/) for various country.

Table 2 represents the obtained values of the parameters by fitting our model with reported data. The parameter p ^⋆ illustrate the degree of social distancing in each country. We notice that the value of the parameter p ^⋆ is more important in Morocco and in Algeria. This reflects the strict lockdown measures implemented in these countries.

Now, Thanks to (3.5), the effective reproduction number is given by

In the next figure, we present the evolution of R t p for each country.

In Figure 4.3, we observe the efficiency of the measure implemented by Moroccan and Algerian governments, that leads to control the outbreak and permit to decrease quickly the R t p reproduction number (recall these two countries started early to implement strict lockdown). Observe that the reproduction number R t p of Germany is close to the one of Morocco and Algeria. Note that Germany has tested for COVID-19 on a larger scale comparing to France, Spain and Italy. Consequently, Germany had a lower death rate comparing to these countries.