Note That the Above Results Allows to Stochastic Processes in Continuous Time

• Journal List
• Springer Open Choice
• PMC8439375

J Math Biol. 2021; 83(4): 34.

Continuous-time stochastic processes for the spread of COVID-19 disease simulated via a Monte Carlo approach and comparison with deterministic models

Fabiana Calleri

Dipartimento di Matematica e Informatica, Università degli Studi di Catania, Viale Andrea Doria 6, 95125 Catania, Italy

Giovanni Nastasi

Dipartimento di Matematica e Informatica, Università degli Studi di Catania, Viale Andrea Doria 6, 95125 Catania, Italy

Vittorio Romano

Dipartimento di Matematica e Informatica, Università degli Studi di Catania, Viale Andrea Doria 6, 95125 Catania, Italy

Received 2021 Mar 16; Revised 2021 Jun 15; Accepted 2021 Aug 29.

Abstract

Two stochastic models are proposed to describe the evolution of the COVID-19 pandemic. In the first model the population is partitioned into four compartments: susceptible S, infected I, removed R and dead people D. In order to have a cross validation, a deterministic version of such a model is also devised which is represented by a system of ordinary differential equations with delays. In the second stochastic model two further compartments are added: the class A of asymptomatic individuals and the class L of isolated infected people. Effects such as social distancing measures are easily included and the consequences are analyzed. Numerical solutions are obtained with Monte Carlo simulations. Quantitative predictions are provided which can be useful for the evaluation of political measures, e.g. the obtained results suggest that strategies based on herd immunity are too risky. Finally, the models are calibrated on data referring to the second wave of infection in Italy.

Keywords: COVID-19, Stochastic process, Epidemic model, Monte Carlo simulation

Introduction

