A stochastic epidemic model of COVID-19 disease

To model the evolution of diseases with extended latency periods and the presence of asymptomatic patients like COVID-19, we define a simple discrete time stochastic SIR-type epidemic model. We include both latent periods as well as the presence of quarantine areas, to capture the evolutionary dynamics of such diseases.


Introduction
There exists a wide class of mathematical models that analyse the spread of epidemic diseases, either deterministic or stochastic, and may involve many factors such as infectious agents, mode of transmission, incubation periods, infectious periods, quarantine periods, etcetera (Allen 2003; Anderson and May 1991;Bailey 1975; Daley and Gani 1999;Diekmann et al. 2013).
A basic model of infectious disease population dynamics, consisting of susceptible (S), infective (I) and recovered (R) individuals were first considered in a deterministic model by Kermack and McKendric (1927). Since then, various epidemic deterministic models have been developed, with or without a time delay (see e.g McCluskey 2009 and Huang et al. 2010). At the same time, many stochastic models have been considered: discrete time models (see e.g Tuckwell and Williams 2007 ; Oli et al. 2006), continuous time Markov chain models and diffusion models (see e.g. Mode and Sleeman 2000). The models obtained in these three categories are of increasing mathematical complexity and allow to study several aspects of the epidemics.
Even if the discrete-time models are the simplest ones, they may help to better define the basic principles of the contagion and to avoid the constraints due to the own definitions of the more sophisticated models.
In this paper we will adapt a simple SIR-type model proposed by Ferrante et al. (2016) and we will divide the population into several classes to better describe the evolution of the COVID epidemic. Here we will need to include latency periods, the presence of asymptomatic patients and different level of isolation. To model the evolution of the epidemic, we will describe the evolution of every single individual in the population, modelling the probability on every day to be infected and, once infected, the exact evolution of its disease until the possible recovery or the death. The construction of the theoretical model is carried out in Section 2, where we are able to compute the probability of contagion and an estimate of the basic reproduction number R 0 . Then, in Section 3 we answer to five research questions by using a simulation of the evolution of the disease. The use of the simulation is justified by the complexity of the model, that prevent to carry out any further exact computation. We are able to see that to stop the epidemic is fundamental to start early with a severe quarantine and that a late starting date or a more soft quarantine makes this procedure almost useless. Moreover, to determine the quarantine it is very important to know the level of infectivity of the asymptomatic, since the more infectious they are, the more important is the quarantine. Finally, as expected, the group immunity plays a very important role to prevent the development of the disease.

A SIR type discrete-time stochastic epidemic model
To model the evolution of epidemics, Tuckwell and Williams (2007) proposed a simple stochastic SIR-type model based on a discrete-time Markovian approach, later generalized by Ferrante et al. (2016) with a SEIHR model. These models, despite their simplicity, are very unrealistic to catch the characteristic of the COVID-19 disease and for this reason in this paper we introduce a more complex system, that we call SEIAHCRD, that better describes this new disease.
Assume that the population size is fixed and equal to n, and that the time is discrete, with the unit for the duration of an epoch one day. Every individual, marked by an integer between 1 and n, belongs to one of the following 8 classes: • the class S includes the individuals susceptible to the disease and never infected before; • the class E includes the individuals in a latency period, i.e. individuals that have been infected but that are still not infectious or sick; • the class I includes the infectious individuals, that are not yet sick, but that will develop later the disease; • the class A includes the infectious individuals, that are not yet sick, but that will NOT develop later the disease, usually referred as Asymptomatic; • the class H includes the infectious individuals, that are sick, but with light symptoms and therefore at home quarantine; • the class C includes the infectious individuals, that are sick and with severe symptoms, and most of the time are hospitalized; • the class D includes the deceased individuals; • the class R includes the recovered individuals.
Assuming that we start at time 0, we will define for any individual i ∈ {1, . . . , n} the family of stochastic processes if the individual i at time t belongs to the class Ξ, where Ξ is equal to S, E, I, A, H, C, D or R, and 0 otherwise. In this way, the total number of individuals in the class Ξ at time t ≥ 0 will be denote by Y Ξ (t) and will be equal to Let us now fix the main assumption on the evolution of the epidemic. At time 0 all the individuals are in S, but one in class I or A and that the evolution of the contagion follows these rules: 1. Daily encounters: each individual i, over (t, t+1], will encounter a number of other individuals equal to N i (t) which we will assume to be a deterministic value; 2. Contagion probability: if an individual who has never been diseased up to and including time t, encounters an individual in (t, t + 1] who belongs to the class I or A, then, independently of the results of other encounters, the encounter results in transmission of the disease with probability q I and q A , respectively.
3. Permanence in the classes: any individual, but one, starts from class S. Once infected he/she moves to class E and so on according to the graph below. The time spent in the classes E, I and A are of r E , r I and r A consecutive days, respectively. These values can be considered deterministic or stochastic. Any individual who enters the class H remains in this class for r HC consecutive days with probability α or for r HR consecutive days with probability 1 − α. Any individual who enters the class C remains in this class for r CD consecutive days with probability λ or for r CR consecutive days with probability 1 − λ. As before, the values of these four numbers of consecutive days can be considered deterministic or stochastic. To conclude, we assume that the individuals once in class R reamin there forever, the same as for the class D.

