An SEIR Model with Time-Varying Coefficients for Analyzing the SARS-CoV-2 Epidemic.

In this study, we propose a time-dependent susceptible-exposed-infected-recovered (SEIR) model for the analysis of the SARS-CoV-2 epidemic outbreak in three different countries, the United States, Italy, and Iceland using public data inherent the numbers of the epidemic wave. Since several types and grades of actions were adopted by the governments, including travel restrictions, social distancing, or limitation of movement, we want to investigate how these measures can affect the epidemic curve of the infectious population. The parameters of interest for the SEIR model were estimated employing a composite likelihood approach. Moreover, standard errors have been corrected for temporal dependence. The adoption of restrictive measures results in flatten epidemic curves, and the future evolution indicated a decrease in the number of cases.

INTRODUCTION

The use of epidemic models permits to simulate disease transmission dynamics to detect emerging outbreaks and to assess public health interventions (Boily et al., 2007; Unkel, Farrington, Garthwaite, Robertson, & Andrews, 2012). With the scope to describe the dynamics of epidemics, standard methods, such as the Susceptible, Infectious, Recovered (SIR) model (Anderson, Anderson, & May, 1992), divide the population into portions of subjects on the basis of their relation concerning the epidemic vector. Here, the focus is on the dynamic such as the depletion of the susceptible portion to the infected one or the possible evolution of the rate of immunization.

However, the standard SIR model, and other extensions as the susceptible‐exposed‐infected‐recovered (SEIR) model, do not take into account the time‐varying nature of epidemics and several attempts were made to overcome this limitation (Boatto, Bonnet, Cazelles, & Mazenc, 2018; Dureau, Kalogeropoulos, & Baguelin, 2013; Kucharski et al., 2020). In particular, most extensions were proposed for adapting the SIR model to specific case studies (Liu & Stechlinski, 2012; Peng, Yang, Zhang, Zhuge, & Hong, 2020) or to include time‐varying coefficients with the scope to estimate epidemic dynamics (Chávez, Götz, Siegmund, & Wijaya, 2017; Fang, Nie, & Penny, 2020; Hooker, Ellner, Roditi, & Earn, 2011).

This article considers a flexible extension of the SEIR model, which incorporates the temporal dynamic connected to the transmission rate parameter, which is one of the most critical indicators for epidemiologists and on the basis of the basic reproduction number . Also, this method allows us to make considerations about both the trend and the prediction of the number of infected cases to evaluate how any possible influencing factors like as the presence of a vaccine or restriction measures taken by the central authorities can affect an epidemic outbreak (Haas, 2020).

The proposed method is applied to the 2019–20 coronavirus pandemic, the coronavirus disease 2019 (COVID‐19), caused by a Severe Acute Respiratory Syndrome CoronaVirus 2 named SARS‐CoV‐2 (World Health Organization, 2020). The World Health Organization declared the outbreak to be a Public Health Emergency of International Concern on January 30, 2020 and recognized it as a pandemic on March 11, 2020 (Eurosurveillance Editorial Team, 2020).

It is worthy to mention that several studies are known in the literature concerning SEIR models with different time‐varying parameter specifications (Hong & Li, 2020; Loli Piccolomini & Zama, 2020; Petropoulos & Makridakis, 2020; Wu, Leung, & Leung, 2020; Zhang et al., 2020). Our proposal differs in a statistical model consistent with counting data, a semiparametric (and therefore more flexible) specification of the time‐varying parameters. We also pay more attention to assessing the uncertainty of estimates.

To study how country‐based mitigation measures influence the course of the SARS‐CoV‐2 epidemic (Anderson, Heesterbeek, Klinkenberg, & Hollingsworth, 2020), we have looked to the ongoing epidemic in the three countries (Italy, Iceland, and the United States) where the adopted mitigation measures have been different (Dong, Du, & Gardner, 2020; Gudbjartsson et al., 2020; Remuzzi & Remuzzi, 2020).

The reference data sets are presented in Section 2, while the proposed model and the statistical inference are illustrated in Sections 3 and 4. In Section 5, we collect our results for the different countries with a proposal for the forecast. We end the article with a brief discussion in Section 6.

COVID‐19 DATA SETS

According to the aim of this article, we used the screening data of daily new cases and the number of the total amount of positive cases of SARS CoV‐2 according to three countries in different phases of an epidemic outbreak. For each country, we considered the data in a time window that starts about 15 days before the adoption of restrictive measures and the sources of the data are described in the Appendix.

At the time of writing, the United States was in the growing phase of the epidemic outbreak with an increasing trend of new cases. In the United States, every single state could decide the need to adopt stay at home measures. In those states that adopted restrictions, we registered different starting dates: the earliest state was Puerto Rico (March 15, 2020), followed by California (March 19, 2020) and New York (March 20, 2020) where the highest increase of new cases was subsequently recorded. Every single state could adopt a different and inhomogeneous panel of restrictive measures and in some states, as Arkansas, Iowa, Nebraska, and North Dakota, the local government never issued stay‐at‐home orders (Lyu & Wehby, 2020). Data were analyzed from March 4, 2020, to April 27, 2020, for a total of 55 days of observation. At the end of the considered period, we reported 820,514 current positive cases, 56,259 deaths, and a total of 988,197 confirmed cases nationwide.

Italy was a country in the middle phase where a first stabilization of the SARS‐CoV‐2 incidence was reported after the restriction measure, and, at the time of writing, we observed the beginning of a decreasing trend. The Italian Government adopted a national home lockdown restriction on March 9, 2020, for all the population followed by more severe measures on March 11, and ordered all nonessential businesses to close on March 22. Data were analyzed from February 23, 2020, to April 28, 2020, for a total of 66 daily observations. At the end of the considered period, we reported 105,813 current positive cases, 26,977 deaths, and a total of 199,414 confirmed cases.

Finally, Iceland was a country in the ending phase, where after stabilization, the incidence of new cases was going down, and the current epidemic outbreak was probably going to the end. The Iceland government adopted stricter measures to slow down the spread of SARS‐CoV‐2 on March 16, 2020, with an active searching strategy of new cases that lead to perform oropharyngeal swabs to about 10% of the entire population. We considered data from February 29, 2020, to April 27, 2020, for a total of 59 days of observation. At the end of the considered period, we reported 158 current positive cases, 10 deaths, and a total of 1,792 confirmed cases.

SEIR MODEL WITH TIME‐VARYING COEFFICIENTS

SEIR Model

We start introducing the SEIR model, which is one of the most used extensions of the standard SIR model, an ordinary differential equation (ODE)‐based epidemiological model (Kermack & McKendrick, 1927). Traditionally, the SEIR model divides a population of hosts into four classes: Susceptible (S), Exposed (E), Infected (I), and Recovered (R). However, in our framework, the last class should collect all the subjects that move outside the (I) status, i.e., recovered and deceased; for this reason, hereafter we denote (R) as Removed status. The model describes how the different portions of the population change over time . In the standard SEIR model, deaths are modeled as flows from the , , , or compartment to outside, because natural deaths are normally not monitored. If , , , and refer to the numbers of individuals in each compartment, then these “state variables” change according to the following system of differential equations:

In Equations (1), is the total population, is the mortality rate, is the transmission rate, is the exposed to infectious rate, and is the removal rate that can broadly assumed to be the sum of , where and are the recovery and the mortality rate, respectively.

In general, is called transmission rate that is the number of people that a positive case infects each day; in our settings, is defined as equal to , where is the contact rate that is the average number of contacts per person in a day, while is the probability of disease transmission in a single contact. However, and cannot be identified on the basis of the current information. The ratio permits to adjust taking to account people who cannot infect each other.

