City-Scale Hybrid Simulation Modeling of COVID-19 Spread

Charles M. Macal, Ph.D., P.E.
Argonne National Laboratory

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

As the COVID-19 pandemic continues to unfold, the epidemiological modeling community continues to support decision makers with forecasts of how the disease could progress under a variety of assumptions: the emergence of COVID variants; the public’s response to resurging COVID-related cases, hospitalizations, and deaths; new vaccination regimes; and potential interventions by decision makers.  In this talk, we describe two complementary Mathematica epidemiological simulation modeling approaches, one an agent-based model (ABM), the other, a compartmental model based on System Dynamics (SD), which have been developed to forecast COVID-19 spread at the city scale, specifically for Chicago. The COVID ABM simulates the movements, contacts, and disease transmission outcomes of every person in a city, installation of facility as they go about their daily activity schedules. Complementing the agent-based model, the System Dynamics COVID model segments the population, compartmentalized by current disease state. Used together, the hybrid system combining the top-down SD model with the bottom-up ABM provides a broader understanding of disease spread than is available with either approach alone. The multiparadigm modeling capability of Mathematica is found to easily and elegantly express essential formalisms of agent-based and System Dynamics modeling. The talk concludes on how agent-based modeling and System Dynamics can work together to best inform public health decision-making.

Key Questions :

How? Can we model disease spread in a city population and the actions and behaviors of people at scale?

Why? Provide forecasts of disease spread with uncertainty ranges to decision makers in real time and understand the possible effects of interventions they might take to mitigate the spread.

Quick Background on Epidemiological Modeling

1. Statistical Models: E.g., The IHME Model

2. Compartmental Models:

Population divided into compartments based on disease state and possibly other characteristics (age)
SEIR differential equation models: Susceptible, Exposed, Infected, Recovered

3. Agent-based (Mechanistic) Models:

Represent all individuals in the population as agents
Agents: characteristics, behaviors, disease state

Modeling A City: People, Places and Activities

Basic Idea of Modeling Disease Transmission

Persons   +   Places   +   Activity Schedules

=>	Mobility (persons to places according to their hourly activity schedules)
	Co-location & Contact
	Infectious-to-Susceptible Disease Transmission Possibility

Modeling the Progression of a Disease (Covid-19)

Compartmental Model:	Disease State of Population Compartment
Agnet-based Model: 	Disease State of Individual People

Compartmental Model

System Dynamics Model with difference equations for disease state updates

State Variable:  state[time, ageGroup, raceEthnicity, zone (zipcode)]

ClearAll["susceptibleIn","susceptibleOut"];

reth=1;
zone="ChicagoV3";

susceptibleIn[0, age_,reth_,zone_] := 0.;
susceptibleIn[t_ /; t > 0, age_,reth_,zone_] := 0;
susceptibleOut[0, age_,reth_,zone_] := 0.;
susceptibleOut[t_ /; t > 0, age_,reth_,zone_] := newcountsMedian[[age,t]]*popPortionRaceEthn[reth](*breakdown newcounts data which is by age, into raceEthnicity portions*);

susceptible[t_/;t>0,age_,reth_,zone_]:=susceptible[t,age,reth,zone]=susceptible[t-1,age,reth,zone]+susceptibleIn[t,age,reth,zone]-susceptibleOut[t,age,reth,zone](*Max[0.,susceptible[t-1,age,reth,zone]+susceptibleIn[t,age,reth,zone]-susceptibleOut[t,age,reth,zone]]*);

ClearAll["exposed","exposedIn","exposedOut"];
exposedIn[0,age_,reth_,zone_]:=susceptibleOut[0,age,reth,zone];
exposedIn[t_/;t>0,age_,reth_,zone_]:=susceptibleOut[t,age,reth,zone];
exposed[0, age_,reth_,zone_] := exposed[0, age,reth,zone ]= 0.(*patient zero for each age group*);exposed[t_/;t>0,age_,reth_,zone_]:=exposed[t,age,reth,zone]=exposed[t-1,age,reth,zone]+exposedIn[t,age,reth,zone]-exposedOut[t,age,reth,zone];
exposedOut[0, age_,reth_,zone_] := exposedOut[0, age,reth,zone] = 0.;
exposedOut[t_ /; t > 0, age_,reth_,zone_] := exposed[t - 1, age,reth,zone]*(1./exposedduration[1,t]);