4.
Transitions between the classes: we assume that any individual can moves between the classes according to the following graph Here β, µ, α and λ denotes the transition probabilities and the transitions occur at the end of the permanence time spent by the individual in the previous class. Note that µ, α, λ are parameters that depends only on the specific nature of the disease, while β depends on this and the number of individuals in the classes I and A.
Any individual, once infected with probability β, follows one of the four paths described here: (a) he/she transits through the states E, I, H, C, D, where he/she remains, respectively, for r E , r I , r HC and r CD days, after which he/she dies and moves to class D.
(b) he/she transits through the states E, I, H, C, R, where he/she remains, respectively, for r E , r I , r HC and r CR days, after which he/she becomes immune and moves to class R.
(c) he/she transits through the states E, I, H, R, where he/she remains, respectively, for r E , r I and r HR days, after which he/she becomes immune and moves to class R.
(d) he/she transits through the states E, A, R, where he/she remains, respectively, for r E and r A days, after which he/she becomes immune and moves to class R.
It is immediate to see that the probability to follow any of these four paths is equal to, respectively, In order to evaluate the probability β of contagion at time t of an individual in class S, we will start defining the probabilities of meeting an individual in the classes S, E, I, A and R as equal. Note that the individuals in classes H and C are in total quarantine and that it is not possible to meet them and that the individuals in class D are removed from the system. So, we will deal with three possibly different encounter probabilities: • p I the probability of meeting an individual belonging to the class I; • p A the probability of meeting an individual belonging to class A; • p S the probability of meeting an individual that is not infectious, that is he/she belongs to the classes S, E or R.
Assuming that the probability of meeting any individual is uniform and independent from the above defined classes, we define these probabilities as In the above formulas, y I , y A , y H , y C and y D denote the number of individuals in classes I, A, H, C and D, respectively. Denoting by j iI , j iA the number of meetings of the i-th individual at time t with individuals in the classes I and A, respectively, the probability to meet this proportion of individuals is where N i (t) denotes the daily encounters of the individual i. We can easily derive the probability of contagion where q I and q A denote the probability of transmission of the specific disease for individuals in classes I and A, respectively, which are usually different. Then the probability of contagion at time t + 1 of a single individual is equal to Substituting (1), we then get when y H + y C + y D < n − 1, while β i = 0 when y H + y C + y D = n − 1. As done by Tuckwell and Williams (2007), we can use these formulas to simulate the spread of an epidemic under these general assumptions. Some results, similar to those presented in Tuckwell and Williams (2007), can be found in Ferrante et al. (2016).
To conclude, let us now consider the basic reproduction number R 0 , i.e. the expected number of secondary cases produced by an infectious individual during its period of infectiousness (see Diekmann et al. 1990). In the present model, this refers to individuals that transits through the classes I or A. Let us recall the threshold value of R 0 , which establishes that an infection persists only if R 0 > 1. As for the SIR-model proposed by Tuckwell and Williams we are not able to derive the exact explicit value of R 0 , but it is possible to extend their results, when N i (t) ≡ N , for all i and t, obtaining that Note that the value of R 0 computed above is the basic reproduction number at the beginning of the disease, when there is one infectious individual that has N contacts in a population with n − 1 susceptible inhabitants. When the disease is in an advanced development, keeping the total number of contacts N some of them won't be with susceptible individuals and the number of cases produced by an infectious individual will be smaller. At time t, let us call S(t) the number of susceptible individuals (class S) and X(t) the number of individuals removed from the population (members of the classes H, C and D). Then, given an individual in class I, set Z i (t) his number of contacts with susceptible individuals during the day t and Q i (t) the number of individuals infected by i at the end of day t. Clearly Z i (t) follows an hypergeomtric distribution with parameters n−X(t)−1, S(t) and N and, given Z i (t) = z, z ∈ {0, . . . , N }, it can be seen that Q i (t) ∼ Bin(z, q I ). Thus Furthermore, using the same ideas in Ferrante et al. (2016), the number of secondary cases corresponding to this individual will be Finally, the number of secondary cases produced by one arbitrary infectious individual (that can be in class I or A) at time t given S(t) and X(t) and that we will call R 0 (t), will be

