MATH564 Homework 7: Infection models

Topics

Screening clinics

Screening clinics are an example of a compartment model

eqs = {A'[t] == 1 - .08 A[t] + 0.01 B[t],B'[t]==2+0.02 A[t]-0.02 B[t]};
icnd = {A[0]==0,B[0]==0};

sol[t_]=Simplify[{A[t],B[t]} /. DSolve[Flatten[{eqs,icnd}],{A[t],B[t]},t][[1]]]

Though this system is linear and should reach a steady state, the above expression and plots below indicate divergence. Explain this behavior. 3 course bonus points! Due date: Oct 26 before class (1:00PM) (email is fine)

Plot[sol[t],{t,0,100}]

SIR model in time

Basic model

SIReqs={S'[t] == -β F[t] S[t], F'[t] == β F[t] S[t] - γ F[t], R'[t] == γ F[t]}

SIeqs = SIReqs[[1;;2]]

SIcnd={S[0]==(10^6),F[0]==1}

DSolve[SIeqs,{S[t],F[t]},t]

sol[t_] = {S[t],F[t]} /.NDSolve[Flatten[{Evaluate[SIeqs /. {β->0.001, γ->0.98}],SIcnd}],{S[t],F[t]},{t,0,10000}][[1]]

Plot[sol[t],{t,0,0.1},PlotRange->All]

Model parameters from data

Get COVID-19 infection data

NCdata=Import["https://covidtracking.com/data/download/north-carolina-history.csv"];

Dimensions[NCdata]

The first row of the spreadsheet data contains data categories

categories=NCdata[[1,;;]]

Dates are given as year-mm-dd, e.g.

NCdata[[2,1]]

Transform into days, starting from earliest available date

t0=AbsoluteTime[NCdata[[-1,1]]]/86400

tNC=Map[(AbsoluteTime[#]/86400-t0)&,NCdata[[2;;,1]]];

Find positions of data entries for infections,  deaths, recovered

iI = Position[categories,"positive"][[1,1]];
iD = Position[categories,"death"][[1,1]];
iR= Position[categories,"recovered"][[1,1]];

Define infectious and recovered populations

INC = NCdata[[2;;,iI]];
RNC = NCdata[[2;;,iD]]+NCdata[[2;;,iR]];

Form time history data

tINC=Reverse[Transpose[{tNC,INC}]];
tRNC=Reverse[Transpose[{tNC,RNC}]];
ListPlot[{tINC,tRNC},Frame->True,GridLines->Automatic,PlotLegends->{"I","R"}]

The initial susceptible population is that of the entire state

S0NC=10.49 (10^6)

Form the susceptible population time series. Plot using decimal logarithm ordinates due to the large scale disparity.

tSNC=Reverse[Transpose[{tNC,S0NC-INC}]];
ListLogPlot[{tINC,tRNC,tSNC},Frame->True,GridLines->Automatic,PlotLegends->{"I","R","S"}]

Define initial conditions for NC

icndNC = {S[0]==S0NC, F[0]== INC[[-1]]}

Define the S-I mode for NC

SINCmodel = Flatten[{SIeqs,icndNC}]

Define a function that returns the S population for some assumed β,γ values

SI[βv_,γv_,tv_]:=Module[{sol,val,SINCmodelv},
SINCmodelv = Evaluate[SINCmodel /. {β->βv,γ->γv}];
sol = {S[t],F[t]} /. NDSolve[SINCmodelv,{S[t],F[t]},{t,0,Max[tv]}][[1]];
val = Table[ sol /. t->tv[[k]],{k,1,Length[tv]}];
Return[Transpose[{tv,val[[;;,1]]}]];
]

Define a list of the times at which we’ll compare the predictions of the above model with the data

tNCv = Reverse[tNC];

Define the error between model prediction and data

ModelError[p_]:=Norm[ SI[p[[1]],p[[2]],tNCv][[;;,2]] - tSNC[[;;,2]]]

Define the gradient of the error

gradModelError[p_,h_,k_] := {ModelError[{p[[1]]+h,p[[2]]}]-ModelError[{p[[1]]-h,p[[2]]}],
ModelError[{p[[1]],p[[2]]+k}]-ModelError[{p[[1]],p[[2]]-k}]}/{2h,2k}

param0={4. (10.^(-9)),(10.^(-12))};
ModelError[param0]

Define a gradient descent iteration

GradDescent[p_,λ_] := {p[[1]],p[[2]]} - λ gradModelError[p,(10.^(-12)),(10.^(-15))]

Determine the scale of the λ step size

{param=GradDescent[param0,(10.^(-24))],ModelError[param]}

First gradient iteration: plot the descent of the function for various λ values

errp=Table[{λ ,{param=GradDescent[param0,λ (10.^(-24))],ModelError[param]}[[2]]},{λ,1,20}];
ListPlot[errp]

Adopt the step size corresponding to the local minimum

param1=GradDescent[param0,13 (10.^(-24))]

Repeat the above: iteration 2

errp=Table[{λ ,{param=GradDescent[param1,λ (10.^(-26))],ModelError[param]}[[2]]},{λ,1,20}];
ListPlot[errp]

param2=GradDescent[param1,12 (10.^(-26))]

Repeat the above: iteration 3

errp=Table[{λ ,{param=GradDescent[param2,λ (10.^(-26))],ModelError[param]}[[2]]},{λ,1,20}];
ListPlot[errp]

param3=GradDescent[param2,14 (10.^(-26))]

Repeat the above: iteration 4

errp=Table[{λ ,{param=GradDescent[param3,λ (10.^(-27))],ModelError[param]}[[2]]},{λ,1,20}];
ListPlot[errp]

No further descent is observed. Adopt the best estimate.

βBest=param3[[1]]; γBest = param3[[2]];
tSbest=SI[βBest,γBest,tNCv];

Compare with data

ListPlot[{tSbest,tSNC},Frame->True,FrameLabel->{"day","S"},PlotLegends->{"model","data"}]

A reasonable approximation is obtained, but clearly the constant parameter SI model does not capture the real data. The above gradient descent procedure can be replaced by a built-in function. First redefine the error to only accept numeric arguments, and also scale the parameters and error to be of order one.

NModelError[βs_ ?NumericQ,γs_ ?NumericQ]:=Norm[ SI[βs (10^(-9)),γs (10^(-12)),tNCv][[;;,2]] - tSNC[[;;,2]]]/(10.^6)

NModelError[4.,1.]

bestfit=FindMinimum[NModelError[βs,γs],{βs,4},{γs,1.}]

Define the best fit and replot

βBestFit=((10 )^(-9))(βs/.bestfit[[2]]); γBestFit =(10^(-12)) (γs/.bestfit[[2]]);
tSBestFit=SI[βBestFit,γBestFit,tNCv];
ListPlot[{tSbest,tSBestFit,tSNC},Frame->True,FrameLabel->{"day","S"},PlotLegends->{"model","best","data"}]

The same conclusion is obtained: the constant parameter does not conform to the measured data

SIR model in space and time

Geographical, population data

Get a list of counties within NC

c=WolframAlpha["Counties in North Carolina",{{"BasicInformation:USCountyData",1},"ComputableData"},PodStates->{"Result__More"}][[1,1]];

counties=AdministrativeDivisionData[LinguisticAssistant, "Subdivisions"];

bordering=EntityValue[counties,"BorderingCounties"];

Length[missingpos=Position[bordering,_Missing]]

counties=Delete[counties,missingpos];
bordering=Delete[bordering,missingpos];

Length[edges=DeleteDuplicates[Sort/@Flatten[MapThread[Thread[UndirectedEdge[##]]&,{counties,bordering}]]]]

Graph[counties,edges,VertexLabels->None]

EntityValue[counties[[1]],"Position"]

EntityValue[counties[[1]],"Population"]

pos=EntityValue[counties,"Position"];
GeoListPlot[pos]

pop=EntityValue[counties,"Population"];

GeoDensityPlot[pos->pop]

Define weighted graph SIR model

Define the weights as proportional to the population  product

α=(10^(-9)); popM=N[QuantityMagnitude[pop]]; nc = Length[pop];

Define the dependent variables: susceptible and infectious populations in each county

Sc=Table[S[n,t],{n,nc}]; Fc=Table[F[n,t],{n,nc}];

Introduce indexing of counties

cidx = Table[ counties[[n]]->n , {n,nc}];

Define indexing for edges

edgidx = Table[edges[[k]] /. cidx , {k,Length[counties]}]

Quarantine NC from adjacent states

NCedges = Select[ edgidx, ( NumericQ[#[[1]]] && NumericQ[#[[2]]] )& ]

Define a weight matrix based on population weights. Start with zero weights

w = Table[0.,{i,nc},{j,nc}];

Assign weights by population product

For[edg=1, edg≤ Length[NCedges], edg++,
  i=NCedges[[edg,1]]; j= NCedges[[edg,2]];
  popi = popM[[i]]; popj = popM[[j]];
w[[i,j]] = popi*popj; w[[j,i]]=w[[i,j]]];

Define the graph Laplacian

Δ[f_,i_]:= Sum[ w[[i,j]](f[[j]]-f[[i]]) , {j,nc}]

Δ[Sc,1]

α=(10.^(-8)); β = 5.33 (10^(-9)); δ=10.α; γ=(10.^(-12));
Seqs = Table[D[S[n,t],t] == -β S[n,t] (F[n,t] + α Δ[Fc,n]) + δ Δ[Sc,n],{n,nc}];
Feqs =  Table[D[F[n,t],t] == β S[n,t] (F[n,t] + α Δ[Fc,n]) -γ F[n,t],{n,nc}];

Sini = Table[ S[n,0] == popM[[n]],{n,nc}];
Fini = Table[ F[n,0] == 0.,{n,nc}];
Fini[[5]] = F[5,0]==100.;
Fini[[11]] = F[11,0]==400.;

ODEs =Flatten[ {Seqs,Feqs,Sini,Fini}];
vars = Flatten[{Sc,Fc}];

sol = NDSolve[ODEs,vars,{t,0.,270.}];
SNC[t_] = Table[S[n,t],{n,nc}] /. sol[[1]];
FNC[t_] = Table[F[n,t],{n,nc}] /. sol[[1]];

Construct a table of plots at 14 day intervals

plts=Table[GeoDensityPlot[pos->FNC[14. k]],{k,0,9}];

Display the space time evolution of the infectious population

GraphicsGrid[ ArrayReshape[plts,{5,2}]]

Problems

Problem 1 (2 points)

The above best fit determination of the infection rate β used data from the entire time history, leading to the conclusion that a constant β value does not conform to the data. Determine β over subintervals of 30 and 60 days, and study whether a variable β model better conforms to the data. Extra credit: Determine β over 7 day intervals, and interpolate this data to obtain β(t). Use the interpolation for β(t) to numerically solve the SI model and compare to measured data.

Problem 2 (3 points)

Choose a reduced set of 10 counties. Determine the evolution of S,I populations in the county subset. Quarantine one of the counties. Compare the evolution when one county is in quarantine with the evolution with no quarantine.  Extra credit: Modify the space time model to include variable in time infectious rates, β(t).