ClearAll["asymptomatic","asymptomaticIn","asymptomaticOut"];
asymptomaticIn[0, age_,reth_,zone_] := asymptomaticIn[0, age,reth,zone] = 0.;
asymptomaticIn[t_ /; t > 0, age_,reth_,zone_] := exposedOut[t, age,reth,zone]*etoaprobabilities[age,t];asymptomaticOut[0, age_,reth_,zone_] := asymptomaticOut[0, age,reth,zone] = 0.;
asymptomaticOut[t_ /; t > 0, age_,reth_,zone_] := asymptomatic[t - 1, age,reth,zone]*(1./asymptomaticduration[1,t]);asymptomatic[0, age_,reth_,zone_] := asymptomatic[0, age,reth,zone] = 0.;asymptomatic[t_/;t>0,age_,reth_,zone_]:=asymptomatic[t,age,reth,zone]=asymptomatic[t-1,age,reth,zone]+asymptomaticIn[t,age,reth,zone]-asymptomaticOut[t,age,reth,zone];

ClearAll["presymptomatic","presymptomaticIn","presymptomaticOut"];
presymptomaticIn[0, age_,reth_,zone_] := presymptomaticIn[0, age,reth,zone] = 0.;
presymptomaticIn[t_ /; t > 0, age_,reth_,zone_] := exposedOut[t, age,reth,zone]*(1. - etoaprobabilities[age,t]);presymptomaticOut[0, age_,reth_,zone_] := presymptomaticOut[0, age,reth,zone] = 0.;presymptomaticOut[t_ /; t > 0, age_,reth_,zone_] := presymptomatic[t - 1, age,reth,zone]*(1./presymptomaticduration[1,t]);presymptomatic[0, age_,reth_,zone_] := presymptomatic[0, age,reth,zone] = 0.;
presymptomatic[t_/;t>0,age_,reth_,zone_]:=presymptomatic[t,age,reth,zone]=presymptomatic[t-1,age,reth,zone]+presymptomaticIn[t,age,reth,zone]-presymptomaticOut[t,age,reth,zone];

ClearAll["infectedNhR","infectedNhRIn","infectedNhROut"];
infectedNhRIn[0, age_,reth_,zone_] := infectedNhRIn[0, age,reth,zone] = 0.;
infectedNhRIn[t_ /; t > 0, age_,reth_,zone_] := presymptomaticOut[t, age,reth,zone]*ptoiprobabilities[age,t][[2]];infectedNhROut[0, age_,reth_,zone_] := infectedNhROut[0, age,reth,zone] = 0.;infectedNhROut[t_ /; t > 0, age_,reth_,zone_] := infectedNhR[t - 1, age,reth,zone]*(1./nhtorduration[1,t]);
infectedNhR[0, age_,reth_,zone_] := infectedNhR[0, age,reth,zone] = 0.;infectedNhR[t_/;t>0,age_,reth_,zone_]:=infectedNhR[t,age,reth,zone]=infectedNhR[t-1,age,reth,zone]+infectedNhRIn[t,age,reth,zone]-infectedNhROut[t,age,reth,zone];

ClearAll["infectedNhD","infectedNhDIn","infectedNhDOut"];
infectedNhDIn[0, age_,reth_,zone_] := infectedNhDIn[0, age,reth,zone] = 0.;infectedNhDIn[t_/;t>0,age_,reth_,zone_]:=presymptomaticOut[t,age,reth,zone] (1.-Total[ptoiprobabilities[age,t]]);infectedNhDOut[0, age_,reth_,zone_] := infectedNhDOut[0, age,reth,zone] = 0.;infectedNhDOut[t_ /; t > 0, age_,reth_,zone_] := infectedNhD[t - 1, age,reth,zone]*(1./nhtodduration[1,t]);infectedNhD[0, age_,reth_,zone_] := infectedNhD[0, age,reth,zone] = 0.;
infectedNhD[t_/;t>0,age_,reth_,zone_]:=infectedNhD[t,age,reth,zone]=infectedNhD[t-1,age,reth,zone]+infectedNhDIn[t,age,reth,zone]-infectedNhDOut[t,age,reth,zone];