Application to a COVID19 type disease
The COVID19 is a highly contagious disease that has appeared at the end of 2019. From all the information that is published every day, often contradictory, in the media, we can extract some properties of the disease. After the contagion, the virus remains in a latent state for 5-7 days before the individual became infectious. Then the individual can begin with symptoms in 3 days or he can continue asymptomatic, but probably infectious during two weeks. It is not known actually the number of asymptomatic people, but it will be probably bigger that the number of people with symptoms. About 75% of persons with symptoms have just light symptoms that last for two weeks and after that they became recovered. The other 25% get sever symptoms after a first period of 7-9 days with light symptoms. A 15% of these patients with severe symptoms die in 4-6 days, while the other 85% will recover after a period of 18-24 days. This disease can be described by the SEIAHRCD model defined before. According to the data above described, we can choose the deterministic values for the permanence in the classes r E = 5, r I = 3, r A = 14, r HR = 14, r HC = 9, r CR = 20, r CD = 4 and for the probabilities µ = 0.5, α = 0.25, and λ = 0.15.
We also have to make an assumption on the ineffectiveness of the individuals when in the classes I and A. We assume that q I = 0.2, and q A = 0.05.
That is, we are assuming that asymptomatic individuals are less infectious than individuals who have symptoms. R 0 , the expected number of secondary cases produced by an infectious individual, is equal with these parameters to We see that the threshold value R 0 = 1 is obtained for N = 1.54 and that R 0 = 3.5 (one of the possible empirical values estimated on the basis of real data) for N = 5.38. We study the spread of the contagion in a closed population of 10000 inhabitants for a period of time of 180 days. We assume that an infectious individual (in class I) arrives to this healthy population and he/she remains there for 3 days. We also assume that N i (t) ≡ N for any i and that therefore β i is constant for any individual. This assumption is strong in the case of a possible quarantine, but it is still reasonable.
In this paper, since the analytical approach to this model is really complicated, we focus on the evolution of the disease by implementing a simulation using environment Maple. All the values presented here are the mean computed over 30 repetitions of the simulation and we report also the confidence intervals.
We know very well that our model cannot explain exactly the COVID19 epidemic, since there are too many unknown aspects about this new disease, but we believe that the study of the behaviour of our model can help to understand the COVID19 epidemic. More precisely, we will answer five Research Questions regarding the dynamics of the disease depending on some of the parameters involved. Particularly, we deal with (1) the importance of the number of contacts, (2)-(3) the effectiveness of a quarantine depending on the moment it begins and on its duration, (4) the role of the asymptomatic depending of their level of infectiveness and (5) what happens with different levels of group immunization.
For each of these situations we study six quantities that we consider of major interest: 1. Class D: The number of deceased, therefore individuals in class D, after 180 days. For all these quantities we give a table with the mean and the 95% confidence interval for the mean in several situations. We also present a plot of one simulation of the number of deaths, the number of individuals in class C and the number of individuals that enter in class H each day for any of these situations along the 180 days. Note that the number of deceased is considered assuming that the health service is able to give the same level of assistance whatever the number of patients is, but probably in the situations where the health service is more stressed the assistance will be worse and the number of deceased may increase.

RQ1: How does the number of contacts N influence the spread of the disease?
We consider the dependence on the number of contacts on the evolution of the disease when N = 10, 5, 4 and 3.

RQ2: How does the beginning date of the quarantine influence the spread of the disease?
Let us check the effects of a quarantine on the evolution of the disease. We assume that N = 10 at the beginning of the spread of the disease and we consider three levels of quarantine, defined by the number of contacts N = 3, N = 2 or N = 0 (the total quarantine). We can also consider what happens depending on the moment that the quarantine starts: • when there are 10 deceased (first 10 individuals in class D), • when there is the first deceased (first individual in class D), • when there is the first individual with severe symptoms (first time the class C in not empty).

