MXB261 Modelling and Simulation Science - Assignment 2 - A Simulation Project - Modelling an Epidemic
MXB261 - Modelling and Simulation Science Assignment 2 - A Simulation Project
Modelling an Epidemic
Epidemics are typically modelled mathematically using so-called compartmental models where a population of size P ∈ N (where N is the set of non-negative integers) is split into sub-populations called compartments that represent groupings of individuals by their state is respect to the given disease. In this project, you will explore the so called susceptible-exposed-infectious-removed (SEIR) epidemic model. Here individuals that are susceptible (S) to infection (i.e., they have no immunity to the disease) can become exposed (E) to infection via interaction with an infectious (I) individual. All exposed individuals eventually become infectious and eventually recover (R) with immunity (i.e., they can’t get infected again). A schematic of the valid transitions of an individual through these different states is given in Figure 1.
Figure 1: Schematic of state transitions and rates for the SEIR epidemic model
Compartmental models can be thought of in exactly the same way as a chemical reaction network, in this sense, the “reactions” are
S+I→α E+I, E→β I, I→ρ R, (1)
where S is the susceptible sub-population, E is the exposed sub-population, I is the infectious sub-population, and R is the recovered sub-population. Note that we assume there is no popu- lation growth during the epidemic so we have the conservation law P = S + E + I + R. The rate parameters are the transmission rate α > 0, the rate of symptom onset β > 0 (note: that 1/β is referred to as the incubation period), and the recovery rate ρ > 0.
Please note that the units of the rate parameters should be carefully considered, in this context you should treat the total population P as the volume of the space for possible interactions (see week 8 Lecture notes on propensity functions).
Task 1: Deterministic modelling of an epidemics
In this task, you will derive a deterministic model of the SEIR system using and ODE, you will then analyse its steady state and perform some simulations using Euler’s method and ode45.
-
(a) Relate the set of “reactions” in Equation (1) to a system of four ODEs describing the change in S(t), E(t), I(t) and R(t) over time t > 0 .
-
(b) Perform an analysis of the equilibrium solutions of the ODE you derived in part (a). That is find the steady state solutions and analyse their stability (i.e., stable or unstable).
-
(c) Derive the explicit Euler’s method of the SEIR ODE system you derived in part (a). Implement the scheme in MATLAB and simulate the system using parameter values α = 1/5, β = 1/14, and ρ = 1/10, initial conditions S(0) = 990, E(0) = 0, I(0) = 10 and R(0) = 0. The time span for the simulation is [0, T ] where T = 360 days, and the step size is h = 0.5 days. Plot the solution and compare with the relevant predicted equilibrium point you identified from part (b).
3
-
(d) Now let α = 10 and keep everything else the same as part (c). Simulate this system using Euler’s method plot against the result of MATLAB’s built in ode45. Compare the two solutions with some form of accuracy measure.
-
(e) Now simulate your model (either using your Euler scheme or ode45) for a range of transmission rate parameters (keeping everything else the same as in part (c)),
α = 0.025, 0.05, 0.075, 0.1, 0.125, 0.15 . . . , 1.0. (2)
For each value of α record the value of S(T) (i.e., the number of individuals who never caught the disease after a year). Plot S(T) against α. Comment on the relationship between α and S(T).
Task 2: Stochastic modelling of an epidemic
In this task, you will explore the stochastic model associated with the SEIR system.
-
(a) ConsiderthenetworkinEquation(1)fromtheperspectiveofadiscrete-statecontinuous- time Markov process. Identify the key components of the stochastic process, including the state vector, propensity functions and state change vectors (i.e., the stoichiometric vectors).
-
(b) Use the process you have defined in part (a) to implement the Gillespie stochastic simulation algorithm (SSA) for the SEIR. Use the same set of parameters and initial conditions as Task 1 part (c) and simulate n = 3 realisations of the stochastic process. Plot and compare the results against the ODE solution in Task 1 part (c).
-
(c) The expected value of a component in the stochastic model at a given time, t, can be estimated though Monte Carlo integration. For example,
X∞ E[S(t)] =
s=0
1 Xn
sPr(S(t) = s | X(0) = x0) ≈ n S(i)(t), (3)
i=1
where S(1)(t), S(1)(t), . . . , S(n)(t) are independent realisations of the susceptible com- ponent at time t. Using the same parameters as in part (b), compute Monte Carlo estimates of E[S(t)], E[E(t)],E[I(t)], and E[R(t)] for t = 0, 0.5, 1.0, . . . , 360.0. Plot these expectations over time and compare with the ODE solution from task 1. Re- peat the process for n = 3, 6, 12, 25, 50, 100, what happens as n increases?
(d) Implement Latin Hypercube Sampling on the 2D space of parameters α, ρ (each in the range [0, 1.0]). For an appropriate mesh size and number of trials, build a population of successful parameter pairs, that reflect the following two conditions (i) S(T) ∈ [S(0) − 20, S(10) − 10] or (ii) S(T ) < S(0) − 100. Plot the successful parameter combinations on as scatter plot colour coded by the condition it satisfied. What is distinctive about the two groups of parameters? Give an interpretation of from the perspective of the epidemic dynamics (you may wish to visualise some realisations of the model under these parameterisations).
Note: For this part You are to write your own code to carry out the Latin Hyper- cube Sampling. You can check that your code is working by comparing results with the built-in MATLAB command lhsdesign, however this comparison should not be included as part of your report.
4
Task 3: Latin Hypercube Sampling and Orthogonal Sampling in 3D
This task requires a Latin Hypercube Sampling to be implemented as for Task 2 part (d), but over a 3D parameter space α,β,ρ (each in the range [0,1.0]). Again, the parameter space is explored to determine regions that correspond to the same conditions as in task 2 part (d).
In this task you must also use a copy of your 3D Latin Hypercube Sampling to then develop an Orthogonal Sampler. Repeat your analysis over 3D parameter space and compare your results for Latin Hypercube Sampling and Orthogonal Sampling.
Visualisation of the successful parameter 3-tuples is required, with appropriate labelling. Discuss the approaches you take to represent the 3rd dimension in your plots in an effective and meaningful way. Write a concise paragraph to explain your results, and compare the two sampling methods. Make sure to then tie in the results of the sampler with those of Task 2 part (d) and Task 3 and discuss what this means for the epidemic dynamics.
Task 4: Spatial agent-based implementation
In this task, we’ll be looking at a spatial stochastic implementation of the SEIR model. A spatial implementation allows exploration of the system dynamics when interactions take place at an individual level rather than at population level. Such an approach also allows us to relax some assumptions (e.g., such as the “well-mixed assumption”)
Construct a grid of cells (101 × 101). This grid represents a city with a “northside” and a “southside” with the grid row of cells at (51,1) to (51,101) representing a river that divides the northside and southside. The river and grid boundary act like reflective boundaries (that is an agent cannot cross them), however, the city has B bridges that are to be placed at random on the river cells (these bridges are treated like any other non-river cell). Populate the grid with S = 990 susceptible agents and I = 10 infectious agents. The susceptible and individuals will be positioned uniformly at random across all city cells, however, the initial infectious agents are to be entirely placed on the northside of the river. There is to be one agent only per cell.
The rules for the simulation are as follows:
-
(a) Each agent, in turn, will attempt to move to a neighbouring cell that is N, S, E or W of its current position;
-
if the agent attempts to move through the boundary, it skips the attempted move and stays in its current location;
-
if the agent attempts to move through the river , it skips the attempted move and stays in its current location;
-
if the new cell is empty, the agent moves;
-
if the new cell is already occupied by an agent, the move does not take place.
-
-
(b) After all moves have occurred, each infectious agent will randomly select a neigh- bouring cell (in N, S, E and W) of its current position, if the selected cell contains a susceptible agent, then the susceptible agent becomes an exposed agent if a uniform random sample (u ∼ U(0,1)) u ≤ α;
-
(c) Each exposed agent will become an infectious agent if a uniform random sample (u ∼ U(0,1)) u ≤ β.
5
(d) Each Infectious agent will become a recovered agent if a uniform random sample (u ∼ U(0,1)) u ≤ ρ.
For a selection of parameter values, and for various bridge numbers B = 10, 20, 30, . . . , 100, simulate the evolution of the system and hence describe the relationship between connec- tivity from the northside and southside of the city (i.e., the number of bridges) and the observed system dynamics. Consider two parameter regimes;
• α=1,β=1/14,ρ=1/20,
• α=0.5,β=1/14,ρ=1/10,
To demonstrate your simulation, track the abundance of susceptible, exposed, infectious and recovered agents after every step for a simulation that lasts 400 steps (or shorter if it reaches equilibrium). Create two figures (one for each parameter regime) that are 2 × 5 subplots where each subfigure is a different numbers of bridges (B = 10,20,30,...,100). These figures should be included in your report and referred to when discussing the effects of the connectivity across the river. Create a movie of the spatial stochastic model above, capturing the output at regular time points as the simulation evolves for each of the eight scenarios. Each movie should not be longer than forty seconds each (ten frames per second). The videos will be uploaded to github.
6
The Group Report
The purpose of the Group Report is to combine the results from the individual tasks, and to compare and contrast the results obtained from the different approaches and strategies.
The report (which should be approximately 10 pages) should have the following structure:
-
an Introduction, setting the scene for the content (max. 1 page);
-
a Methods (by Task) section, where you describe the approaches taken, for each Task; include any relevant mathematical details (e.g., analysis of equilibrium, algorithm pseudo- code etc) (2-3 pages).
-
a Results and Discussion section, where you combine, compare and contrast the results for parameter regions and the characteristics of the system dynamics; include here plots that justify your answers; also include results and figures from Task 4 (the spatial agent- based simulation) (5-6 pages);
-
a Conclusions section, where you summarise the main results in no more than a few sentences (max. 1 page).
Note: Page limits are a guide. We only ask that you roughly aim to follow these limits. The Oral Presentation
The purpose of the oral presentation will be to present your methods and key findings in an interesting and informative way to an audience of experts in your field (i.e. you can use technical jargon).
The presentation is intended to be a group effort and will go for a maximum of 6 minutes (aim for 5 mins, hard cut off at 6 mins) where each group member will present one slide each discussing one of the tasks, and then a few extra slides for the group to present a brief introduction and their final conclusions.
The second part of the presentation will then consist of 4-5 minutes of questions from the audi- ence.
7
Guide to the Marking Schedule
Marks Breakdown Individual component
Structure 2 Solution 6
Methods (by task) 2
Results 2
Github 3
Oral component 5
Subtotal 20
2 Code is well structured and well documented.
0 Code is not structured and/or not documented.
6 Solution is correct, all aspects of the task have been con-
sidered.
3 Parts of the solution may be incorrect, most aspects of
the task have been considered.
0 Little/no solution.
2 Methods for the task have been clearly and correctly
outlined for reproducibility.
1 Methods for the task have been outlined, but may not
be correct, clear or reproducible.
0 Little/no methods.
2 Relevant and concise figures that clearly demonstrate the
solution. Figures are clear, correctly and completely la-
belled.
0 Irrelevant, excessive or no figures that do not clearly
demonstrate the solution.
3 Code has been committed to github at least once with
commit comments that reflect the changes made to the
code.
0 No github commits or no reasonable commit comments
reflecting changes made.
5 Conveys a depth of understanding of the programming
and issues for the Individual Task.
2.5 Conveys some understanding of the programming and
issues for the Individual Task.
0 No oral discussion, or poor level of understanding con-
veyed.
8
Group component
Structure and presentation
Introduction
Discussion and comparison
Conclusions
1 1
Report has a clear and logical structure. Presen- tation is of a professional standard. No major spelling/grammatical errors.
1 1
0 Report does not have a clear or logical structure.
Presentation is not professional. Major spelling or
grammatical errors.
Introduction provides a short relevant background
for the exercises performed in the report. Motivation
for each task and an outline of the structure of the
report is provided.
0.5 Introduction missing key components, but outlines the main features of the report.
0 Little/no introduction.
4 4
Discussion shows concisely presented, but extensive, detailed and relevant insight. Results from individual tasks have effectively been compared. Plots, figures and tables are clear, concise and relevant.
2 2
Conclusions accurately and concisely summarise
2 Discussion shows some relevant insight. Results from individual tasks have been compared. Plots, figures and tables may not be clear, concise and relevant.
0 Little/no discussion or comparison.
main findings of the project.
1 Conclusions summarise some main findings of the
project.
0 Little/no conclusions.
Subtotal 8 Total 28
9
Reflection (individual) x4
Structure and presentation 1
Structure and presentation 3
Subtotal 3 x 4 Total 40
1 Reflection clearly articulates the progress made by the individual and group over the last week. The mark award to group members is clearly motivated.
0 Reflection is not submitted or the mark submitted for group members is not clearly motivated.
3 Average of marks awarded by group members for an individuals contribution
10