The pandemic of COVID-19 has scourged the world since the beginning of 2020.¹ 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 ∈ [ 0 , + ∞ [ 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 β > 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 α ∈ [ 0 , 1 ] to die and 1 - α 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. τ 1 is the delay time for getting infected and Δ t ∼ β - 1 is the frequency of the contacts between the individuals of S and I

Moreover, let τ 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 τ 2 and τ 3 be the heal and dead time respectively, counted after the incubation period. First, we will assume that τ 1 , τ 2 , τ 3 are constant; in a second step this assumption will be relaxed to simulated effects as those related to the lockdown.

Let us consider S , I , R , D as four random variables which take values in 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

with I 0 ≪ 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 S , I , corresponding to the infection mechanism. We think of the two random variables as unsynchronized in time. In particular, we set

and we consider the new set S , I ~ and suppose that at time t the random variables take values s , i , that is

After a small period of time Δ t > 0 , Δ t ≪ min i = 1 , … , 3 τ i the state of system changes in

S ( t + Δ t ) , I ~ ( t + Δ t ) = s + m , i + n ,

2

being m , n ∈ - 1 , 0 , + 1 .

The transition probability is defined as follows

p ( s , i ) → ( s + m , i + n ) ( Δ t ) = P ( S ( t + Δ t ) , I ~ ( t + Δ t ) = s + m , i + n | S ( t ) , I ~ ( t ) = s , i ) .

In this way we define a continuous time Markov chain and the transition probability can be written as

p ( s , i ) → ( s + m , i + n ) ( Δ t ) = β s i N Δ t + o ( Δ t ) , ( m , n ) = ( - 1 , + 1 ) 1 - β s i N Δ t + o ( Δ t ) , ( m , n ) = ( 0 , 0 ) o ( Δ t ) , otherwise .

The second set of random variables is I , R , coupled to the first one. We set

I ~ ~ ( t ) = I ( t + τ 1 + τ 2 ) , R ~ ~ ( t ) = R ( t + τ 1 + τ 2 )

and consider the couple I ~ ~ , R ~ ~ . Let us suppose that at time t the random variables I ~ ~ , R ~ ~ take values (i,r), that is

After a small period of time Δ t > 0 the state of the system changes in

I ~ ~ ( t + Δ t ) , R ~ ~ ( t + Δ t ) = ( i + n , r + u ) ,

4

being n , u ∈ - 1 , 0 , + 1 . Moreover, since a healing at time t + τ 1 + τ 2 is related to an infection in the past interval [ t + τ 1 , t + τ 1 + Δ t ] , we need to know the values assumed by the random variable I ~ at t and t + Δ t . In this case the transition probability is defined as follows

p ( i , r ) → ( i + n , r + u ) ( Δ t ) = P ( I ~ ~ ( t + Δ t ) , R ~ ~ ( t + Δ t ) = i + n , r + u | I ~ ~ ( t ) , R ~ ~ ( t ) = i , r , I ~ ( t ) = j , I ~ ( t + Δ t ) = j + 1 ) .

Therefore, a non-Markovian continuous time stochastic process is defined with the transition probability

p ( i , r ) → ( i + n , r + u ) ( Δ t ) = ( 1 - α ) Δ t + o ( Δ t ) , ( n , u ) = ( - 1 , + 1 ) 1 - ( 1 - α ) Δ t + o ( Δ t ) , ( n , u ) = ( 0 , 0 ) o ( Δ t ) , otherwise .

Finally, the third set of random variables is I , D , which is also coupled to the first one. Now we set

I ^ ( t ) = I ( t + τ 1 + τ 3 ) , D ^ ( t ) = D ( t + τ 1 + τ 3 )

and consider the couple I ^ , D ^ . Let us suppose that at time t the random variables take values (i,d), that is

After a small period of time Δ t > 0 the state of the system changes in

I ^ ( t + Δ t ) , D ^ ( t + Δ t ) = ( i + n , d + v ) ,

6

being n , v ∈ - 1 , 0 , + 1 . Moreover, since a death at time t + τ 1 + τ 3 is related to an infection in the past interval [ t + τ 1 , t + τ 1 + Δ t ] , we need to know the values assumed by the random variable I ~ at t and t + Δ t . In this case the transition probability is defined as follows

p ( i , d ) → ( i + n , d + v ) ( Δ t ) = P ( I ^ ( t + Δ t ) , D ^ ( t + Δ t ) = i + n , d + v | I ^ ( t ) , D ^ ( t ) = i , d , I ~ ( t ) = j , I ~ ( t + Δ t ) = j + 1 ) .

Therefore, a non-Markovian continuous time stochastic process is defined with the transition probability

p ( i , d ) → ( i + n , d + v ) ( Δ t ) = α Δ t + o ( Δ t ) , ( n , v ) = ( - 1 , + 1 ) 1 - α Δ t + o ( Δ t ) , ( n , v ) = ( 0 , 0 ) o ( Δ t ) , otherwise .

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

1. 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.

2. We suppose there exists a probability η ∈ [ 0 , 1 ] to be asymptomatic and, consequently, 1 - η is the probability to show symptoms once an individual has been infected. Moreover, we assume that, after a certain time τ 1 , a symptomatic individual is recognized and isolated into a subclass, called L, of lonely individuals of I.

3. Lonely individuals are not infectious anymore. Asymptomatic people heal after a time τ 2 . Lonely symptomatic individuals can die with probability α in a time τ 4 or heal with probability 1 - α in a time τ 3 .

4. Finally, removed individuals become immune to COVID-19 for a short period of time or forever. Let λ ∈ [ 0 , 1 ] be the probability that COVID-19 confers a short immunity of time length τ 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 S , A , I , L , R , D as a six dimensional random variable which can assume values in N 0 6 and depends on time t ≥ 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

  ( S ( t ) , A ( t ) , I ( t ) ) = ( s , a , i ) .

  7

  After a small period of time Δ t > 0 , Δ t ≪ min i = 1 , … , 5 τ i , the state of the system changes in

  ( S ( t + Δ t ) , A ( t + Δ t ) , I ( t + Δ t ) ) = ( s + m , a + k , i + n ) ,

  8

  with m , k , n ∈ - 1 , 0 , 1 .

  The transition probability is given by

  p ( s , a , i ) → ( s + m , a + k , i + n ) ( Δ t ) = P ( S ( t + Δ t ) , A ( t + Δ t ) , I ( t + Δ t ) = s + m , a + k , i + n | S ( t ) , A ( t ) , I ( t ) = s , a , i ) .

  In this way we define a time continuous Markov chain with transition probabilities which can be written as

  p ( s , a , i ) → ( s + m , a + k , i + n ) ( Δ t ) = η β s i N Δ t + o ( Δ t ) , if ( m , k , n ) = ( - 1 , + 1 , 0 ) ( 1 - η ) β s i N Δ t + o ( Δ t ) , if ( m , k , n ) = ( - 1 , 0 , + 1 ) 1 - β s i N Δ t + o ( Δ t ) , if ( m , k , n ) = ( 0 , 0 , 0 ) o ( Δ t ) , otherwise .

• In our model, we suppose that all the infected individuals will be isolated after a time τ 1 . Let us consider (I,L). We set

  I ~ ( t ) = I ( t + τ 1 ) , L ~ ( t ) = L ( t + τ 1 ) .

  In this case a lone infected individual at time t + Δ t + τ 1 is related to a symptomatic infection in the past interval [ t , t + Δ t ] . Regarding the transition probability

  p ( i , l ) → ( i + n , l + j ) ( Δ t ) = P ( I ~ ( t + Δ t ) , L ~ ( t + Δ t ) = i + n , l + j | I ~ ( t ) , L ~ ( t ) = i , l , I ( t ) = f , I ( t + Δ t ) = f + 1 )

  for some nonnegative integer f, we have

  p ( i , l ) → ( i + n , l + j ) ( Δ t ) = 1 , if ( n , j ) = ( - 1 , + 1 ) 0 , otherwise .

• Now we consider the pair of random variables (L,R). To define the healing process we set

  L ~ ~ ( t ) = L ( t + τ 1 + τ 3 ) , R ~ ~ ( t ) = R ( t + τ 1 + τ 3 )

  and suppose that

  After a small period of time Δ t > 0 the state of the system changes in

  ( L ~ ~ ( t + Δ t ) , R ~ ~ ( t + Δ t ) ) = ( l + j , r + u )

  with j , u ∈ - 1 , 0 , 1 . Moreover, since in this case a healing at time t + Δ t + τ 1 + τ 3 is related to an infection with symptoms in the past interval [ t , t + Δ t ] , we need to know the values assumed by the random variable I at t and t + Δ t . This is also equivalent to knowing the values assumed by the random variable L at t + τ 1 and t + Δ t + τ 1 . In this case the transition probability is defined as

  p ( l , r ) → ( l + j , r + u ) ( Δ t ) = P ( L ~ ~ ( t + Δ t ) , R ~ ~ ( t + Δ t ) = l + j , r + u | L ~ ~ ( t ) , R ~ ~ ( t ) = l , r , L ( t + τ 1 ) = f , L ( t + Δ t + τ 1 ) = f + 1 ) ,

  for some nonnegative integer f. In this way a time continuous non-Markovian stochastic process is defined whose transition probability can be written as

  p ( l , r ) → ( l + j , r + u ) ( Δ t ) = ( 1 - α ) Δ t + o ( Δ t ) , if ( j , u ) = ( - 1 , + 1 ) 1 - ( 1 - α ) Δ t + o ( Δ t ) , if ( j , u ) = ( 0 , 0 ) o ( Δ t ) , otherwise .

• With the same arguments, a dead process is described by ( L ^ , D ^ ) where

  L ^ ( t ) = L ( t + τ 1 + τ 4 ) , D ^ ( t ) = D ( t + τ 1 + τ 4 ) .

  In this case a death at time t + Δ t + τ 1 + τ 4 is related to an infection with symptoms in the past interval [ t , t + Δ t ] and thus the gain of one unit to the variable L in [ t + τ 1 , t + Δ t + τ 1 ] . The transition probability

  p ( l , d ) → ( l + j , d + v ) ( Δ t ) = P ( L ^ ( t + Δ t ) , D ^ ( t + Δ t ) = l + j , d + v | L ^ ( t ) , D ^ ( t ) = l , d , L ( t + τ 1 ) = f , L ( t + Δ t + τ 1 ) = f + 1 ) ,

  for some nonnegative integer f, is given by

  p ( l , d ) → ( l + j , d + v ) ( Δ t ) = α Δ t + o ( Δ t ) , if ( j , v ) = ( - 1 , + 1 ) 1 - α Δ t + o ( Δ t ) , if ( j , v ) = ( 0 , 0 ) o ( Δ t ) , otherwise .

• A further process we introduce is the healing of an asymptomatic individual. Let us consider (A,R). We set

  A ^ ( t ) = A ( t + τ 2 ) , R ^ ( t ) = R ( t + τ 2 ) .

  In this case a healing at time t + τ 2 + Δ t is related to an asymptomatic infection in the past interval [ t , t + Δ t ] . The transition probability

  p ( a , r ) → ( a + k , r + u ) ( Δ t ) = P ( A ^ ( t + Δ t ) , R ^ ( t + Δ t ) = a + k , r + u | A ^ ( t ) , R ^ ( t ) = a , r , A ( t ) = f , A ( t + Δ t ) = f + 1 ) ,

  for some nonnegative integer f, reads

  p ( a , r ) → ( a + k , r + u ) ( Δ t ) = 1 , if ( k , u ) = ( - 1 , + 1 ) 0 , otherwise .

• 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

  S ^ ^ ( t ) = S ( t + τ 5 ) , R ^ ^ ( t ) = R ( t + τ 5 )

  Let us suppose that at time t the two-dimensional random variable ( S ^ ^ , R ^ ^ ) takes the value (s,r). After a small period of time Δ t > 0 the state of the system changes in

  ( S ^ ^ ( t + Δ t ) , R ^ ^ ( t + Δ t ) ) = ( s + m , r + u ) ,

  with m , u ∈ - 1 , 0 , 1 . Moreover, since in this case a healed individual can come back to the class of susceptible people after a certain time τ 5 , it is needed to know the values assumed by the random variable R at time t and t + Δ t . In this case the transition probability

  p ( s , r ) → ( s + m , r + u ) ( Δ t ) = P ( S ^ ^ ( t + Δ t ) , R ^ ^ ( t + Δ t ) = s + m , r + u | S ^ ^ ( t ) , R ^ ^ ( t ) = s , r , R ( t ) = f , R ( t + Δ t ) = f + 1 ) ,

  for some nonnegative integer f, is given by

  p ( s , r ) → ( s + m , r + u ) ( Δ t ) = λ + o ( Δ t ) , if ( m , u ) = ( + 1 , - 1 ) 1 - λ + o ( Δ t ) , if ( m , u ) = ( 0 , 0 ) o ( Δ t ) , otherwise .

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 τ 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

X ( t ) = X 1 ( t ) , X 2 ( t ) , … , X N ( t ) ∈ D N ,

9

for t ≥ 0 , where D = - 1 , 0 , 1 , 2 and N is the population size. We indicate by X i ( t ) ∈ D for i = 1 , … , 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 β , that is

ξ 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 τ 1 , the individual i changes its state in infectious,

At time t j , 1 + τ 1 , the destiny d ∈ - 1 , 2 of the new infectious is established too accordingly to a Bernoulli distribution with probability α , that is d ∼ B ( 1 , α ) . If the destiny is to heal then after a time τ 2 the individual state changes from infectious to recovered; if the destiny is to die then after a time τ 3 the state of the i-th individual changes from infectious to dead:

X i ( t j , 1 + τ 1 + τ 2 ) = 2 if d = 2 , X i ( t j , 1 + τ 1 + τ 2 ) = - 1 if d = - 1 .

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

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 , … , k , be the time at which the algorithm ends at the rth simulation. We set

the average of the final process times. After introducing the error as

as stopping criterion we adopt

Here 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

where E = - 1 , 0 , 1 , 2 , 3 , 4 . 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 ∈ 1 , 3 by a Bernoulli distribution having probability η ∈ [ 0 , 1 ] . If s = 3 the individual is asymptomatic and, after the time τ 2 , she/he will heal, and therefore

If s = 1 the person is symptomatic infected and, after the time τ 1 , she/he will be isolated, that is

Moreover, the destiny d ∈ - 1 , 2 of such an individual is established according to a Bernoulli distribution with probability α . If the destiny is to heal then after a time τ 3 the individual state changes from infectious to recovered; if the destiny is to die then after a time τ 4 the individual state changes from infectious to dead. That evolution can be described as follows

X i ( t j , 1 + τ 1 + τ 3 ) = 2 if d = 2 , X i ( t j , 1 + τ 1 + τ 4 ) = - 1 if d = - 1 .

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 λ ∈ ( 0 , 1 ) . Furthermore, if the immunity is lost, one has two possibilities:

if the individual suffered from an asymptomatic infection we set

if she/he suffered from a symptomatic infection we set

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

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 τ 1 , the number of susceptible people decreases by a quantity depending on the amount of infected at a previous time t - τ 1 ; the amount of recovered and dead people after a time τ 2 and τ 3 respectively is proportional to the amount of people who have been infected at the previous time t - τ 2 and t - τ 3 respectively. From those considerations, we propose the following model

where H ( · ) 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

and some simple algebraic manipulations, the system can be written more explicitly as

The system must be augmented assigning the functions

S ( t ) = Φ 1 ( t ) t ∈ [ - max ( τ 2 , τ 3 ) , 0 ] , I ( t ) = Φ 2 ( t ) t ∈ [ - max ( τ 1 + τ 2 , τ 1 + τ 3 ) , 0 ] .

As customary we assume that Φ 1 and Φ 2 are continuous in the considered intervals. Regarding the other variables it is realistic to take R ( 0 ) = D ( 0 ) = 0 . Specifically we assume Φ 1 ( t ) = S 0 ∈ ] 0 , 1 [ and Φ 2 ( t ) = I 0 = 1 - S 0 ∈ ] 0 , 1 [ with I 0 ≪ 1 . Therefore at t = 0 we have

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)

holds. In fact, summing up the equations (22a)-(22d) one has

where Y = S + I + R + D . Since Y ( 0 ) = 1 , it follows that Y ( t ) = 1 ∀ t > 0 .

The equilibrium points of the system (22) are the following.

• Endemic solutions:

  with I ∗ , R ∗ , D ∗ ∈ [ 0 , 1 ] satisfying I ∗ + R ∗ + D ∗ = 1 .

• Disease-free solutions:

  with S ∗ , R ∗ , D ∗ ∈ [ 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

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

λ λ - β S ∗ e - τ 1 λ + ( 1 - α ) e - ( τ 1 + τ 2 ) λ + α e - ( τ 1 + τ 3 ) λ = 0 ,

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 < … < t M = T max of constant time step Δ t . We introduce the numerical approximations

S k ≈ S ( t k ) , I k ≈ I ( t k ) , R k ≈ R ( t k ) , D k ≈ D ( t k ) ,

for k = 0 , 1 , … , M , and discretize the system (22) as follows

S k + 1 = S k - Δ t β S k I k 1 N H ( t k - τ 1 ) , I k + 1 = I k + Δ t β S k I k 1 N H ( t k - τ 1 ) - ( 1 - α ) β S k 2 I k 12 N H ( t k - τ 1 - τ 2 ) - α β S k 3 I k 13 N H ( t k - τ 1 - τ 3 ) , R k + 1 = R k + Δ t ( 1 - α ) β S k 2 I k 12 N H ( t k - τ 1 - τ 2 ) , D k + 1 = D k + Δ t α β S k 3 I k 13 N H ( t k - τ 1 - τ 3 ) .

where the indexes k 1 , k 2 , k 12 , k 3 and k 13 are given by

k j = max 0 , t k - τ j Δ t , j = 1 , 2 , 3 , k 1 m = max 0 , t k - τ 1 - τ m Δ t , m = 2 , 3 ,

with ⌊ · ⌋ the floor function.

Note that at each time step the condition

S k + 1 + I k + 1 + R k + 1 + D k + 1 = S k + I k + R k + D k

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 α 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 τ 1 is on average 5–6 days, but it can be as long as 14 days. For the healing and dead time τ 2 and τ 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 Table1.

Table 1

Parameters adopted for the simulation of the SIRD model

Parameter	Value
α	0.006
τ 1	6 d
τ 2	12 d
τ 3	13 d

The τ 's are expressed in days (d)

More controversial is to fix the contact frequency β . It should be around β ≃ 1 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: β = 1.2 , 1 , 0.8 , 0.6 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 β 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 β 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 β changes in β ′ whether the fraction of infected individuals reaches 10% of the total population size. Remember that 1 / β is the average contact time. So, if social distancing measures are adopted, they can be modeled as a reduction of β , 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 β ′ . 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 β ′ : number of infected individuals versus days. Note that the peak lowers by decreasing β ′

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 ∗ ∈ [ 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 β = 1 . Similar results are obtained with other values of β

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, η and τ 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 τ 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 λ to have temporary immunity equal to 0.1 and the duration 90 days. The list of the adopted values is reported in Table2. In the plots L is included in I.

Table 2

Parameters adopted for the simulation of the SAI(L)RD model

Parameter	Value
α	0.006
η	0.4
λ	0.1
τ 1	6 d
τ 2	14 d
τ 3	12 d
τ 4	13 d
τ 5	90 d

In Fig.7 we show the numerical solutions of the SAI(L)RD model in the case of N = 1000 and β = 1 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 β 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 β = 1 d - 1 and then we set β = 1 / 10 d - 1 when the fraction of infected individuals reaches 0.1. After 60 days we switch to β = 1 / 3 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 β . 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 λ 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 λ = 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 × 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.² 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 β in the SIRD and SAI(L)RD models as a piece-wise constant function of time, as follows

β = β 0 , if t ≤ t 0 , β 1 , if t < t 0 ≤ t 1 , β 2 , if t > t 1 ,

25

where β 0 , β 1 , β 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

J = I - I d L 2 ( [ T 1 , T 2 ] ) + R - R d L 2 ( [ T 1 , T 2 ] ) + D - D d L 2 ( [ T 1 , T 2 ] ) ,

26

being I, R and D the numerical values obtained with the model and I d , R d and D 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 Table3. 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 Table3 (continuous lines) and real Italian data (dashed lines)

Table 3

Fitting parameters for the SIRD model

Parameter	Value
α	0.0405
τ 1	7.62 d
τ 2	28.9 d
τ 3	30.0 d
β 0	0.0493 d - 1
β 1	0.1164 d - 1
β 2	0.0422 d - 1
t 0	73 d
t 1	100 d

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 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 Table4. 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 η , λ and τ 5 that are the same of Table2 instead.

Comparison between SAI(L)RD model with parameters of Table4 (continuous lines) and real Italian data (dashed lines)

Table 4

Fitting parameters for the SAI(L)RD model

Parameter	Value
α	0.0517
τ 1	7.66 d
τ 2	31.0 d
τ 3	31.1 d
τ 4	30.1 d
β 0	0.0768 d - 1
β 1	0.3446 d - 1
β 2	0.1026 d - 1
t 0	73 d
t 1	102 d

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 Tables1 and ​ 2 , in particular the time delays τ 2 and τ 3 in the SIRD model and τ 2 , τ 3 and τ 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.

Acknowledgements

The authors G.N. and V.R. acknowledge the support from INdAM (GNFM).

Funding

Open access funding provided by Università degli Studi di Catania within the CRUI-CARE Agreement.

Footnotes

Contributor Information

Fabiana Calleri, Email: ti.oohay@irellacanaibaf.

Giovanni Nastasi, Email: ti.tcinu@isatsan.g.

Vittorio Romano, Email: ti.tcinu.imd@onamor.

References

• Allen LJS. An introduction to stochastic processes with applications to biology. 2. Boca Raton: CRC Press, Taylor & Francis Group; 2010. [Google Scholar]
• Ansumali S, Kaushal S, Kumar A, Prakash MK, Vidyasagar M. Modelling a pandemic with asymptomatic patients, impact of lockdown and herd immunity, with applications to SARS-CoV-2. Annu Rev Control. 2020;50:432–447. doi: 10.1016/j.arcontrol.2020.10.003. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
• Bagarello F, Gargano F, Roccati F. Modeling epidemics through ladder operators. Chaos Solitons Fractals. 2020;140:110193. doi: 10.1016/j.chaos.2020.110193. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
• Calafiore GC, Novara C, Possieri C. A time-varying SIRD model for the COVID-19 contagion in Italy. Annu Rev Control. 2020;50:361–372. doi: 10.1016/j.arcontrol.2020.10.005. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
• Capasso V, Bakstein D. An introduction to continuous-time stochastic processes: theory, models, and applications to finance, biology, and medicine. Basel: Birkhäuser; 2015. [Google Scholar]
• Faranda D, Alberti T. Modeling the second wave of COVID-19 infections in France and Italy via a stochastic SEIR model. Chaos. 2020;30:111101. doi: 10.1063/5.0015943. [PubMed] [CrossRef] [Google Scholar]
• Feng S, Chen M, Zhan N, Fränzle M, Xue B. Taming delays in dynamical systems. Unbounded verification of delay differential equations. In: Dillig I, Tasiran S, editors. "Computer aided verification", CAV 2019. Lecture notes in computer science. Cham: Springer; 2019. [Google Scholar]
• Fine P, Eames K, Heymann DL. "Herd immunity": a rough guide. Clin Infect Dis. 2011;52(7):911–916. doi: 10.1093/cid/cir007. [PubMed] [CrossRef] [Google Scholar]
• Giordano G, Blanchini F, Bruno R, Colaneri P, Di Filippo A, Di Matteo A, Colaneri M. Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy. Nat Med. 2020;26:855–860. doi: 10.1038/s41591-020-0883-7. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
• Guidotti E, Ardia D. COVID-19 data hub. J Open Sour Softw. 2020 doi: 10.21105/joss.02376. [CrossRef] [Google Scholar]
• Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proc Roy Soc Lond Ser A. 1927;115(772):700–721. doi: 10.1098/rspa.1927.0118. [CrossRef] [Google Scholar]
• Murray JD. Mathematical biology I. An introduction. New York: Springer-Verlag; 2002. [Google Scholar]
• Oran DP, Topol EJ. Prevalence of asymptomatic SARS-CoV-2 infection. Ann Intern Med. 2020;173:362–367. doi: 10.7326/M20-3012. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
• Peng L, Yang W, Zhang D, Zhuge C, Hong L (2020) Epidemic analysis of COVID-19 in China by dynamical modeling. arXiv:2002.06563
• Rihan FA, Alsakaji HJ, Rajivganthi C. Stochastic SIRC epidemic model with time-delay for COVID-19. Adv Differ Equ. 2020;2020:502. doi: 10.1186/s13662-020-02964-8. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
• Seow J, Graham C, Merrick B, et al. Longitudinal observation and decline of neutralizing antibody responses in the three months following SARS-CoV-2 infection in humans. Nat Microbiol. 2020;5:1598–1607. doi: 10.1038/s41564-020-00813-8. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
• Whittle P. The outcome of a stochastic epidemic-a note on Bailey's paper. Biometrika. 1955;42(1/2):116–122. doi: 10.2307/2333427. [CrossRef] [Google Scholar]
• World Health Organization (2020a) Estimating mortality from COVID-19, scientific brief, 4 August 2020. WHO Reference Number: WHO-2019-nCoV-Sci\_Brief-Mortality-2020.1
• World Health Organization (2020b) Transmission of SARS-CoV-2: implications for infection prevention precautions, scientific brief, 9 July 2020. WHO Reference Number: WHO/2019-nCoV/Sci\_Brief/Transmission\_modes/2020.3
• World Health Organization (2020c) Criteria for releasing COVID-19 patients from isolation, scientific brief, 17 June 2020. WHO Reference Number: WHO/2019-nCoV/Sci\_Brief/Discharge\_From\_Isolation/2020.1
• Zhang Y, You C, Cai Z, Sun J, Hu W, Zhou X-H. Prediction of the COVID-19 outbreak in China based on a new stochastic dynamic model. Sci Rep. 2020;10:21522. doi: 10.1038/s41598-020-76630-0. [PMC free article] [PubMed] [CrossRef] [Google Scholar]
• Zhou P, Yang XL, Wang XG, et al. A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature. 2020;579:270–273. doi: 10.1038/s41586-020-2012-7. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

fortunehimsen.blogspot.com

Source: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8439375/

Share this post