ClearAll["infectedSevere","infectedSevereIn","infectedSevereOut"];
infectedSevereIn[0, age_,reth_,zone_] := infectedSevereIn[0, age,reth,zone] = 0.;
infectedSevereIn[t_ /; t > 0, age_,reth_,zone_] := presymptomaticOut[t, age,reth,zone]*ptoiprobabilities[age,t][[1]];infectedSevereOut[0, age_,reth_,zone_] := infectedSevereOut[0, age,reth,zone] = 0.;
infectedSevereOut[t_ /; t > 0, age_,reth_,zone_] := infectedSevere[t - 1, age,reth,zone]*(1./severeduration[1,t]);infectedSevere[0, age_,reth_,zone_] := infectedSevere[0, age,reth,zone] = 0.;
infectedSevere[t_/;t>0,age_,reth_,zone_]:=infectedSevere[t,age,reth,zone]=infectedSevere[t-1,age,reth,zone]+infectedSevereIn[t,age,reth,zone]-infectedSevereOut[t,age,reth,zone];

ClearAll["hospD","hospDIn","hospDOut"];
hospDIn[0, age_,reth_,zone_] := hospDIn[0, age,reth,zone] = 0.;
hospDIn[t_/;t>0,age_,reth_,zone_]:=infectedSevereOut[t,age,reth,zone] (1.-Total[stohprobabilities[age,t]]);hospDOut[0, age_,reth_,zone_] := hospDOut[0, age,reth,zone] = 0.;hospDOut[t_ /; t > 0, age_,reth_,zone_] := hospD[t - 1, age,reth,zone]*(1./hospitaldeadduration[1,t]);hospD[0, age_,reth_,zone_] := hospD[0, age,reth,zone] = 0.;
hospD[t_/;t>0,age_,reth_,zone_]:=hospD[t,age,reth,zone]=hospD[t-1,age,reth,zone]+hospDIn[t,age,reth,zone]-hospDOut[t,age,reth,zone];

ClearAll["hospIcuD","hospIcuDIn","hospIcuDOut"];
hospIcuDIn[0, age_,reth_,zone_] := hospIcuDIn[0, age,reth,zone] = 0.;
hospIcuDIn[t_ /; t > 0, age_,reth_,zone_] := infectedSevereOut[t, age,reth,zone]*stohprobabilities[age,t][[3]];
hospIcuDOut[0, age_,reth_,zone_] := hospIcuDOut[0, age,reth,zone] = 0.;
hospIcuDOut[t_ /; t > 0, age_,reth_,zone_] := hospIcuD[t - 1, age,reth,zone]*(1./hospitalicudeadduration[1,t]);hospIcuD[0, age_,reth_,zone_] := hospIcuD[0, age,reth,zone] = 0.;
hospIcuD[t_/;t>0,age_,reth_,zone_]:=hospIcuD[t,age,reth,zone]=hospIcuD[t-1,age,reth,zone]+hospIcuDIn[t,age,reth,zone]-hospIcuDOut[t,age,reth,zone];

ClearAll["hospIcuR","hospIcuRIn","hospIcuROut"];
hospIcuRIn[0, age_,reth_,zone_] := hospIcuRIn[0, age,reth,zone] = 0.;
hospIcuRIn[t_ /; t > 0, age_,reth_,zone_] := infectedSevereOut[t, age,reth,zone]*stohprobabilities[age,t][[2]];hospIcuROut[0, age_,reth_,zone_] := hospIcuROut[0, age,reth,zone] = 0.;
hospIcuROut[t_ /; t > 0, age_,reth_,zone_] := hospIcuR[t - 1, age,reth,zone]*(1./hospitalicurecoveredduration[1,t]);hospIcuR[0, age_,reth_,zone_] := hospIcuR[0, age,reth,zone] = 0.;hospIcuR[t_/;t>0,age_,reth_,zone_]:=hospIcuR[t,age,reth,zone]=hospIcuR[t-1,age,reth,zone]+hospIcuRIn[t,age,reth,zone]-hospIcuROut[t,age,reth,zone];

