Message Boards

GROUPS:

Staff Picks Biological Sciences Data Science Mathematics Medical Sciences Social Science Graphics and Visualization Wolfram Language Modeling Statistics and Probability

Arnoud Buzing

[NB] The SIR Model for Spread of Disease

Arnoud Buzing, Wolfram Research

Posted 2 years ago

19487 Views

21 Replies

14 Total Likes

Follow this post

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus Introduction The SIR (Susceptible-Infected-Recovered) model for the spread of infectious diseases is a very simple model of three linear differential equations. It examines how an infected population spreads a disease to a susceptible population which transforms into a recovered population. This model assumes no one dies from the disease and that the population is also isolated and stable (no new people arrive by birth or migration during the model simulation). The independent & dependent variables In this model, the independent variable is t (time, in days): In[]:= Clear[t] The dependent variables are S (number of susceptible individuals, at time t), I (number of infected individuals, at time t), and R (number of recovered individuals, at time t). To simplify the differential equations we will use the scaled versions of the S, I, and R variables. The scaled variables are s, i, and r. The scale runs from 0 (representing a zero population) to 1 (representing the full population). In[]:= Clear[s,i,r] The susceptible population starts near 1, with almost everyone being susceptible. It drops over time. The infected population starts with a positive small value near 0. The recovered population starts at 0, with no one having recovered. It increases over time. At any given time s+i+r=1 , because the sum represents the full scaled population which is 1. Each infected individual is assumed to spread its disease to b new individuals each day. If b is high the disease will spread faster. In[]:= Clear[b] The fraction k represents the fraction of the infected population that recovers each day. For example, if the disease lasts 7 days, then k is 1/7. In[]:= Clear[k] The differential equations The first equation describes how the susceptible population changes over time. Because we start with a susceptible population of 1 which decreases over time, we can assume the rate of change will be negative. Also the rate of change is proportional to the number of contacts b (more contacts decreases the susceptible population faster). Additionally the rate of change is proportionally related to the current susceptible and infected population (but not the recovered population): equationS=s'[t]-bs[t]i[t] Out[]= ′ s [t]-bi[t]s[t] The second equation describes how the infected population changes over time. It has two terms. The first term is the same as the right hand side of the previous equation, except the sign is now positive. It represents the increase in infected population. The second term is negative because it represents the decrease of infected cases (people who recover): equationI=i'[t]bs[t]i[t]-ki[t] Out[]= ′ i [t]-ki[t]+bi[t]s[t] The final equation describes how the recovered population changes over time. The right hand side is the same as the last term of the previous equation. It represents the transition of infected cases to recovered cases: equationR=r'[t]ki[t] Out[]= ′ r [t]ki[t] Before solving the differential equations we need to set the parameters b and k. Let' s assume that, on average, each infected person spreads their disease to one person every two days: b=0.5 Out[]= 0.5 Here we assume that the average recovery time is 14 days, which means that, on average, 1/14th of the infected population recovers: k=1/14.0 Out[]= 0.0714286 We also need initial values at t=0 for the susceptible, infected, and recovered population. We need a small non-zero value for the initial infected population, so will will set it to 0.0001. We will set the initial value for the susceptible population to 0.9999 (everyone else, who is not infected). And we will set the initial recovered population to 0. We pass the differential equations and the initial values to NDSolve and request a solution for t=100 days: In[]:= solution=NDSolve[{equationS,equationI,equationR,s[0]0.9999,i[0]0.0001,r[0]0.0000},{s,r,i},{t,100}]; Next we get the individual solutions for S, I, and R: In[]:= solutionS=First[s/.solution];solutionI=First[i/.solution];solutionR=First[r/.solution]; We can plot the results. Here, the blue line represents the susceptible population over time, red is infected, and green is recovered: Plot[{solutionS[t],solutionI[t],solutionR[t]},{t,0,100},PlotRange{0,1.01},PlotStyle{Blue,Red,Green}] Out[]= Alternatively we can stack these plots to highlight their sum is always 1: Plot[{solutionS[t],solutionS[t]+solutionI[t],solutionS[t]+solutionI[t]+solutionR[t]},{t,0,100},PlotRange{0,1.01},Filling{1Axis,2{1},3{2}},PlotStyle{Blue,Red,Green}] Out[]=

POSTED BY: Arnoud Buzing

21 Replies

Sort By:

Luc Hens

Luc Hens, Vrije Universiteit Brussel

Posted 2 years ago

Why, if the average recovery time is 14 days, does this imply that on average one third of the infected population recovers?

POSTED BY: Luc Hens

Arnoud Buzing

Arnoud Buzing, Wolfram Research

Posted 2 years ago

Sorry, that was a typo. I will see if I can fix the notebook text.

POSTED BY: Arnoud Buzing

Arnoud Buzing

Arnoud Buzing, Wolfram Research

Posted 2 years ago

I think I fixed it.

POSTED BY: Arnoud Buzing

Robert Nachbar

Robert Nachbar, Wolfram Research Inc.

Posted 2 years ago

The statement is still incorrect. In the model with these parameter values, all the susceptible individuals get infected and all recover. If the recovery rate is fast enough, R0 = b N / k = b / k will be less than 1 and the infection will die out before everyone gets infected. The exact value can be worked out from the equations in Brauer.

