Message Boards

GROUPS:

Staff Picks Biological Sciences Data Science Medical Sciences Social Science Graphics and Visualization Wolfram Language Modeling Events & Media

Daniel Lichtblau

[NB] A simple COVID-19 spread model

Daniel Lichtblau, Wolfram Research

Posted 2 years ago

4641 Views

3 Replies

15 Total Likes

Follow this post

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

A simple COVID-19 spread model

Introduction

Over the past two to three weeks I’ve been working on some simple modeling of COVID-19 spread, using assumptions that seem appropriate for localized regions with shelter-in-place restrictions. A first stab was reported in this Wolfram Community post (I will refer to this as “the prior post” in a few places below). I have since refined some of the code, but the basic idea is this. After the initial set of infections is noted in any given city, a partial quarantine takes place with social distancing and shelter-in-place or related restrictions taking effect. In such circumstances I believe the SEIR-type models are not going to be terribly useful because they use mass-action kinetics to give spread rates. This simply does not apply when the population mingling is minimal.

As in the SEIR family of models, we denote those currently infected, as a function of time, by

i(t)

. Those who have recovered or died are denoted by

r(t)

. The cumulative number ever infected includes both of these groups, and I will refer to that function as

p(t)

(the

being for “ever tested positive”). What I noticed is that, based on reporting of US cases at the county level, the quotient of new COVID-19-positive cases to current ones serves as a good surrogate for the

from the epidemiology literature (the average numbe of people infected by an infected individual). In math terminology, the quotient I utilize is the logarithmic derivative of

i(t)

, that is, it is

i'(t)/i(t)

It turns out I was not the only person at the time to focus on such derivatives for analyzing COVID-19. At about the same time, Niayesh Afshordi, a physicist at the University of Waterloo and the Perimeter Institute for Theoretical Physics, wrote a nice blog with the title "Physical Measure of a Pandemic". He had the useful idea that while we do not know actual numbers of infected individuals, we do have hard numbers for deaths. (Actually this has been called into question recently. But we do know that the difference is likely to be well under a factor of two.) As do not have any such knowledge for infected-vs.-reported, mortalities can give accurate information in ways that infected counts cannot. I will say a little more about this later. (Another difference between his work and mine is that he does modelling whereas I do modeling. Canadians. What else need be said?) We (and a couple of colleagues) recently met up via Zoom. So I know yet another difference is that he is developing more extensive models than what I present here, treating both location and time as variables. (This is ambitious and I cannot do anything like it, but he can. His field is General Relativity-- space-time is their playground.)

Back to

thatlogarithmicderivativei'/i

. We can use this quotient, for example, to estimate the effective

value, if we have some data on average length of infectivity, average percentage of asymptomatic cases, etc. Even in the absence of such data, a number of plots across many US cities indicated that there is a narrow band for the "critical value" of this quotient. This value, for a given city, is the crossing point between rising and falling numbers of new cases. It seems that when it is under 0.06 the value will be falling, and when it exceeds 0.1 or so the new cases will be rising. The point is that one can gauge how a city is faring by charting a value that is very simple to compute, and moreover free from model assumptions. I will remark that for

the critical value is 1. When infected people themselves infect on average more than one other person, there is exponential growth in infection counts. When they infect on average less than one person it becomes exponential decline. (When it is exactly one there is linear growth. But in practice it is never exactly one.)

While this might be of some use it is not really enough to help predict the near-term case load, and of course that is a very important piece of information for health care facilities and professionals within any given municipal area. In this notebook I show a small refinement to the model and then use it to provide scenarios over the next few weeks. These look, to me, to be mostly optimistic. Which makes me pessimistic as to whether the model is really capturing the dynamics correctly. The nice thing is that the next few weeks will provide that information. If validated by upcoming data, this then can become a useful tool for future outbreaks. If it should be invalidated, well, back to school for its author (loud sob).

Acknowledgements

