# Message Boards

Answer
(Unmark)

GROUPS:

4

# A SEIRD Model For COVID-19 Using DDEs

Posted 1 year ago

A SEIRD Model For COVID-19 Using DDEs
Introduction I use as starting point the post from Robert Nachbart about Epidemiological Models and I follow his notation to some extent. It has an excellent, pedagogical introduction to the subject and a great read for non-experts like myself. There, it is shown that the SEIR models, using ordinary differential equations (ODEs), capture the qualitative behavior of the COVID-19 pandemic but fail to provide good fits to actual data, even when considering more realistic models that take into account important features like the large number of asymptomatic cases and a quarantine compartment is added. In this post the use of Delayed Differential Equations (DDEs) is considered instead. For that purpose, I chose data from 7 countries where the epidemic has already past its peak and the evolution (at first glance) seems uniform, i.e., without apparent new outbreaks. The model and procedures are first presented showing data from Germany, followed by summary results for Austria, Iceland, Israel, New Zealand, South Korea, and Switzerland. The main drawback of using DDEs is that they are prone to oscillations and overshooting. This is particularly problematic for the fitting/minimizing algorithms, that like to wander freely the parameter space in search for an optimal solution. Therefore, they are not suited like the ODEs models to build up complexity: a few free parameters proved already challenging to handle. Having said that, they appear to produce reasonably good fits for the few cases considered here, although some model parameters seem to take unrealistic values. For a discussion of the results look for the Discussion at the end. For this presentation purposes, I moved the Initialization section to the bottom of the notebook.
General Procedure (Germany)
Getting Data Germany lockdown measures on Mar 20 day 54 Population = 83 millions country="Germany";region="";{confirmedData,deadData,recoveredData}=getData2[country,region];(*interpolatingrecovereddata*)recoveredData=interF@recoveredData;infectedData=confirmedData-deadData-recoveredData;countryData={infectedData,recoveredData,deadData};fitData=addTime/@countryData;cdates; Last@cdates(*latestupdate*) Day: Sun 24 May 2020 (*lockdowndate*){{day0}}=Position[cdates,DateObject["20Mar2020"]] {{54}} ndays=Length@fitData[[1]];imax=Max@fitData[[1,All,2]];pop=83000000.;
SEIRD Model The model uses 5 compartments: the susceptible , exposed ℰ, infectious ℐ, recovered ℛ, and dead populations. The transition rates between the compartments are the transmission rate β, the inverse of the latent (incubation) period ζ, the recovery rate γ, and the mortality rate of infectious individuals μ. Time delays τ r τ d SEIRD Model with 2 delays can be simplified neglecting the dead removed from the infectious population in the 3rd Eq., taking into account that always γ >> μ, that greatly improves the equation solver stability:
refresh:=(eqns={ ′ ′ ℰ ′ ℐ ′ ℛ ′
Manual Fit First we must explore manually the parameter space looking for a good fit, since to use e.g. NMinimize available procedures, we need to find relatively good, reasonable values to start with. {βo,ζo,γo,τro,μo,τdo}={1.456,0.11,0.0566,6.56,0.0024,3.6}; =170000.; rg=2;(*controler'srange*)refresh;manualFit
idata(*savedparameters{βo,ζo,γo,τro,μo,τdo}*){1/βo,1/ζo,1/γo,1/μo} {0.763836,0.0803526,0.064121,4.3,0.00328095,3.6} {1.30918,12.4451,15.5955,304.789}
NMinimize Sum of Squares Estimates (SSE) Since γ>>μ, we can fit the dead population data separately and minimize first the sum of squares estimates (SSE) for the infected and recovered models. Note that only the data after the lockdown is taken into account. ref0:=(Mean@Table[(recoveredData[[t]]-recoveredModel[iβ,iζ,iγ,iτr,μo,τdo][t])^2+(infectedData[[t]]-infectedModel[iβ,iζ,iγ,iτr,μo,τdo][t])^2,{t,day0+1,ndays}]) {βo,ζo,γo,τro,μo,τdo}={0.76,0.08,0.064,4.3,0.0033,3.6}; =168000;refresh; Timing[minsol=NMinimize[ref0,{{iβ,0.5,1.0},{iζ,0.06,0.09},{iγ,0.05,0.08},{iτr,3.,7.}},Method"NelderMead"]]Beep[] {94.8005,{6.4148× 6 10 {βo,ζo,γo,τro}={iβ,iζ,iγ,iτr}/.minsol[[2]] {0.873186,0.0682186,0.0624517,4.74667} 2nd step: the same for the dead population model. ref1:=(Mean@Table[(deadData[[t]]-deadModel[βo,ζo,γo,τro,iμ,iτd][t])^2,{t,day0+1,ndays}]) =168000;refresh; Timing[minsol2=NMinimize[ref1,{{iμ,0.002,0.004},{iτd,2.,4.}},Method"NelderMead"]]Beep[] {3.37586,{7135.3,{iμ0.00332575,iτd5.25082}}} {μo,τdo}={iμ,iτd}/.minsol2[[2]] {0.00332575,5.25082}
Search for optimum population size Since the population size parameter is part of the initial conditions of the DDEs, one cannot use NMinimize directly to obtain the best fit value (my guess is that it tries to compile the SSE quality function that needs to be “refreshed” for each different value of ). Since this is not critical, I just brute-force a simple search: {μo,τdo}={0.00332575,5.25082}; solset={};o=173000.;Timing[Do[=o+(i-1)1000.;refresh;minsol=NMinimize[ref0,{{iβ,0.6,1.2},{iζ,0.06,0.09},{iγ,0.05,0.08},{iτr,3.,7.}},Method"NelderMead"];Print["i = ",i," = ",];Print[minsol];AppendTo[solset,Flatten@{,minsol}];,{i,1,5}]]Beep[] i = 1 = 173000. {4.99283× 6 10 i = 2 = 174000. {4.96759× 6 10 i = 3 = 175000. {5.01199× 6 10 =solset[[2,1]]{βo,ζo,γo,τro}={iβ,iζ,iγ,iτr}/.solset[[2,3;;6]] 174000. {0.987554,0.0584478,0.0632589,5.24472} Final search for μ and τ d =174000;refresh; Timing[minsol2=NMinimize[ref1,{{iμ,0.002,0.004},{iτd,2.,4.}},Method"NelderMead"]]Beep[] {2.26067,{8960.1,{iμ0.00331479,iτd5.27758}}} {μo,τdo}={iμ,iτd}/.minsol2[[2]] {0.00331479,5.27758}
Results (*BestFitParameters*)=174000;{βo,ζo,γo,τro,μo,τdo}={iβ,iζ,iγ,iτr,iμ,iτd}/.{iβ0.987,iζ0.0584,iγ0.063,iτr5.24,iμ0.0033,iτd5.27};
Using another scale to better appreciate the dead population model:
We modeled only the behavior under the strict quarantine imposed after lockdown (day 54). Note that the initial data points are poorly fitted, in part also because the chosen quality function (SSE) puts more weight on larger values.
Fatality Rate Evolution (*peakofinfectiouscases*)msol=NMaximize[{infectedModel[βo,ζo,γo,τro,μo,τdo][t],1<t<150.},t] {67270.,{t72.123}} (*atthepeak&limitt>>1*)tpk=t/.msol[[2]];{frpk,frlim}=(100.deadModel[βo,ζo,γo,τro,μo,τdo][t]/(recoveredModel[βo,ζo,γo,τro,μo,τdo][t]+infectedModel[βo,ζo,γo,τro,μo,τdo][t]+deadModel[βo,ζo,γo,τro,μo,τdo][t])/.t#)&/@{tpk,200} {1.88843,4.97511}
Basic Reproduction Number R 0 One of the many drawbacks of using DDEs is that the analysis of its dynamics is not straightforward. On the other hand, although it poorly fits the data, from the standard SEIR model is easy to obtain an estimate for R 0 (*peakofinfectedcases*)ipeak=msol[[1]]{ipeak/,ipeak/pop} 67270. {0.386609,0.00112117} NSolve[(1+Log[x])x(1-ipeak/)&&x>1,x]//Quiet {{x3.81174}} NSolve[(1+Log[x])x(1-ipeak/pop)&&x>1,x]//Quiet {{x1.04929}} One might interpret the effective R 0
Showing Fits for 6 Other Studied Countries
Austria =15050;{βo,ζo,γo,τro,μo,τdo}={iβ,iζ,iγ,iτr,iμ,iτd}/.{iβ1.52201,iζ0.105816,iγ0.0565933,iτr6.44976,iμ0.00242167,iτd3.52964};
Detailed view of deadModel fit:
Iceland =1800;{βo,ζo,γo,τro,μo}={iβ,iζ,iγ,iτr,iμ}/.{iβ1.24757,iζ0.102812,iγ0.0642217,iτr7.57979,iμ0.000363};
Detailed view of deadModel fit:
New Zealand =1460;{βo,ζo,γo,τro,μo,τdo}={iβ,iζ,iγ,iτr,iμ,iτd}/.{iβ0.45527,iζ0.270821,iγ0.0613337,iτr5.49906,iμ0.00094351,iτd8.07569};
Detailed view of deadModel fit:
Switzerland =29700;{βo,ζo,γo,τro,μo,τdo}={iβ,iζ,iγ,iτr,iμ,iτd}/.{iβ2.741,iζ0.06468,iγ0.0613,iτr6.207,iμ0.00416,iτd1.715};
Detailed view of deadModel fit:
Israel =17100;{βo,ζo,γo,τro,μo,τdo}={iβ,iζ,iγ,iτr,iμ,iτd}/.{iβ1.24671,iζ0.0647559,iγ0.0399033,iτr9.42441,iμ0.0007575,iτd0.};
South Korea =10200;{βo,ζo,γo,τro,μo,τdo}={iβ,iζ,iγ,iτr,iμ,iτd}/.{iβ0.659621,iζ0.13471,iγ0.0325641,iτr8.60056,iμ0.000859913,iτd3.22902};
Summary
Data Table Description: pop is the total population; , β , ζ , γ , τ r τ d R 0 R e
Linear Correlation Matrix
Discussion The first thing to notice is that best-fits for the presented model are obtained taking an initial susceptible population that is considerably smaller than the country total population. This seems at odds with the fact that this is a novel virus and therefore everyone is susceptible and likely to spread it, even if having mild or no symptoms at all. Nevertheless, one of the main assumptions of these models, that were borrowed from the analysis of chemical reactions, is that there should be an homogenous mixing of the reactants, or equivalently between the infectious and the susceptible populations. This can be hardly the case under the lockdown measures imposed in all these countries. Thus, one possible interpretation is that the limited mixing results in an effective population of reduced size and a fast β rate of transmission within it. The obtained incubation periods (i.e, time between infection and being able to transmit it), given by 1/ζ, show a large spread, and in most cases much longer than a mean of about 5 days reported by WHO based on China patients. On the other hand, recovery times, given by 1/γ, are more similar and consistent with reports of about 2-3 weeks. The impact of the disease wildly varies between countries. There are many reasons to expect a large spread in parameter values, a priori that should not be a concern, since there are no uniform international criteria for testing, time taken to report recoveries, and not even dead counts. Differences in age distribution, health care systems, fraction of the population with comorbidities (e.g., diabetes, hypertension), etc. To conclude, although this model can be taken simply as a tool to fit the data, it would need to be refined to obtain more realistic time scales (particularly the incubation time). However, it is hard to build up on it due to stability problems. Also, I find unsatisfying that the initial population almost acts as a free parameter, representing very different fraction of the total population in each case.
Initialization
Useful Functions addTime[set_]:=Transpose[{Range[Length@set],set}] fdate[x_]:={2000+#[[3]],#[[1]],#[[2]]}&@ReadList[StringToStream[x],Number,RecordSeparators{"/"}] interF[set_]:=Module[{f,n,aux},n=Length@set;aux=addTime@set;f=Interpolation[First/@Gather[aux,#1[[2]]#2[[2]]&],InterpolationOrder1];Round@f@Range[n]]
Plotting Functions pstyles4={ColorData["Crayola"]["Green"],ColorData["Crayola"]["Sunglow"],ColorData["Crayola"]["Scarlet"],ColorData["Crayola"]["NavyBlue"]};styles3={ColorData["Crayola"]["Scarlet"],ColorDat |