Quarantine beginning after first 10 deceased.
Let us recall that we assume that before the quarantine N = 10. We get It is clear that the efficacy of these quarantine, including the total quarantine, is almost null, since at this level of the disease the 99% of the population has been infected (case N = 0).   In this case we can note the efficacy of a quarantine if it is rigorous, with a significant difference between N = 3 and N = 2. Note that in the total quarantine, there will be more that 100 deceased since at this level of the disease more of the 50% of the population has been infected (case N = 0).

Quarantine beginning after the first individual enters in class C
In this case the quarantine is effective. For N = 2 we reduce the mortality of about 50% and for N = 0 the number of deceased is just 20. Note that when the first individual enters in C there are about 1000 individuals infected.  3.3 RQ3: How does the duration of the quarantine influence the spread of the disease?
Let us consider now the duration of the quarantine. We consider quarantines, beginning at first individual in C, of duration 20, 60 and 120 days. We also deal with two levels of quarantine, beginning from a number of contacts of N = 10 we pass to a quarantine with N = 3 or N = 1.

3.3.1
The case with N = 3  When the quarantine is done with N = 3, its effects in the number of deceased are not very important although we do a long quarantine of 120 days. Nevertheless, we are able to reduce the maximum in class C. To value the importance of the quarantine we need a more restricted quarantine.

3.3.2
The case with N = 1  We can see here that with a strong quarantine (N = 1) and with a long quarantine (k = 120), the disease can be stopped. In the two other cases, k = 20 or k = 60 the final numbers of deceased are very similar, since after a short stop the disease returns to growth up arriving at the same levels of the initial outbreak and having two clear peaks.

RQ4: How does the infectivity of asymptomatic individuals influence the spread of the disease?
One of the main problems in the study of the COVID19 disease is the role of the asymptomatic. In our general framework we have supposed that the infectivity of an asymptomatic individual (class A) is one fourth of the infectivity of an infective individual that will have symptoms (class I), that is, q A = 0.05 and q I = 0.2. Here we consider two other cases: in the first both probabilities are equal (q A = q I = 0.2) and in the second case the asymptomatic individuals are not infectious (q A = 0). We also consider the case without quarantine and the case with quarantine (from N = 10 to N = 3) after the first individual in class C.  We can see that the behaviour of the disease is very similar in the three cases and the main difference just consists in the velocity of its spread.

The asymptomatics' role with quarantine
We can see the importance of knowing the role of asymptomatics in implementing quarantine. The more infectious the asymptomatic are, the more effective the quarantine is.
3.5 RQ5: How does the group immunity influence the spread of the disease?
We consider finally the group immunity. We consider three cases, where the immunity group is of 50%, 70% and 90%, respectively. We can see that an immunity of the 50% is not enough to control the disease since almost all the individuals without immunity become infected after the 180 days.  On the other hand, with an immunity of the 70% the disease is well controlled. Moreover if the immunity regards the 70% of the population, in 11 over 30 trials there aren't any deaths and in one of the trials the disease does not infect any person. If we consider an immunity of the 90% only in 8 of the 30 trials there have been infected people with a maximum of 12 individuals infected and no one deceased.

Conclusions
Like any mathematical model, the model presented in this paper does not exactly describe COVID19 disease, even if it shares with this most of its characteristic. Furthermore, for this new disease many aspects are still unknown, such as the level of infectivity of asymptomatic patients.
Our goal is to understand how a disease similar to the COVID19 spreads over in a closed population and to answer to some specific research questions. What we have obtained can be summarized as follows: 1. Reducing the number of contacts of each individual the spread of the disease slows down, the pressure on the health system reduces, but we end up with a similar number of deceased.
2. It is basic to start the quarantine at least when the first severe patient is detected. Waiting for the first deceased leads to many more additional deceased.
3. To raise the quarantine it is very important to know the level of infectivity of the asymptomatic. The more infectious they are, the more important is the quarantine.
4. To stop the disease you must perform a strict and long quarantine and you must start this as soon as possible.
5. Group immunity is very important to prevent the development of the disease.
Even if most of the answers are the expected ones, we have obtained these through a sound stochastic epidemic model, that despite of its simplicity probably is able to catch most of the peculiarity of this disease.
We believe that a more sophisticated version of this model and more elaborated simulations can allow us to answer to more complex questions, but probably it will be better to wait when a deeper knowledge on this disease will be available.