FYS4150 - Project 5 Disease Modeling  - Solved

Theoretical background and description of the system 
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 NN people which are divided into three separate groups: 

Susceptible (S): those without immunity to the disease,
Infected (I): those who are currently infected with the disease,
Recovered (R): those who have been infected in the past and have developed an immunity to the disease.
A person can move from one group to another only in the cyclic order suggested by the name of the model S→I→R→SS→I→R→S. The rate of transmission aa, the rate of recovery bb and the rate of immunity loss cc help describe the flow of people moving between the three groups. The population is assumed to mix homogeneously and total population is assumed to remain constant so that 

N=S(t)+I(t)+R(t).(1)(1)N=S(t)+I(t)+R(t).

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: 

S′I′R′=cR−aSIN=aSIN−bI=bI−cR(2)(2)S′=cR−aSINI′=aSIN−bIR′=bI−cR

Note that for small populations, one may choose to write aSIaSI instead of aSINaSIN 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′R′ can just be omitted: 

S′I′=c(N−S−I)−aSIN=aSIN−bI(3)(3)S′=c(N−S−I)−aSINI′=aSIN−bI

The steady state is found by setting both equations in (3) equal to zero. Let ss, ii, and rr denote the fractions of people in SS, II, and RR, respectively. Then the fractions of people in each group at equilibrium are: 

s∗i∗r∗=ba,=1−ba1+bc,=bc1−ba1+bc.(4)(4)s∗=ba,i∗=1−ba1+bc,r∗=bc1−ba1+bc.

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 bb must be less than the rate of transmission aa 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<ab<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. 

Part a) setting up the differential equations 
We will create four different populations to investigate the effect of increasing the rate of recovery bb 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 aa and the rate of immunity loss cc between populations, but let the rate of recovery bb 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 aa, recovery by bb, and immunity loss cc for populations AA, BB, CC, and DD. 

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. 

Part b), moving to a Monte Carlo simulation 
The set of differential equations in (2) inherently assume that SS, II, and RR 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→IS→I, I→RI→R, and R→SR→S. To obtain the specific values for each probability, we first notice from (2) that in a small time step ΔtΔt, the number of people moving from SS to II is approximately aSINΔtaSINΔt. Likewise, about bIΔtbIΔt move from II to RR and cRΔtcRΔt move from RR to SS. 

Now we let ΔtΔt be small enough so that at most one person moves from a given group to another. Since 

max{aSINΔt}=aN(N2)2Δt=aN4Δt,max{bIΔt}=bNΔt,max{cRΔt}=cNΔt,(5)(5)max{aSINΔt}=aN(N2)2Δt=aN4Δt,max{bIΔt}=bNΔt,max{cRΔt}=cNΔt,

the time step is then given by 

Δt=min{4aN,1bN,1cN}.(6)(6)Δt=min{4aN,1bN,1cN}.

This construction allows us to reinterpret the values aSINΔtaSINΔt, bIΔtbIΔt, and cRΔtcRΔt as transition probabilties: 

P(S→I)P(I→R)P(R→S)=aSINΔt,=bIΔt,=cRΔt.(7)(7)P(S→I)=aSINΔt,P(I→R)=bIΔt,P(R→S)=cRΔt.

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 NN, aa, bb, and cc 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. 

Part c) Improvements, vital dynamics 
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 ee is the birth rate, dd is the death rate, and dIdI is the death rate of infected people due to the disease, then the modified differential equations are given by: 

S′I′R′=cR−aSIN−dS+eN=aSIN−bI−dI−dII=bI−cR−dR(8)(8)S′=cR−aSIN−dS+eNI′=aSIN−bI−dI−dIIR′=bI−cR−dR

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. 

Part d) Seasonal Variation 
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 aa be given by 

a(t)=Acos(ωt)+a0,(9)(9)a(t)=Acos(ωt)+a0,

where a0a0 is the average transmission rate, AA is the maximum deviation from a0a0, and ωω 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. 

Part e) Vaccination 
Diseases with available vaccinations allow people to move directly from SS to RR, 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 ff 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 

S′I′R′=cR−aSIN−f=aSIN−bI=bI−cR+f(10)(10)S′=cR−aSIN−fI′=aSIN−bIR′=bI−cR+f

Add now the effect of vaccination via the parameter ff and discuss again your results. Play around with different values of the parameter ff. 

Introduction to numerical projects 
Here follows a brief recipe and recommendation on how to write a report for each project. 

Give a short description of the nature of the problem and the eventual numerical methods you have used.