[MUSIC]

>> Modelling the Ebola epidemic, we have to make some

decisions as to what it is we want to use and what it is we don't want to use.

This is dominated by the available data.

So the available data gives us something about how many cases occurred and

how many people were involved.

An appropriate choice for the modeling paradigm is called SIR models.

That simply stands for susceptible infected and removed.

It's a standard technique in disease modeling.

Now we go on to use Julia for a purpose.

The purpose here is to model the West African Ebola epidemic.

And by this we just mean a mathematical model that will predict for

us the numbers that are observed.

The particular kind of model we use will be SIR.

And it takes a little while to explain it and then we still we still also have to

talk little bit about how the Julia code works but then we can plug the results.

So SRI, the name of the model stands for

three kinds of people involved in the disease.

The first kind are the susceptibles.

These are people who might still get ill.

Initially there are S at node, that time equals zero,

but then S time goes on, this changes.

And so we write S(t) for

the number of people who are susceptible to the disease at time t.

At any given time as the disease progresses,

there maybe some infected people, which we symbolize with i.

And of course initially S(t) equal zero.

There are some people who are ill because otherwise the epidemic would not start.

And then finally people previously infected but

no longer infected called removed.

We assume that they remain non infectious from that time onward.

So this is not a disease like flu where people can get re-infected all the time.

Unfortunately, remove does not mean that they have recovered.

The illness isn't a joke.

So it's important to notice that S, I and R do not have to be whole numbers.

All we want is for S, I, and R to be close to what is observed.

And by close we just mean in some few percent.

If we can get a model that gets within a few percent of the values

that are observed we have done brilliantly well.

Epidemiology is not an exact science.

Getting the numbers within a few percent unfortunately is only possible with

hindsight usually.

We will do this with the SIR model for our able epidemic.

Okay so let's get the equations.

What we do is that we assume that time goes into discrete steps.

And the distance from one time to another,

the duration of time form one step to the next will always be dt.

That means that T sub I which is the time of the I steps is t sub zero.

Which is zero, plus I times dt, so

t sub I is just I times dt.

So code from t sub I to t sub I plus one,

we want to compute values for S of t I plus one.

I at D sub I plus one and RD sub I plus one.

To predict these we will only use the values at the previous time step so

that's s and t sub I and S t sub I and R t sub I plus one.

And they are what we use, so we need a formula to compute S and

t sub I plus one using these.

And similarly for I.

To compute I t sub I plus one using those three.

And so our model is a set of three equations, and they all have this form.

The new value is the old value, plus the gains, minus the losses.

So that's a very obvious form in which models are derived.

And all of the effort of modeling goes into working out how you might compute

at every time step gains and the losses.

And for susceptibles,

of course, the only thing that can happen is that a susceptible becomes infected.

The simplest SIR models use the law of mass action.

And what it says is that all susceptibles are exactly equivalent and

they all have an equal chance of meeting an infected.

Similarly all infecteds are exactly equal,

and they all have an exactly equal chance of meeting a susceptible.

And so the number of meetings between infected and

susceptibles is just proportional to the product S times I.

That's the goal of mass.

Now the number of infections that will result from meetings is not equal to SMI.

It's not equal equal to S times I.

It's just proportional to it.

We will say that this is a constant,

that the probability that it's acceptable is infected when

an infected person is exactly the same for all people at all times.

And so then, lambda SI is the rate at which susceptibles become infected.

This is a rate per day, because we are using days as our time variable.

And so the actual loss over a time step dt in terms of day.

So say dt is one tenth of a day.

Then actually only a tenth of lambda SI should get infected.

And so lambda SI times dt is the rate at which

the lambda SI times dt is the number of new infections that have generated.

During this time step dt.

The new value is the old value minus the losses.

Now we go on to infected people.

So infected people stop being infected by being removed or being recovered.

And again we say that all infected people have the same chance of recovery

at every day.

So some constant fraction is removed every day.

But this is the fraction gamma.

And so therefore gamma I is the number of infecteds that get removed per day.

And so gamma I dt, that's the loss of infecteds.

The kind of infecteds we already know, because every bit of

population that gets lost from susceptible gets added to infected.

And so here is the new infecteds, the previous level of infected

plus all the ones that are newly infected minus the ones that are now being removed.

Because it's an SIR model, once you're removed you stay there, so

all that happens is that there's a gain term.

And this gain term is of course the same as the last term from the infecteds, and

so that's our third equation.

And our SIR model is just three equations for S, I, and R.

So now we go on to writing the Julia code.

What I choose to do is to write a single function, which takes S, I and

R all at once, and gives back the new values of S, I and R all at once.

So these functions of course have to know the values of the parameters.

And the way we're going to pause these values to Julia is not to write them

as part of the argument of the function.

But just to declare them at the level where the function is created,

and then the function will see them on the inside of the function.

This is quite surprising.

Many experienced programmers would not expect to see this.

The way it works is the following.

