The Dynamics of Distributions in Continuoustime Stochastic Models
Introduction
The pandemic of COVID-19 has scourged the world since the beginning of 2020. Footnote 1 The responsible virus is the SARS-CoV-2, identified in China at the end of 2019 (Zhou et al. 2020). Governments are constantly looking for ways to predict and contain the spread of such an illness in order to monitor the public health and to prevent economic and social issues.
Epidemic models constitute a branch of interest in applied mathematics since several years. They are adopted not only to study epidemics properly but also to predict social phenomenon or the behavior of biological systems. The simplest epidemic model is called SIR model. It looks at a population split out into three compartments: susceptible, infected and removed. The SIR model was introduced the first time in 1927 (Kermack and McKendrick 1927) and many variations have been proposed to study diseases with complex behaviors and other phenomena (Murray 2002). Epidemic models can also be formulated by means of the theory of stochastic processes. The first application of stochastic processes to epidemics was presented in 1955 (Whittle 1955) and more recently several applications have been proposed (Capasso and Bakstein 2015). In some cases there is an equivalence between the two approaches (Allen 2010). Concerning the new pandemic, many mathematical models have been proposed. In Ansumali et al. (2020), Calafiore et al. (2020), Giordano et al. (2020) some deterministic epidemic models for COVID-19 based on ordinary differential equations have been proposed. In Zhang et al. (2020) a stochastic dynamic model has been introduced while in Faranda and Alberti (2020), Rihan et al. (2020) the authors propose mathematical models based on stochastic differential equations. A novel model based on an operatorial approach as in quantum mechanics can be found in Bagarello et al. (2020).
In this paper we would like to introduce some epidemic models based on stochastic processes, taking into account peculiarities of the COVID-19 disease. We proposed two models. In the first one we consider that COVID-19 has an incubation period in which people are apparently healthy and after that they become infected and are also able to infect other people. Since COVID-19 has a quite high fatality rate, the removed people have been split in two sub-classes: healed and dead. We suppose that an individual recovers or dies after a fixed time from infection. In the second model we would like to include asymptomatic people, i.e. infectious individuals without severe or identifiable symptoms. They seem to play an important role in the diffusion of the virus because usually they don't know to be infectious. Since it is not clear whether during the incubation an individual is infectious or not, in this model we make the assumption that it is possible. Moreover, it is not ascertained so far whether and for how long people preserve the immunity to the virus. Therefore, we consider the possibility for a healed individual to lose immunity and to become susceptible again.
The plan of the paper is as follows. In Sects. 2 and 3 the two stochastic models are introduced; in Sect. 4 we present the Monte Carlo algorithm adopted for simulations; in Sect. 5 we propose a deterministic model to assess the validity of the one introduced in 2; in Sect. 6 we show and comment the numerical simulations, in particular we test the proposed models by considering the second wave of COVID-19 in Italy.
A SIRD model for COVID-19 disease
Let us consider a fixed (no births and no deaths) population of N individuals split out into four compartments: susceptible S, infected I, removed R and dead people D. In principle, it is also possible to include the so-called vital dynamics by introducing birth and death rates but in the typical time scale of the pandemic spread the effects can be considered negligible. We suppose that the number of individuals in each class evolves in time \(t\in [0,+\infty [\) because of two mechanisms: susceptible individuals become infected and infected individuals recover or die. We call the introduced model SIRD.
To describe the infection mechanism, we suppose that the rate of new infectious cases is proportional to the number of susceptible individuals S(t) times the fraction of infected people I(t)/N. The proportionality factor is denoted by \(\beta >0\) which represents the average number of contacts of a person per unit time (the day in our case) that result in an infection of a susceptible individual.
In relation to the recovery mechanism, we suppose that an infected individual has a probability \(\alpha \in [0,1]\) to die and \(1-\alpha \) to heal. The situation is schematized in Fig. 1.
Flow chart of the SIRD model. Above the arrows between the compartment I and the compartments R, D we indicate the times for the transition while below the arrows the probabilities of changing compartment. \(\tau _1\) is the delay time for getting infected and \(\varDelta t \sim \beta ^{-1}\) is the frequency of the contacts between the individuals of S and I
Moreover, let \(\tau _1\) be the incubation time which, in this model, is the time after that an individual who contracted the virus is infected and becomes infectious. Finally, let \(\tau _2\) and \(\tau _3\) be the heal and dead time respectively, counted after the incubation period. First, we will assume that \(\tau _1, \tau _2, \tau _3\) are constant; in a second step this assumption will be relaxed to simulated effects as those related to the lockdown.
Let us consider \(\left( S,I,R,D \right) \) as four random variables which take values in \({\mathbb {N}}_0^4\) and depend on time t. The epidemic is triggered by a small number of infected individuals \(I_0\), that is at time \(t=0\) we have
$$\begin{aligned} (S,I,R,D)=(N-I_0,I_0,0,0), \end{aligned}$$
with \(I_0 \ll S_0\).
We split the set (S,I,R,D) into three subsets in each of which there will be defined a stochastic process. The first set is \(\left( S,I \right) \), corresponding to the infection mechanism. We think of the two random variables as unsynchronized in time. In particular, we set
$$\begin{aligned} {\tilde{I}}(t)&= I(t+\tau _1) \end{aligned}$$
and we consider the new set \(\left( S,{\tilde{I}} \right) \) and suppose that at time t the random variables take values \(\left( s,i \right) \), that is
$$\begin{aligned} \left( S(t),{\tilde{I}}(t) \right) = \left( s,i \right) . \end{aligned}$$
(1)
After a small period of time \(\varDelta t>0\), \(\varDelta t \ll \displaystyle {\min _{i=1, \ldots ,3}} \tau _i\) the state of system changes in
$$\begin{aligned} \left( S(t+\varDelta t),{\tilde{I}}(t+\varDelta t) \right) = \left( s+m,i+n \right) , \end{aligned}$$
(2)
being \(m,n\in \left\{ -1, 0, +1 \right\} \).
The transition probability is defined as follows
$$\begin{aligned} p_{(s,i)\rightarrow (s+m,i+n)}(\varDelta t)&= P\Big ( \left( S(t+\varDelta t),{\tilde{I}}(t+\varDelta t) \right) = \left( s+m,i+n \right) \\&\quad \Big \vert \, \left( S(t),{\tilde{I}}(t) \right) = \left( s,i \right) \Big ). \end{aligned}$$
In this way we define a continuous time Markov chain and the transition probability can be written as
$$\begin{aligned} p_{(s,i)\rightarrow (s+m,i+n)}(\varDelta t) = \left\{ \begin{aligned}&\beta s \frac{i}{N} \varDelta t + o(\varDelta t),&\qquad (m,n)=(-1,+1)\\&1-\left( \beta s \frac{i}{N} \right) \varDelta t + o(\varDelta t),&\qquad (m,n)=(0,0)\\&o(\varDelta t),&\qquad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
The second set of random variables is \(\left( I,R \right) \), coupled to the first one. We set
$$\begin{aligned} \tilde{{\tilde{I}}}(t)&= I(t+\tau _1+\tau _2),\\ \tilde{{\tilde{R}}}(t)&= R(t+\tau _1+\tau _2) \end{aligned}$$
and consider the couple \(\left( \tilde{{\tilde{I}}},\tilde{{\tilde{R}}} \right) \). Let us suppose that at time t the random variables \(\left( \tilde{{\tilde{I}}},\tilde{{\tilde{R}}} \right) \) take values (i,r), that is
$$\begin{aligned} \left( \tilde{{\tilde{I}}}(t),\tilde{{\tilde{R}}}(t) \right) = (i,r). \end{aligned}$$
(3)
After a small period of time \(\varDelta t>0\) the state of the system changes in
$$\begin{aligned} \left( \tilde{{\tilde{I}}}(t+\varDelta t),\tilde{{\tilde{R}}}(t+\varDelta t) \right) = (i+n,r+u), \end{aligned}$$
(4)
being \(n,u\in \left\{ -1, 0, +1 \right\} \). Moreover, since a healing at time \(t+\tau _1+\tau _2\) is related to an infection in the past interval \([t+\tau _1,t+\tau _1+\varDelta t]\), we need to know the values assumed by the random variable \({\tilde{I}}\) at t and \(t+\varDelta t\). In this case the transition probability is defined as follows
$$\begin{aligned} p_{(i,r)\rightarrow (i+n,r+u)}(\varDelta t)&= P\Big ( \left( \tilde{{\tilde{I}}}(t+\varDelta t),\tilde{{\tilde{R}}}(t+\varDelta t) \right) = \left( i+n,r+u \right) \\&\quad \Big \vert \, \left( \tilde{{\tilde{I}}}(t),\tilde{{\tilde{R}}}(t) \right) = \left( i,r \right) , {\tilde{I}}(t)=j, {\tilde{I}}(t+\varDelta t)=j+1 \Big ). \end{aligned}$$
Therefore, a non-Markovian continuous time stochastic process is defined with the transition probability
$$\begin{aligned} p_{(i,r)\rightarrow (i+n,r+u)}(\varDelta t) = \left\{ \begin{aligned}&(1-\alpha ) \varDelta t + o(\varDelta t),&\qquad (n,u)=(-1,+1)\\&1-(1-\alpha ) \varDelta t + o(\varDelta t),&\qquad (n,u)=(0,0)\\&o(\varDelta t),&\qquad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
Finally, the third set of random variables is \(\left( I,D \right) \), which is also coupled to the first one. Now we set
$$\begin{aligned} {\hat{I}}(t)&= I(t+\tau _1+\tau _3),\\ {\hat{D}}(t)&= D(t+\tau _1+\tau _3) \end{aligned}$$
and consider the couple \(\left( {\hat{I}}, {\hat{D}} \right) \). Let us suppose that at time t the random variables take values (i,d), that is
$$\begin{aligned} \left( {\hat{I}}(t), {\hat{D}}(t) \right) = (i,d). \end{aligned}$$
(5)
After a small period of time \(\varDelta t>0\) the state of the system changes in
$$\begin{aligned} \left( {\hat{I}}(t+\varDelta t), {\hat{D}}(t+\varDelta t) \right) = (i+n,d+v), \end{aligned}$$
(6)
being \(n,v\in \left\{ -1, 0, +1 \right\} \). Moreover, since a death at time \(t+\tau _1+\tau _3\) is related to an infection in the past interval \([t+\tau _1,t+\tau _1+\varDelta t]\), we need to know the values assumed by the random variable \({\tilde{I}}\) at t and \(t+\varDelta t\). In this case the transition probability is defined as follows
$$\begin{aligned} p_{(i,d)\rightarrow (i+n,d+v)}(\varDelta t)&= P\Big ( \left( {\hat{I}}(t+\varDelta t), {\hat{D}}(t+\varDelta t) \right) = \left( i+n,d+v \right) \\&\quad \Big \vert \, \left( {\hat{I}}(t), {\hat{D}}(t) \right) = \left( i,d \right) , {\tilde{I}}(t)=j, {\tilde{I}}(t+\varDelta t)=j+1 \Big ). \end{aligned}$$
Therefore, a non-Markovian continuous time stochastic process is defined with the transition probability
$$\begin{aligned} p_{(i,d)\rightarrow (i+n,d+v)}(\varDelta t) = \left\{ \begin{aligned}&\alpha \varDelta t + o(\varDelta t),&\qquad (n,v)=(-1,+1)\\&1-\alpha \varDelta t + o(\varDelta t),&\qquad (n,v)=(0,0)\\&o(\varDelta t),&\qquad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
A SAI(L)RD model for COVID-19
At variance with the SIRD, the second model we are going to introduce also contemplates another two compartments: the class A of asymptomatic individuals and the class L of isolated infected people. We call it SAI(L)RD model. Its detailed features are summarized below
- A1.
We suppose that if an effective contact occurs between an infected (symptomatic or not) and a susceptible individual then the latter becomes infected, with or without symptoms, and infectious at the same time.
- A2.
We suppose there exists a probability \(\eta \in [0,1]\) to be asymptomatic and, consequently, \(1-\eta \) is the probability to show symptoms once an individual has been infected. Moreover, we assume that, after a certain time \(\tau _1\), a symptomatic individual is recognized and isolated into a subclass, called L, of lonely individuals of I.
- A3.
Lonely individuals are not infectious anymore. Asymptomatic people heal after a time \(\tau _2\). Lonely symptomatic individuals can die with probability \(\alpha \) in a time \(\tau _4\) or heal with probability \(1-\alpha \) in a time \(\tau _3\).
- A4.
Finally, removed individuals become immune to COVID-19 for a short period of time or forever. Let \(\lambda \in [0,1]\) be the probability that COVID-19 confers a short immunity of time length \(\tau _5\), after which recovered individuals come back to the class of susceptible people and, in principle, can suffer a reinfection.
The situation is represented in Fig. 2. Note that both people belonging to the class A and I can infect. Moreover, a person who has become infected on account of a contact with an asymptomatic individual can go to the class A or I.
Flow chart of the SAI(L)RD model. Above the arrows we indicate the times for the transition while below the arrows the probabilities of changing compartment
We consider \(\left( S,A,I,L,R,D \right) \) as a six dimensional random variable which can assume values in \({\mathbb {N}}_0^6\) and depends on time \(t\ge 0\). Now we split the set of random variables into several coupled sub-systems.
-
The first set is (S,A,I) corresponding to the infection mechanism. We suppose that at time t the random variables take values (s,a,i), that is
$$\begin{aligned} \Big ( S(t),A(t),I(t) \Big ) = (s,a,i). \end{aligned}$$
(7)
After a small period of time \(\varDelta t>0\), \(\varDelta t \ll \displaystyle {\min _{i=1, \ldots ,5}} \tau _i\), the state of the system changes in
$$\begin{aligned} \Big ( S(t+\varDelta t),A(t+\varDelta t),I(t+\varDelta t) \Big ) = (s+m,a+k,i+n), \end{aligned}$$
(8)
with \(m,k,n\in \left\{ -1,0,1 \right\} \).
The transition probability is given by
$$\begin{aligned}&p_{(s,a,i)\rightarrow (s+m,a+k,i+n)}(\varDelta t) \\&\qquad = P\Big ( \left( S(t+\varDelta t),A(t+\varDelta t),I(t+\varDelta t) \right) = \left( s+m,a+k,i+n \right) \\&\qquad \quad \, \Big \vert \, \left( S(t),A(t),I(t) \right) = \left( s,a,i \right) \Big ). \end{aligned}$$
In this way we define a time continuous Markov chain with transition probabilities which can be written as
$$\begin{aligned}&p_{(s,a,i)\rightarrow (s+m,a+k,i+n)}(\varDelta t) \\&\qquad \qquad = \left\{ \begin{aligned}&\eta \beta s \frac{i}{N} \varDelta t + o(\varDelta t),&\quad \text{ if } \quad (m,k,n)=(-1,+1,0)\\&(1-\eta )\beta s \frac{i}{N} \varDelta t + o(\varDelta t),&\quad \text{ if } \quad (m,k,n)=(-1,0,+1)\\&1-\left( \beta s \frac{i}{N} \right) \varDelta t + o(\varDelta t),&\quad \text{ if } \quad (m,k,n)=(0,0,0)\\&o(\varDelta t),&\quad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
-
In our model, we suppose that all the infected individuals will be isolated after a time \(\tau _1\). Let us consider (I,L). We set
$$\begin{aligned}&{\tilde{I}}(t) = I(t+\tau _1),\\&{\tilde{L}}(t) = L(t+\tau _1). \end{aligned}$$
In this case a lone infected individual at time \(t + \varDelta t + \tau _1\) is related to a symptomatic infection in the past interval \([t,t+\varDelta t]\). Regarding the transition probability
$$\begin{aligned}&p_{(i,l)\rightarrow (i+n,l+j)}(\varDelta t) = P\Big ( \left( {\tilde{I}}(t+\varDelta t),{\tilde{L}}(t+\varDelta t) \right) = \left( i+n,l+j \right) \\&\qquad \qquad \, \Big \vert \, \left( {\tilde{I}}(t),{\tilde{L}}(t) \right) = \left( i,l \right) , I(t)=f,I(t+\varDelta t)=f+1 \Big ) \end{aligned}$$
for some nonnegative integer f, we have
$$\begin{aligned} p_{(i,l)\rightarrow (i+n,l+j)}(\varDelta t) = \left\{ \begin{aligned}&1 ,&\quad \text{ if } \quad (n,j)=(-1,+1)\\&0,&\quad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
-
Now we consider the pair of random variables (L,R). To define the healing process we set
$$\begin{aligned} \tilde{{\tilde{L}}}(t) = L(t+\tau _1+\tau _3),\\ \tilde{{\tilde{R}}}(t) = R(t+\tau _1+\tau _3) \end{aligned}$$
and suppose that
$$\begin{aligned} \Big ( \tilde{{\tilde{L}}}(t),\tilde{{\tilde{R}}}(t) \Big ) = (l,r). \end{aligned}$$
After a small period of time \(\varDelta t >0\) the state of the system changes in
$$\begin{aligned} \Big ( \tilde{{\tilde{L}}}(t+\varDelta t),\tilde{{\tilde{R}}}(t+\varDelta t) \Big ) = (l+j,r+u) \end{aligned}$$
with \(j,u\in \left\{ -1,0,1 \right\} \). Moreover, since in this case a healing at time \(t + \varDelta t + \tau _1+\tau _3\) is related to an infection with symptoms in the past interval \([t,t+\varDelta t]\), we need to know the values assumed by the random variable I at t and \(t+\varDelta t\). This is also equivalent to knowing the values assumed by the random variable L at \(t+\tau _1\) and \(t + \varDelta t + \tau _1\). In this case the transition probability is defined as
$$\begin{aligned}&p_{(l,r)\rightarrow (l+j,r+u)}(\varDelta t) = P\Big ( \left( \tilde{{\tilde{L}}}(t+\varDelta t), \tilde{{\tilde{R}}}(t+\varDelta t) \right) = \left( l+j,r+u \right) \\&\qquad \qquad \, \Big \vert \, \left( \tilde{{\tilde{L}}}(t),\tilde{{\tilde{R}}}(t) \right) = \left( l,r \right) , L(t+\tau _1)=f,L(t + \varDelta t + \tau _1) = f+1 \Big ), \end{aligned}$$
for some nonnegative integer f. In this way a time continuous non-Markovian stochastic process is defined whose transition probability can be written as
$$\begin{aligned} p_{(l,r)\rightarrow (l+j,r+u)}(\varDelta t) = \left\{ \begin{aligned}&(1-\alpha )\varDelta t + o(\varDelta t),&\quad \text{ if } \quad (j,u)=(-1,+1)\\&1-(1-\alpha )\varDelta t + o(\varDelta t),&\quad \text{ if } \quad (j,u)=(0,0)\\&o(\varDelta t),&\quad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
-
With the same arguments, a dead process is described by \(({\hat{L}},{\hat{D}})\) where
$$\begin{aligned} {\hat{L}}(t) = L(t+\tau _1+\tau _4),\\ {\hat{D}}(t) = D(t+\tau _1+\tau _4). \end{aligned}$$
In this case a death at time \(t+ \varDelta t + \tau _1+\tau _4\) is related to an infection with symptoms in the past interval \([t,t+\varDelta t]\) and thus the gain of one unit to the variable L in \([t+\tau _1,t+\varDelta t + \tau _1]\). The transition probability
$$\begin{aligned}&p_{(l,d)\rightarrow (l+j,d+v)}(\varDelta t) = P\Big ( \left( {\hat{L}}(t+\varDelta t),{\hat{D}}(t+\varDelta t) \right) = \left( l+j,d+v \right) \\&\qquad \qquad \, \Big \vert \, \left( {\hat{L}}(t),{\hat{D}}(t) \right) = \left( l,d \right) , L(t+\tau _1)=f,L(t+\varDelta t +\tau _1)=f+1 \Big ), \end{aligned}$$
for some nonnegative integer f, is given by
$$\begin{aligned} p_{(l,d)\rightarrow (l+j,d+v)}(\varDelta t) = \left\{ \begin{aligned}&\alpha \varDelta t + o(\varDelta t),&\quad \text{ if } \quad (j,v)=(-1,+1)\\&1-\alpha \varDelta t + o(\varDelta t),&\quad \text{ if } \quad (j,v)=(0,0)\\&o(\varDelta t),&\quad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
-
A further process we introduce is the healing of an asymptomatic individual. Let us consider (A,R). We set
$$\begin{aligned}&{\hat{A}}(t) = A(t+\tau _2),\\&{\hat{R}}(t) = R(t+\tau _2). \end{aligned}$$
In this case a healing at time \(t+\tau _2 + \varDelta t\) is related to an asymptomatic infection in the past interval \([t,t+\varDelta t]\). The transition probability
$$\begin{aligned}&p_{(a,r)\rightarrow (a+k,r+u)}(\varDelta t) = P\Big ( \left( {\hat{A}}(t+\varDelta t),{\hat{R}}(t+\varDelta t) \right) = \left( a+k,r+u \right) \\&\qquad \qquad \, \Big \vert \, \left( {\hat{A}}(t),{\hat{R}}(t) \right) = \left( a,r \right) , A(t)=f,A(t+\varDelta t)=f+1 \Big ), \end{aligned}$$
for some nonnegative integer f, reads
$$\begin{aligned} p_{(a,r)\rightarrow (a+k,r+u)}(\varDelta t) = \left\{ \begin{aligned}&1,&\quad \text{ if } \quad (k,u)=(-1,+1)\\&0,&\quad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
-
The last process we are going to introduce is the one involving a removed individual who comes back to the class of susceptible people after a certain time. To define the process, we set
$$\begin{aligned}&\hat{{\hat{S}}}(t)=S(t+\tau _5),\\&\hat{{\hat{R}}}(t)=R(t+\tau _5) \end{aligned}$$
Let us suppose that at time t the two-dimensional random variable \((\hat{{\hat{S}}},\hat{{\hat{R}}})\) takes the value (s,r). After a small period of time \(\varDelta t>0\) the state of the system changes in
$$\begin{aligned} (\hat{{\hat{S}}}(t+\varDelta t),\hat{{\hat{R}}}(t+\varDelta t))=(s+m,r+u), \end{aligned}$$
with \(m,u\in \left\{ -1,0,1 \right\} \). Moreover, since in this case a healed individual can come back to the class of susceptible people after a certain time \(\tau _5\), it is needed to know the values assumed by the random variable R at time t and \(t+\varDelta t\). In this case the transition probability
$$\begin{aligned}&p_{(s,r)\rightarrow (s+m,r+u)}(\varDelta t) = P\Big ( \left( \hat{{\hat{S}}}(t+\varDelta t),\hat{{\hat{R}}}(t+\varDelta t) \right) = \left( s+m,r+u \right) \\&\qquad \qquad \, \Big \vert \, \left( \hat{{\hat{S}}}(t),\hat{{\hat{R}}}(t) \right) = \left( s,r \right) , R(t)=f,R(t+\varDelta t)=f+1 \Big ), \end{aligned}$$
for some nonnegative integer f, is given by
$$\begin{aligned} p_{(s,r)\rightarrow (s+m,r+u)}(\varDelta t) = \left\{ \begin{aligned}&\lambda + o(\varDelta t),&\quad \text{ if } \quad (m,u)=(+1,-1)\\&1-\lambda + o(\varDelta t),&\quad \text{ if } \quad (m,u)=(0,0)\\&o(\varDelta t),&\quad \text{ otherwise }. \end{aligned} \right. \end{aligned}$$
The major advantage to adopt a stochastic model is the possibility to easily add further more sophisticated features. Indeed, the delays are assumed constants but it is possible to consider in turn the times \(\tau _i\) as random variables obeying suitable probability distributions. However, in average we get the same results.
The Monte Carlo method for the simulations
An efficient simulation of both the SIRD and SAI(L)RD models can be performed by a Monte Carlo approach. The details are outlined in the next subsections.
SIRD model
Firstly we describe the method adopted for SIRD model. The state of the system is represented by a time-dependent random vector
$$\begin{aligned} {\mathbf {X}}(t) = \left( X_1 (t), X_2 (t), \ldots , X_N(t) \right) \in D^N, \end{aligned}$$
(9)
for \(t\ge 0\), where \(D=\left\{ -1, 0, 1, 2 \right\} \) and N is the population size. We indicate by \(X_i(t)\in D\) for \(i=1,\ldots ,N\) the trajectory of an individual in time, i.e. the time-evolution of the states assumed by the i-th person. D is a set of labels where 0 represents a susceptible individual, 1 an infected one, 2 a healed person and \(-1\) a dead individual.
To fix the initial infective people, at time \(t=0\) a number \(I_0\) of individuals are labeled by 1 randomly, all the others are susceptible thus labeled by 0. We note that the pure process of the encounters is Markovian, no matter it leads to an infection or not. Therefore, for each infected individual a contact time t is determined according to the exponential distribution of scale parameter \(\beta \), that is
$$\begin{aligned} t = -\frac{1}{\beta }\log \xi , \end{aligned}$$
(10)
\(\xi \) being a random number uniformly distributed in [0, 1].
Let us suppose that the smallest contact time is that of the j-th individual, \( t_{j,1}\) (the second index indicates the first temporal step of the individual j). At this point another individual i is chosen randomly. If it belongs to the susceptible class then the contact is effective and after the incubation, i.e. a period of time \(\tau _1\), the individual i changes its state in infectious,
$$\begin{aligned} X_i(t_{j,1}+\tau _1)=1. \end{aligned}$$
(11)
At time \(t_{j,1} +\tau _1\), the destiny \(d\in \left\{ -1,2\right\} \) of the new infectious is established too accordingly to a Bernoulli distribution with probability \(\alpha \), that is \(d\sim {\mathcal {B}}(1,\alpha )\). If the destiny is to heal then after a time \(\tau _2\) the individual state changes from infectious to recovered; if the destiny is to die then after a time \(\tau _3\) the state of the i-th individual changes from infectious to dead:
$$\begin{aligned} \begin{aligned}&X_i(t_{j,1} +\tau _1+\tau _2)=2&\qquad \text{ if }\quad d=2,\\&X_i(t_{j,1} +\tau _1+\tau _2)=-1&\qquad \text{ if }\quad d=-1. \end{aligned} \end{aligned}$$
(12)
After the choice of the individual i, the individual j still continues to infect unless in the meantime he has recovered or passed away. Another random infection time \(t_{j,2}\) is generated according to (10) and we set
$$\begin{aligned} t_j= t_{j,1} + t_{j,2}. \end{aligned}$$
Once again we determine the infected individual having associated the minimum time and iterate the procedure. The algorithm ends when there are not susceptible individuals any more.
In order to record the time evolution of the system, a time grid is fixed and at each time of such a grid we count the number of individuals in the several classes. Moreover, to reduce the statistical noise an averaging procedure is applied as follows. We perform the entire simulation k times. Let \(T_r\), \(r = 1, 2, \ldots , k\), be the time at which the algorithm ends at the rth simulation. We set
$$\begin{aligned} m_k = \frac{T_1+T_2+\cdots +T_k}{k}, \end{aligned}$$
(13)
the average of the final process times. After introducing the error as
$$\begin{aligned} \varepsilon _k = \vert m_{k+1}-m_k \vert , \end{aligned}$$
(14)
as stopping criterion we adopt
$$\begin{aligned} \varepsilon _k<\text{ tol }, \qquad k>N_{min}. \end{aligned}$$
(15)
Here \(\text{ tol }\) is a numerical tolerance and \(N_{min}\) is a minimum number of iterations which are required to prevent early stops of the procedure.
SAI(L)RD model
In a similar way the SAI(L)RD model can be simulated. Now we have
$$\begin{aligned} {\mathbf {X}}(t)\in E^N, \end{aligned}$$
where \(E=\left\{ -1,0,1,2,3,4 \right\} \). The label \(-1\) represents a dead person, 0 a susceptible individual, 1 an infected one, 2 a healed person, 3 an asymptomatic individual, 4 a lone infected one.
Even in this case, at time \(t=0\) we randomly select \(I_0\) individuals we label as infected, i.e. by 1. A contact time t is determined according to (10) for each infected individual. Let us suppose that the minimum contact time is that of the j-th individual, \(t_{j,1}\). At this point another individual i is chosen randomly. Now at variance with the SIRD model, if the latter is susceptible her/his state changes as follows: we determine the symptomaticity \(s\in \left\{ 1,3\right\} \) by a Bernoulli distribution having probability \(\eta \in [0,1]\). If \(s=3\) the individual is asymptomatic and, after the time \(\tau _2\), she/he will heal, and therefore
$$\begin{aligned} X_i(t_{j,1}+\tau _2) = 2. \end{aligned}$$
(16)
If \(s=1\) the person is symptomatic infected and, after the time \(\tau _1\), she/he will be isolated, that is
$$\begin{aligned} X_i(t_{j,1}+\tau _1) = 4. \end{aligned}$$
(17)
Moreover, the destiny \(d\in \left\{ -1,2\right\} \) of such an individual is established according to a Bernoulli distribution with probability \(\alpha \). If the destiny is to heal then after a time \(\tau _3\) the individual state changes from infectious to recovered; if the destiny is to die then after a time \(\tau _4\) the individual state changes from infectious to dead. That evolution can be described as follows
$$\begin{aligned} \begin{alignedat}{2}&X_i(t_{j,1}+\tau _1+\tau _3)=2&\qquad \text{ if }\quad d=2,\\&X_i(t_{j,1}+\tau _1+\tau _4)=-1&\qquad \text{ if }\quad d=-1. \end{alignedat} \end{aligned}$$
(18)
Finally there is also the possibility that a healed individual loses the immunity, coming back to the susceptible class. We take that into account by generating a random number according to a Bernoulli distribution of parameter \(\lambda \in (0,1)\). Furthermore, if the immunity is lost, one has two possibilities:
if the individual suffered from an asymptomatic infection we set
$$\begin{aligned} X_i(t_{j,1} + \tau _2+\tau _5) = 0; \end{aligned}$$
(19)
if she/he suffered from a symptomatic infection we set
$$\begin{aligned} X_i(t_{j,1} +\tau _1+\tau _3+\tau _5) = 0. \end{aligned}$$
(20)
After the choice of the individual i, the individual j still continues to infect unless in the meantime she/he has recovered without losing the immunity or passed away. Another random infection time \(t_{j,2}\) is generated according to (10) and we set
$$\begin{aligned} t_j= t_{j,1} + t_{j,2}. \end{aligned}$$
Again we determine the infected individual having associated the minimum time and iterate the procedure. The algorithm ends when all individuals cannot change their state any longer.
To reduce the statistical noise we have adopted the same technique as the SIRD model presented above.
A deterministic SIRD model with delays
In order to check the validity of the SIRD stochastic model proposed above, a deterministic delayed SIRD model is devised as well. Since the disease has an incubation time \(\tau _1\), the number of susceptible people decreases by a quantity depending on the amount of infected at a previous time \(t-\tau _1\); the amount of recovered and dead people after a time \(\tau _2\) and \(\tau _3\) respectively is proportional to the amount of people who have been infected at the previous time \(t-\tau _2\) and \(t-\tau _3\) respectively. From those considerations, we propose the following model
where \(H(\cdot )\) represents the Heaviside step function. In that way we include the effects of an incubation time and healing or death times.
For further analysis it is convenient to work with proportions. After the substitution
$$\begin{aligned} S \mapsto \frac{S}{N}, I \mapsto \frac{I}{N}, R \mapsto \frac{R}{N}, D \mapsto \frac{D}{N} \end{aligned}$$
and some simple algebraic manipulations, the system can be written more explicitly as
The system must be augmented assigning the functions
$$\begin{aligned} \begin{aligned} S(t)&=\varPhi _1(t)&\qquad t\in [- \max (\tau _2,\tau _3),0],\\ I(t)&=\varPhi _2(t)&\qquad t\in [- \max (\tau _1+ \tau _2, \tau _1 + \tau _3),0]. \end{aligned} \end{aligned}$$
As customary we assume that \(\varPhi _1\) and \(\varPhi _2\) are continuous in the considered intervals. Regarding the other variables it is realistic to take \(R(0) = D(0) = 0\). Specifically we assume \(\varPhi _1 (t) = S_0 \in ]0,1[\) and \(\varPhi _2 (t) = I_0 = 1 - S_0 \in ]0,1[\) with \(I_0 \ll 1\). Therefore at \(t=0\) we have
$$\begin{aligned} S(0) + I(0) + R(0) + D(0) = 1. \end{aligned}$$
The presence of the delays makes the qualitative analysis of the system rather cumbersome, so a complete phase portrait is a daunting task. However, some insights can be deduced anyway.
Along the solution of the system (22)
$$\begin{aligned} S(t) + I(t) + R(t) + D(t) = 1 \end{aligned}$$
holds. In fact, summing up the equations (22a)-(22d) one has
$$\begin{aligned} \dot{Y}(t) = 0, \end{aligned}$$
where \(Y=S+I+R+D\). Since \(Y(0) =1\), it follows that \(Y(t) = 1\) \(\forall t > 0\).
The equilibrium points of the system (22) are the following.
-
Endemic solutions:
$$\begin{aligned} S=0, I = I^*, R = R^*, D = D^*. \end{aligned}$$
(23)
with \(I^*,R^*,D^* \in [0,1]\) satisfying \(I^*+ R^*+ D^* = 1\).
-
Disease-free solutions:
$$\begin{aligned} S=S^*, I = 0, R = R^*, D = D^*. \end{aligned}$$
(24)
with \(S^*,R^*,D^* \in [0,1]\) satisfying \(S^*+ R^*+ D^* = 1\).
If we linearize around the generic endemic critical point, one gets the following characteristic equation for the eigenvalues
$$\begin{aligned} \lambda \left( \lambda + \beta I^* \right) =0 \end{aligned}$$
which shows that the endemic stationary states are linearly stable for any \(I^* > 0\).
The linearization around the disease-free stationary points leads to a much more complex characteristic equation
$$\begin{aligned} \lambda \left[ \lambda - \beta S^* \left( e^{- \tau _1 \lambda } + (1 - \alpha ) e^{- (\tau _1 + \tau _2) \lambda } + \alpha e^{- (\tau _1 + \tau _3) \lambda } \right) \right] =0, \end{aligned}$$
which in general admits infinite solutions in the complex plane due to the functional nature of the evolution equations. The only viable way to get the eigenvalues is to resort to a numerical procedure (Feng et al. 2019). Therefore, since the primary goal is to have a comparison with the stochastic model, we look directly at the numerical solutions of the system (22). To this aim, we adopt a first order finite differences scheme.
Let us fix a temporal grid \(0=t_0<t_1<\ldots <t_M=T_{max}\) of constant time step \(\varDelta t\). We introduce the numerical approximations
$$\begin{aligned} S_k \approx S(t_k),\quad I_k \approx I(t_k),\quad R_k \approx R(t_k),\quad D_k \approx D(t_k), \end{aligned}$$
for \(k=0,1,\ldots ,M\), and discretize the system (22) as follows
$$\begin{aligned} \begin{aligned} S_{k+1} =&S_k - \varDelta t \beta S_k \frac{I_{k_1}}{N}H(t_k-\tau _1),\\ I_{k+1} =&I_k + \varDelta t \left\{ \beta S_k \frac{I_{k_1}}{N}H(t_k-\tau _1) - (1-\alpha )\beta S_{k_2}\frac{I_{k_{12}}}{N}H(t_k-\tau _1-\tau _2)\right. \\&\left. - \alpha \beta S_{k_3}\frac{I_{k_{13}}}{N}H(t_k-\tau _1-\tau _3)\right\} ,\\ R_{k+1} =&R_k + \varDelta t (1-\alpha )\beta S_{k_2}\frac{I_{k_{12}}}{N}H(t_k-\tau _1-\tau _2),\\ D_{k+1} =&D_k + \varDelta t \alpha \beta S_{k_3}\frac{I_{k_{13}}}{N}H(t_k-\tau _1-\tau _3). \end{aligned} \end{aligned}$$
where the indexes \(k_1\), \(k_2\), \(k_{12}\), \(k_3\) and \(k_{13}\) are given by
$$\begin{aligned} \begin{aligned} k_j&= \max \left\{ 0, \left\lfloor \frac{t_k-\tau _j}{\varDelta t} \right\rfloor \right\} ,&\qquad j=1,2,3,\\ k_{1m}&= \max \left\{ 0, \left\lfloor \frac{t_k-\tau _1-\tau _m}{\varDelta t} \right\rfloor \right\} ,&\qquad m=2,3, \end{aligned} \end{aligned}$$
with \(\lfloor \cdot \rfloor \) the floor function.
Note that at each time step the condition
$$\begin{aligned} S_{k+1} + I_{k+1} + R_{k+1} + D_{k+1} = S_{k} + I_{k} + R_{k} + D_{k} \end{aligned}$$
is satisfied.
Results of the simulations
SIRD model
Concerning the SIRD model introduced in Sect. 2, we perform some numerical simulations by adopting the algorithm of Sect. 4. A crucial point is to fix the parameters entering the model. About the mortality \(\alpha \) we consider the infection fatality ratio (IFR), i.e. the ratio between the number of deaths from disease and the number of infected individuals, whose value is reported in the range 0.5–1% (World Health Organization 2020a). According to World Health Organization (2020b), the incubation period \(\tau _1\) is on average 5–6 days, but it can be as long as 14 days. For the healing and dead time \(\tau _2\) and \(\tau _3\) respectively, we remark that the commonly adopted criteria for discharging patients from isolation are the following: for symptomatic patients, 10 days after the symptom onset, plus at least 3 additional days without symptoms; for asymptomatic individuals, 10 days after positive test (World Health Organization 2020c). In all the simulations of the present paper we have assumed \(I_0=1\). The adopted values are reported in Table 1.
More controversial is to fix the contact frequency \(\beta \). It should be around \(\beta \simeq 1\) \(\hbox {d}^{-1}\) according to some estimations Peng et al. (2020) but it varies with time and as consequence of measures of social restriction by the authorities. For such a reason we have performed the simulations of the stochastic SIRD model for several values of this parameter: \(\beta =1.2,1,0.8,0.6\) \(\hbox {d}^{-1}\). The results are shown in Fig. 3 in the case of a population of 1000 individuals. The same cases have been also simulated when \(N=10,000\). The qualitative behavior is essentially the same, in particular the value of the maximum percentage of infected, but with a temporal dilation (Fig. 4).
Behavior of SIRD model with a population of 1000 individuals and \(\beta \) of 1.2 (top-left), 1 (top-right), 0.8 (bottom-left) and 0.6 (bottom-right). In the inset the values for large times are magnified
Behavior of SIRD model with a population of 10000 individuals and \(\beta \) of 1.2 (top-left), 1 (top-right), 0.8 (bottom-left) and 0.6 (bottom-right). In the inset the values for large times are magnified
Typical features of all the simulations are the reduction of susceptible people and a non monotone behaviour of the number of infected ones. The latter first increases and then tends to zero. The peak of the infected people is considered [see Fine et al. (2011)] as the value which represents the state when the herd immunity is reached without any actions by the authorities in charge for the health issues, also according to other models, e.g. see Ansumali et al. (2020). From a quantitative point of view our results indicate that the herd immunity is reached when about 70–80% of the population is infected which is a quite pessimistic foresight. Therefore, political strategies based on a pure herd immunity appear too risky because they could lead to a 0.6% of dead people due to the COVID-19 (see the inset in Fig. 3).
To avoid the drawbacks mentioned above, worldwide governments are assuming restrictions on free movement of people, the so-called lockdown, to contain the spread of the pandemic. In order to simulate the effect of a lockdown in our stochastic SIRD model, we suppose that the parameter \(\beta \) changes in \(\beta '\) whether the fraction of infected individuals reaches 10% of the total population size. Remember that \(1/\beta \) is the average contact time. So, if social distancing measures are adopted, they can be modeled as a reduction of \(\beta \), that is by extending the average contact time among the individuals. In Fig. 5 we show the curve of infected individuals for several values of \(\beta '\). It is evident that the peak of infected people lowers even if we have longer tails. This is quite realistic because the disease still remains but the number of recovered people increases in a slower way. Of course, if the aim is to alleviate the burden of hospitalized patients, the presence of a longer time to get the disappearance of the disease is a minor matter. Apparently with the lockdown it seems that the herd immunity is reached with a lower percentage of infected people than the case without lockdown.
Behavior of the SIRD model with a population of 1000 individuals and several values of \(\beta '\): number of infected individuals versus days. Note that the peak lowers by decreasing \(\beta '\)
To assess the validity of the model, in Fig. 6 we show a comparison between numerical results obtained by the stochastic SIRD model of Sect. 2 and the deterministic one presented in Sect. 5. There is a good agreement by obtaining a cross validation of both models. It is noteworthy that the numerical solution tends to the disease-free stationary critical points \(S^* = I^* =0\), \(R^*, D^* \in [0,1]\) with \(R^* + D^* = 1\).
Comparison between the stochastic (continuous lines) and deterministic (dashed lines) SIRD model in the case of \(N=1000\) and \(\beta =1\). Similar results are obtained with other values of \(\beta \)
SAI(L)RD model
As further improvements, here we also include the presence of asymptomatic and infected isolated individuals by presenting the results of the simulations obtained with the SAI(L)RD model. Regarding the parameters for the asymptomatic infections, \(\eta \) and \(\tau _2\), the literature reports that the proportion of people who become infected and remain asymptomatic throughout infection seems to be in the range 40–45% and they can transmit the virus for a period of about 14 days (Oran and Topol 2020). Concerning the parameter \(\tau _1\), we assume that symptomatic individuals can transmit the virus and they are isolated after an average period of 5–6 days because the illness is detected by tests. During the isolation they can heal or die with the same arguments of the SIRD model.
Since it is not known how long antibody responses will be maintained or whether they will provide protection from reinfection (Seow et al. 2020), we suppose that healed individuals may get a temporary or permanent immunity. We set the probability \(\lambda \) to have temporary immunity equal to 0.1 and the duration 90 days. The list of the adopted values is reported in Table 2. In the plots L is included in I.
In Fig. 7 we show the numerical solutions of the SAI(L)RD model in the case of \(N=1000\) and \(\beta =1\) \(\hbox {d}^{-1}\). The main distinctive feature with respect to the results obtained by the SIRD model is that after about 100 days we observe a new availability of susceptible individuals due to the loss of immunity. This along with a second wave of infection which, however, has a lower peak. Asymptotically we get again a disease-free situation with about 6% of dead people. We remark that by isolating the infected people the herd immunity is guaranteed by a peak of infected of about 30% to which about 20% of asymptomatic individuals must be added with a total of about 50% of people with disease. Again the strategy based on reaching the herd immunity can be deemed as to avoid because too costly in terms of hazard for the life of the population.
Behavior of SAI(L)RD model with a population of 1000 individuals and \(\beta \) set as 1
The above findings strongly support the need of the restrictive measures from a quantitative point of view. In order to analyze the effect of a lockdown we have also adopted the SAI(L)RD model. We start by taking \(\beta =1\) \(\hbox {d}^{-1}\) and then we set \(\beta =1/10\) \(\hbox {d}^{-1}\) when the fraction of infected individuals reaches 0.1. After 60 days we switch to \(\beta =1/3\) \(\hbox {d}^{-1}\) considering some restrictions still valid after the lockdown. The obtained numerical solutions are shown in Fig. 8.
Behavior of SAI(L)RD model with a population of 1000 individuals in the case of a lockdown
The values of the maximum for both infected and asymptomatic individuals are lower. Moreover, we observe that a second wave of infection is present which is less intense than the case of constant \(\beta \). However, the asymptotic number of the dead people with the lockdown is only slightly improved. The main effect of the lockdown is to alleviate the congestion in the intensive care because the infections are spread over a longer time. Finally, we would like to remark that if the probability of immunity loss \(\lambda \) is very high then periodic waves of persistent infection will show up for some years, as indicated in Fig. 9. Note that the asymptotic number of the dead people is about 1.7%.
Behavior based on the SAI(L)RD model with a population of 1000 individuals in the case of a probability of reinfection \(\lambda =\) 0.8
As last remark, if the size of the population is greater a dilation of the time is observed but the main features remain the same as the case of 1000 people.
Modeling the second wave of the spread of COVID-19 in Italy
To assess the validity of the proposed models, we will calibrate them on real data available in the literature. In particular we consider the daily collected values of infected, recovered and dead people in Italy during the second epidemic wave reported in Guidotti and Ardia (2020). The reason of the choice is that such data are provided with a stable number of tests per day and they are homogeneously collected in all the country. We consider the data in the temporal range: July 22, 2020–February 14, 2021. The raw data are shown in Fig. 10. The infected people at day zero of infection are about \(1.23\times 10^4\).
Raw data of COVID-19 spreading in Italy from February 24, 2020 until May 20, 2021. The dashed lines delimit the period time of the second infection wave whose data have been used for the calibration of the models
In order to have a more precise fitting between our models and data, in the investigated period of time we recognize three phases of diffusion. The first one is characterized by a low diffusion period starting from the end of July and lasting a couple of months. Subsequently, a higher diffusion is observed. On October 14, 2020 and, with more stringent actions, on November 6, 2020 the Italian government adopted contrast measures to contain the COVID-19 diffusion. Footnote 2 Finally, at the beginning of 2021 such measures have been reduced due to a weakening of the spread. For such reasons we consider the parameter \(\beta \) in the SIRD and SAI(L)RD models as a piece-wise constant function of time, as follows
$$\begin{aligned} \beta = \left\{ \begin{alignedat}{2}&\beta _0,&\qquad \text{ if }\quad t\le t_0,\\&\beta _1,&\qquad \text{ if }\quad t<t_0\le t_1,\\&\beta _2,&\qquad \text{ if }\quad t>t_1, \end{alignedat} \right. \end{aligned}$$
(25)
where \(\beta _0, \beta _1, \beta _2\) and \(t_0, t_1\) are parameters to be estimated.
Upon these considerations, we would like to find the parameters needed in the proposed models by a non-linear optimization procedure based on the Nelder–Mead method. The target is the minimization of the functional
$$\begin{aligned} J = \left\| I-I_{\text{ d }} \right\| _{L^2([T_1,T_2])} + \left\| R-R_{\text{ d }} \right\| _{L^2([T_1,T_2])} + \left\| D-D_{\text{ d }} \right\| _{L^2([T_1,T_2])}, \end{aligned}$$
(26)
being I, R and D the numerical values obtained with the model and \(I_{\text{ d }}\), \(R_{\text{ d }}\) and \(D_{\text{ d }}\) the same quantities given by data. \(T_1\) and \(T_2\) represent the first and the last day respectively of the considered time period of infection.
Regarding the SIRD model, the results are shown in Fig. 11 and the fitting parameters are reported in Table 3. The parameter \(t_0\) corresponds to October 3, 2020 while the parameter \(t_1\) to October 30, 2020. It is a reasonable time period for the high spreading before strong restrictions were adopted. There is a qualitative and quantitative good agreement between the simulation results and the data.
Comparison between SIRD model with parameters of Table 3 (continuous lines) and real Italian data (dashed lines)
Analogously, the optimization procedure described above is applied to the SAI(L)RD model. In particular, we adopt the model (25) and the quantities I and \(I_{\text{ d }}\) in the functional (26) have to be intended in this case as the total amount of symptomatic and asymptomatic infected individuals. The simulation results are shown in Fig. 12 and the fitting parameters are reported in Table 4. The good agreement with real data is searched for the total amount of infected people \(A+I\); the single curves A and I represent the projection due to the model. The optimization is not performed on parameters \(\eta \), \(\lambda \) and \(\tau _5\) that are the same of Table 2 instead.
Comparison between SAI(L)RD model with parameters of Table 4 (continuous lines) and real Italian data (dashed lines)
Again the model gives results which are in a qualitative and quantitative good agreement with the data.
Note that the parameters obtained with the calibration on the specific problem of this subsection are quite different with respect to those in Tables 1 and 2 , in particular the time delays \(\tau _2\) and \(\tau _3\) in the SIRD model and \(\tau _2\), \(\tau _3\) and \(\tau _4\) for the SAI(L)RD model. This can be ascribed to the fact that in the previous subsections the values of the parameters have been inferred from the worldwide literature. Here a specific situation has been analyzed; moreover, different criteria of defining the healing have been adopted in different countries.
Conclusions
Two stochastic models for simulating the evolution of the pandemic SARS-CoV-2 have been proposed. By using a Monte Carlo method, realistic situations have been investigated, obtaining insights about the possibility to get the herd immunity and the effects of measures as social distancing.
The models and the numerical approach have been tested by considering a deterministic version. The good agreement between the stochastic and deterministic results provides a cross validation.
The models are quite flexible and allow us an easy inclusion of the effects of a lockdown. As a specific test-case, we have successfully reproduced the second wave of the spread of COVID-19 in Italy during the period July 22, 2020–February 14, 2021.
The evolution we have considered does not take into account a campaign of vaccination but this can be included with a moderate additional effort.
Source: https://link.springer.com/article/10.1007/s00285-021-01657-4
0 Response to "The Dynamics of Distributions in Continuoustime Stochastic Models"
Post a Comment