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

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.

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}]]​​ ];

In[]:=

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

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

In[]:=

{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∞}

Out[]=

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)

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

2

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

In[]:=

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

In[]:=

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

Out[]=

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.

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];

In[]:=

For example:

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

In[]:=

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

Out[]=

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.

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]]​​ ];

In[]:=

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

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

In[]:=

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

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

In[]:=

p1-

-

timeStep

patience

2



Out[]=

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:

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

In[]:=

For example:

countStates[testPopulation]

In[]:=

Infected1,Vulnerable9

Out[]=

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:

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.]

In[]:=

{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}

Out[]=

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

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"​​ ]]

In[]:=

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.

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

In[]:=

Infected Immune Vulnerable

Out[]=

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

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

In[]:=

Infected Immune Vulnerable

Out[]=

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:

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

In[]:=

Infected Immune Vulnerable

Out[]=

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

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

In[]:=

Infected Immune Vulnerable

Out[]=

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

In[]:=

Infected Immune Vulnerable

Out[]=

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

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

In[]:=

Infected Immune Vulnerable

Out[]=

Interactive simulation

We can add a simple user interface using Manipulate

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]

In[]:=

Out[]=

​

Visualization of agents

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

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];

In[]:=

We can now view our testPopulation data set :

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

In[]:=

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:

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]

In[]:=

Out[]=