The parameters and are strictly dependent on the specific disease causing the epidemic and on the fraction of susceptible population. The parameter is set equal to where is the incubation period that may be higher for asymptomatic subjects; is the recovery rate calculated as where is the average duration of the disease in days.

Moreover, unlike the full specification, we do not consider the effect of births in model (1), and therefore, represents the number of new infected.

Based on this parametrization, we can define the reproduction number, , as The index conveys the strength of contagious in an epidemic outbreak. In the case of both and , can be approximated by .

Time‐Varying Parameter Specification

The standard SEIR model does that the parameters , and are time‐invariant. However, the characteristics of an epidemic suggest us that these parameters can vary. In particular, the overall mortality rate may increase if the number of deaths in a population directly or indirectly attributable to the disease (i.e., the insufficient capacity of health services) rises. The rate may also vary according to social distancing policies or, even, the isolation of infected people.

We aim to evaluate as the actions taken by the governments, and how a different degree of travel restrictions, social distancing, or limitation of the people movement can affect the epidemic curve of the infectious population. Our working hypothesis is that if there is an effect of the actions, they only affect the transmission rate of the epidemic, . For this reason, we propose to modify this parameter over time, namely, Since the function takes positive values, in the estimation step, we consider the following log‐linear specification: where , , are natural cubic spline basis functions evaluated at equally spaced knots in addition to the boundary knots.

The representation in (3) has the advantage that the estimation of reduces to the estimation of the coefficients . We refer to the next subsection for a short discussion about the number of knots and their positions. The time‐dependent transmission rate allows us to define a time‐dependent version (the basic reproduction number) as follows: This index permits to evaluate the strength of contagious over a temporal window comparing with the removal rate . To constraint between 0 and 1, in the estimation process, we reparametrize as =.

The system (2) is a system of nonlinear ODEs, which must be solved numerically. In this article, we use the ODE solver lsode (Hindmarsh, 1983) as it has implemented in the R package deSolve. If we suppose that and are known parameters, the (numerical) solutions , , and , an depends on the (vector of) parameters .

STATISTICAL INFERENCE

Different agencies in the world that daily update and publish data sets of epidemic data that contain at least three time series: the total number of infected, the number of dead, and the number of recovered (see Section 2 for more details). We derive from these time series the daily number of current positive cases and the daily number of new positive cases , recorded at day , .

Usually, the time series are supposed to be realizations of a stochastic version of the compartmental models. The different versions can be broadly classified into continuous models and discrete models. In the first group, fall the continuous‐time Markov chains (CTMCs) and the stochastic differential equations (SDEs) (Allen, 2008). In the second group, a discrete‐time approximation to the stochastic continuous‐time model is considered (Lekone & Finkenstädt, 2006). There exists an extensive literature on calibrating the stochastic models against time series with different inferential approaches (Andersson & Britton, 2012; Dureau et al., 2013; Finkenstädt & Grenfell, 2000; Hooker et al., 2011; Ionides, Bretó, & King, 2006).

Instead in this article, we follow the simplest idea that the solutions of the system (2) are actually the expectations at days of as many counting random variables. More precisely, we model the observed counts as Then the estimate of the parameter already defined is obtained by maximizing the independence log‐likelihood (Chandler & Bate, 2007)

Note that is not a “true” log‐likelihood but an instance of a composite likelihood (Lindsay, 1988) since it does not seem reasonable to assume that and are mutually and temporally independent. However, even though the model is not correctly specified, the maximum composite likelihood estimator, , is still a consistent and asymptotically Gaussian estimator with asymptotic variance under mild conditions (Chandler & Bate, 2007; Jacod & Sørensen, 2018).

The variance can be estimated by the sandwich estimator . The “bread” matrix is given by with . In the presence of time dependence, the “meat” matrix is given by the heteroskedasticity and autocorrelation consistent (HAC) estimator where is a vector of weights (Andrews, 1991).

