The goal of this project is to develop a Monte Carlo simulation of the spread of an infectious disease. We will use the classical SIRS model of to develop the necessary transition probabilities for the simulation. Interpreting the simulation requires a thorough understanding of the SIRS model, so a deterministic approach is used in conjunction with Monte Carlo methods.
The main purpose of creating such a simulation is to investigate how a disease spreads throughout a given population over time. Then we can make predictions about whether or not a certain disease has the capacity to establish itself within the population, i.e. a fraction of the population remains infected after the system reaches equilibrium.
The SIRS model considers an isolated population of \( N \) people which are divided into three separate groups:
For the simple case, we will assume that the dynamics of the epidemic occur during a time scale much smaller than the average person's lifetime. Hence the effect of the birth and death rate of the population is ignored.
From these assumptions, a set of coupled differential equations can be constructed to form the classical SIRS model: $$ \begin{equation} \begin{split} S'&=cR-\frac{aSI}{N}\\ I'&=\frac{aSI}{N}-bI\\ R'&=bI-cR\\ \end{split} \label{_auto2} \end{equation} $$
Note that for small populations, one may choose to write \( aSI \) instead of \( \frac{aSI}{N} \) since the number of susceptibles which become infected depend more on the absolute number of infected people rather than the infected fraction of the population.
Though this set does not have analytic solutions like the closely-related SIR model, the equilibrium solutions are simple to obtain. The constraint in (1) reduces this three dimensional system into a two dimensional one so that the equation for \( R' \) can just be omitted: $$ \begin{equation} \begin{split} S'&=c(N-S-I)-\frac{aSI}{N}\\ I'&=\frac{aSI}{N}-bI\\ \end{split} \label{_auto3} \end{equation} $$
The steady state is found by setting both equations in (3) equal to zero. Let \( s \), \( i \), and \( r \) denote the fractions of people in \( S \), \( I \), and \( R \), respectively. Then the fractions of people in each group at equilibrium are: $$ \begin{equation} \begin{split} s^* &= \frac{b}{a},\\ i^* &= \frac{1-\frac{b}{a}}{1+\frac{b}{c}},\\ r^* &= \frac{b}{c}\frac{1-\frac{b}{a}}{1+\frac{b}{c}}. \end{split} \label{_auto4} \end{equation} $$
Here, the asterisk (\( ^* \)) signifies that these fractions are at equilibrium. Notice that each fraction must be a number between 0 and 1, and that the three fractions must add up to 1. Thus the equations in (4) suggest that the rate of recovery \( b \) must be less than the rate of transmission \( a \) for the number of infected people at equilibrium to be greater than zero. In other words, the disease establishes itself in the population only if \( b < a \).
For this project you can collaborate with fellow students and you can hand in a common report. This project (together with projects 3 and 4) counts 1/3 of the final mark.
We will create four different populations to investigate the effect of increasing the rate of recovery \( b \) and crossing the "threshold". Each population consists of 400 people, 100 of which were initially infected and 300 of which were initially susceptible. We fix the rate of transmission \( a \) and the rate of immunity loss \( c \) between populations, but let the rate of recovery \( b \) be varied because it is the one parameter which can be reasonably controlled by the actions of human society (access to health care, development of new medicines, etc). Note that the rates have units of inverse time. The table here lists a possible set of parameters for the four different populations.
| Rate | A | B | C | D |
|---|---|---|---|---|
| a | 4 | 4 | 4 | 4 |
| b | 1 | 2 | 3 | 4 |
| c | 0.5 | 0.5 | 0.5 | 0.5 |
The rates of transmission are given by the parameter \( a \), recovery by \( b \), and immunity loss \( c \) for populations \( A \), \( B \), \( C \), and \( D \).
We treat the population as a continuous variable. Develop a Runge-Kutta fourth order code to solve the above equations.
Discuss your results and present your interpreations.
The set of differential equations in (2) inherently assume that \( S \), \( I \), and \( R \) are continuous variables, when they are discrete in reality. We can use the idea of randomness to remedy this issue by defining a set of transition probabilities for the possible moves a person can take from one state to another: \( S\rightarrow I \), \( I \rightarrow R \), and \( R \rightarrow S \). To obtain the specific values for each probability, we first notice from (2) that in a small time step \( \Delta t \), the number of people moving from \( S \) to \( I \) is approximately \( \frac{aSI}{N}\Delta t \). Likewise, about \( bI\Delta t \) move from \( I \) to \( R \) and \( cR\Delta t \) move from \( R \) to \( S \).
Now we let \( \Delta t \) be small enough so that at most one person moves from a given group to another. Since $$ \begin{equation} \begin{split} &\max \Big\{ \frac{aSI}{N}\Delta t \Big\} = \frac{a}{N}\left(\frac{N}{2}\right)^2\Delta t=\frac{aN}{4}\Delta t,\\ &\max \Big\{ bI\Delta t \Big\} = bN\Delta t,\\ &\max \Big\{ cR\Delta t \Big\} = cN\Delta t,\\ \end{split} \label{_auto5} \end{equation} $$ the time step is then given by $$ \begin{equation} \Delta t = \min \Big\{ \frac{4}{aN}, \frac{1}{bN}, \frac{1}{cN} \Big\}. \label{_auto6} \end{equation} $$ This construction allows us to reinterpret the values \( \frac{aSI}{N}\Delta t \), \( bI\Delta t \), and \( cR\Delta t \) as transition probabilties: $$ \begin{equation} \begin{split} P(S\rightarrow I) &= \frac{aSI}{N}\Delta t,\\ P(I \rightarrow R) &= bI\Delta t,\\ P(R \rightarrow S) &= cR\Delta t.\\ \end{split} \label{_auto7} \end{equation} $$ For each possible move, a random number between 0 and 1 is generated. If the number is less than the probability for the move, the move is taken.
Develop now a Monte Carlo algorithm which implements the last equations using the same parameters \( N \), \( a \), \( b \), and \( c \) for each population. Compare your results to those obtained with the deterministic differential equation solver and discuss your findings.
You will need to find the fraction of people when the system has been equilibrated. You will need to find the equlibrium expectation values and corresonding standard deviations.
The same principles used in this simple model can be extended to include more details about the population and disease.
For instance, vital dynamics can be easily added to the system so that the model can describe the spread of diseases which occur over longer stretches of time. If \( e \) is the birth rate, \( d \) is the death rate, and \( d_I \) is the death rate of infected people due to the disease, then the modified differential equations are given by: $$ \begin{equation} \begin{split} S'&=cR-\frac{aSI}{N}-dS+eN\\ I'&=\frac{aSI}{N}-bI-dI-d_I I\\ R'&=bI-cR-dR\\ \end{split} \label{_auto8} \end{equation} $$ Here, we have assumed that all the babies born into the population are initially susceptible.
Add now these additions to your ODE solver and Monte Carlo solver and discuss the results.
For diseases such as influenza, the rate of transmission depends largely on the time of year. During the colder months, individuals are more likely to spend time in closer proximity to one another, resulting in a rate of transmission which oscillates. So we can let \( a \) be given by $$ \begin{equation} a(t)=Acos(\omega t) + a_0, \label{_auto9} \end{equation} $$
where \( a_0 \) is the average transmission rate, \( A \) is the maximum deviation from \( a_0 \), and \( \omega \) is the frequency of oscillation. Study the the effect of seasonal variations as well, with both the ODE and Monte Carlo solver. Try different values of the above new parameters and discuss your results.
Diseases with available vaccinations allow people to move directly from \( S \) to \( R \), breaking the cyclic structure of the SIRS model. We assume that a susceptible individual's choice to become vaccinated does not depend on how many other susceptibles are vaccinated. We may, however, assume that the rate of vaccination \( f \) can depend on the time, since this rate may oscillate during the course of a year and/or increase as awareness and medical research increases. Then the system of differential equations become $$ \begin{equation} \begin{split} S'&=cR-\frac{aSI}{N}-f\\ I'&=\frac{aSI}{N}-bI\\ R'&=bI-cR+f\\ \end{split} \label{_auto10} \end{equation} $$
Add now the effect of vaccination via the parameter \( f \) and discuss again your results. Play around with different values of the parameter \( f \).
Here follows a brief recipe and recommendation on how to write a report for each project.
The preferred format for the report is a PDF file. You can also use DOC or postscript formats or as an ipython notebook file. As programming language we prefer that you choose between C/C++, Fortran2008 or Python. The following prescription should be followed when preparing the report: