Multinational modeling of SARS-CoV-2 spreading dynamics: Insights on the heterogeneity of COVID-19 transmission and its potential healthcare burden

Data sources

We used the time series of the COVID-19 reported by the Johns Hopkins Coronavirus Resource Center (hereafter JHCRC), starting from January 21 to April 14, 2020.^17,18 We selected 18 countries corresponding to ∼10% of the total number of countries currently reported in the JHCRC database. Such a particular choice was made based on the consistency of the data with the SEIR model as evaluated by R², which is used as a metric of the goodness of fit. We only take into account countries with a number of confirmed cases larger than one per 100k and a number of deaths exceeding 200. Country-level data on population size was extracted from the Department of Economic and Social Affairs of the United Nations.¹⁹ To evaluate the burden that the pandemic lays on the country-level healthcare system, we used the number of hospital beds reported by the World Bank Database.²⁰

The SEIR model

The total population N is considered to be constant and is partitioned into four disjoint groups, namely: susceptible, S(t), exposed, E(t), infected, I(t), and recovered, R(t), which are all functions of time. Individuals who have had contact with an infected person are not themselves instantaneously infectious, but rather move to the exposed compartment E(t) where they develop infectious symptoms at an incubation rate σ.²¹ A set of coupled ordinary differential equations (ODEs) describing the spreading of the disease and the evolution of the partitions can define the SEIR model as follows: subject to the restriction where s(t) = S(t)/N, e(t) = E(t)/N, i(t) = I(t)/N, r(t) = R(t)/N, the coefficients to be determined β ≥ 0, γ ≥ 0 and σ ≥ 0 are the effective contagious contact rate of spread, the removal (including deaths and recovered patients) rate, and the incubation rate, respectively. The time t is taken in days. A local epidemic scenario occurs when the number of the first infected patients increases, namely,. The basic reproductive number, R = β/γ, commonly gives the number of secondary infections generated by the first infected individual.

We solve the coupled system using SciPy’s implementation of the explicit Runge-Kutta method of order 8 (DOP853).²² This methodology enables us to estimate β and γ, as well as the initially exposed population E₀ ≡ E(t = 0) by fitting the progression of R(t) and I(t). We use the trust-reflective non-linear least-squares method to minimize the residuals between the SEIR model and the data from the time series. The values R(t) and I(t) were taken as linear combinations of the number of confirmed cases C(t), confirmed recoveries R_r(t) and confirmed deaths R_d(t) with R(t) = R_r(t) + R_d(t)and I(t) = C(t) − R(t). For purposes of following the progression at different stages of the epidemic, the values I_n = I(t = t_n), R_n = R(t = t_n) and E_n = E(t = t_n) can be used as inputs of the SEIR model to determine the evolution of the system.

Fitting SEIR parameters

We estimate the contagious contact rate (β) by fitting the data at three different stages of the COVID-19 spread. The parameter β is difficult to estimate accurately because the transmission rate is susceptible to changes related to government intervention. Hence, β can range several orders of magnitude, depending on the government response to contain the outbreak locally. Early reports from Wuhan estimated a value of 0·064 at the initial stage of the spread,²³ but after clinical support was adopted by the Chinese government and considered in the models, it increased to 0·322. On the other hand, the removal and incubation rates (γ and σ, respectively) are not directed by population controls; according to empirical data, it was shown that these rates could range in the order of days (from 2 to 14 days).²⁴ Because of this, we impose the following bounds on the possible values that these fitted parameters can take: 1/β ∈ [0,10], 1/γ ∈ [1 day, 1 year], E(t = 0) ∈ [0, N], and finally 1/σ= 5·1 days which is a robust estimation of the average incubation time.²⁵

Piecewise SEIR fitting and projection

We fitted I(t) and R(t) piecewise and optimized the interval ranges until the value of R²was higher than 0·9 (0·8 for the particular case of Mexico which is discussed in the following sections of this work). We defined the ratio of transmission rates from subsequent intervals (Q_n = β_n+1/ β_n), for n = 1,2,3 being the interval number. This estimate is useful to identify the changes in the COVID-19 spread dynamics upon changes in public health policies. Functional prediction bands pb(t) = f(x) ± δ(x) were obtained using δ(x) =(variance)×(inverse of the t-student cumulative distribution).

Model projections on healthcare burden

A fundamental question for decision-makers is whether the healthcare system will be able to cope with the increasing number of patients infected with SARS-CoV-2.^26–28 To add this practical perspective into our analysis, we considered the number of hospital beds per capita.²⁰ We split the projected number of COVID-19 into three categories according to their severity: (a) mild cases that do not require hospitalization accounting for 80% of the projected infections, (b) severe disease cases which do require hospitalization, accounting for 15% of the projected infections and (c) acute disease cases, which require admission to the intensive care unit, accounting for 5% of the number of projected infections.²⁹ We have performed projections for four different scenarios: one assuming the best fit of the SEIR model at the latest stage, and three scenarios where we assume a 10%, 50%, and 90% decrease in β. As a metric of the level of healthcare burden under each scenario, we define the healthcare stress as the ratio between the number of infections that require hospitalization at the peak of the disease and the number of available beds.²⁰ The lethality rates were used to assess the status of the COVID-19 outbreak.

Confirmed cases and lethality rates in different countries

The current officially reported lethality rate for COVID-19 is approximately 3·4%.¹ Our fitted data indicates a global lethality rate of about 6·87%. Since this estimation considers only confirmed cases, we argue that the former value could be an overestimation introduced by underreporting the actual number of infections. A cluster of European countries (UK, Italy, Spain and France) whose calculated lethality rate range from 10·3 to 12·7% and their number of confirmed cases exceeds 100 cases per 100k encompasses the most affected countries (figure 1). Germany, Austria, Switzerland and the US reported >100 cases per 100k and their lethality rate is similar or less than that of China (4·02%). South Korea, Singapore, Australia and Canada also can be clustered together, being countries with a lower lethality rate (< 3%) and confirmed cases less than 100 per 100k. Brazil, Ecuador and Iran have lethality rates higher than China and their number of confirmed cases is currently increasing (>10 per 100k). Finally, Japan and Mexico show a number of confirmed cases less than 10 per 100k. The lethality rate for Mexico (6·47%) is higher than that of China, however, in this study, Mexico shows the lowest number of confirmed cases per 100k habitants. Such a comparison is only a geographical indication that the epidemic has different dynamics of transmissibility.

• Download figure
• Open in new tab

Figure 1.

Number of confirmed COVID-19 cases per 100k habitants versus the lethality rate (in %) for different countries. The lethality rate was calculated as the quotient of deaths by COVID-19 divided by the confirmed cases according to the data published from JHCRC, on April 14, 2020.

Differences in SEIR model fitting across countries

In figure 2b, the time series of the active infections and removed individuals shows an effective change in both contagious contact and removal rates, leading to a crossover in February 15 with a transitory decline in the number of active infections. A second cross-point around March 21 is observed, when the number of active infections takes over the number of recoveries. These features can be interpreted to be due to a lag between the pandemic’s onsets across different regions of the world. Such timing differences introduce an additional heterogeneity due to local reporting and containment policies. To address the heterogeneity of the COVID-19 outbreak, the SEIR model is used to model the progression of COVID-19 in 18 countries (see Supplementary Material). The crossing point between the infected and removed curves is currently missing on all of the countries studied, with a clear exception for China, South Korea and Austria. Because of the inherently high effectiveness of the transmission of SARS-CoV-2, it is expected that the epidemic will extend for several months in some countries. According to the current trend, we project that the peak in the number of infections worldwide will take place after the third quarter of 2020 (figure 2c). We estimate that a 90% reduction in the contagious contact rate (0.1β) will be required to bring the peak of on the number of infections to a date within the second quarter of 2020 (figure 2d).

• Download figure
• Open in new tab

Figure 2.

Fit of the SEIR model to worldwide data. In all cases, the solid lines represent to our data fitted, while markers correspond to official data. The x-axis is the period of time in format mm/dd. Upper left: total infections, recoveries, deaths and infected cases as per official reports from JHCRC. Figure 2a shows the global picture of the progress of COVID-19 epidemic as reported by JHCRC, including the SEIR model fit and confidence intervals for both infected and removed (shadow intervals indicate the piecewise fitting for obtaining the effective infectious contact rate of spread for each period). In solid orange circles is the total number of confirmed cases per day C(t), the number of confirmed recoveries R_r(t) and deaths R_d(t) are indicated in open green squares and open red triangles respectively. The active number of infections at time t estimated as I(t) = C(t) − R_r(t) − R_d(r) is indicated by open blue circles. Figure 2b shows the selection of intervals used to fit the mode globally: (a) from the date of the first detected case to the point 1 indicated on the infections curve, (b) on the range indicated between the arrows labeled as 2 on the infections curve and (c) from the point labeled as 3 on the infections curve to April 14, 2020. Figure 2c indicates the projections accordingly the best-fit model for infections, removed and recovery cases. Lethality rate is shown for comparison. Figure 2d indicates the severity of the spread, in which solid lines represent splitting the projected curve into severity categories, whereas dashed lines represent the projection with 90% reduction of current contagious contact rate.

