Message Boards

GROUPS:

Francisco Rodríguez

[NB] Epidemic simulation with a polygon container

Francisco Rodríguez, WOLFRAM|Alpha

Posted 1 year ago

6212 Views

4 Replies

21 Total Likes

Follow this post

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

enter image description here

Introduction

Inspired by this Washington Post article: https://www.washingtonpost.com/graphics/2020/world/corona-simulator/ ; I tried to replicate the simulation adding a few more parameters, and the possibility of using a different bounding box, to illustrate it with a known form and create more awareness of the current situation on how different decisions can change the outcome.

This is an agent based model, using random walk to model how agents move, and using a given polygon as “container”.

An important difference with Washington Post article model is that here I’m including also the possibility of dying besides recovering, and to be a bit more closer to the fact that there is a limited health care capacity, there is also a different state for cases that need care and don’t get it, and for cases that can be attended.

In this notebook I only run once for each scenario, to actually get statistical results we would need to run it many times.

Simulation specification

States

The following states are used in this model:

◼

0: Healthy

◼

1: Infected

◼

2: Infected with symptoms

◼

3: Severe case without proper care

◼

4: Severe case with care

◼

10: Recovered

◼

11: Dead

Parameters

◼

“ProbabilityContagion”: probability to get the infection when “touching” an infected agent.

◼

“DaysAsymptomaticMin”: minimum number of days to start showing symptoms

◼

“DaysAsymptomaticMax”: maximum number of days to start showing symptoms

◼

“ProbabilitySevere”: probability of having a severe case

◼

“ProbabilitySevereDies”: probability a severe case (with proper care) dies

◼

“ProbabilitySevereWithoutCareDies”: probability a severe case without any care dies

◼

“HealthCapacity”: total capacity to attend severe cases

◼

“DaysOutcomeMin”: minimum number of days to get an outcome (recovered or dead)

◼

“DaysOutcomeMax”: maximum number of days to get an outcome (recovered or dead)

◼

“TimeStep”: time for each movement step

◼

“ReinfectionPropability”: probability to get the infection again (for a recovered case)

◼

“StatesInfecting”: Which states infect others, that is, if all infected cases will infect others, or if only infected without symptoms will, or only cases under care won’t. Possible recommended values: 1; 1 | 2; 1 | 2 | 3; 1 | 2 | 3 | 4.

Min/max parameters are used to compute a normal distribution around them with a mean centered between min and max, and standard deviation of

1

max-min

Simulation code

Data Structures

I will use a list of integer positions, and a list of states of the form {state, days infected, movement}, with the following values:

◼

state: one of the state values (0, 1, 2, 3, 4, 10, 11}

◼

days infected: time counter used only on infected cases

◼

movement: 1 if agent moves; 0 if it doesn’t; and -1 if it doesn’t move (while severe infected), but used to move (so, when it gets cured, it will move again).

At the end the full state of the situation will be given by {positions, states}.

Constants

Let’s associate a color with each state:

In[]:=

statecolor=Association[0Gray,1Pink,2Red,3Darker@Red,4Brown,10Green,11Black]

Out[]=

0

,1

,2

,3

,4

,10

,11



Pattern of infected states:

In[]:=

$infectedstates=1|2|3|4;

Create simulation case

We want to generate the initial conditions for a given container (defined by a region, entity, rectangle, polygon). We will rescale the region in case we want to fit a given number of entities with a target density.

Given a region/entity and number of elements, this returns the rescaled region, a region member function, the intial position and states.

About Region

We need a function to modify the input region and rescale it such way that we can fit integer coordinates for the agents, and ensure there is enough room for them. This function will return the modified region, the corresponding region member function and the rescale transform used to change the region. That rescale transform function is needed later to also transform the coordinates used as infection starting point.

In[]:=

Options[prepareRegion]={};prepareRegion[r:_Region|_Entity|_Polygon,targetarea_]:=Module[{region,area,newregion},If[MatchQ[r,_Entity],region=EntityValue[r,"Polygon"]/.g_GeoPositiong["LongitudeLatitude"];,region=r;];area=Area[region];If[!NumericQ[area],Return[{$Failed,$Failed,$Failed}];];newregion=RegionResize[region,Abs[Subtract@@RegionBounds[region][[1]]]Sqrt[targetarea/area]];{newregion,RegionMember[newregion],RescalingTransform[RegionBounds[region],RegionBounds[newregion]]}];

Initial conditions

The initial conditions for the simulation are the position of the agents and the state of each one. To generate the random points inside the region we can use RandomPoint, but because we are rounding the values, we need to use Select to verify the are still inside afeter rounding and remove duplicates, after that, we need to add more iteratively, to ensure we get the requested number of elements (and remember to always add a limit counter to avoid an infinite loop!).

createPositions[n_,region_,regionq_:Automatic]:=Module[{memberq=Replace[regionq,AutomaticRegionMember[region]],positions,i=0},positions={};While[Length[positions]<n&&i++<n,positions=DeleteDuplicates[Join[positions,Select[Round@RandomPoint[region,n-Length[positions]],memberq]]];];If[Length[positions]<n,$Failed,positions]];

To initialize the state of each point, we use a random 1/0 initialization using the isolation proportion to define if they move or not, and then we use different strategies to define which ones are infected. It can be a given number (chosen randomly), or we can provide a fixed list (in case that index is computed externally), or use UpTo as head to say we want the first index to be the ones infected.

In[]:=

createStates[n_,infected:{__Integer}|_Integer|UpTo[_Integer],isolationProp_]:=Module[{states},states=Table[{0,0,RandomChoice[{isolationProp,1-isolationProp}{0,1}]},n];If[!TrueQ[And@@MatchQ[infected,_List]],Return[$Failed]];Switch[infected,_List,states[[infected]]=ConstantArray[{1,0,1},Length[infected]];,_UpTo,states[[;;infected[[1]]]]=ConstantArray[{1,0,1},infected[[1]]];,_,states[[RandomSample[Range[n],infected]]]=ConstantArray[{1,0,1},infected];];states];

In[]:=

searchInfectionFocus[posisitions_,point_,transformation_,i_]:=Nearest[posisitions"Index",transformation[point],i];

Simulation case

The main function to generate the initial state for the simulation calls the previous initialization functions. This function is the one that returns the region, the region member function, the initial positions and states.

In[]:=

Options[createSimulation]:={"Infected"1,"Density"Automatic,"InfectionPoint"None,"Isolation"0.1}createSimulation[n_,r:_Region|_Entity|_Polygon,OptionsPattern[]]:=Module[{region,regionq,positions,transform,infected,density=OptionValue["Density"]/.Automatic1/10,ninfected,isolation,states},{region,regionq,transform}=prepareRegion[r,n/density];positions=createPositions[n,region];ninfected=OptionValue["Infected"];Which[MatchQ[OptionValue["InfectionPoint"],_Entity],infected=searchInfectionFocus[positions,EntityValue[OptionValue["InfectionPoint"],"Position"]["LongitudeLatitude"],transform,ninfected];,MatchQ[OptionValue["InfectionPoint"],{_?NumericQ,_?NumericQ}],infected=searchInfectionFocus[positions,OptionValue["InfectionPoint"],transform,ninfected];,True,infected=ninfected;];isolation=OptionValue["Isolation"];states=createStates[n,infected,isolation];{region,regionq,positions,states}];

Evolution

Main routines

The main rutines are the ones to run the whole simulation and compute the evolution of each step.

In[]:=

runSimulation[positions_,states_,totalsteps_,regionq_,parameters_]:=NestList[evolutionStep[#,regionq,parameters]&,{positions,states},totalsteps];

For each step we do 3 things:

◼

First, compute the movement of agents, avoiding overlapping and ensuring all remain inside the region (using the region member function).

◼

Second, compute new infections and increase the infection evolution of each infected agent. In this step it is computed if one gets worse, recovers or dies.

◼

Third, compute which ones get health care, which is based on the health care capacity. This is done by taking a look on the total health care capacity, and how many cases are in the state “infected with care”, and if there are available free slots, some random sever cases without care are changed to the care state.

In[]:=

evolutionStep[{pos_,state_},memberq_,params_]:=Module[{newpos,newposi,newstates,contacts,infected,contacted,severe,usedcapacity,n=Length[pos]},(*Computemovement*)newpos=pos;Do[If[!MatchQ[state[[i]],{_Integer,_?NumericQ,_Integer}],Print[state[[i]]];Abort[]];newposi=If[state[[i,3]]>0,RandomChoice[{{-1,0},{1,0},{0,-1},{0,1}}]+newpos[[i]],newpos[[i]]];If[!MemberQ[newpos,newposi]&&memberq[newposi],newpos[[i]]=newposi;];,{i,1,n}];(*Infectionevolution/transmision*)contacts=getContacts[newpos];infected=Flatten@Position[state[[All,1]],params[["StatesInfecting"]]];contacted=Complement[Cases[contacts,{x_,Alternatives@@infected}|{Alternatives@@infected,x_}x],infected];newstates=Table[Switch[state[[i,1]],0/;MemberQ[contacted,i]&&RandomReal[]<params[["ProbabilityContagion"]],ReplacePart[state[[i]],11],10/;MemberQ[contacted,i]&&RandomReal[]<params[["ReinfectionPropability"]],{1,0,state[[i,3]]},$infectedstates,infectedEvolution[state[[i]]+{0,params[["TimeStep"]],0},params],_,state[[i]]],{i,1,n}];(*Distributehealthcapacity*)usedcapacity=Count[newstates,{4,_,_}];If[usedcapacity<params[["HealthCapacity"]],severe=Flatten@Position[newstates[[All,1]],3];If[Length[severe]>0,severe=RandomSample[severe,Min[Length[severe],params[["HealthCapacity"]]-usedcapacity]];newstates[[severe,1]]=4;]];{newpos,newstates}];

Helper functions

Helper functions are called from the main routines.

To compute which agent is neighbour of which one, we first sort them by x coordinate, then by y coordinate, and after grouping them, we have a smaller subset to test the actual neighbours.

In[]:=

getContacts[pos_]:=With[{r=Range[Length[pos]]},Select[Intersection[Sort/@(Sequence@@Subsets[#,{2}]&/@Split[SortBy[r,pos[[#,1]]&],Abs[pos[[#1,1]]-pos[[#2,1]]]≤1&]),Sort/@(Sequence@@Subsets[#,{2}]&/@Split[SortBy[r,pos[[#,2]]&],Abs[pos[[#1,2]]-pos[[#2,2]]]≤1&])],closeq[pos[[#[[1]]]],pos[[#[[2]]]]]&]];

In[]:=

closeq[p1_,p2_]:=MatchQ[Abs[p1-p2],{0|1,0|1}]

This is one of the core routines, how an infected agent evolves. It takes into consideration all the time parameters and the probabilities to jump between the different infected states. Depending on the random values, an infected person can become symptomatic at some point, but be a light case or severe, or if very very lucky recover directly. If it is a light case, it will eventually recover, and only severe cases may recover or die (depending of care state it has a different probability for each outcome).

In[]:=

infectedEvolution[state_,params_]:=With[{day=state[[2]],dayssymptoms:=randomNormal[params[["DaysAsymptomaticMin"]],params[["DaysAsymptomaticMax"]]],daysoutcome:=randomNormal[params[["DaysOutcomeMin"]],params[["DaysOutcomeMax"]]]},Switch[state[[1]],1,If[day>dayssymptoms,If[day>daysoutcome,{10,day,Abs@state[[3]]},{If[RandomReal[]<params[["ProbabilitySevere"]],3,2],day,-Abs@state[[3]]}],If[day>daysoutcome,{10,day,Abs@state[[3]]},state]],2,If[day>daysoutcome,{10,day,Abs@state[[3]]},state],3,If[day>daysoutcome,If[RandomReal[]<params[["ProbabilitySevereWithoutCareDies"]],{11,day,0},{10,day,Abs@state[[3]]}],state],4,If[day>daysoutcome,If[RandomReal[]<params[["ProbabilitySevereDies"]],{11,day,0},{10,day,Abs@state[[3]]}],state],_,state]]

In[]:=

randomNormal[min_,max_]:=RandomVariate[NormalDistribution[(min+max)/2,1/Sqrt[max-min]]]

Other helper functions

Helper functions to determine when a simulation has reached a stable point: no more infected cases.

In[]:=

stableQ[allstates_]:=Total[Lookup[Rule@@@Tally[allstates[[-1,2,All,1]]],List@@$infectedstates,0]]===0;

In[]:=

stableTime[allstates_]:=Module[{data},data=Rule@@@#&/@(Tally/@allstates[[All,2,All,1]]);data=Total[Lookup[#,List@@$infectedstates,0]]&/@data;First[Flatten[Position[data,0,1,1]],None]];

In[]:=

stableTimeDays[allstates_]:=With[{length=Length[allstates],te=stableTime[allstates]},If[IntegerQ[te],Min[length,1+Ceiling[te,24]],length]];

Display functions

We build a main display function which takes all the steps done in the simulation, and the region, and can plot it for any intermediate step.

In[]:=

showAll//Clear;Options[showAll]={"Language""es",Alignment"Horizontal","MapSize"500,"Title"None};showAll[allstates_,regionbase_,index_Integer,maxsteps_,params_,OptionsPattern[]]:=Module[{time,max,g,plot,summary,signature,title,i},i=Clip[index,{1,Length[allstates]}];max=Replace[maxsteps,AutomaticLength[allstates]];title=OptionValue["Title"];If[title===None,title=Sequence[];,title=Text@Style[title,Large];];time=Text@Row[{"t=",showTime[i,params[["TimeStep"]]]}];g=Show[diplayRegion[regionbase,allstates[[i]]],ImageSizeOptionValue["MapSize"]];plot=Show[miniplot[allstates,i,max,"paso"params[["TimeStep"]]],ImageSize250];summary=displaySummary[allstates[[i,2]],"Language"OptionValue["Language"]];signature=Text@Style[Replace[OptionValue["Language"],{"en""Created in Wolfram Language by FJRA","es""Hecho en Wolfram Language por FJRA"}],Small,Gray];Switch[OptionValue[Alignment],"Horizontal",Column[{title,Row[{g,Column[{time,plot,summary}]},AlignmentTop],signature},AlignmentCenter],"Vertical2",Column[{title,time,Row[{plot,summary}],g,signature},AlignmentCenter],_,Column[{title,time,g,Row[{plot,summary}],signature},AlignmentCenter]]];showAll[___]:=$Failed;

Display helper functions


Running a single simulation case

Simulation with Peru background

We need to define the parameters for the simulation in an Association. We are setting a mortality probability of 10% of the severe cases with care (which are already a 20%, giving a global 1%), and 80% for the sever cases without care. Also the time step is given as 1h. Another important parameter is the health care capacity, set to 3 as example, for this total population of 200.

In[]:=

parameters=Association["ProbabilityContagion"0.9,"DaysAsymptomaticMin"5,"DaysAsymptomaticMax"14,"ProbabilitySevere"0.2,"ProbabilitySevereDies"0.1,"ProbabilitySevereWithoutCareDies"0.8,"HealthCapacity"3,"DaysOutcomeMin"12,"DaysOutcomeMax"21,"TimeStep"1/24,"ReinfectionPropability"0.0001,"StatesInfecting"1|2|3];

Let’s run it for 10% isolation

To initialize the simulation, we can use Peru country entity, and set that the first infect started close to Lima. We should remember this is a model using the polygon as container, and it should not be confused as a predictive model using actual geographic places.

{region,regionq,initialpos,initialstates10}=createSimulation200,

Peru

COUNTRY

,"Infected"1,"Density"0.075,"InfectionPoint"

Lima

CITY

,"Isolation"0.1;

We run the simulation over 90 days.

In[]:=

simulationperu10=runSimulation[initialpos,initialstates10,24*90,regionq,parameters];//AbsoluteTiming

Out[]=

{73.1277,Null}

To create a lighter manipulate with each step, we can use the time to get all infected gone as ending step.

In[]:=

With[{max=stableTimeDays[simulationperu10]},Manipulate[showAll[simulationperu10,region,i,max,parameters,"Title""Isolation 10%","Language""en",Alignment"Vertical2"],{i,1,max,1}]]

Out[]=

Now with 90% isolation

To simulate again with a different isolation parameter, we don’t need to return the region and regionq function, because we already had them for the same place. Then we run and plot with manipulate this new case.

{initialpos90,initialstates90}=createSimulation200,

Peru

COUNTRY

,"Infected"1,"Density"0.075,"InfectionPoint"

Lima

CITY

,"Isolation"0.90[[{3,4}]]

In[]:=

simulationperu90=runSimulation[initialpos90,initialstates90,24*90,regionq,parameters];//AbsoluteTiming

Out[]=

{53.4722,Null}

In[]:=

With[{max=stableTimeDays[simulationperu90]},Manipulate[showAll[simulationperu90,region,i,max,parameters,"Title""Isolation 90%","Language""en",Alignment"Vertical2"],{i,1,max,1}]]

Out[]=

Visualize both together

We can now put both simulations with only one different parameter in the same manipulate using the maximum time to get the stable state.

In[]:=

With[{max=Max[stableTimeDays[simulationperu10],stableTimeDays[simulationperu90]]},Manipulate[Row[{showAll[simulationperu10,region,i,max,parameters,"Title""Isolation 10%","Language""en",Alignment"Vertical2"],showAll[simulationperu90,region,i,max,parameters,"Title""Isolation 90%","Language""en",Alignment"Vertical2"]}],{i,1,max,1}]]

Out[]=

Simulation with Great Britain island background

For the Great Britain island, let’s use the same parameters, but increase the health capacity.

In[]:=

parametersb=parameters=Association["ProbabilityContagion"0.9,"DaysAsymptomaticMin"5,"DaysAsymptomaticMax"14,"ProbabilitySevere"0.2,"ProbabilitySevereDies"0.1,"ProbabilitySevereWithoutCareDies"0.8,"HealthCapacity"5,"DaysOutcomeMin"12,"DaysOutcomeMax"21,"TimeStep"1/24,"ReinfectionPropability"0.0001,"StatesInfecting"1|2|3];

Great Britain island polygon

If we want to use a nice projection, we cannot use the automatic polygon extraction from the code, but we can project it manually:

In[]:=

projection={"LambertAzimuthal","Centering"GeoPosition[{54.31580165`,-2.2320619`}]};

In[]:=

britain=OuterPolygon

Great Britain

ISLAND

["Polygon"]/.g_GeoPositionGeoGridPosition[g,projection]["GridXY"]

Out[]=

Polygon

Number of points: 1246

Embedding dimension: 2



In[]:=

london=GeoGridPosition

London

CITY

,projection["GridXY"]

Out[]=

{0.0229868,-0.0487984}