On a Coupled Time-Dependent SIR Models Fitting with New York and New-Jersey States COVID-19 Data.

This article describes a simple Susceptible Infected Recovered (SIR) model fitting with COVID-19 data for the month of March 2020 in New York (NY) state. The model is a classical SIR, but is non-autonomous; the rate of susceptible people becoming infected is adjusted over time in order to fit the available data. The death rate is also secondarily adjusted. Our fitting is made under the assumption that due to limiting number of tests, a large part of the infected population has not been tested positive. In the last part, we extend the model to take into account the daily fluxes between New Jersey (NJ) and NY states and fit the data for both states. Our simple model fits the available data, and illustrates typical dynamics of the disease: exponential increase, apex and decrease. The model highlights a decrease in the transmission rate over the period which gives a quantitative illustration about how lockdown policies reduce the spread of the pandemic. The coupled model with NY and NJ states shows a wave in NJ following the NY wave, illustrating the mechanism of spread from one attractive hot spot to its neighbor.

1. Introduction and Model

From a scientific perspective, the COVID-19 pandemic has highlighted the crucial role of mathematical and statistical models in providing guidance for health policies. Expressions such as “flatten the curve”, the “apex”, the “plateau” have been widely heard in media and employed by decision makers to explain their choices regarding rules and policies during this critical period. In this short article, we first introduce a simple Susceptible Infected Recovered (SIR) model, in which we adjust a key parameter k standing for a control on the Susceptible-Infected rate, and secondarily the death rate, in order to fit the data of the pandemic in New York (NY) state in March 2020, and provide predictions for a near future. Then, we add a node in the model to take into account the daily fluxes between NY and New Jersey (NJ) states. Note that these two close states, are, up to the day of publication of this article, the most severely hit by the pandemic in United States. Of course, the coupling may be extended to other states. However, in this article, we restrict ourselves to NY and NJ. Accordingly, the main key points of this article are that, (1) it highlights the dynamics and epidemiological characteristics which have been discussed in press and health policies; it highlights qualitatively how lockdown policies have decreased the spread of the virus and provides prediction and explanation of an upcoming apex, (2) it fits real data provided for the New York state and (3) it fits the data of NJ state by considering coupled equations taking into account the daily fluxes between NY and NJ. This provides a quantitative visualization of how the virus may spread from an attractive hot spot (New York City in NY state) towards close states trough the daily fluxes of commuters.

We especially focus on fitting the total number of cases tested positive for COVID-19 as well as the number of deaths in both NY and NJ states. We also give insights in prediction of the number of people needing hospitalization in NY state.

SIR models are very classical in literature. For some reader’s convenience, we mention here a few contextual elements and references. The simplest classical SIR model is the Kermack–McKendrick (KMcK) model which goes back to 1927, see [1,2]. It is written

In this original KMcK model, the population splits into three classes. The class S stands for susceptible, who can catch the disease, I stands for infective, who have the disease and can transmit it and R stands for the removed, namely, those who have or have had the disease but do not transmit it anymore. Note that our terminology is slightly different as explained below. In (1), the dynamics follow the scheme with respective transfer rates between classes of k and r.

SIR type models, and more generally mathematical models of epidemics have in fact a significant history. We refer to [2] for a textbook on these models and a brief history of epidemics. Models have become more sophisticated and may include more compartments such as exposed, infective asymptomatic, infective with symptoms, and also reservoirs such as bats, and include stochastic dynamics. Recently, SIR type models have been widely used in the context of the COVID-19 pandemic, to model the spread all around the world, see for example [3,4,5,6,7,8,9] and reference therein cited. Here are also some examples of references for SIR models in other epidemic diseases: Dengue [10], Chikungunya [11,12] and Ebola [13]. See also references therein cited.

In the present article, we first consider the following model

This simple model has classically three classes: S for susceptible, I for infected and R for recovered. Specifically, the class I is intended to represent all the people who bear effectively the virus at a given time, and can transmit it if in contact with other people. It includes all infected people with or without symptoms, reported or not. There are some differences with Equation (1). First, it includes a death rate . Even though the number of deaths does not appear explicitly as a variable, it is simply given by the integral . Note that in expression the rate of contamination from S to I is proportional to the proportion of susceptible (S) in the whole population (). This is a classical expression standing for the fact that the probability for each individual in the I class to spread the virus among the class S is proportional to the portion of S in the whole population, see for example [14,15] and references therein cited. This rate is corrected by a crucial coefficient which is intended to fit the real transfer rate and which contains the effects inherent to the properties of the virus (for example change of propagation rate due to genetic mutation of the virus) or to specific policies (like quarantine, social distancing, lockdown). This time dependence allows us to adjust the dynamics to fit the data. This is a specificity of our model and turns it into a non-autonomous equation. This time dependence of k is obviously relevant in our model since the rate of transfer from S to R is the main target of health policies and is subsequently subject to vary over time. Secondarily, we also allow the death rate to vary. Many internal or external factors may affect the death rate among which are concomitant lethal disease, temperature, hospital conditions. More significantly, one has to note that the rate transfers considered here are instantaneous transfer between compartments, and the function is different from the Case Fatality Rate (CFR). Recall that the CFR is the death rate per confirmed case over a given period of time, and is a typical indicator for death rate. In South Korea, the country which led the highest number of tests, it has been reported to be of 1 percent, see [16]. In China as of 20 February, this rate varied between 3.8 in the region of Wuhan and 0.7 in others regions, see [17] and also [18]. Using Table 1 and Table 2 gives a CFR in NY from 1 March to 1 April of per cent. Since our model fits the data, by definition, it fits the CFR. At the end of the epidemic, if the whole number of people that contracted the virus would be tested positive the CFR would provide the probability to die for an individual who catches the virus. However, during a growing phase such as the month of March considered here, the CFR has large variations. Moreover, there is a delay between the time a person is tested positive and the time he dies. The time between symptom onset and death has been reported to range from about 2 to 8 weeks, see [18], the typical average being 23 days according to [16]. One could for example introduce a time-integral death expression to take into account these informations. However, the time-window considered here is short and corresponds to the beginning of reported cases in NY. Furthermore, in this short article, we wanted to focus on a simple model able to fit data, highlight relevant dynamics and provide estimations. Since a person in the I compartment will either recover or die, above remarks on the function d hold for the coefficient r. In the present work, for sake of simplicity, we have set the r coefficient to the constant value . This is a simplification which is classical in SIR models. Note that setting the coefficient r to a constant value is equivalent to assume that people recovering between times and t, which is given by would also write according to Equation (2). A more meaningful expression for would be standing to the fact that people recovering between time and t have been infected from a period ranging from 0 to t, with a transfer rate given by . This would lead to the equation . Once r is set to a constant value, to fit data over a reasonable period, numerical tries provide a unique choice for constants k and d. It was still possible to use different values of r to fit the data. However, different values of r result in different dynamics over time, and notably different time for apex. Our choice of was made to provide dynamics that seem relevant to us regrading the timely dynamics beyond the data. In particular, too smaller of values for r would provide dynamics with a late apex and less relevant regarding the timely effects of the disease. Other studies have considered models with a constant r, with different values. See for example [19] and references therein cited. Remarkably, varying r, and making it depend on time, provides some freedom to later fit data over a longer period of time, taking in account the end of the first wave in NY.

Table 1

Total number of cases reported in NY state from 1 March to 1 April. See [20].

Day	3/1	3/2	3/3	3/4	3/5	3/6	3/7	3/8	3/9	3/10	3/11
Number of Cases	1	1	2	11	22	44	89	106	142	173	217
Day	3/12	3/13	3/14	3/15	3/16	3/17	3/18	3/19	3/20	3/21	3/22
Number of Cases	326	421	610	732	950	1374	2382	4152	7102	10356	15168
Day	3/23	3/24	3/25	3/26	3/27	3/28	3/29	3/30	3/31	4/1
Number of Cases	20875	25665	33066	38987	44635	53363	59568	67174	75832	83804

Table 2

Total number of deaths reported in NY state from 1 March to 1 April. See [20].

Day	3/1	3/2	3/3	3/4	3/5	3/6	3/7	3/8	3/9	3/10	3/11
Number of Deaths	0	0	0	0	0	0	0	0	0	0	0
Day	3/12	3/13	3/14	3/15	3/16	3/17	3/18	3/19	3/20	3/21	3/22
Number of Deaths	0	0	2	6	10	17	27	30	57	80	122
Day	3/23	3/24	3/25	3/26	3/27	3/28	3/29	3/30	3/31	4/1
Number of Deaths	159	218	325	432	535	782	965	1224	1550	1941

Upon the above discussions, our strategy is rather simple: set r to a constant value. Then choose a constant k which fits the data of positive cases during a period of time. Next, choose a constant d to fit deaths data for the same period of time. Note that since the parameter k has much more effect than the small parameter d, this procedure is possible and efficient. After, repeat the procedure over a subsequent period. The overall procedure results in a constant function r and two piecewise constant functions and . For our model and the given data, the procedure was efficient allowing us to provide these choices by successive tries. Note that this could be done automatically by cooking an algorithm to set the parameters following these guidelines. Note finally that our assumptions do not a include birth rate for the susceptible population but rather focus on the short time effects of the disease.

2. Numerical Simulations, Data and Dynamics

2.1. Fitting the Total Number of Infected People and the Number of Deaths

Data from total infected people and total number of deaths in New York state are available at the New York times web site, see [20]. We downloaded the data from there from 1 March to 1 April, which makes 32 days. We reported the number of total cases in Table 1 and the number of deaths in Table 2.

For numerical simulations of Equation (2), the parameter r was set to . Then, in order to fit the data of the total number of infected people, we chose: and to fit the total number of deaths, we chose:

Initial conditions were set to 19453556 being the number of people living in NY state in 2019 according to USA Census bureau. Regarding reported data in Table 1 and Table 2, as a matter of fact, not all infected people in March 2020 were tested. To take into account this fact in our model, we adjusted the parameters to fit the quantity from the model, with the total number of infected people (data in Table 1) minus the total number of deaths (data in Table 2), i.e., we assumed here that only 20 percent of living people that have contracted the virus has been tested positive. Note that shortly after the first submission of this work, NY state reported the result of a random antibody test over a sample of 3000 individuals. Statewide, 13.9 percent of people tested positive, ranging form 21.2 percent in New York City (NYC), to 3.6 percent in upstate counties. See [21,22]. Comparing with data from [20], this gave a ratio of approximatively 10; only about 10 percent of people having caught the virus in NY would have tested positive. Statistical estimates of related ratio, depending on the number of estimates provide a range of , see [23]. Simulation of system (2), and comparison with data resulting from Table 1 and Table 2 are plotted in Figure 1a. In Figure 1b, we plotted the total number of deaths from the model and compared it to the data. In Figure 1c, we have plotted four curves corresponding to various for different simulations of (2).

Figure 1

This figure illustrates the simulation of system (2) and how it fits the data. In (a), we have plotted the quantity as a function of time in red. The blue dots correspond to the data retrieved from [20]. Analogously, in (b), we have plotted the quantity , which represents the total number of deaths according with the model, as a function of time in red. The blue dots correspond to the data retrieved from [20]. In (c), we have illustrated the quantity corresponding to different values of : the curve in red corresponds to the simulation of (2) with and for all time. The curve in green corresponds to the simulation of (2) with and for , and and for . The curve in pink corresponds to the simulation of (2) with as given in (3), i.e., and , , with . It illustrates how the health policies flatten the curve. In (d), we have again plotted the solution for as in (3), for a longer period.

The curve in red corresponds to the simulation of (2) with and for all time.

The curve in green corresponds to the simulation of (2) with and for , and and for .

The curve in pink corresponds to the simulation of (2) with as given in (3), i.e., and , , with .

For reader’s convenience, these values are summarized in Table 3. This panel illustrates how the health policies flatten the curve. In Figure 1d, we have again plotted the solution for as in (3), for a longer period.

Table 3

Summary of the values used in Figure 1c to obtain the curves I1(t), I2(t) and I3(t). Recall that k1 = 1.057, k2 = 0.9, k3 = 0.67, k4 = 0.71, d1 = 0.0016, d2 = 0.00232, d3 = 0.00232, d4 = 0.0068, t0 = 0, t1 = 21, t2 = 24, t3 = 27 and t4 = 32.

	k(t)	d(t)
I_1(t)	k(t)=1.057	d(t)=0.0016
I_2(t)	1.057 if t<21; 0.9 otherwise	0.0016 if t<21; 0.00232 otherwise
I_3(t)	k_i if t_i−1≤t<t_i,i∈{1,…,4}	d_i if t_i−1≤t<t_i,i∈{1,…,4}

2.2. Fitting the Total Number of People at Hospital

Next, we have retrieved data corresponding to the number of people at the hospital between 16 March and 1 April, see Table 4. During this period, and after, NY state officially reported daily useful charts and statistics, on the local spread. We then computed an estimation of people being effectively at hospital at a given date. We denote the total number of hospitalizations at a given date. In order to fit the solution of (2) with this data, we performed a linear regression between and . The coefficients a and b of the linear regression were determined by the least-square method, leading to the formulas: with

Table 4

Total number of daily current hospitalizations reported in NY state from 16 March to 1 April. See [23,24].

Day	3/16	3/17	3/18	3/19	3/20	3/21	3/22	3/23	3/24	3/25
Total Number of Hospitalized	326	496	617	1042	1496	2043	2629	3343	4079	5327
Day	3/26	3/27	3/28	3/29	3/30	3/31	4/1
Total Number of Hospitalized	6481	7328	8503	9517	10929	12226	13383

The result is plotted in Figure 2a. It clearly shows a good approximation by two distinct lines, corresponding to and