So, if we assign to the variable b the value a, that is no, I'll be here.

Then the function f of x is created so when I run this, that executes.

Then this line executes and it creates a function x.

Now, this point the value of b.

Is actually known.

So when we run this, we get 56 as the outcome.

If we change that to six, of course that becomes 48.

But if we change that to four, then that becomes 24.

So makes our function simple to write.

We have a population vector with three elements.

One, two and three and they will be the susceptibles, the infected, and

the removeds.

And then they'll create the new s is the old susceptibles minus lambda

times the susceptibles times the infected times dt.

So we need a value of dt as well.

And the new I, is that?

Here I used gam rather than gamma because Julia has a built in function, gamma,

that I don't want to interfere with.

So I use different variable name and of the function name.

Eventually I have to pass it back.

Now this makes a row.

If you put commas in here this would have been a comma so

new is comma newI comma newR would be a column.

But by putting in spaces there, this is an array with one row and three parts.

Now we compiles and we set our parameter dt, which is the time state length.

And lambda and

gamma which are the epidemic parameters which specify some sort of input argument.

And recall that this is a way of using one assignment to set three different values.

So we can set S, I, and R.

S becomes 1000, I ten, R becomes 20.

Notice that only 1000 is a float, the other two are integers.

In fact, what will happen is that Julia will convert ten and 20 to integers.

It's not the best form of programming, but Julia forgives us that, and

that's part of why it's easy and quick to write you the app code.

If we want to optimize the running of the Julia code then,

we might find in the end that, this is one of the reasons why our code is running

slower than we might wanted, right.

So, what does it do?

It gives us an output of 975,

this is susceptible as minor lambda times infected.

So, the number ten that become 10,000

then gam is one tenth that means 1000.

Okay so the function works and so we can go back to the top and

now we discuss how we look through all the times.

We want that function to run from start time to the end time and

that can't be infinite so we define an end time t final.

We will use 610 day to fit data from Wikipedia.

So what that means is that there are n steps.

Tf divided by dt, but now you've gotta watch out.

Computer arithmetic is not infinitely precise.

So if we compute n times dt, so we've defined n this way.

And now we calculate Tf as n times dt,

this should be exactly true mathematically, algebraically.

But, in a computer we cannot be sure that this is the case.

So if we wanted to, we tested this is an equality might not work exactly.

So what I would like to do instead is to use time as a variable,

and to keep explicit track of this number I.

Then I know exactly how many steps I want to take, at the end of it,

the time I reach might not be exact the Tf that I started with.

But it will be extremely close to it.

Actually, the take home message is that you should not rely on exact equality,

even if you can reasonably expect it.

Computer code that relies on exact equality is much more prone to error,

and it really requires professional programmers to use.

We know exactly how many values we are going to need,

then we can create an array to hold all of them.

If we take n steps, we need n plus one, that is of course we need a value for

T equals zero and then the values n steps on from there.

So total number of n plus one values.

So here's the concept code, now this code doesn't actually do anything.

The reason it doesn't do anything is that I only have comment lines inside the loop.

This is just to say, okay final 610 days, so

the n steps is an integer, that's why I call the round function.

And round with type integer means that it's going to round it to integer types.

So that gives me an integer number of steps.

And in the results I will hold in a variable, three columns n plus one rows.

And also a vector of times, which will just be n plus one float values.

We need to specify our initial values and our initial result values.

Right.

Now, we are ready to run the model.

So I just cut and paste this code here, initializing everything.

Set all the values we need.

We need to set lambda, we need to set gam, we need to set dt, we need to set tfinal,

we need to set S not, we need to set I not, we need to set R not.

None of these can be taken for granted.

And then we initialize.

We calculate the number of steps.

We initialize a vect.

We initialize on the right, or results and for the time vector,

we initialize this array.

So we've created the array here and then we initialize its first row.

We're not going to compute this row.

We're also not going to compute this time.

Having done that, we can now just take n steps,

now if we want to take a step so we have one.

The next value is going to go to two.

And then if the step is two the next value goes to three and

in the end we end up with the n plus one steps that we put in our array.

That's it now we have all of our numbers.

Now it's just a question of plotting the results.

And again we use pyplot as the back end of plots.

And we just say timvec versus result values so as you remember.

Result values is an array with three columns.

So, this will plot, again,

timevec, one curve for each of the columns of result values.

And that's what it looks like.

And of course, we should put up labels.

For now, we know that blue is the susceptibles.

They crash down the red spikes very quickly, and

that comes down and the green goes out very quickly.

So certainly within a hundred days and

probably within 30 days, the epidemic has spiked and gone.

And everybody has been ill.

It was a day over here in which half of the 10,000 people were infected.

So that establishes the basics of an SIR model for

the Ebola virus epidemic in West Africa.

And we are now ready to go on in the next lecture.

To do some interrogation of the model.

SIR Models of Disease Dynamics

Julia Scientific Programming

Meet the Instructors