ClearAll["hospR","hospRIn","hospROut"];
hospRIn[0, age_,reth_,zone_] := hospRIn[0, age,reth,zone] = 0.;
hospRIn[t_ /; t > 0, age_,reth_,zone_] := infectedSevereOut[t, age,reth,zone]*stohprobabilities[age,t][[1]];hospROut[0, age_,reth_,zone_] := hospROut[0, age,reth,zone] = 0.;hospROut[t_ /; t > 0, age_,reth_,zone_] := hospR[t - 1, age,reth,zone]*(1./hospitalrecoveredduration[1,t]);hospR[0, age_,reth_,zone_] := hospR[0, age,reth,zone] = 0.;
hospR[t_/;t>0,age_,reth_,zone_]:=hospR[t,age,reth,zone]=hospR[t-1,age,reth,zone]+hospRIn[t,age,reth,zone]-hospROut[t,age,reth,zone];

ClearAll["icuR","icuRIn","icuROut"];
icuRIn[0, age_,reth_,zone_] := icuRIn[0, age,reth,zone] = 0.;
icuRIn[t_ /; t > 0, age_,reth_,zone_] := hospIcuROut[t, age,reth,zone];
icuROut[0, age_,reth_,zone_] := icuROut[0, age,reth,zone] = 0.;
icuROut[t_ /; t > 0, age_,reth_,zone_] := icuR[t - 1, age,reth,zone]*(1./icurecoveredduration[1,t]);icuR[0, age_,reth_,zone_] := icuR[0, age,reth,zone] = 0.;
icuR[t_ /; t > 0, age_,reth_,zone_] := icuR[t, age,reth,zone] = icuR[t - 1, age,reth,zone] + icuRIn[t, age,reth,zone] - icuROut[t, age,reth,zone];

ClearAll["icuD","icuDIn","icuDOut"];
icuDIn[0, age_,reth_,zone_] := icuDIn[0, age,reth,zone] = 0.;
icuDIn[t_ /; t > 0, age_,reth_,zone_] := hospIcuDOut[t, age,reth,zone];
icuDOut[0, age_,reth_,zone_] := icuDOut[0, age,reth,zone] = 0.;
icuDOut[t_ /; t > 0, age_,reth_,zone_] := icuD[t - 1, age,reth,zone]*(1./icudeadduration[1,t]);icuD[0, age_,reth_,zone_] := icuD[0, age,reth,zone] = 0.;
icuD[t_ /; t > 0, age_,reth_,zone_] := icuD[t, age,reth,zone] = icuD[t - 1, age,reth,zone] + icuDIn[t, age,reth,zone] - icuDOut[t, age,reth,zone];

ClearAll["recovered","recoveredIn"];recoveredIn[0, age_,reth_,zone_] := recoveredIn[0, age,reth,zone] = 0.;recoveredIn[t_ /; t > 0, age_,reth_,zone_] := asymptomaticOut[t, age,reth,zone] + infectedNhROut[t, age,reth,zone] +  icuROut[t, age,reth,zone] + hospROut[t, age,reth,zone];
recovered[0, age_,reth_,zone_] := recovered[0, age,reth,zone] = 0.;
recovered[t_ /; t > 0, age_,reth_,zone_] := recovered[t, age,reth,zone] = recovered[t - 1, age,reth,zone] + recoveredIn[t, age,reth,zone];