POSTED BY: Robert Nachbar

Arnoud Buzing

Arnoud Buzing, Wolfram Research

Posted 2 years ago

Sure, but why are the values I am using incorrect? I am just using example values (which anyone can vary to experiment with this model?)

POSTED BY: Arnoud Buzing

Robert Nachbar

Robert Nachbar, Wolfram Research Inc.

Posted 2 years ago

I see the difficulty now. I was interpreting the statement "Here we assume that the average recovery time is 14 days, which means that, on average, 1/14th of the infected population recovers” too broadly, that is, I was assuming you meant the total number of infected and not the current number of infected. What you have is indeed correct. Sorry for the confusion.

POSTED BY: Robert Nachbar

Luc Hens

Luc Hens, Vrije Universiteit Brussel

Posted 2 years ago

The sentence would be clearer if it would say: "... 1/14th of the infected population recovers every day:"

POSTED BY: Luc Hens

Robert Nachbar

Robert Nachbar, Wolfram Research Inc.

Posted 2 years ago

No, this is an instantaneous rate, not a daily rate.

POSTED BY: Robert Nachbar

Daniel Lichtblau

Daniel Lichtblau, Wolfram Research

Posted 2 years ago

It has to have units of recoveries/time, and by design the time units are days.

POSTED BY: Daniel Lichtblau

Robert Nachbar

Robert Nachbar, Wolfram Research Inc.

Posted 2 years ago

Yes, the units are day^-1, but the fraction of infected that recover in a given day is the integral Integrate[k i[t], {t, t1, t2}] / i[t1] where t2 - t1 = 1 and i[t] is the solution from NDSolve, is not necessarily exactly 1/k.

POSTED BY: Robert Nachbar

Christos Papahristodoulou

Christos Papahristodoulou, Mälardalen University

Posted 2 years ago

What about if the Recovery equation changes to something like: equationR = r'[t] == (1 - w) k i[t], where w = a "death" parameter, such as w = 0.005? Then, DSolve[equationR, {r[t]}, {t}] gives: r[t] -> C[1] + Inactive[Integrate][(k i[K[1]] - k w i[K[1]]), {K[1], 1, t}]

POSTED BY: Christos Papahristodoulou

Arnoud Buzing

Arnoud Buzing, Wolfram Research

Posted 2 years ago

Yes, I think this basic set of differential equations can be extended in many ways to make it more useful and precise. It would be interesting to see what you can do with it!

POSTED BY: Arnoud Buzing

Moderation Team

Moderation Team, WOLFRAM

Posted 2 years ago

-- you have earned *Featured Contributor Badge* Your exceptional post has been selected for our editorial column *Staff Picks* http://wolfr.am/StaffPicks and Your Profile is now distinguished by a *Featured Contributor Badge* and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: Moderation Team

Updating Name

Posted 2 years ago

How do I calculate b from real data? Here are the cases the begin of infecteds in São Paulo Brazilian State data = {1,1,1,2,2,2,2,3,6,10,13,16,19,30,42,56,65,136,152,164,240,286,396,459,631,745,810}; Clear[a, k]; NonlinearModelFit[data, a Exp[k t], {a, k}, t]["BestFitParameters"] {a -> 2.35641, k -> 0.219352} It's somehow related to k?

POSTED BY: Updating Name

Arnoud Buzing

Arnoud Buzing, Wolfram Research

Posted 2 years ago

I am not sure how to best calculate b and k from actual data. This model is very simple and the real world is quite complex (different spreading rates for different people and places, and time-dependent spreading rates due to changing policies). It may be possible to get a fit for a current situation, but this fit probably has limited predictive value. This model does have value to understand how changing the parameters affects the spreading of infectious diseases.

POSTED BY: Arnoud Buzing

David G

Posted 2 years ago

I was wondering if you could change the WindowSize option to Automatic for this notebook. I downloaded it from WolframCloud and it displayed in Full Screen mode, making it impossible (AFAIK) to resize it. I was finally about to open options dialog and reset the WindowSize option. Would others have the same issue?

POSTED BY: David G

Y Getachew

Posted 2 years ago

Can this be replicated with discrete time?

POSTED BY: Y Getachew

Arnoud Buzing

Arnoud Buzing, Wolfram Research

Posted 2 years ago

I think so. But I have not tried it.

POSTED BY: Arnoud Buzing

Y Getachew

Posted 2 years ago

II tried to replicate it in discrete time but RSolve/RSikvevakye couldn't solve it. Attachments: SIR_Dicrete_time.nb

POSTED BY: Y Getachew

Robert Nachbar

Robert Nachbar, Wolfram Research Inc.

Posted 2 years ago

I've just finished making a stochastic solver for compartmental models, and will be making a post to Community later this week.

POSTED BY: Robert Nachbar

Daniel Lichtblau

Daniel Lichtblau, Wolfram Research

Posted 2 years ago

I would suggest either `RecurrenceTable` or an explicit loop for obtaining a numerical solution in discrete time steps. The system might not be amenable to exact analytic results such as would be provided by `RSolve`.

POSTED BY: Daniel Lichtblau

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback

[NB] The SIR Model for Spread of Disease

Introduction

The independent & dependent variables

The differential equations