So far, European countries still lead in the number of confirmed cases and COVID-19 deaths per capita, followed by the United States (see Table 2). With the exception of Japan and Mexico, it seems that other countries studied have managed to slow down the spreading of the epidemic as revealed by lower β₃ and Q_1,2 < 1 values (figure 3). In the case of Mexico, the goodness of fit is limited (R² > 0.8) due to a different sampling methodology. Namely, as a response to AH1N1 2009 outbreak in Mexico, the healthcare system implemented a proper sentinel surveillance model, which samples only a few of all SARS-CoV-2 infections;³⁰ therefore, the effect of public policy measures for this country might not be accurately estimated using our predictions. Substantial declines in spreading rates were observed for Austria, China and South Korea. Notably, Austria was an early responder of the epidemic spread by enforcing more stringent restrictions, and thus it seems that is preparing to lift lockdown.³¹ Removal rates (associated to both: deaths and recoveries in figure 2c) appear to be increasing in all countries, which goes along with a rise in the estimations of lethality rates, and it might also be related with an effective increase in the recovery rates as available data on the progress of the epidemic becomes available.

• Download figure
• Open in new tab

Figure 3.

Comparison of the fitted values for the SEIR model for 18 countries. Figure 3a shows the total confirmed cases, recoveries and deaths per 100k, colored in orange, green and red, respectively. Figure 3b shows the fitted contagious contact rate (β) for each of the different intervals. Figure 3c shows the fitted removal rate (γ) for each of the different intervals. Figure 3d shows the ratio Q_n = β_n+1/β_n between the different intervals, which is proposed as a measure of the effectiveness of the containment measures. For panels (b-d) the missing values indicate that only two intervals were identified in the corresponding dataset. The scale of the vertical axis in (d) was adjusted to a maximum of only 100. Nonetheless, Mexico overpasses this limit drastically in Q₂, which could be attributed to a low early detection rate due to the implementation of a sentinel model for epidemiological surveillance (see Discussion and Table 1).

Two distinct values of Q are defined per country to reflect changes in their contagious contact rate for the three fitted intervals (figure 3d). Values of Q above the dashed line indicate an increase of β from interval n to interval n + 1, whereas values below that threshold indicate a reduction in the transmission rate. Q₁ gives a metric of how susceptible the contagious contact rate is to updates in the reported number of cases, as more information becomes available since the early stages of COVID-19 spread for each country. Except for Switzerland, Austria, Iran China and Japan, all other studied countries show Q₁ > 1. Overall, the projections show that Q₂ < 1 with the exception for Japan that recently reported a surge in new cases; and Mexico, which has recently entered the exponential phase of the epidemic.

Projections on the healthcare burden and effects of containment measures

Assuming that treatment and prevention of the disease do not have an immediate impact on the parameters of the model, our projections indicate that tighter measures should be taken to effectively reduce the spread of COVID-19. Our projections also suggest that stringent measures are required in most countries; with the only exceptions of Australia, Austria, China, Iran and South Korea, which could either continue with their ongoing policies or, up to some safe extent, relax their local interventions (figure 4). As expected, a monotonic decrease in the number of cases at the peak of the infection I_peak is observed with a decrease in the value of β (see figure 4b). For 12 out of the 18 countries analyzed (figure 4c), the current trend indicates that the number of severe infections will surpass the healthcare capacity (see details in Table 2). According to our projections, Germany, South Korea, Spain and Switzerland would avoid a saturation of the healthcare capacity if they will apply only with a 10% reduction in their current transmission rate (0.9β_latest). We observe that for four out of the 18 countries analyzed, a 50% reduction (0.5β_latest) is enough to avoid the saturation of the available health care infrastructure. We also see that a reduction of 90% (0.1β_latest) will be a failsafe intervention for all the analyzed countries. For such a level of intervention, it is expected that 11 out of 18 countries (not including Australia, Austria, China and South Korea, which seem to have already passed their infection peaks) will reach their peak before the third quarter of 2020. Additionally, it is observed that five out of the 18 analyzed countries will significantly delay the peak of the infection even at 90% reduction. For such cases, the course of the epidemic will be extended further beyond the 365-day simulation window. The countries that follow the former trend are Brazil, Ecuador, Japan, the United Kingdom and the United States.

• Download figure
• Open in new tab

Figure 4.

Projected infections for 18 countries using different values of the contagious contact rate (β) as related to the effectiveness of the containment. For all figures, cyan color is for the current value of β as estimated by our methodology, blue (0.9□), purple (0.5□), and magenta (0.1□). The period estimated until the local peak of the epidemic has been reached is shown in Figure 4a, in which the dashed line indicates the maximum simulation time. Figure 4b shows the total number of cases at the local peak of the epidemic. Figure 4c shows the healthcare stress, defined as the number of cases at the peak requiring hospitalization (20% of the total number of cases) divided by the number of available hospital beds as reported in the World Bank database.