Message Boards

GROUPS:

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

Jon McLoone

[NB] Agent based epidemic simulation

Jon McLoone, Wolfram Research

Posted 2 years ago

11214 Views

5 Replies

13 Total Likes

Follow this post

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus

enter image description here

Inspired by the model at Washington Post the following computational essay creates and visualizes an agent-based simulation of epidemic spread in a small community. The rules that control the agent's behaviour can be easily modified or extended.

The basic assumptions are that agents move freely and without purpose within a square space. If a vulnerable agent is touched by an infected agent, they immediately become infectious for a fixed period of time and before becoming permanently immune to the disease. Agents do not interact with each other in any other way.

Agents

Each agent in the model represents a person with their own individual attributes. For the simplest version of this model we give these people the following attributes:

Position: Their location in 2D space, given as an {x,y} pair.
Direction: Their direction of travel in space, given as a number from 0 to 2π
Status: Their disease status, given as either “Immune”, “Vulnerable” or a number which indicates that they are infected. The value of the number describes how long they will remain infected.
Isolated: A number to indicate their isolation status. Positive number indicates the time remaining until they self-isolate. Negative indicates they are isolated. Isolation means that the person no longer moves around, though they can still be visited by other people who have not self-isolated.

We start by creating a function to generate an initial random population of n agents. The agents start at random positions within a square space with a random initial direction.

In[]:=

initialize[n_(*Sizeofpopulation*),diseased_,immune_,recoveryTime_,isolationRate_,delay_,space_]:= Map[Join[#,<| "Position"RandomReal[space,{2}], "Direction"RandomReal[2Pi], "Isolated"RandomChoice[{isolationRate,1-isolationRate}{delay,Infinity}]|> ]&, Join[ Table[<|"Status"recoveryTime|>,{diseased}], Table[<|"Status""Immune"|>,{immune}], Table[<|"Status""Vulnerable"|>,{n-diseased-immune}]] ];

For example, here is a population of 10 people with 1 diseased and 25% set to self-isolate in 5 days.

In[]:=

testPopulation=initialize[10(*populationsize*),1(*diseased*),0(*immune*),20.(*recoverytime*),0.25(*isolationrate*),5(*isolationdelay*),20(*spaceavailable*)]

Out[]=

{Status20.,Position{10.3648,18.5176},Direction1.56859,Isolated5,StatusVulnerable,Position{15.6157,19.0584},Direction0.412972,Isolated∞,StatusVulnerable,Position{12.22,3.20619},Direction5.67977,Isolated∞,StatusVulnerable,Position{1.33157,3.86196},Direction2.94968,Isolated∞,StatusVulnerable,Position{4.86843,11.209},Direction1.63861,Isolated∞,StatusVulnerable,Position{6.13558,17.4433},Direction6.17746,Isolated∞,StatusVulnerable,Position{17.8775,2.29511},Direction3.57282,Isolated∞,StatusVulnerable,Position{8.11507,19.5828},Direction4.36879,Isolated∞,StatusVulnerable,Position{6.82111,19.8707},Direction3.26842,Isolated∞,StatusVulnerable,Position{0.232099,11.1059},Direction2.73981,Isolated∞}

Updating the agents

At each time step we must update each agent’s properties according to some rules:

1) The agent must change position according to their Direction.
2) If the agent is infected, the time remaining for their infection (Status) must count down.
3) The agents Direction is updated. For this model we will assume small random changes of direction - a purposeless meander.
4) The agent’s Isolation status is updated. If the agent is waiting for self-isolation the timer must count down, but in this model we also assume that if they are already self-isolated there is a chance that they abandon self-isolation depending on a Patience parameter.(See below)

In[]:=

updateAgent[agent_Association,timeStep_Real,patience_Real]:= Joinagent, <|"Position"->If[agent["Isolated"]≤0,agent["Position"],agent["Position"]+timeStep{Cos[agent["Direction"]],Sin[agent["Direction"]]}], "Status"Which[NumberQ[agent["Status"]]&&agent["Status"]<=0,"Immune",NumberQ[agent["Status"]],agent["Status"]-timeStep,True,agent["Status"]], "Direction"agent["Direction"]+RandomReal[timeStep{-0.1,0.1}]; "Isolated"Ifagent["Isolated"]≤0&&RandomReal[]<1-

timeStep

patience

,Infinity,agent["Isolated"]-timeStep |>

In[]:=

updateAgent["Status"20,"Position"{10.3,11.2},"Direction"1.973`,"Isolated"10,0.1(*timeStep*),10.(*patience*)]

Out[]=

Status19.9,Position{10.2609,11.292},Direction1.973,Isolated9.9

We assume that the agents are in a constrained space, so after updating the position must check if they are leaving this space. If the agent is outside of the box, change their direction so that they return to it on the next update.

In[]:=

bounceAgent[agent_,space_]:= Which[ agent["Position"][[1]]>space,Append[agent,"Direction"RandomReal[{Pi/2,3Pi/2}]], agent["Position"][[1]]<0,Append[agent,"Direction"RandomReal[{-Pi/2,Pi/2}]], agent["Position"][[2]]>space,Append[agent,"Direction"RandomReal[{Pi,2Pi}]], agent["Position"][[2]]<0,Append[agent,"Direction"RandomReal[{0,Pi}]], True,agent];

For example:

In[]:=

bounceAgent[<|"Position"{20.01,0.5},"Direction"0|>,20]

Out[]=

Position{20.01,0.5},Direction1.86937

All these updates can be performed on each agent in isolation. The final update must consider the population as a whole. We assume that if a vulnerable agent touches an infected agent (defined as a EuclideanDistance <= 1) then the vulnerable agent is immediately infected for a fixed period.

In[]:=

infectPopulation[population_List,recoveryTime_]:= Block[infected, (*Forefficiencysplitthepopulationintothreegroups*) vulnerable=Select[population,#Status=="Vulnerable"&], immune=Select[population,#Status=="Immune"&], infected=Complement[population,immune,vulnerable]; Join[ immune, infected, Map[(*Foreachvulnerableperson*) Function[{v},(*passoninfection*) With[{pos=v["Position"]}, If[Min[EuclideanDistance[pos,#]&/@infected[[All,"Position"]](*Smallestdistancetoaninfectedagent*)]<1, Append[v,"Status"recoveryTime], v] ]], vulnerable]] ];

We now create a function to perform all of these updates for the entire population:

In[]:=

updatePopulation[population_,timeStep_,recoveryTime_,patience_,space_]:= infectPopulation[ bounceAgent[ updateAgent[#,timeStep,patience], space]&/@population, recoveryTime];

Agent patience

We assume that patience is randomly distributed with a decision made at each time-step on whether to abandon isolation. The half-life of the population remaining in isolation is set with the patience parameter. In order to calculate the individual decision probabilities we must solve the following equation

In[]:=

Quiet@FullSimplify[Assuming[{timeStep>0,patience>0,0<p<1},Solve[Probability[t>0,tBinomialDistribution[patience/timeStep(*Numberofdecisions*),p]]1/2,p]]]

Out[]=

p1-

timeStep

patience



Analyzing behaviour

Now that we can simulate the behaviour of our agents, we need to collect data from experiments for each time-step we are interested in counting the number of agents in each state:

In[]:=

countStates[population_List]:=Counts[ReplaceAll[population[[All,"Status"]],_Real"Infected"]]

For example:

In[]:=

countStates[testPopulation]

Out[]=

Infected1,Vulnerable9

We now build an experiment by generating a population and then repeatedly applying the updatePopulation function to it and counting the states until there are no infected agents remaining:

In[]:=

runExperiment[n_(*Sizeofpopulation*),diseased_Integer,immune_Integer,recoveryTime_,isolationRate_,delay_Real,patience_Real,timeStep_Real,space_]:= Block[{population=initialize[n,diseased,immune,recoveryTime,isolationRate,delay,space],history={}}, While[Count[population[[All,"Status"]],_Real]>0, population=updatePopulation[population,timeStep,recoveryTime,patience,space]; AppendTo[history,countStates[population]]]; history]

In[]:=

experiment=runExperiment[10,5,0,20.,0.,0.,0.,1.,20.]

Out[]=

{Infected5,Vulnerable5,Infected5,Vulnerable5,Infected5,Vulnerable5,Infected5,Vulnerable5,Infected5,Vulnerable5,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Infected6,Vulnerable4,Immune5,Infected1,Vulnerable4,Immune5,Infected1,Vulnerable4,Immune5,Infected1,Vulnerable4,Immune5,Infected1,Vulnerable4,Immune5,Infected2,Vulnerable3,Immune5,Infected2,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune6,Infected1,Vulnerable3,Immune7,Vulnerable3}

This simulation data can be analyzed in many ways, but here we will simply visualize it:

In[]:=

historyPlot[history_]:= If[ Length[history]<2, Graphics[{},ImageSize400], ListLinePlot[ history[[All,"Infected"]], history[[All,"Immune"]], history[[All,"Vulnerable"]] , ImageSize450, PlotStyleNone, PlotLayout"Stacked", Filling{1{Bottom,Red},2{{1},Green},3{{2},Gray}}, PlotLegendsSwatchLegend[{Red,Green,Gray},{"Infected","Immune","Vulnerable"}], PerformanceGoal"Quality" ]]

Example experiments

For all experiments, I will use 200 agents in a 25 x 25 space and 2 initial infections.

With no isolation strategy infections rapidly rise to 100% following an approximately logistic curve.

In[]:=

historyPlot[runExperiment[200(*agents*),2(*infected*),0(*immune*),40.(*recoverytime*),0.(*isolationrate*),0.(*isolationdelay*),0.(*isolationpatience*),1.(*timestep*),25(*space*)]]

Out[]=

	Infected
	Immune
	Vulnerable

With the rapid implementation of a 100% isolation policy the infections can be halted.

In[]:=

historyPlot[runExperiment[200(*agents*),2(*infected*),0(*immune*),40.(*recoverytime*),1(*isolationrate*),10.(*isolationdelay*),2000.(*isolationpatience*),1.(*timestep*),25(*space*)]]

Out[]=

	Infected
	Immune
	Vulnerable

Such a policy is unrealistic. Let us first consider a very strong policy but a less patient population who gradually tires of the policy. The spread is briefly halted but then continues. The whole population is still infected but the delay has been sufficient to reduce the peak infection rate:

In[]:=

historyPlot[runExperiment[200(*agents*),2(*infected*),0(*immune*),40.(*recoverytime*),1(*isolationrate*),10.(*isolationdelay*),100.(*isolationpatience*),1.(*timestep*),25(*space*)]]

Out[]=

	Infected
	Immune
	Vulnerable

Lower initial uptake of isolation but greater patience produces a similar effect but prolongs the outbreak.

In[]:=

historyPlot[runExperiment[200(*agents*),2(*infected*),0(*immune*),40.(*recoverytime*),.95(*isolationrate*),10.(*isolationdelay*),200.(*isolationpatience*),1.(*timestep*),25(*space*)]]

Out[]=

	Infected
	Immune
	Vulnerable

In[]:=

historyPlot[runExperiment[200(*agents*),2(*infected*),0(*immune*),40.(*recoverytime*),.99(*isolationrate*),10.(*isolationdelay*),1000.(*isolationpatience*),1.(*timestep*),25(*space*)]]

Out[]=

	Infected
	Immune
	Vulnerable

For very strong isolation policies the outcome is better but can cause waves of infection.

In[]:=

historyPlot[runExperiment[200(*agents*),2(*infected*),0(*immune*),40.(*recoverytime*),0.99(*isolationrate*),10.(*isolationdelay*),1000.(*isolationpatience*),1.(*timestep*),25(*space*)]]

Out[]=

	Infected
	Immune
	Vulnerable

Interactive simulation

We can add a simple user interface using Manipulate

In[]:=

Manipulate[historyPlot[runExperiment[200(*agents*),2(*infected*),0(*immune*),40.(*recoverytime*),isolationRate,isolationDelay,isolationPatience,1.(*timestep*),25(*space*)]],{isolationRate,0.,1,Appearance"Labeled"},{isolationDelay,0.,80,Appearance"Labeled"},{{isolationPatience,5000.},50.,5000},SaveDefinitionsTrue]

Out[]=

Visualization of agents

To give a more detailed animated view of the spread we first need to create a visualization of the population:

In[]:=

diseasePlot[data_,recoveryTime_,space_]:= Graphics[ PointSize[0.05], Map[ Which[#Status==="Vulnerable",Gray, #Status==="Immune",Green, True,Blend[{Green,Red,Red},#Status/recoveryTime]], Point[#Position]&, data] , PlotRange{{0,space},{0,space}}, FrameTrue, ImageSize300, FrameTicksFalse];

We can now view our testPopulation data set :

In[]:=

diseasePlot[testPopulation,20(*recoverytime*),20(*space*)]

Out[]=

User interface

We can now add a user interface to this. For a smooth animation we need to reduce the timeStep between updates, but for performance reasons we do not need to accumulate all of these, now more frequent, simulation results:

In[]:=

DynamicModule[{timeStep=0.04(*Controlsanimationspeed*),recoveryTime=40.,space=25.,n=80,immune=0,population={},history={},time=0},Manipulate[population=initialize[n,diseased,immune,recoveryTime,isolationRate,delay,space];Grid[Button["Run simulation",history={};time=0;While[Count[population[[All,"Status"]],_Real]>0, Pause[0.01];(*AllowDynamicstorefreshforsmootherdisplay*) population=updatePopulation[population,timeStep,recoveryTime,patience,space]; If[Round[Mod[time,1/timeStep]]0,AppendTo[history,countStates[population]]];(*Donotrecordeverystep*) time++],Method"Queued"],SpanFromLeft,{Dynamic[diseasePlot[population,recoveryTime,space],TrackedSymbols{population}],Dynamic[historyPlot[history],TrackedSymbols{history}]}],{diseased,1,10},{{patience,2000.},50.,2000.},{isolationRate,0,1},{{delay,10},0,50}],SaveDefinitionsTrue]

Out[]=

POSTED BY: Jon McLoone

5 Replies

Sort By:

Michael Kelly

Michael Kelly, Wolfram Research Inc.

Posted 2 years ago

Jon I really like the clean coding style you have adopted here. The use of agents identified as Associations with updated properties and an updating rule is one that can be modified for a wide range of stochastic problems with a large number of parameters. I used agent models for economics and finance back in the 90's. Whenever I tried to model neoliberal economic models as advocated by the Chicago School of economics, I found that they only ever worked if you assumed everything worked perfectly like your experiment No 2 with immediate 100% isolation and that everyone in the market had perfect knowledge of prices. Such situations only arose about one in every 100,000 experiments. Journals were not interested in such experiments, especially since corporations owned the journals. Strange how that happens. Michael

POSTED BY: Michael Kelly

David Gurarie

Posted 1 year ago

Human interactions and ensuing transmission, are not described by randomly colliding balls. They involved scheduled and random social contacts, and their spatial and temporal organization are paramount.

POSTED BY: David Gurarie

Jon McLoone

Jon McLoone, Wolfram Research

Posted 1 year ago

Of course. This was meant to be a simple model and the illustration of an approach. By adding many more property keys (each persons social calendar, spatial environment, type of current activity etc) and more complex update rules (going to places in the caledar, calculating the type of contact based on each agents current location and activity) one could model that. Though performance will likely become a consideration.

POSTED BY: Jon McLoone

Michael Kelly

Michael Kelly, Wolfram Research Inc.

Posted 1 year ago

David The agent based simulations described above are not random. They take into account the known data to represent features such as the rate of infection, isolation, recovery and immunity. See for instance the variables included in the function runExperiment. Michael

POSTED BY: Michael Kelly

David Gurarie

Posted 1 year ago

Those features are relevant to disease progress, of course, and for some interventions. But my key reservation concerns the nature (patterns, types) of social mixing. Randomness could play some role (people may bump into each other 'randomly', once in a while). But large part of interactions are structured, along the lines mentioned by Jon. One example is Christopher Wolfram model, based on social network structure. We are working on IBM with structured social mixing that goes beyond networks, and we are planning to post it soon.

POSTED BY: David Gurarie

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback