When we worked with the SIR models, we made the assumption that the population of S$,$ I,

and R individuals was always large enough to justify a continuous$-$time, deterministic model.

In practice, this is rarely the case at the beginning of an outbreak, when there are very few

infected individuals; some diseases will take off, whereas others will fizzle out if the I$($t$=0)$

don$$t manage to infect anyone before they recover. It this problem, we will approach the SIR

model through stochastic simulation. In this model,

i$.$ The total population size is Ntot, which does not change over the time of the simulation.

ii$.$ At any given moment there will be S susceptible, I infected, and R $=$ Ntot $$ S $$ I

recovered individuals, where S$,$ I, R are the number $($rather than the proportions$)$ of

individuals.

iii. The total rate for any transition $($S $->$ I or I $->$ R$)$ is $\backslash$lambda tot $=$ IS$($R$0/$Ntot$)+$ I. Hence,

waiting times can be drawn from an Exp$(\backslash$lambda tot$)$ distribution.

iv$.$ Given an event, the probability that one of the infectives recovers is I$/\backslash$lambda tot; otherwise, a

susceptible gets infected.

v$.$ Starting from initial populations I$($t$=0)=$ Iinit, S$($t$=0)=$ Ntot $$ I$($t$=0),$ R$($t$=0)=0,$

we repeat steps i$$iv until the system hits a fixed point.

Given this,

$($a$)$ State the criterion for ending the simulation. That is$,$ what is the fixed point?

$($b$)$ In step $($iii$),$ we asserted that the total rate for any transition is $\backslash$lambda tot $=$ IS$($R$0/$Ntot$)+$ I.

What is the reasoning for this?

$($c$)$ Implement the above using Gillespie$$s algorithm. Your function should take arguments

Iinit, Ntot and R$0.$ It should return an array of four columns: time t for the step, S$,$

I, and R $=$ Ntot $$ S $$ I. $($Note: R$0$ is the reproduction number for the disease, not

R$($t$=0)$; be sure that you don$$t accidentally conflate them in your code.$)$

$($d$)$ Run the stochastic SIR model five or more times with Ntot $=3000,$ I$($t$=0)=3,$ R$0=2\mathrm{.}2.$

Plot I vs t for each one $($or all of them on the same axis with different colours$).$ What

do you observe?

$($e$)$ Modify your function from part $($c$)$ to return only Ipeak, the maximum number infected

over the course of the whole simulation and tf $,$ the duration of the simulation $($i$.$e$.,$ time

until it hits a fixed point$).$ Note that you need not keep thw whole timeseries for this,

just Ipeak and tf $.$ Then, write another $$wrapper$$ function that will call that function

Nsim times. The output of your wrapper function should be a two$-$column array with

Nsim rows, where the first column is Ipeak and the second column is tf for each run.

$($f$)$ Using the above function, plot a histogram of tf when Nsim$=100.$ Comment on what you

observe; what does it imply for the distribution of outbreak durations?