Epidemic Modelling: Simulations using Stochastic Methods

By Angeline Xiao

Branching Processes are a method of modelling the processes of populations that evolve independently with chance. These stochastic models can be used to simulate epidemics. Say we have an infected individual, they will inflict more cases with a certain probability, and the new infected individuals will independently continue to infect at the same rate. [1] The offspring distribution, and mean daily offspring are very useful in showing how epidemics may spread stochastically, given different parameters and models. Over summer, as part of the UoA summer research programme, we conducted epidemic simulations for a variety of models and infection rates.

We can determine if a branching process will die or not based on its mean daily offspring, popularly known as an Reff value. If each individual has on average less than 1 offspring per generation then it will almost surely die. We call this process the subcritical process. If the mean offspring is 1, this is the critical process, and the population will still, almost surely die. If Reff is greater than 1 then we call this a supercritical process, and it is not certain that the process will die out, although it is still possible. [2]

We can have branching processes in both discrete and continuous time. Both have their uses, with discrete time being a simpler model to use and understand. Continuous time branching processes are still useful as viruses do not adhere to human set intervals such as days, and a more accurate picture can be painted.

For a discrete time branching process model, the offspring of a population can be expressed by the sum of all the offspring that each member of the population will have, independent of each other. For a population in generation n, Zn, the population in the next time generation, Zn+1, can be expressed by the sum of independent identically distributed random variables which is same offspring distribution, Xn,i , where Xn,i represents the offspring per individual per time generation.

That is,

To simulate total active cases, we define that when an individual has n infected offspring in a new day, n+1 of those offspring are new infections, and 1 offspring is the individual surviving. If an individual has 0 offspring then they have 'recovered'.

Perhaps the most intuitive offspring distribution is the Poisson distribution. Every day each infected individual infects more individuals at a rate λ, where λ-1 is the mean number of individuals they infect (as they are not infecting themselves again). 100 iterations of this simulation were performed with λ rates of 0.8, 1 & 1.2 to demonstrate a subcritical, critical & supercritical process. We can observe the trend towards death of the subcritical process population while the supercritical process population explodes (Fig. 1).

Figure 1

An extremely important but less intuitive offspring distribution we can use is the Geometric distribution. We define a Geometric distribution to represent the probability of the number of failures before the first success in a sequence of Bernoulli trials, where we can set the parameter p to be the probability of success. In context, this would signify the number of individuals infected before an individual stops infecting for the time unit.

For X ~ Geom(p),

The Geometric distribution has special properties that allow us to explicitly calculate the branching process with a Geometric offspring distribution using generating functions. This allows us to conduct sanity checks on our simulations before we dive into more complex, non-calculable simulations. To compare subcritical, critical, and supercritical processes, p values of 5/9, 1/2 & 5/11 were used so that we have a mean daily infection rate of 0.8, 1, and 1.2 respectively. We can see that the simulation run for our 3 mean daily offspring rates looks similar to the Poisson model (Fig. 2).

Figure 2

However what we are really interested in is the extinction proportion. On an arbitrary day (take day 10) we can see the distribution of extinction for each other mean daily offspring rates. (Fig. 3, Fig. 4, and Fig. 5).

Figure 3

Figure 4

Figure 5

By comparing the histogram on day 10 to the explicit survival probability on day 10, we can confirm the accuracy of our simulations.

Using the generating function in the geometric case, we can find the extinction probability on day n,

, for an initial population of 1 explicitly. [3] Given the offspring distribution:

Let us assume p does not equal q with μ = p/q for n > 1

Notice

Hence for:

We can apply this formula to find πn for our simulations.

For p of 5/9, π10 = 0.9765

For p of 1/2, π10 = 10/11

For p of 5/11, π10 = 0.8074

To properly compare with the simulations run, these results need to be adjusted for 10 initial population, which is equivalent to having 10 independent branching processes with population 1 all dying out.

Therefore with an initial population of 10:

For p of 5/9, π10 = 0.9765^10 = 0.7884 (simulated: 0.789).

For p of1/2, π10 = (10/11)^10 = 0.3855(simulated: 0.365).

For p of 5/11, π10 = 0.8074^10 = 0.1177 (simulated: 0.111).

We can see that the theoretical values line up with our simulated extinction proportions.

From these basic models we are able to adjust and add variables, such as shortening the generation time, or adding more complex offspring distributions (such as piecewise distributions) to accurately model epidemics.

As previously mentioned, branching processes can also operate in continuous time, where after a certain time calculated by an exponential distribution, an individual infects a number of individuals represented by an offspring distribution. Specifically, in a birth death process, an individual infects a new individual at rate α recovers at rate β independently. These birth and death times are exponentially distributed. [4]

Let the population at time t be Zt.

Z0 = 1

If α/β < 1 then the population will die out, and for α/β > 1,

, and the population will not necessarily die out.

Two independent exponential distributions with α = population • infection rate, and β = population • recovery rate, were simulated. The smaller of the two results is taken as the time before the first event, with the events either being infection or recovery. A critical process occurs with infection rate = recovery rate. A one-type continuous branching process at critical was simulated below (Fig. 6). We can introduce different strains or types of viruses in the same model in what is called a multi-type branching process. (Note: This is possible for both a continuous and discrete time model, we are simply choosing to demonstrate it with a continuous model) A multi-type model of a branching process is when there are more than one type of individual in the population. For the different types of virus that exist, they each have independent infection and recovery rate, as well as a rate to change from one type to another.

Figure 6

A specific example could be when a virus has an incubation period and an infectious period like Covid -19 did. When people are first infected, they are in the incubation period and are not infectious. After an exponentially distributed amount of time, they become infectious, and then will recover after that following an exponential distribution. The graph below shows a simulation of a two type branching process with an incubation and infectious period, with a critical infection and recovery rate (Fig. 7).

Figure 7

The simulations shown for branching processes in continuous time are all birth-death processes, which is a subset of branching processes. This simulation can also be run for different offspring distribution, such as a geometric offspring distribution, and for multi-types with more than two types.

Branching process models can be adapted and used to model epidemics, with a variety of offspring distributions, and using both continuous and discrete time. This report showed the different ways that stochastic models can be used to model epidemics. Simulating stochastic models can be very useful for results that cannot be explicitly calculated. The validity of the results in this report can be confirmed by comparing them to the explicit calculations. This research can serve as a base for more sophisticated stochastic population models combining more types, time inhomogeneity, as well as population size dependence branching.

This project was conducted under the UoA Summer Research Scholarship. I would like to thank the UoA Department of Statistics for making this possible, and especially to my supervisor Simon Harris, who has very patiently guided me through this project.

References

1. Theodore Edward Harris et al. The theory of branching processes, volume 6. Springer Berlin, 1963.

2. Geoffrey Grimmett and David Stirzaker. Probability and random processes. Oxford university press, 2020.

3. S imon Harris. Stats 225 Coursebook. University of Auckland, 2020.

4. Rinaldo B Schinazi. Continuous time branching processes. In Classical and Spatial Stochastic Processes, pages 151–173. Springer, 2014.