With the aim of forecasting the spread of the epidemic outside the observed period, the number and the positions of knots in (3) play a crucial role. The higher the number of nodes, the less smooth the function . In this way, however, there is a risk of overfitting the data. On the other hand, the trend of the outside the observation time interval is mainly determined by the basis functions corresponding to the boundary knots.

We select the number of knots by maximizing the composite likelihood information criterion (CLIC) (Varin & Vidoni, 2005) The criterion has a strong analogy with the Akaike information criterion (AIC). In fact, measures the goodness‐of‐fit similarly to the log‐likelihood and the penalty tr() reduces to if the model (4) is correctly specified, i.e., if Equation (5) is the “true” log‐likelihood.

We could locate the internal knots to reflect policy interventions. However, it is very difficult to hypothesize the immediate effects of these policies and a simpler choice has been to place temporally equally spaced nodes. As for the boundary knots, it was chosen to place them at the beginning of the period and one week after the last observation available to obtain more stable estimates in the forecast period.

RESULTS

The epidemic outbreak showed different patterns in the selected time window: the reported SARS‐CoV‐2 cases in the United States were rapidly increasing, reaching a peak and subsequent stabilization of the number of daily new cases with an incidence of about 10 cases 100.000 inhabitants; in Italy, a drop up to 2.5 daily new SARS‐CoV‐2 cases 100.000 was reported, after an initial growth that reached a peak of incidence of about 10 cases 100.000 people similar to the United States; in Iceland, the incidence of new SARS‐CoV‐2 cases knew a huge peak ( 25 cases 100.000 people), then a decreasing trend and finally a limited number of new cases in the last considered day (Fig. 1). At the end of the temporal window, prevalence of the disease was quite different among the three considered countries: we registered 250, 180, and 45 current SARS‐CoV‐2 cases 100.000 people in the United States, Italy, and Iceland, respectively.

Fig 1

Daily new (solid black line) and current (dashed red line) cases of SARS‐CoV‐2 in (a) the United States, (b) Italy, and (c) Iceland.

In the literature, the incubation duration of the SARS‐CoV‐2 was estimated as (Wang et al., 2020), and therefore, we set the specific parameter .

The overall mortality rate was calculated as where the life expectancy (LE) is 78.5 years in the United States, 83.2 years in Italy, and 82.2 years in Iceland, respectively. The total population (N) in 2020 was 329.23 million (United States), 60.32 million (Italy), and 0.36 million (Iceland) of inhabitants. The starting values , , , and for the numerical resolution of the system (2) were set as follows:

, i.e., the number of currently infected on the first day of the data set (United States: 142, Italy: 155, Iceland: 1);

equal to the number of currently recovered on the first day of the data set (United States: 7, Italy: 0, Iceland: 0);

where is the number of new infected on the first day of the data set (United States: 68, Italy: 66, Iceland: 2);

.

We have tried several values for the number of basis function , i.e., from 3 to 8, and we found that and minimize the value of CLIC for Italy/United States and Iceland, respectively.

The estimate of (see Fig. 2) showed an overall decreasing pattern of the transmission rate across the selected countries. In particular, in the United States, the estimate of reached a peak close to 0.8 after about 10 days by the beginning of the epidemic outbreak, denoting an uncontrolled situation, moving to values of approximately 0.1 after about 45 days of the epidemic, with a predicted scenario of a slightly decreasing trend and a great amount of uncertainty. In Italy, the estimate of was moving from initial values of 0.75 to value close to 0.05 after about 50 days of observations. The estimate of was lower than those reported for the United States. The estimate of in Iceland showed a fast decreasing trend from values a bit over 1.0 to about 0 at day 40. The 30‐days prediction for is practically zero, denoting the end of the current phase of the epidemic.

Fig 2

The estimate of curves and 95% confidence bands (in gray color) for the United States, Italy, and Iceland. The dashed line represents the 30‐day predicted evolution.