ClearAll["dead","deadIn"];
deadIn[0, age_,reth_,zone_] := deadIn[0, age,reth,zone] = 0.;
deadIn[t_ /; t > 0, age_,reth_,zone_] := infectedNhDOut[t, age,reth,zone] + icuDOut[t, age,reth,zone] + hospDOut[t, age,reth,zone];
dead[0, age_,reth_,zone_] := dead[0, age,reth,zone] = 0.;
dead[t_ /; t > 0, age_,reth_,zone_] := dead[t, age,reth,zone] = dead[t - 1, age,reth,zone] + deadIn[t, age,reth,zone];

Compartmental Model

Running the System Dynamics Model

reth={"W",0};
zone="ChicagoV3";
simResultsSD=Table[({
susceptible[t,#,reth,zone],
exposed[t,#1,reth,zone],
asymptomatic[t,#1,reth,zone],
presymptomatic[t,#1,reth,zone],
infectedNhR[t,#1,reth,zone],
infectedNhD[t,#1,reth,zone],
infectedSevere[t,#1,reth,zone],
hospR[t,#1,reth,zone],
hospIcuR[t,#1,reth,zone],
hospIcuD[t,#1,reth,zone],
hospD[t,#1,reth,zone],
icuR[t,#1,reth,zone],
icuD[t,#1,reth,zone],
recovered[t,#1,reth,zone],
dead[t,#1,reth,zone]}&)/@ages,{t,simDays}];

Compartmental Model Results

Agent-based Model

Develop Synthetic Populations 1-2-3

1. Persons
 	2. Places
 	3. Activities

Simulate the synthetic population over a 4-week period as people move from place to place (every hour) according to their daily activity schedules
 	1. Co-location by place
 	2. Disease transmission among people at each place

Synthetic Populations 1-2-3 for the ABM (1 of 3)

person Association

personPattern

AssociationQ[person]
KeyValueMap[Rule,person][[1]]
Print[Style["There are persons: ",Bold],Length@person];
pkey=randomKey[person];
showPerson@pkey

Synthetic Populations 1-2-3 for the ABM (2 of 3)

place Association

placePattern

AssociationQ[place]
KeyValueMap[Rule,place][[1]]
Print[Style["There are places: ",Bold],Length@place];
plkey=randomKey[place];
showPlace@plkey

Synthetic Populations 1-2-3 for the ABM (3 of 3)

schedule Association

schedulePattern

AssociationQ[schedule]
KeyValueMap[Rule,schedule][[1]]
Print[Style["There are schedules: ",Bold],Length@schedule];
skey=randomKey[schedule]
showSchedule[{skey,schedule@skey}]

Summary Synthetic Population

Simulating Disease Spread Over Time

Synthetic populations (person, place, schedule) are static

Simulation updates person states over time as people:

move and engage in activities at various places (location state)
become infected and transition through ensuing disease states (disease state)
adopt their daily schedules (schedule state)
become vaccinated (vaccination state)

State			state (composite)	{covid19State, locState, schedState, vaccState}
Disease state		covid19State 	 	covidS[time_,covidState {“S”,”E”,”A”,..}]
Location state	 	locState			locS[time_, locState {places}]
Daily schedule	 	schedState 		schedS[time_,  schedState: {1224..}]
Vaccination state	vaccState 		vaccS[time_, vaccState {“Moderna 1”,”Pfizer Bi-valent Booster”..}]

Agent-based Model

Events:

An event occurs to a person at a time.  Events have either (1) been previously scheduled (and come from the agent's event queue) or (2) are spontaneously generated.

The effect of an event is to (1) update the agent state via the event, (2) log the former state in the agent's history, and (3) schedule an associated next event (if any) in the person's event queue.

Events types: 
	Disease state: covid19State = covidS[timeEnteredState_,covidState {"S","E","A",..}]
	Location state: locState = locS[timeEnteredState_,locState {places}]
	Vaccination state: vaccState = vaccS[timeEnteredState_,vaccState {"Moderna 1",Pfizer Bi-valent Booster"..}]
	Daily schedule: schedState = schedS[timeEnteredState_,schedState {1224..}]

Next Events Table:

updatePersonNowEvent: mapped over person Association at time tNow for an event.

ClearAll[updatePersonNowEvent];

calls["updatePersonNowEvent"]=Sort[Names["*duration"]~Join~
Names["*probabilities"]~Join~{"getSchedulePersonNow"}];
called["updatePersonNowEvent"]={"MAIN"};

(* usage *)updatePersonNowEvent::usage="updatePersonNowEvent[key_String|key_Integer,now:{tNow_?NonNegative,weekday_?BooleanQ,_},event:eventHead_[tEvent_?NonNegative,iEvent_Integer|iEvent_String,___]]
updates person person@key state at time tNow for event eventHead_[tEvent_?NonNegative,iEvent_Integer|iEvent_String,___] occurence. Returns next event scheduled in event queue evQ."

updatePersonNowEvent::argumentsDontMatch="Arguments do not match. Pattern should match:
updatePersonNowEvent[key_String|key_Integer,now:{tNow_?NonNegative,weekday_?BooleanQ,_},event:eventHead_[tEvent_?NonNegative,iEvent_Integer|iEvent_String,___]].";

updatePersonNowEvent::wrongNumberOfArguments="Wrong number of arguments. Function call pattern should match: 
 updatePersonNowEvent[p:personPattern,now:{tBegin_?NonNegative,weekday_?BooleanQ,_}].";

(* don't know whether we need this *)updatePersonNowEvent[key_String|key_Integer,now:{tNow_?NonNegative,weekday_?BooleanQ,_},{}]:=Print[Style[#,14,Bold,Red]&@"MESSAGE: Event is {}."];

(* Main function *)
updatePersonNowEvent[key_String|key_Integer,now:{tNow_?NonNegative,weekday_?BooleanQ,_},event:eventHead_[tEvent_?NonNegative,iEvent_Integer|iEvent_String,___]]:=Module[{ageG,ageP,nextEvent,nextEventLocS,nextEventSchedS,nextEventTestS,nextEventVaccS,newlocS,newschedS,appendevQlocS},

If[updatePersonNowEventprintQ,
Print[Style[#,18,Bold,Blue]&@"ENTER updatePersonNowEvent::event: ",Style[#,18,Bold,Red]&@event];
Print["key: ",key]
];

ageP=person[key][[$ageSlot]](* assumes person ageG in personPattern is in $ageSlot *);
ageG=Min[9,Floor[ageP/10]+1];
Switch[eventHead,

covid19S,
(
(nextEvent=
Switch[iEvent,

"E",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[exposedduration[1,tNow],exposedduration[2,tNow]]],If[RandomReal[]<etoaprobabilities[ageG,tNow],"A","P"]],

"A",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[asymptomaticduration[1,tNow],asymptomaticduration[2,tNow]]],"R"],

"P",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[presymptomaticduration[1,tNow],presymptomaticduration[2,tNow]]],r=RandomReal[];Which[r<ptoiprobabilities[ageG,tNow][[1]],
"SH",r<ptoiprobabilities[ageG,tNow][[1]]+ptoiprobabilities[ageG,tNow][[2]],
"INhR",True,"INhD"]],

"INhR",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[nhtorduration[1,tNow],nhtorduration[2,tNow]]],"R"],

"INhD",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[nhtodduration[1,tNow],nhtodduration[2,tNow]]],"D"],

"SH",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[severeduration[1,tNow],severeduration[2,tNow]]],r=RandomReal[];Which[r<stohprobabilities[ageG,tNow][[1]],"HR",r<stohprobabilities[ageG,tNow][[1]]+stohprobabilities[ageG,tNow][[2]],"HIR",r<stohprobabilities[ageG,tNow][[1]]+stohprobabilities[ageG,tNow][[2]]+stohprobabilities[ageG,tNow][[3]],"HID",True,"HD"]],

"HR",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[hospitalrecoveredduration[1,tNow],hospitalrecoveredduration[2,tNow]]],"R"],

"HD",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[hospitaldeadduration[1,tNow],hospitaldeadduration[2,tNow]]],"D"],

"HIR",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[hospitalicurecoveredduration[1,tNow],hospitalicurecoveredduration[2,tNow]]],
"IR"],

"HID",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[hospitalicudeadduration[1,tNow],hospitalicudeadduration[2,tNow]]],"ID"],

"IR",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[icurecoveredduration[1,tNow],icurecoveredduration[2,tNow]]],"R"],

"ID",
covid19S[tEvent+24.` RandomVariate[GammaDistribution[icudeadduration[1,tNow],icudeadduration[2,tNow]]],"D"],

"D",Nothing,
"R",Nothing,
_,Nothing
];

If[updatePersonNowEventprintQ,
Print["\nbefore evQ@key: ",evQ@key];
Print["before"," covid19State@key: ",covid19State@key];
Print["before evH@key: ",evH@key];
Print["before state@key: ",state@key]
];
(*update evQ, state, evH*)
evH@key=Prepend[evH@key,covid19State@key](* add current state to history *);
evQ@key=Prepend[DeleteCases[evQ@key,event],nextEvent](* remove event from evQ. schedule next covid event(s) to add to evQ *);
covid19State@key=event(* update current state *);
If[updatePersonNowEventprintQ,
Print["\nupdate evQ@key: ",evQ@key];
Print["update covid19State@key: ",covid19State@key];
Print["update evH@key: ",evH@key];
Print["update state@key: ",state@key]
])
),


locS,
If[updatePersonNowEventprintQ,
Print["\nbefore evQ@key: ",evQ@key];
Print["before locState@key: ",locState@key];
Print["before evH@key: ",evH@key];
Print["before state@key: ",state@key]
];

evH@key=Prepend[evH@key,locState@key](* append current state to history *);
nextEvent=Nothing;
evQ@key=DeleteCases[evQ@key,event](* remove event from evQ. locS event does not schedule any new events, so nothing to add to evQ *);
locState@key=event(* update current state *);

If[updatePersonNowEventprintQ,
Print["\nupdate evQ@key: ",evQ@key];
Print["update locState@key: ",locState@key];
Print["update evH@key: ",evH@key];
Print["update state@key: ",state@key]
],


schedS,
If[updatePersonNowEventprintQ,
Print["\nbefore update evQ@key: ",evQ@key];
Print["before locState@key: ",locState@key];
Print["before schedState@key: ",schedState@key];
Print["before evH@key: ",evH@key];
Print["before state@key: ",state@key]
];

evH@key=Join[{locState@key,schedState@key},evH@key](* append current schedS and locaS states to history *);
{newlocS,newschedS,appendevQlocS}=getSchedulePersonNow[{key,person@key},{tEvent,weekday,_}];
evQ@key=Prepend[
Join[
DeleteCases[evQ@key,event],
appendevQlocS],
nextEvent=schedS[tEvent+24,0]
](* remove event locS. (schedule) add next schedS event to evQ*);
locState@key=newlocS(*update location state*);
schedState@key=newschedS(*update schedule state*);

If[updatePersonNowEventprintQ,
Print["\nupdate evQ@key: ",evQ@key];
Print["update locState@key: ",locState@key];
Print["update schedState@key: ",schedState@key];
Print["update evH@key: ",evH@key];
Print["update state@key: ",state@key]
],


testS,
If[updatePersonNowEventprintQ,
Print["\nbefore update evQ@key: ",evQ@key];
Print["before testState@key: ",testState@key];
Print["before evH@key: ",evH@key];
Print["before state@key: ",state@key]
];

evH@key=Prepend[evH@key,testState@key](*add current state to history*);
nextEvent=Nothing;
evQ@key=DeleteCases[evQ@key,event];
testState@key=event(* update current state *);

If[updatePersonNowEventprintQ,
Print["\nupdate evQ@key: ",evQ@key];
Print["update testState@key: ",testState@key];
Print["update evH@key: ",evH@key];
Print["update state@key: ",state@key]
],


vaccS,
If[updatePersonNowEventprintQ,
Print["\nbefore update evQ@key: ",evQ@key];
Print["before vaccState@key: ",vaccState@key];
Print["before evH@key: ",evH@key];
Print["before state: ",state@key]
];

evH@key=Prepend[evH@key,vaccState@key](* add current state to history *);
nextEvent=Nothing;
evQ@key=DeleteCases[evQ@key,event];
vaccState@key=event(* update current state *);

If[updatePersonNowEventprintQ,
Print["\nupdate evQ@key: ",evQ@key];
Print["update vaccState@key: ",vaccState@key];
Print["update evH@key: ",evH@key];
Print["update state: ",state@key]
],


(* eventHead *)_,Nothing(* default is current state is unchanged *)
];
nextEvent
];


(* Message::argumentsDontMatch: *)
updatePersonNowEvent[key_,{tNow_,weekday_,_},eventHead_[tEvent_,iEvent_Integer|iEvent_String,___]]:=Module[{trueFalseRule={True->Style[True,Blue,Bold,14],False->Style[False,Red,Bold,14]}},
Message[updatePersonNowEvent::argumentsDontMatch,argumentsDontMatch];
Print@TableForm[{
{Style[#,Bold,14]&@"key",key,Head@key,"String|Integer",MatchQ[key,_String|_Integer]/.trueFalseRule},
{Style[#,Bold,14]&@"tNow",tNow,Head@tNow,"?NonNegative",MatchQ[tNow,_?NonNegative]/.trueFalseRule},
{Style[#,Bold,14]&@"weekday",weekday,Head@weekday,"?BooleanQ",MatchQ[weekday,_?BooleanQ]/.trueFalseRule},
{Style[#,Bold,14]&@"tEvent",tEvent,Head@tEvent,"?NonNegative",MatchQ[tEvent,_?NonNegative]/.trueFalseRule},
{Style[#,Bold,14]&@"iEvent",iEvent,Head@iEvent,"String",MatchQ[iEvent,_String]/.trueFalseRule}
},TableHeadings->{Automatic,Style[#,GrayLevel[.3],Bold,14]&/@{"Variable","Value","Head","Condition","Test"}}
];
Abort[](* End processing if reportIteration fails *)
];

(* Message::wrongNumberOfArguments: *)
updatePersonNowEvent[args___]:=(
Message[updatePersonNowEvent::wrongNumberOfArguments,wrongNumberOfArguments];
Print[Style[#,Bold,14]&@"Arguments: ",Short[args,12]];
Abort[](* End processing if fails *)
);

Running the ABM

Main Program:

Table[
weekday=weekdayQ@tNow;

placeToPersons=GroupBy[Normal[locState[[All,2]]],Last->First];
timingContagionModel=First@Timing[(contagionModel[#1,nContactsMax=8,prStoE=.5,{tNow,weekday,_}]&)/@KeyValueMap[List,placeToPersons]];


(key=#1;updatePersonNowEvent[key,{tNow,weekday,_},#]&/@Cases[evQ@key,_[_?((#<tNow)&),__](*_[t_/;t<tNow,__]*)])&/@Keys[person];

If[Mod[tNow,printEveryNHours]<0.01,
Print[Style["Day "<>ToString[Ceiling[tNow/24]],12,Bold,"Courier"]," ",Style["tNow "<>ToString[tNow],12,Bold,"Courier"]," ",
DateObject[DatePlus[DateObject[dateStart,"Hour"],Quantity[tNow,"Hours"]],"Day"]," ",DateObject[DatePlus[DateObject[dateStart,"Hour"],Quantity[tNow,"Hours"]],"Week"]]
](*print simulation progress*);

logPersonsCovid[](*Log person disease states*)

,{tNow,tInitial=1,tFinal=simDays*24,dT=1}]

ABM Results

evQ
evH
state

pkey=randomKey[person]
showPerson[pkey]
showPerson[pkey,"State"]

Conclusions

Enormous benefits from incorporating two modeling paradigms together in one modeling system (Mathematica notebook)

Challenge:
Simulating 2.7M+ persons, moving hourly between 1.3M places and  interacting at those places for a period of 4 or more weeks
Computation Time An Issue: Esimtaed compute time is 15 hours.