I thank my Wolfram colleague Mohammad Bahrami for a code modification that provided a big improvement over my clumsy attempt to portray infection curves. I also thank Robert Nachbar and Michael Trott for several helpful suggestions.

Obtaining data

Get case data

In[]:=

AbsoluteTiming[timeseries=ResourceFunction["NYTimesCOVID19Data"]["USCountiesTimeSeries"];]

Out[]=

{38.8901,Null}

Here I use number of major US cities. I found the counties and located their various positions in the time series dataset using a rather pedestrian approach detailed in the prior post. Note that new reports (reports from counties that have not previously had cases) get appended at the end, so the numbers below do not change over time.

In[]:=

sanfran=timeseries[[8]];chicago=timeseries[[2]];la=timeseries[[5]];nyc=timeseries[26];dallas=timeseries[159];miami=timeseries[[170]];orlando=timeseries[[266]];detroit=timeseries[[194]];seattle=timeseries[[22]];houston=timeseries[[49]];nola=timeseries[[146]];atlanta=timeseries[[30]];bergenC=timeseries[[36]];essexC=timeseries[[230]];hudsonC=timeseries[[125]];philly=timeseries[[153]];phoenix=timeseries[[4]];ftlau=timeseries[[55]];memphis=timeseries[[112]];denver=timeseries[[52]];vegas=timeseries[[45]];baltimore=timeseries[[180]];boston=timeseries[[7]];

Recapturing existing data

Here is a way to compute an approximation of

p'(t)/p(t)

from existing data. It is essentially a simplified and arguably better version of code from the previous post.

In[]:=

covid19RatioFunction[ent_,t_]:=Module[{caseTS=ent["Cases"],totals,perday,quots},totals=caseTS[[2,1,1]];perday=Differences[totals];quots=perday/ListConvolve[{1,1}/2.,totals];Interpolation[Transpose[{Range[Length[quots]],quots}],t,InterpolationOrder3]]

Again I will remark that this is ignoring the effect of

r(t)

during the time span in question.

Worked examples that recover past 21 days of behavior using an ODE

We now set up an ordinary differential equation (ODE) to simulate the recent past data. The code was a bit fiddly and what I found to be generally reliable was to set the initial value as an average of days 0 and 1; this is largely due to how I estimated the derivatives, as difference quotients, in covid19RatioFunction above.

In[]:=