The estimates of ranged from a low rate in the United States and Italy, 0.012 and 0.025, respectively, to a higher rate in Iceland, 0.080 (Table I). Then the removal duration, i.e., the reciprocal of , was estimated at 85.0 days (95% CI: 65.5–110.5) for the United States, at 40.2 days (95% CI: 37.6–43.1) for Italy, and at 12.5 days (95% CI: 9.9–15.9) for Iceland.

The model fitting was deemed satisfactory (Figs. 3 and 4) both with respect to the number of new cases and to the cumulative positive cases. Our major findings were: in the United States, the current positive cases were going to increase, reaching a probable maximum after the window of the next 30 days. In Italy, the epidemic outbreak had known its maximum in the number of positive cases around April 20th, and the tendency was for a slight decline. In Iceland, the peak of positive cases was registered on April 10th, associated with a rapid decreasing phase and a low number of new cases in the last observed days; in this case, the SARS‐CoV‐2 epidemic was going to be overcome approximately at the end of May.

Fig 3

The expected number of exposed, infected, and recovered subjects for the United States, Italy, and Iceland, based on the model parameter estimates. The dotted points indicate the observed number of infected cases.

Fig 4

The expected number of new infected cases for the United States, Italy, and Iceland, based on the model parameter estimates. The dotted points indicate the observed number of new infected cases. The dashed line indicates the last day observed.

The estimated trend for appears quite different among the selected countries (see Fig. 5). The value was reached on different dates in Iceland (March 28th) and in Italy (April 16th), while in the United States is expected to be achieved only in the end of May.

Fig 5

Estimated values for the United States, Italy, and Iceland and predicted evolution. The Y‐axis is in the log‐scale. The dashed lines indicate and the last day observed, respectively.

DISCUSSION

International health organizations recommend to implement public health and social measures to slow or stop the spread of SARS‐CoV‐2, reaching the full engagement of all members of society (World Health Organization, 2020). Countries have adopted different public health and social measures depending on the local specific historical evolution of the SARS‐CoV‐2 pandemic and on their health system capacity. Our analysis considers data of three countries, the United States, Italy, and Iceland, which have, on one side, different geographic and demographic characteristics and, on the other one, as many dissimilar approaches in terms of public health policies and restrictive measures concerning the ongoing epidemic.

Our proposal allows us to estimate an epidemiological SEIR model with a time‐varying transmission rate () with the scope to assess the timeline and the strength of the effects produced by the adopted restrictive measures. The removal rate was estimated, considering both two different time series (daily new and current positive cases). We avoided considering estimates of clinical SARS‐CoV‐2 recovery rate and specific mortality rate calculated by others, with the scope to comment on the knowledge provided by the analyzed data on the removal rate estimated by our model.

In the United States, the transmission rate of SARS‐CoV‐2 was very high at the beginning of the considered temporal window, and its reduction appears to be later and slower in comparison with those observed in Italy and even more in Iceland. This difference may be viewed as a result of the approach adopted by the United States in the epidemic onset based only on a limited and nonhomogeneous containment measures (Parodi & Liu, 2020).

The adoption of this strategy in the United States was a decision of particular importance since the COVID‐19 onset began a few days later than Italy, where, on March 9, 2020, the Italian Government set Europe's first nationwide restriction on movement due to the incoming SARS‐CoV‐2 epidemic. The estimates on the Italian transmission rate confirm the control of the epidemic wave after approximately 20 days of home restrictions, but with a high mortality toll in comparison with the preceding Chinese epidemic (Rubino, Kelvin, Bermejo‐Martin, & Kelvin, 2020). A social distancing and passive testing of symptomatic cases was the Italian strategy to contain the epidemic. Positive cases with few symptoms were confined in home isolation. However, there was a consistent amount of asymptomatic, which remained undetected, contributing to spread the epidemic (Flaxman et al., 2020; Lavezzo et al., 2020).

