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.

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.

Forecasting based on past results

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 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.

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).

We take a look at New York City.

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.

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

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.

Here are similarly dated predictions for Chicago.

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

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

Here again are curves for three West Coast cities.

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).

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.