A SEIRD Model For COVID-19 Using DDEs
Introduction
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.
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)
General Procedure (Germany)
Getting Data
Getting Data
Germany lockdown measures on Mar 20 day 54
Population = 83 millions
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*)
(*lockdowndate*){{day0}}=Position[cdates,DateObject["20Mar2020"]]
{{54}}
ndays=Length@fitData[[1]];imax=Max@fitData[[1,All,2]];pop=83000000.;
SEIRD Model
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 and are used for the recoveries and deaths, respectively. As usual, given the short period of the pandemic, it is assumed that the births and natural deaths have a null net effect and the total population is constant.
τ
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:
′ βℐ[t][t] |
′ ℰ βℐ[t][t] |
′ ℐ |
′ ℛ |
′ |
refresh:=(eqns={[t]-β[t]ℐ[t]/,[t]-ζℰ[t]+β[t]ℐ[t]/,[t]ζℰ[t]-γℐ[t-τr],[t]γℐ[t-τr],[t]μℐ[t-τd]};iconds:={[t/;t≤0]==,ℰ[t/;t≤0]==0,ℐ[t/;t≤0]==1,ℛ[t/;t≤0]==0,[t/;t≤0]==0};{susceptibleModel,infectedModel,exposedModel,recoveredModel,deadModel}={,ℐ,ℰ,ℛ,}/.ParametricNDSolve[Join[eqns,iconds],{,ℐ,ℰ,ℛ,},{t,0,200},{β,ζ,γ,τr,μ,τd}];)
′
′
ℰ
′
ℐ
′
ℛ
′
Manual Fit
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)
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×,{iβ0.873186,iζ0.0682186,iγ0.0624517,iτr4.74667}}}
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
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[]
Results
Results
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
Fatality Rate Evolution
Showing Fits for 6 Other Studied Countries
Showing Fits for 6 Other Studied Countries
Austria
Austria
Detailed view of deadModel fit:
Iceland
Iceland
Detailed view of deadModel fit:
New Zealand
New Zealand
Detailed view of deadModel fit:
Switzerland
Switzerland
Detailed view of deadModel fit:
Israel
Israel
South Korea
South Korea
Summary
Summary
Data Table
Data Table
Linear Correlation Matrix
Linear Correlation Matrix
Discussion
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.
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
Initialization
Useful Functions
Useful Functions
Plotting Functions
Plotting Functions
Manual fit - SEIRD with 2 Delays (τ_r & τ_d)
Manual fit - SEIRD with 2 Delays (τ_r & τ_d)
This routine is a slightly modified version of the code published by Robert Nachbar
Fetch Data From Johns Hopkins University
Fetch Data From Johns Hopkins University
Import COVID19 Data
Import COVID19 Data