Iceland has the advantage that the epidemic outbreak started later than Italy; we observed that the Icelandic transmission rate quickly moved to values close to 0 after only 15 days of the restrictive measures (Gudbjartsson et al., 2020). These results were mainly attributable to an active searching strategy of asymptomatic positive cases organized by the National Health Service, which lead to be tested about 6% of the Iceland population at the date of April 2, 2020. However, the presence of a free voluntary private screening program estimated that the fraction of undetected infections by the Icelandic health service ranged from 88.7% to 93.6% (Stock, Aspelund, Droste, & Walker, 2020).

The comparison between the United States, Italy, and Iceland was certainly affected by a regional/state variation to COVID‐19 response. However, this administrative‐level granularity plays a role in the diffuseness of nonpharmacological health measures in the first phase of an epidemic outbreak. A faster and more aggressive response of each local administrative unit helps to contain the contagious spread and to achieve a relative control of the disease (Lancet, 2020). Despite fears of the negative consequences on their economy, Italy and Iceland experienced a contagion control in a relative short period. With a rapid and structured stay‐at‐home order and an assertive infection control measures, Iceland reduced the required time to flatten the curve.

The estimated values for reflect both the locally adopted swab policy and the specific phase of the epidemic wave: in fact, the active monitoring in Iceland provides a reliable value for the removal rate that is of about 12 days in line with that measured in China (14 days, Wang et al., 2020). In Italy, the controlling strategy implies that after a first positive swab test, a control swab will be repeated after a period of home isolation and this fact implies a longer time to obtain the healing confirmation.

In the United States, the low number of removed subject in comparison with the high increase in the incidence makes challenging a reliable estimation of the removal rate.

The proposed model has become a standard approach to estimate the transmission rate in a dynamic context (Godio, Pace, & Vergnano, 2020; Hong & Li, 2020; Loli Piccolomini & Zama, 2020; Petropoulos & Makridakis, 2020; Wu et al., 2020; Zhang et al., 2020). The model has the double scope of having a real‐time monitoring and of supplying possible evolution scenarios. The estimated fluctuations of were driven by gradual changes in the behavior of the population at risk as a consequence of the adopted restrictions. With respect to other specifications, our approach has the advantage of employing a basis of splines that allow us a high grade of flexibility for the estimation.

We estimated the parameters for and through a composite likelihood considering the information provided by both the occurrence of new SARS‐CoV‐2 cases and the current positive cases; in order to cope with the possible presence of heteroscedasticity and autocorrelation in the data, we estimated consistent standard errors combining a sandwich variance estimator and an HAC correction. Even though the model formulation has some ancestors, our proposal differs in two aspects from the aforementioned literature: we allow the data to indicate the shape of the function using a semiparametric approach. This feature could help to identify the best health intervention policy in a country; HAC variance estimators permit to reduce the bias allows to correct the underestimation of the variance of the estimator and therefore to produce future scenarios with a more appropriate margin of uncertainty.

There are several other limitations to our analysis. We used plausible biological SARS‐CoV‐2 parameters for the SEIR model based on updated numbers (i.e., ), but these values may be refined as more comprehensive data become available. The predicted values for are valid only in the absence of future changes to the restrictions, which is not likely to happen if an intermittent social distancing measures will be adopted (Ferguson et al., 2020).

Our results point out that the transmission rate in the United States, Italy, and Iceland showed a decline after the introduction of restriction measures. Despite this common trend, some differences in terms of timeline and impact are present. In particular, U.S. experts argue that more helpful tools are needed in order to reach the control of the epidemic wave (Parmet & Sinha, 2020).

The adoption of restrictive measures results in flatten epidemic curves and thus the distribution of the SARS‐CoV‐2 cases in a more extended period, with respect to an uncontrolled epidemic outbreak. In the absence of a specific vaccine, the high number of susceptibles and the relaxation of restrictions taken represent a cause of future outbreaks.

Table I

Estimated Values of the Parameter and Its Relative 95% Confidence Interval

Country	γ	(95% CI)
United States	0.012	[0.009–0.015]
Italy	0.025	[0.023–0.027]
Iceland	0.080	[0.063–0.101]