key: cord-0949597-xqy5tkqq
authors: Hachtel, G. D.
title: The Futures of the Pandemic in the USA : A Timed Intervention Model
date: 2020-08-25
journal: nan
DOI: 10.1101/2020.08.23.20180174
sha: 14f9e2c2b8b69596bce4f02fa28257d597c64e1b
doc_id: 949597
cord_uid: xqy5tkqq

We propose a novel Timed Intervention SPEIQRD extended SEIR model for predicting the evolution of the Covid 19 Pandemic in the USA. The model can, by parameter, assignment, be reduced to the model of Hong et al, which appeared in February, 2020. Novel aspects of the proposed model include 1. Formulation of a "Protected" population P, which can be viewed as a "Sheltered in Place", unexposed population which, starting at time t = Tau_P, builds up and stores a reservoir of unexposed Population; 2. This "Protection Intervention" provides the basis for a second, "Timed Release" Intervention: on receiving a "reopening signal" at time t = tau_R$, this second intervention initiates a release of the stored population P back into the general population S. 3. Selection of model parameters to "optimize" the approximation of the model up to the present, and then projecting simulation of the model based on the value thus obtained 100,200, 365, and 730 days into the future. This model shows excellent qualitative results and quantitative results that are good, considering the chosen nationwide scope. These results compare favorably to University of Washington projections, as published Daily on the Worldometers.info website. The qualitative and quantitative behavior of all 7}state variables S,P,E,I,Q,R,D will be illustrated and discussed, in this and possible future further Intervention cyclesfootnote{"How long, I wondered, will this thing, last?} Lyric from "A Foggy Day" by George and Ira Gershwin}

The pandemic has led to world wide burst of data science, modeling and prediction research, much of which focuses on extensions of the dynamic SIR model of Reference [8] , developed for modeling the Spanish Flu pandemic of 1917-18, which was exacerbated by the trench warfare of WW1. Some recent relevant academic research has been reported in References [2] , [6] , and [7] . Here we propose an extension of Reference [2] .

When the Pandemic began it was hard to avoid following the data reported daily by the CDC, Johns Hopkins University (arcGIS world maps) [5] (tn the USA Covid-19 first took hold in Seattle). Later, the www.Worldometers.info, and University of Washington IHME (as quoted on Worldometer) became go-to sources of data. The data in this paper is recorded daily into a cumulative CSV file with 3 fields: Total Cases CCases, Active cases ACases, and Deceased cases CDeaths. Recovered Cases RCases are computed from the subtraction RCases = CCases − ACases − CDeaths.

In Figure 1 2 is shown the daily record of statistics often emphasized in these data sources. At the left is shown the daily Confirmed Case Count Per Day, and Deaths per Day scaled up by a factor of 10. Annotations are never scaled. On the right is shown the daily Death Counts per Day. Note "per Day" data must be computed from the difference between consecutive values, thus introducing noise. These noisy Counts per Day data have been smoothed by fitting 15knot Natural B-Splines 3 . The data covers the interval March 3 (Day 0) to today, August 20 (Day 165). The smooth blue curve is an example of Gaussian Modeling, as discussed below in Section 3.4 below. Figure 2 demonstrates the primary motivation for developing the Timed Intervention Model. At the left, the Confirmed Cases per Day data show the alarming rise observed between Day 0 and Day 43, as well as the equally alarming second upsurge, which began around June 15 (Day 105).

The waveforms appear to have 3 separate phases. The first for t < τ P ≡ 43 could be called a "Denial" Phase. Then, in the second phase, for τ P < t ≤ τ R ≡ 105, the Deaths per day declined gradually. Finally the data turned sharply upward again. To borrow a term from the financial sector, this could be called a "Correction Phase". The third Phase, for t > τ R , looks like another "Denial" phase. The classical Predator-Prey models of the Lotka-Volterra Equations do not show such nuanced three-phase behavior. So, the search began for a model that does. After a few failed attempts, the model of Equation 1 and Equation 4 was settled upon.

To be especially noted is the similarity in the shape of this curve to that of the spline approximation to the noisy data on the right of Figure 1 .

Note also in Figure 2 that a fourth phase is also indicated for t > τ P 2 . This phase will be investigated in a subsequent report.

Our dynamical pandemic model is represented by the following state equations:

(1) This definition is completed by the specification of the first and second intervention parameters α(t) and φ(t), given in Equation ??.

A similar definition applies to the Second Intervention Parameter φ:

In this paper, N is regarded as fixed and γ = beta/RN 0 as derived, leaving just 14 assignable parameters in Equation 4 . For any given simulation run, the model is defined by the values assigned to these 14 free parameters and the 7 Initial Conditions, so there are 14+7 = 21 free parameters in principle. However, in this paper, the Initial Conditions of Equation 6 are used exclusively. The values shown above are the default values. For reproducibility, any set of results show in this paper are defined by defining multipliers of these default values which are stored in the data source. 

3 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. . Figure 2 that in the 3rd phase, the Cases per Day has peaked and is currently falling, resemblmg a repeat of the first phase

This model is conservative in the sense that it assumes the total population is constant. That is, S +Ṗ +Ė +İ +Q +Ṙ +Ḋ = 0,

Differences between this model and that of Reference [2] are apparent in the definitions of andṠ andṖ , and inṘ, andḊ. These changes enabled the model to obtain excellent qualitative and, in some respects, good quantitative results for the USA as a whole. It should be noted that each individual term in the right hand sides of Equation 1has an equal and opposite term in a subsequent or preceding equality. Consequently, we have Equation 5, and the total population is fixed by Construction.While fixed total population is not strictly true in the USA, this is an important aspect of the model.

Of special interest are The proposed model assumes that re-infection cannot occur. Also, we note that this model includes novel direct transitions. The (ν,and ρ) parameters control transitions from the Infected state I to the Recovered state R and Deceased state D. These transitions and parameters are not found in Reference [2] , Since Equation 4 defines 14 free parameters, that the parameter space is is 14-dimensional (or 120 if we include the initial conditions for the 7-dimensional State Space). The essence of modelling is to solve the Parameter Identification Problem: Find the point in the 14 (or 20) dimensional parameter space for which the model results best fit real world data. To solve this rigorously, a second round of research is proposed: find the optimal 13(21) parameter assignments using the approach of Refreence [1] Equation 1 specifies a set of 7 ordinary differential equations In which thė E andṠ equations have a bilinear form which first appeared in Equation 29 of Reference [8] . More recently, they are also known as Lorenz Equations [3] , who used them to model compartmental layers of the atmosphere. Such equations are known to exhibit unstable, and even chaotic, behavior. Fortunately, in our limited exploration of the 11-dimensional parameter space, no exotic behavior has been observed. Perhaps, the last 4 linear equalities constitute a dissipative damping effect which inhibits oscillatory behavior.

The odeint function from the python3 SCIP Y distribution was used to solve these 7 ordinary ordinary differential equations over the closed time interval t ∈ [0, T ], where t represents time in days, and T is the Final time in days. Note that since we extend the model 100 days into the future, the current time, denoted T oday, is near the middle of the interval t ∈ [1, T ], In this simulation 5 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020 . . only unit time steps are considered in displayed data from the model Of course th odeint function takes arbitrarily smaller steps as required. Fortunately, for the interesting part of the parameter space, the solutions of these equations are well-behaved, and have been used frequently as epidemiology models for over a century, including multiple times this millenium. Computationally, this simulation is easy, even trivial.

The above model has evolved significantly from the Covid-19 model described by Peng et. al in Reference [2] (February, 2020). The chosen scope is the USA as a whole, in which our form of government inhibits the imposition of a uniform approach to controlling the pandemic. It is proposed to treat individual states of interest, e.g. Colorado, California, and New York, etc, in a third phase of research in which we will apply an instance of the above model for each hot spot state, as well as the USA as a whole. This will be referred to as the distributed [7] model. Computational experience leads to the expectation that even a distributed model including all 50 states and the District of Columbia as well as the entire USA is still computationally practical.

reading the papers cited above, research commenced,and after getting some initial results, a search of the literature for other work on Covid-19 modeling was begun. Early discoveries included Reference [4] , and Reference [6] The first of these came out before the Pandemic started in China. It is focussed on control theoretic and Chaotic analysis in the phase space. As discussed above such exotic behavior does not arise in our numerical work. After Covid-19 first broke out in December 2019, the work of Reference [2] came too light and published data from Johns Hopkins University, the CDC, and Worldometers was followed.

The title of Reference' [6] is worthy of note : "The unreasonable effectiveness of simple models". Its authors state categorically that modeling the Pandemic must be some form of SEIR approach. And yet early results contradicted the direct applicability of this approach, which leads to more or less symmetrical results well modeled by Gaussian and/or Logistic (integral ofGaussian) functicns. With this approach important populations reach their peak, and then fall toward extinction, far too rapidly to fit measured data. The "Dr. Chris Murray Model" [5] doesn't have this problem, but despite being used and quoted widely, its details have not been published.

It was felt that additional numerical work was needed, with timely dissemination of the modeling approach. In particular new and publushed projections of Mortality were still needed. Alarmed by the early data, the project felt urgent, even more urgent due to alarming recent data. In Reference [7] , A Distributed, Difference equation approach Model of Community Spread was used that was sufficiently general to include re-opening and mask rules. Results were given for 6 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020 . . total deaths nationwide as well as on a state by state basis. This work has been ahead of the curve in predicting the sharp uptick in National and State Total cases.

Numerical studies required the choice of the 7 initial Conditions and the 14 free parameters. In the sequel only the initial conditions S = N, P = 0, E = 50000, I = 2500, Q = 0, R = 0, D = 0 (6) were employed. 

For these chosen values, the results are illustrated in Figure 3 .

In order to compare these results to source data, the simulated Recovered and Deaths per Day populations for and just the first 154 days up to the current date, and for the chosen parameters is plotted in Figure 3 (1) that prior to the Protection Intervention, the susceptible population is drained by both the linear term and the powerful bilinear term in the equation forṠ. However, after the first Intervention, it is drained much more slowly due the intervention term φ 0 P (t).

Until the release of the Protected, Population P builds up correspondingly, and consequently drains S. After the release, S is resupplied from P , so P is 7 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. drained until t = T On the right is shown the corresponding waveform for the Protected Population P . Note the interesting twin inflections in P .

Qualitatively it can be seen that S and P are (roughly) complementary, that is, S ∼ N − P . Consequently in the sequel focus is mainly on P .

Also, note that there is still interesting behavior beyond the end of a 100 Day Projection. Hence 200 Day Projections are adopted in the sequel.

Early in the 6-month course of research it was noted that many derived population statistics to date can be accurately represented by Gaussian Approximations, for example, the no intervention simulated waveforms for E and I.

When 100-and 200-day Projections were added, the simulated solution of the above 7 differential equations became distinctly tri-modal in accordance with the initial Shelter in Place Orders, which effectively ended the Denial Phase, followed by apparently premature Re-opening and Easing of mask restrictions. Apparently, this leads to a second Denial Phase. It is conjectured that the future of the Pandemic, could be Cyclical, that is, a convergent series of Denial, Correction phases. On this somber view the Pandemic could be with us for years. The conditions for which this cyclical behavior can occur are discussed in Section 4.

For large populations, the law of large numbers declares that for any random variable representing a property of a single member of the population, no matter how it is distuributed, the population as a whole for that random variable will have a normal distribution of that property. As we showed in Figure 3 , E and I appear to be normally distributed, that is, show (almost) Gaussian waveforms with parameters mu and sigma.. However it is to be noted that after the maximum, the value of sigma appears to have increased. This is a salient feature of Figure 3 discussed above, and in Section'3.2 below.

For source populatioins R and D (not yet shown) a logistic (error) function (Integral of Gaussian) is to be expected, since these states and their dis-9 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. tributions, are cumulative. From a probabilistic perspective, the parameters α, β, γ, δ, λ, κ can be interpreted as conditional probabilities.

When we look at the overall response of the model, including the last 100 or 200 days of Projection, we see that the Timed Release effect adds a second mode in a bimodal response, as shown at the left of Figure 5 , At the left, projected data for E and I,from late in our mode simulation is shown in solid lines, whereas the Gaussian model prediction is shown in dashed lines. The vertical bars show the fitted values for the Gaussian means muE and muI. The fit is excellent.

At the right of the figure, the Deaths per Day day up to the current date only (so no simulation is involved) is considered. For this case a double, or split, Gaussian is employed. This means one Gaussian is fit to the first phase, or rising part of the curve. For the second phase, a second Gaussian with same maximum value, but a larger standard deviation sigma was chosen to match the slow decline of this phase. This gives an excellent fit until the onset of the third phase, after which the Deaths per Day data turns up again and diverges sharply after the Timed Release phase is entered. However, since reported Deaths is a lagging indicator, we expect the difference to increase significantly. We observe that the Exposed Population is constant or even a decreasing Gaussian tail.

Shown in the blue waveform at the right of Figure 5 Part (b), is an example of modeling with a split Gaussian. This means that a rising Gaussian is fit to the smoothed data up to the maximum to date (day 43, or April 19), and also a separate, falling Gaussian to the falling data This produces an excellent fit, at least prior to the Timed Release. After that , the falling Gaussian diverges markedly from the source data which has turned upward, suggesting that a second Protection intervention might be called for.

Finally, note the similarities and differences between Figure 5 Part (b) and Figure 3 Part (b) (bottom row). (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. 12 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. 

Similarly, to obtain the data on the right of Figure 8 , φ has been reset to the default value φ G Both columns of this figure show some remarkable features-instead of falling after the Timed Release, P rises on the left and is rendered almost constant on the right. This results in P being only slightly reduced at the final time t F . The implication of this shift is evident in the Waveforms for E,I, D, andḊ.

The effects of the final two variations are illustrated in Figure 9 . The results are similar overall but have some subtleties which warrant further study.l

It is also interesting to consider the question: "How many Survived in this Simulation? Assuming that re-infection is impossible, the R population has definitely survived. To answer this question it was decided that the S and P Populations should be considered to have survived, because of the likely, but eventual, emergence of a successful vaccine. Thus the survivors are defined by the sum Ω = R + S + P + Q. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. 15 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. 

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. Figure 11 Shows similar results. These last two figures suggests that the survivors are maximized by maximizing the Second Intervention parameter φ

We have proposed and studied a Timed Intervention model and applied it to the population of the USA as a whole. It was shown that this model can predicts 17 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020 . . a 3 (or more) phase dynamic response over a wide range of conditions and parameter sets The dynamic behavior of our SPEIQRD model is explored in detail and found to be quite rich in producing waveforms of interest.

A salient finding of this study is its relatively conservative projection of the cumulative current death toll in the USA-about 168000, to date, for a parameter set which leads to Death per Day data consistent with reported values to date. However, the data in the left column of Figure 7 , second figure from the bottom, predicts 5.5 million deaths, 200 days out. The USA evidently needs coordinated intervention to avoid these catastrophic scenarios.

The data suggest that the calamitous rapid rise and fall of Predator-Prey models cannot be eliminated, but can, with the proper intervention cycle, be pushed into the future. Unfortunately, the proper intervention might well be a much more rigorous lockdown than yet attempted.

The sensitivity of the results to variations in the Reproduction number R 0 and the Second Intervention parameter φ were analyzed. This study produced some striking 3-Phase results.

A definition of the Survivors of a given Pandemic was given and analyzed. It was shown that Survivor Count is maximized if φ is maximized.

Further model development should include: 1. A fully optimized solution to the Parameter Identification problem defined above is particularly needed-until this is done, results cannot be said to be characteristic of the model, rather than of the parameter sets chosen; 2. Inclusion of data for an extended projection of 730 days so the implied long term maxima would be fully revealed;

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 25, 2020. . https://doi.org/10.1101/2020.08. 23.20180174 doi: medRxiv preprint

2 100 day Projeiction Results with No Interventions

Effect of Timed Protection and Release Interventions on Susceptible and Protected populations

200 day Projection Results with Both Interventions

Long Term effects of Second Intervention Parameter

4 2 Noisy and Smoothed Cases per day Data showing the 3 Phases induced by the "Interventions

Day Projections;Ḋ no Projection,with No Intervention and RN 0 = 10

Both Interventions RN 0 = 20, 15

Survivors Ω(t) and Recovered R(t) for φ = φ G * 10 * * 0.5 (left)

Survivors Ω(t) and Recovered R(t) for φ = φ G /10 * * 0.5 (left) and φ = φ G /10

The sparse tableau to network analysis and design

Epidemic analysis of covid-19 in china by dynamical modeling

The nature and theory of the general circulation of the atmosphere

Analysis and control of an seir epidemic system with nonlinear transmission rate. Mathematical and Computer Modeling

The unreasonable effectiveness of simple models

Projection of covid-19 cases and deaths in the us as individual states re-open. wordpress.com Columbia University

A contribution to the mathematical theory of epidemics