Figure 2

This figure illustrates how to provide an estimation for people needing hospitalization thanks to Equation (2) and statistical methods. Panel (a) illustrates an approximation of by a vector and another vector where and , are the coefficients obtained thanks to the least-square method. Panel (b) then provides a prediction of people in need of hospitalization by plotting the quantity if , otherwise.

Then, a prevision of number of people needing hospitalization can be provided by the formula:

The result is plotted in Figure 2b.

2.3. Dynamics

Summing up the equations in (2) and looking for stationary solutions yield the following theorem.

We assume that initial condition of system (2) satisfies and

Non-negative stationary solutions of the system are given by

It is worth noting that this theorem, whose proof is relatively straightforward, provides two simple but relevant interpretations from the applicative point of view. The first thing is that the stationary solutions are given by

3. Two Coupled SIR Systems Fitting COVID-19 for NY and NJ States

It is estimated that around 400,000 people used to commute from NJ to NY before the lockdown policies due to the pandemic. It is natural to integrate these effects in the model, since those commuters played the role of vector for virus before the lockdown. We therefore built a small network by coupling two nodes, representing NY (node 1) and NJ (node 2). To this end, we coupled two copies of the local model (2) to take into account the daily fluxes between those two states. This led to the following model:

The assumptions are analogous to those given Section 2, but we aimed here to take into account the daily fluxes between NJ and NY. For sake of simplicity, we considered that there was a flux of people coming from NJ to NY in the morning and returning home at night. The coupling functions, which are non-autonomous here are given by: for node 1 and for node 2, where functions stand for the densities of population coming from mode j and going into node i. In the remaining of the manuscript we assume that these functions do not depend on and R and drop the superscripts. Furthermore, we assume and to be periodic of period 1 (one day), with a Gaussian profile and respective apex at 8:30 am and 6:30 pm. Note that the amplitude is multiplied by a coefficient greater than 1, in comparison with to take into account the amount of population int NY and NJ. We want to emphasize here the attractivity of NYC, and therefore assume that the fluxes are mainly from NJ to NY in the morning and form NY to NJ in the evening. The function integrates policies of lockdown: after March 23, in the model, the daily fluxes are divided by 1000. Therefore, the function is a piecewise constant function given by:

Initial conditions were set to for node 1 and for node 2. Note there are no infected cases at initial time in NJ. This means that in our model, initial spread in NJ follows from infection in NY. The same methods as in Section 2 were used to fit the data for both NY and NJ for the network model (6). Illustrations are provided in Figure 3. It shows how the model fits the data. In Figure 3a, we have plotted the quantity as a function of time in red. Recall that, in the model , represents the number of people in the population which has been infected by the virus and are still alive. The blue dots correspond to the data retrieved from [20] and plots the total number of infected minus the number of total deaths. Analogously, in Figure 3b we have plotted the quantity as a function of time, in red. The blue dots correspond to the data retrieved from [20]. Finally, Figure 3c illustrates and , which represent respectively the infected population in NY and NJ. It shows how the curve of NJ follows the curve in NY, with a small attenuation.

Figure 3

This figure illustrates the simulation of system (6) and how it fits the data. In (a), we have plotted the quantity as a function of time in red. Recall that, in the model , represents the number of people in the population which has been infected by the virus and are still alive. The blue dots correspond to the data retrieved from [20] and plots the total number of infected minus the number of total deaths. Analogously, in (b) we have plotted the quantity as a function of time in red. The blue dots correspond to the total number of deaths in NJ according with data retrieved from [20]. Panel (c) illustrates and , which represent respectively the infected in NY and NJ.

An analog theoretical result as in Section 2 holds for solutions of (6).

We assume that initial condition of system (6) satisfies and

For any non-negative are solutions of (6). Furthermore,

As in Remark 1, we want to point out some interpretation of this theorem in the context of the pandemic. Here, we have described some solutions with free epidemic component (

4. Conclusions

In this article, we have considered a simple non-autonomous SIR model to fit the data of COVID-19 in New York state. The model illustrates and quantifies how acting on the control allows us to flatten the curve of infected people over the time. From the model, using classical statistical methods, it is then possible to provide predictions of the number of people in needs of hospitalization. Lastly, we have fitted data from NJ state thanks to a coupled SIR model taking into account the daily fluxes between NJ and NY. It allows us to predict similar dynamics in NY and NJ, with a delay and small attenuation. Note that, despite its simplicity, our model fits, as good as more sophisticated models, the available data during the growing phase of March 2020. In a forthcoming work, we aim to fit data with the model over a longer period. This would provide an accurate estimation of the different parameters.