city=nyc;ratiocity=covid19RatioFunction[city,t];upper=Length[city["Cases"][[2,1,1]]]-1;gap=21;lower=upper-gap;ndsoln=NDSolveValue[{p'[t]ratiocity*p[t],p[lower]Mean@city["Cases"][[2,1,1,lower;;lower+1]]},p[t],{t,lower,upper}];Show[{Plot[ndsoln,{t,lower,upper}],ListPlot[With[{vals=city["Cases"][[2,1,1,lower;;upper]]},Transpose[{Range[lower,upper],vals}]]]},PlotRangeAll]

Out[]=

In[]:=

city=sanfran;ratiocity=covid19RatioFunction[city,t];upper=Length[city["Cases"][[2,1,1]]]-1;lower=upper-gap;ndsoln=NDSolveValue[{i'[t]ratiocity*i[t],i[lower]Mean@city["Cases"][[2,1,1,lower;;lower+1]]},i[t],{t,lower,upper}];Show[{Plot[ndsoln,{t,lower,upper}],ListPlot[With[{vals=city["Cases"][[2,1,1,lower;;upper]]},Transpose[{Range[lower,upper],vals}]]]},PlotRangeAll]

Out[]=

Estimating actual cases

Before proceeding to the forecast model I want to indicate some crude ways to estimate actual cases given only numbers of those who have tested positive.

Method 1: Guess actual case count to be five times higher than stated on a particular “initial” date.

This is actually useless, but since it is likely to be the first thing that comes to mind it should be mentioned. It is silly because it simply bakes the result right into the guess. That is to say, since we have a linear ODE, five times higher on day 0 means five times higher today.

In[]:=

city=nyc;ratiocity=covid19RatioFunction[city,t];upper=Length[city["Cases"][[2,1,1]]]-1;gap=30;lower=upper-gap;ndsoln=NDSolveValue[{p'[t]ratiocity*p[t],p[lower]5Mean@city["Cases"][[2,1,1,lower;;lower+1]]},p[t],{t,lower,upper}];Show[{Plot[ndsoln,{t,lower,upper}],ListPlot[With[{vals=city["Cases"][[2,1,1,lower;;upper]]},Transpose[{Range[lower,upper],vals}]]]},PlotRangeAll]

Out[]=

Method 2: We guess by how many days prior we had actually hit a given count

Now at least we can get realistic scenarios based on assumption that the

p'/p

vs time curve approximately hold not just for tested-positives but for all positives. First I show the curve if we assume that the count from the initial day (or the average of the first two days, per notes above) was actually hit three-and-a-half days prior.

In[]:=

back=3.5;ndsoln=NDSolveValue[{p'[t]ratiocity*p[t],p[lower-back]Mean@city["Cases"][[2,1,1,lower;;lower+1]]},p[t],{t,lower-back,upper}];Show[{Plot[ndsoln,{t,lower-back,upper}],ListPlot[With[{vals=city["Cases"][[2,1,1,lower;;upper]]},Transpose[{Range[lower,upper],vals}]]]},PlotRangeAll]

Out[]=

Now suppose it was six days prior.

In[]:=

back=6;ndsoln=NDSolveValue[{p'[t]ratiocity*p[t],p[lower-back]Mean@city["Cases"][[2,1,1,lower;;lower+1]]},p[t],{t,lower-back,upper}];Show[{Plot[ndsoln,{t,lower-back,upper}],ListPlot[With[{vals=city["Cases"][[2,1,1,lower;;upper]]},Transpose[{Range[Length[vals]]+lower-1,vals}]]]},PlotRangeAll]

Out[]=

A problem with this method, which I will describe but not show, is that it is quite sensitive to how many days we use to model prior history (this is the value of the variable gap used above, not the variable back). I believe this is to some extent from the fact that earlier data involves smaller numbers and hence tends to be more “jumpy” relative to later. But this is only a conjecture.

A third way to gauge actual cases might be as follows. Take the existing mortality data for a given city and use the range of factors for mortalities computed in places that have done more accurate testing of infected individuals e.g. South Korea. This will still be imperfect, since mortality rates are not known exactly and also can vary by location (based on overall health and age factors for any given city). But it might be a useful way to do the estimation.

Modeling near-term scenarios

We only have data for reported cases, and there are estimates suggesting that the number of cases not reported might be many times larger. As in the previous post, I make the assumption that the day-to-day ratios for the unreported cases are similar to those for the reported cases, in any given location.

I mistakenly believed that the data showed currently infected patients in any given county. Closer reading of the New York Times’ Github repository page indicates that what is actually reported is the total that have ever tested positive in a given county, that is,

p(t)

rather than

i(t)

. So what I was approximating as a ratio was not

i'(t)/i(t)

but rather

i'(t)/p(t)

. While it is unfortunate for the patients and hospitals, these do not differ by much because to date the numbers that have recovered are small compared to the numbers that tested positive in any given county. So I use these values, charted day by day, as an approximation to

i'(t)/i(t)

Going forward we do need to account for the recoveries. In particular, we have the equality

i'(t)=p'(t)-r'(t)

. I will use this in the model. We will proceed as follows. We assume some average value, call it

, for how many days, after testing positive, it takes to recover. For most examples I will take this to be 21. We then go back that many days in the current data and solve a simple differential equation of the form

i'(t)=k(t)i(t)

where

k(t)

is just a function that interpolates those

i'/i

ratios estimated from existing data. In other words, we are forcing a differential equation to recapitulate existing infection counts (in “other-other” words, we are cheating). Since some approximations are used along the way this does not perfectly recapture existing data. But plots will indicate that it comes close in practice.

The point of doing this is that we now have a continuous function of time that closely approximates the past

days. We now use this data to model the upcoming

days. Specifically, it gives us

r(t)

and thus

r'(t)

. Readers knowledgeable about differential equations will recognize that I might instead model the system with delay differential equations, using

r(t)=i(t-b)

. But I will be using the NDSolve function in Mathematica, and I’ll just say that NDSolve does not like me very much (“barely tolerates me” would be more accurate). So I will go about this in a more pedestrian manner.

This is not the only use to be made of

r(t)

. People who may have read the prior post will have noted that, pretty much everywhere I showed, that important ratio function

k(t)

. So we really want to model that function as we move forward. We do not want to extrapolate from the prior data because, well, that simply makes too many assumptions that are difficult to justify (mostly because they are, invariable, false). A better way to estimate the ratio is to use values from a day prior as we move forward. Which means I actually do use a delay approach in the code. But in a way that NDSolve seems to find forgivable.

Forecasting based on past results

I start with a function that takes a location, an average

for recovery time lag, and a string name for the location. It also has an optional fourth argument

for capturing history only to

days ago. Its intended usage will be explained later. The first part simply encapsulates code from the examples that recover the past b days of data. The next part uses this, as mentioned in the description of the model, to move forward

days. I should mention that it is possible, at the expense of more complicated code, to go back more days than we go forward (and vice versa, if the delay differential solver is not throwing a tantrum).

The code for this second part is straightforward except for the part that approximates

k(t)

(that all-important ratio of newly infected to currently infected, day by day) as we move forward. The idea is to average values from recent days. This is done in order to “smooth” the behavior and thus suffer less from vagaries of single data points. For those values that are at times prior to the start, we use the result from the first computation of past-to-present behavior. This is sufficient daya for NDSolve to use as it integrates forward (readers unfamiliar with the workings of a numeric differential equation solver should just ignore this). One caveat is that this code uses several approximations. Another is that it is far enough from my expertise that it could well not be doing quite what I expect. But the plots seem to give reasonable ranges of future outlooks based on the average recovery time. This can be useful even though the model is imperfect and the implementation might also have flaws.

In[]:=

currentAndProjected[city_,gap_,name_String,back_:0]:=Module[{cases,ratiocity,upper,lower,ndsoln,p,r,firstpval,lastRatio,ndsolnNextPeriod,ndsolnprime,listplot,t,from,to,incr},ratiocity=covid19RatioFunction[city,t];cases=city["Cases"][[2,1,1]];upper=Length[cases]-1-back;lower=upper-gap;firstpval=Mean@cases[[lower;;lower+1]];{ndsoln,ndsolnprime}=NDSolveValue[{p'[t](ratiocity/.tt+lower)*p[t],p[0]firstpval},{p[t],p'[t]},{t,0,gap}];lastRatio=Mean[ratiocity/.tRange[upper-1,upper]];{from,to,incr}={9/2,1/4,1/10};ndsolnNextPeriod=NDSolveValue[{p'[t]((Mean[Table[p'[j],{j,t-from,t-to,incr}]]+ndsolnprime)/(Mean[Table[p[j],{j,t-from,t-to,incr}]]+ndsoln))*p[t]-ndsolnprime,p[t/;t≤0](ndsoln/.tt+gap)},p,{t,0,gap}];listplot=With[{vals=city["Cases"][[2,1,1,lower;;upper]]},Transpose[{Range[Length[vals]]-1,vals}]];Plot[{If[t>gap,None,ndsoln],If[t>gap,ndsolnNextPeriod[t-gap],None]},{t,0,2gap},Epilog{Red,Opacity[.6],PointSize[Medium],Point[listplot]},AxesOrigin{0,0},PlotStyle{Automatic,Dashed},PlotLegends{"retrodiction","forecast"},PlotLabelStringJoin[name," caseload "]]]

Some examples of past-to-present with future forecast

I will start with a solution for recent behavior in Bergen County (this is the county in northern New Jersey that has had the most cases in that state to date).

In[]:=

currentAndProjected[bergenC,21,"Bergen County"]

Out[]=

	retrodiction
	forecast

In this scenario the number of active reported cases peaks in a week or so at just above 13,000 (it is currently just over 11,400).

Now we posit a longer recovery time of 28 days.

In[]:=

currentAndProjected[bergenC,28,"Bergen County"]

Out[]=

	retrodiction
	forecast

For this longer recovery time the peak will not happen for two weeks, and it will exceed 17,000. Note that this can happen even though the per-day seems to be dropping over these past several days. The more prolonged the illness, the more time for the overall case load to rise and the more delay before we get a drop. One will also notice that the first plots, that recap the behavior of the past few weeks, are not quite identical in these two scenarios. As I mentioned above, the code for this is a bit intricate and results are not perfect. That stated, both plots show solutions that come close to matching the reference values (shown as red dots).

Here are scenarios for Hudson and Essex counties (containing Jersey City and Newark, respectively).

In[]:=

currentAndProjected[hudsonC,21,"Hudson County"]

Out[]=

	retrodiction
	forecast

In[]:=

currentAndProjected[essexC,21,"Essex County"]

Out[]=

	retrodiction
	forecast

For 21 days recovery times, Hudson fares relatively worse than Essex and Bergen. Indeed, Bergen county has nearly twice the population of Hudson but only around 20-25% more cases at present. Moreover the

i'/i

ratio for Bergen is notably lower than for Hudson, and this serves to explain why the future modeling has these different outcomes. This difference might be from any number of reasons, with likely culprits being socioeconomic differences (correlating with health and environmental issues) and population density (which is quite high in Hudson County).

We take a look at New York City.

In[]:=

currentAndProjected[nyc,20,"NYC"]

Out[]=

	retrodiction
	forecast

One can only hope this model has merit. Certainly I do not know whether this will be so. We get to find out over the next few weeks.

Below are plots for a number of other cities. Most use recovery times of 21 days but, for variety, a few do not.

In[]:=

Apply[currentAndProjected[#1,21,#2]&,Transpose[{{sanfran,seattle,nola,atlanta,memphis,denver,vegas,dallas,detroit,ftlau,baltimore,boston,miami},{"San Francisco","Seattle","New Orleans","Atlanta","Memphis","Denver","Las Vegas","Dallas","Detroit","Fort Lauderdale","Baltimore","Boston","Miami"}}],{1}]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



Below are plots for the five largest cities after NYC. Here I use a recovery average time of 24 days.

In[]:=

Apply[currentAndProjected[#1,24,#2]&,Transpose[{{la,chicago,houston,philly,phoenix},{"Los Angeles","Chicago","Houston","Philadelphia","Phoenix"}}],{1}]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



Backing up in time for comparison

One question to ask is whether this approach has any level of self-consistency. That is to say, are predictions for a given county changing wildly over time? It should be clear that they will generally shift, since they are based on noisy data and forcasting the future is like, well, forcasting the future. Here is where the optional argument to currentAndProjected comes into use. For a given city we can do our past-and-future plot, then compare it with what we would have obtained had we run it on data from several days ago. Below I show some examples. The findings are good in some ways. For one, predicted day of peak caseloads do not seem to change by much in most cities I have tested. Curve shapes do not change too much. Predicted maximal caseloads do change but even that does not seem drastic in most cases.

Here are examples using New York City, with data from 6, 3, and 0 days ago respectively.

In[]:=

Map[currentAndProjected[nyc,21,"New York City",#]&,Range[6,0,-3]]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



Here are similarly dated predictions for Chicago.

In[]:=

Map[currentAndProjected[chicago,21,"Chicago",#]&,Range[6,0,-3]]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



We repeat the curves for the North Jersey counties, and below show similar curves using data from three days ago.

In[]:=

Apply[currentAndProjected[#1,21,#2,0]&,Transpose[{{bergenC,hudsonC,essexC},{"Bergen County","Hudson County","Essex County"}}],{1}]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



In[]:=

Apply[currentAndProjected[#1,21,#2,3]&,Transpose[{{bergenC,hudsonC,essexC},{"Bergen County","Hudson County","Essex County"}}],{1}]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



Here are (again) curves for some of the larger cities, and below are similar curves using data from three days ago.

In[]:=

Apply[currentAndProjected[#1,21,#2,0]&,Transpose[{{phoenix,philly,houston},{"Phoenix","Philadelphia","Houston"}}],{1}]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



In[]:=

Apply[currentAndProjected[#1,21,#2,3]&,Transpose[{{phoenix,philly,houston},{"Phoenix","Philadelphia","Houston"}}],{1}]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



Here again are curves for three West Coast cities.

In[]:=

Apply[currentAndProjected[#1,21,#2,0]&,Transpose[{{sanfran,seattle,la},{"San Francisco","Seattle","Los Angeles"}}],{1}]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



These have reported cases going back two months or more. So we can backtrack by a week without too much risk (for areas with relatively shorter data periods, one cannot go back too far without losing accuracy).

In[]:=

Apply[currentAndProjected[#1,21,#2,7]&,Transpose[{{sanfran,seattle,la},{"San Francisco","Seattle","Los Angeles"}}],{1}]

Out[]=



	retrodiction
	forecast

	retrodiction
	forecast

	retrodiction
	forecast



While consistency is not perfect, we do see basic fidelity in most cases, at least if we do not go back too far. This gives some confidence that the method might apply whenever one has sufficient past data. Of course we will know soon enough whether it really captures the behavior with any reasonable accuracy.

Summary

I began with the premise that approximating the logarithmic derivative of the infection function might give useful information for assessing how well a location is getting its infection rate under control. This did appear to provide useful information, insofar as there seemed to be a small band in which critical values all lie. In cities with values above that band the rates are rising, and for those below it they are falling. So far so good.

I took this as an indication that the infection rate could be modeled using this ratio, along with a few assumptions about how to assess it via shifted values moving forward. This is on the speculative side, and it might turn out that I’ve simply had delusions of adequacy, as they say (alas, it will be far from the first time). If this method, or modest adjustments thereof, do turn out to have short-term predictive value, that can be useful in future scenarios where local transmission rates are affected by some level of overall quarantine. This all remains to be seen.

POSTED BY: Daniel Lichtblau

3 Replies

Sort By:

Daniel Lichtblau

Daniel Lichtblau, Wolfram Research

Posted 2 years ago

I have updated the model to better account for recoveries. I actually now use two different approximations to form a delay differential equation. While there remains some amount of hand-waving in the justifications, the results seem to be in accordance with current data for many of the locales I checked. Also some of the miss-prediction might be from increases in testing rather than underlying changes in the disease spread dynamics. Since the notebook contains a lot of graphics and is therefor on the longish side I am attaching it rather than as a Wolfram Cloud notebook. Edit: It seems the notebook is not appearing. Not sure why but I will make inquiries, as they say. Edit 2: It was added though not by me. I need to get a magic touch upgrade or something. Edit 3: I was informed that the attachment is difficult for some to download (maybe a size limitation on some computers?), so I have now posted directly as a Cloud-based notebook. Attachments: COVID19_spreadra...nb

POSTED BY: Daniel Lichtblau

Moderation Team

Moderation Team, WOLFRAM

Posted 2 years ago

-- you have earned *Featured Contributor Badge* Your exceptional post has been selected for our editorial column *Staff Picks* http://wolfr.am/StaffPicks and Your Profile is now distinguished by a *Featured Contributor Badge* and is displayed on the Featured Contributor Board. Thank you!