Abstract
As the COVID-19 epidemic rages, the total number of cases and deaths is raising across the world. This is certainly a challenging situation particularly for countries in a perilous socio-economic situation. With the rapid evolution of the disease, once introduced, making decisions becomes challenging particularly when past experience or data are not so relevant. This is when mathematical models could come handy to predict the evaluations of the epidemic. We designed a modified SEIR-type mathematical model which we assessed with up-to-date data of COVID-19 cases and deaths around the world, and in which the rate of transmission is derived from the observed delayed fatality rate. We used Brazil as an example for middle income countries and applied our model to the Brazilian context. We implement age stratification and a horizontal confinement followed by a “vertical confinement” exit strategy.
Introduction
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) outbreak has been ongoing for 5 months now [1]. Since it was first reported in Dec 2019 in China [2], the virus rapidly made its way to other parts of the world taking pandemic proportions[3]. The number of cases and deaths exponentially incremented reaching a total of 1.5 million confirmed cases and 88 thousand deaths in early April 2020. Recent diseases outbreaks that spilled over from animals such as Ebola [4, 5] or Avian Influenza [6] have been described as specific to developing countries. COVID-19 has been breaking this myth as the virus has been particularly exceptional at breaching in inside all developing countries and challenging their health system. In Europe, Italy has been particularly affected. With 140 thousand cases, the Italian national health system’s has been struggling to effectively respond to the exponentially increasing flow of patients in need for intensive care [7]. The United States recently surpassed China in total number of cases (420 thousand), becoming a particular hot bed in this phase of the pandemics [8]. By the time this article is published, there will not be a place on earth where the virus did not cause any damage. West African countries such as Sierra Leone just reported their first cases [9] and catastrophic scenario similar to the 2016 Ebola outbreak is possible.
The threats of COVID-19 on countries that started to count cases prompted us to develop a model to describe the evolution of the epidemic and its effects the health care system. Mathematical models are a powerful tool that proved their importance in previous epidemiological disasters such as the Ebola virus [10, 11], smallpox [12], or influenza [13], contributing to the understanding of the dynamics of disease and providing useful predictions about the potential transmission of a disease and the effectiveness of possible control measures, which can provide valuable information for public health policy makers [14]. SIR-type models, also known as Kermack-McKendrick model [15], consists of a set of differential equations and has been applied to a variety of infectious diseases. Although considered as simple, SIR models have been of great help to stop epidemics in the past such as putting in place effective vaccination protocols [16].
Here we develop an SIR type compartmental models for COVID-19 including both symptomatic and asymptomatic, quarantined, and hospitalized while taking into consideration differences by age groups. We also analysed at the effect of confinement during a specific period of time. Contrary to similar epidemiological models, the proposed SIR model is initiated by the first confirmed COVID-19 death. Numerical simulations of the deterministic models are compared with real numbers of the ongoing outbreak in different countries. Moreover, the deterministic framework in which we operate greatly simplifies model analysis and allows a more thorough comparison of the various intervention strategies.
We focus on this paper on the case of Brazil, where the pandemics counts 16 000 confirmed and cases and 800 fatalities at the time of writing. The country has 35 682 ICU beds according to government data of Feb 2020 [17]. The first official B
SARS-CoV-2 case in Brazil was reported in in São Paulo on February 26th and the first official SARS-CoV-2 case was confirmed on March 19th.Shortly after, the government ordered a general lockdown first on Rio de Janeiro March 22nd, then on other regional urban centers (such as?). There are no clear estimate on the percentage population that is actually in confinement. However, it is estimated to be around is 56% according to satellite data 1,
Given socio-economic consequences of a lockdown, particularly on a middle income country, decision makers are considering a vertical confinement as an exit strategy to the regular lockdown. Vertical confinement is understood as reducing contact to a specific age group that is more at risk of contracting and developing SARS-CoV-2 [18], as opposed to horizontal (or general) confinement that does not discriminate between age groups. Throughout this manuscript, we will present the model which we validated with data on other countries. We then applied the model to the specific SARS-CoV-2 scenario in Brazil. We finally tested the effect of both general confinement and the vertical one on the epidemic curve.
Needs a global conclusion on application of the model to other low middle income countries and how this could help decision makers make efficient decisions to stop causalities.
The model
We used a modified version of an SIR-type deterministic compartmental model to trace COVID-19 epidemic evolution in an isolated population of N individuals2. We assumed that a population could be subdivided into the following compartments:
Susceptible (S): COVID-19 naive individuals,
Confined (C): subset of susceptibles removed from the epidemics (by e.g. social distancing).
Exposed (E): Susceptible that have been exposed to infective individuals,
Asymptomatic (A): Infected and infective but showing mild or no symptoms
Symptomatic (I): Infected and infective but showing symptoms described in the literature,
Quarantined (Q): Symptomatic that are not infective,
Hospitalized (H) Symptomatic, not infective, who are being treated,
Removed (R) People removed from the epidemic dynamics by recovering or passing away.
We split that the population in subcategories by age (range, 0-10, 10-20, 20-30, 30-40, 40-50, 50-60, 60-70, 70-80, and 80+ years old) and that rates should vary with age [18].
Taking into consideration the 8 compartments and the 9 age groups, the model is described by a set of 72 coupled non-linear equations:
For each compartment X the age sub-bins add up to X = ∑i Xi and compartments are such that S + C + E + A + I + Q + H + R = N, with N = ∑i Xi being the total population; Ni is the population in each age bin.
Eqs. (1)–(8) describe a compartmentalization of the population and the flow between the compartments. Contact with infected individuals removes a fraction of the susceptible (S) population at a rate given by λ, referred to as infection force, making them exposed (E) to (SARS-CoV-2). Exposed (E) becomes infectious at the rate σ; a fraction p of them becoming symptomatic (I) and (1 – p) asymptomatic (A). The symptomatic (I) are removed from the infective force and become quarantined (Q)at a rate γ, the asymptomatic (A) are removed at a rate θ, a fraction w of them going in remission (a fraction (1 – w) becomes symptomatic. A fraction qi of the quarantined (Q) are hospitalized at a rate ξ. The hospitalized (H) are removed at a rate η. The average fatality rate is µi.
The timescales corresponding to σ, γ, θ, ξ, and η are the incubation time tσ = σ-1 the infectious interval tγ = γ-1, the remission time tθ = θ-1, the time to hospitalization tξ = ξ-1, and the average length of hospital stay tη = η-1.
The infection force is driven by the infected, both symptomatic (I) and asymptomatic (A) where we use the shorthand notation and β is the infection rate, related to the reproduction number ℛ(t) via
Lock-down consists of having a fraction of the susceptible population removed from the epidemic dynamic by moving them from Si to Ci at a rate ψi. Similarly, lifting the lock-down is done by placing Ci into Si at the rate ϕi. We consider these functions to be Dirac deltas where tlock and tlift are the time (in days) of lock-down and of lifting of the lock-down, respectively. To allow for partial demographic lock-downs (e.g., 80% lock-down of the population are in complete lock-down), ai and bi are allowed to vary by age (e.g., 80% lock-down of the 40’s age group population are in complete lock-down). The flow chart between compartments is shown in Fig. 1.
Other diagnostic quantities are the numbers Ui of people in need of an intensive care unit (ICU) bed where ζi is the fraction of hospitalized patients that need critical care. Both ζi and the hospitalization fraction qi are age-stratified.
For integration, we use a standard Runge-Kutta algorithm, with timesteps
Model validation
In this section we present details about the model validation and strategies to verify characteristic timescales and other parameters.
Model fit to the 2020 COVID-19 epidemic
We consider the susceptible population (S) as the total population of a country since at the onset of outbreak no one is immune to the virus yet. Model parameters, shown in Table 1, were based on previous knowledge of Coronaviruses and early reports and research on COVID-19 [19]. The age-dependent parameters fatality rate µi, fraction of infectious that are hospitalized qi, and fraction of hospitalized that need critical care are shown in Table 2.
Table notes: Values taken from [18].
Because all these timescales are much smaller than a human lifetime, aging of the population is ignored and no upward flow between the age sub-compartments (i → i + 1) is considered. Population pyramids are taken from UN data 3, and split into the pre-defined age bins.
We derive ℛ (t) from the available statistics since knowledge on the real number of infected is not clear. The most reliable indicator in this situation is the number of deaths. Given a fatality rate µ and an average time τ between exposure and death, the number of dead at a time t + τ will equal the fatality rate times the number of people that got exposed at time t. Assuming that confinement dynamics do not play a role (although it is trivial to include it), the equation is the following:
Taking the continuous limit and substituting Eq. (1) where we also write tr = t + τ for the retarded time. Summing over all age bins D = ∑ i Di we have the cumulative death rate on the LHS, which is an observable and ⟨µS⟩ = ∑ i µiSi. We can then substitute Eq. (9) and solve for (t) as a function of time
Since death occurs an average of τ days after infections, we setup the integration to start τ days before the first reported COVID-19 death, i.e., t = 0 means tr = τ. Finally, to start the integration we need to define the initial number of exposed individuals.
This should be where t0 is the time of the first death and is the age-weighted fatality rate. According to current knowledge of the epidemics, τ ≈ 14 days [18].
We compared our model predictions with official data on cases and deaths for multiple countries, as tracked by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University 4. We plot in the left panel of Fig. 3 the fatality rate for a number of countries, which corresponds to the left hand side of Eq. (18). We apply Eq. (19) to convert this data into ℛ (t), feeding this value into Eq. (1)-Eq. (8) to start the SEIR evolution. The populations I(t) and S(t) that enter in Eq. (19) are then calculated to update ℛ (t). The resulting values are plotted in the right-hand-side of Fig. 3.
The timescales σ, γ, θ, and ξ, as well as the fractions p and w, are found by Markov chain Monte Carlo (MCMC) fitting, with priors as given in Table 1 and explained in the Supporting Information (Markov Chain Monte Carlo).
Finally, we compare the cumulative number of hospitalizations calculated from our model with the number of confirmed COVID-19 cases. For a country that is not doing massive testing and only reporting COVID-19 as acute cases reach the hospital, these curves should match reasonably well.
Results and Discussion
Brazil epidemic scenario
Fig. 4 represents the modeled epidemic scenario in Brazil over 400 days time (checked a posteriori as the time when the number of symptomatic cases ends with a 70% horizontal lockdown). Parameters determined by the MCMC modeling are shown in Fig. 5, being , and . ℛ (t) at present is hovering around 2.
Fig. 4a shows the evolution of the compartments of exposed (E), asymptomatic (A), symptomatic (I), and hospitalized (H), in linear scale. Fig. 4b shows the same curve of H but also the fraction of hospitalizations needing ICU (U), in log scale. The epidemic is starting at March 1st and number of symptomatic is predicted to end at July 1st.
Table notes: Values taken from [18].
The peak of symptomatics is predicted for May 17th with 20 million symptomatics. Consequently, there is a predicted rise in number of hospitalized, reaching saturation on May 3rd and peaking at May 22nd with 106 hospitalized. ICU beds will reach saturation on May 3rd, when the ≈ 35 thousand ICU beds in Brazil are occupied (since the estimate assumes that all ICU beds should be occupied with coronavirus patients, which is not realistic, the collapse should in fact happen sooner). Demand for ICU will get higher until reaching a peak on May 22nd with 300 000 patients. The cumulative number of deaths June 1st is 106.
Fig. 4 contrasts the predicted cumulative numbers of infected persons (orange line), hospitalized persons (blue line), and deaths (black line). In this figure, the cumulative number of confirmed cases (yellow dots) and deaths (black dots). The cumulative number of hospitalized is very close to the actual confirmed cases. This is expected as Brazil is not doing testing on a massive scale.
Horizontal Lockdown
In Fig. 6 we check the effect of horizontal confinement, defined as equal percentage of the population confined at any age bin. There is a change in the epidemic dynamic when horizontal confinement is applied in different rates Fig. 6 The plots show (a) the number of hospitalizations, (b) the number of ICU cases, and (c) the number of fatalities, as a function of the degree of social distancing. Confinement was implemented at time t = 0 corresponding to March 22 when the first measurement of social distancing was implemented. To not overwhelm the health care system capacity (≈ 4 × 104) ICU beds, the level of social distancing should be over 70%. As mentioned in the introduction, estimates are that Brazil is maintaining 56% (with state to state variation from a maximum of 64.7% to a minimum of 53.7%). At this low level of distancing, capacity should be reach in less than 50 days, which is in agreement with the dynamical R(t) model in Fig. 4.
Vertical Lockdown
We vary now the degree of confinement by age bin, characterizing the vertical confinement. Fig. 7 shows the number of hospitalizations in a model where confinement was implemented, broken down by age bins. The upper plots show horizontal confinement with different proportions of the population (same as Fig. 6 but broken down by age and in logarithmic scale). Confinement was implemented at the same time as in Fig. 6. The other rows explore vertical confinement. In the second column 60% of the population under 40 is confined, but the population older than 40 is confined to a higher degree, at 90% (solid blue line) and 99% (dashed blue line). The cyan line marks the same model as the upper plots, where 60% of the population is confined, irrespective of age. The 3rd, 4th, and 5th row of plots show the same analysis but confining 60% up 50, 60, and 70 years old, respectively. As seen in the cyan line, the number of hospitalized rises from 30-60 and falls for 70 onwards. That is because even though 70+ are more likely to be hospitalized, the number of 30-60 is much higher in the population.
Fig. 8 shows the same results for the fraction of hospitalized that needs ICU. Fig. 9 shows results from the same suite of models but for the number of fatalities. For the number of ICU cases, there is no significant difference past age 60, with only a minor uptick at the 70-80 age range. Collapse of health care system can be avoided if vertical confinement is instored on people who are 60 or older, but at the expense of a significant number of extra ICU cases for the 50-60 age bin. At 60% confinement, hundred of thousands of deaths are seen in the 60-70, 70-80, and 80+ age bins. The number lowers to 50 000 in the 90% confinement. As noted before, vertical confinement for 60 years old and older, leads to a significant number of deaths for the 50-60 age bin (over 50 000). Vertical confinement at 50 years old leads to much lower death rate for this age segment.
Finally, we look at vertical confinement as an exit strategy. In Fig. 10 we model a release from lockdown on May 1st, according to two scenarios: full release for the population under 50 (dashed line) and full release for the population under 60 (solid line). The population past this age is kept at 90% confinement. The upper plots show the susceptible (S) and confined compartments (C), normalized by the number of individuals in the respective age bin. The second row from top to bottom shows the number of hospitalizations, the third row the number of ICU cases, and the last row the cumulative number of fatalities. As the population is released from global confinement, the number of H/U /D peaks at 400 000/50 000/120 000 in the 50-60 age bin alone, that bears the lion’s share of morbidity. Keeping the 50-60 age population in 90% confinement lowers the statistics significantly, with the health care system at capacity, and the number of deaths per age bin about 25 000, with 60+ years olds having the same fatalities as 40-50 years old.
Limitations
As in any setting,the outbreak response strategy plays a crucial role in quality of the outputs models can give.Since the identification of first case, In Brazilian response strategy has been changing over time. At first, only international travelers admitted to hospitals had access to SARS-CoV-2 testing. Now there are diagnostic clinics and universities involved in COVID-19 testing, but there is no national massive testing strategy in place. Besides, each Brazilian state has the authority to put in place their own strategy to address the epidemic. The states of São Paulo and Rio de Janeiro, containing the largest metropoles in the country, adopted larger strategies of isolation with schools and stores closed early on while similar strategies had not yet been adopted in other states. Bottom line, the resulting morbidity and mortality rates can change significantly, resulting in dramatically different output numbers as the number of infected people or the number of hospital beds needed. It is necessary to have massive testing strategy in place to have better prediction accuracy of the models.
Our model estimate hundreds of thousands of infected people in Brazil at April 1st. This is more than the number of expected cases in the country while we write this article, considering the estimated under notification of cases [20] and do nothing to control the infection. It is possible that the actual number be less than that although it is also important to notice Brazil has not done a real lockdown so far.
The model also ignores mobility, in the sense that it does consider travel to and from the country. Given that right now we are the stage of community tranmission, this limitation should not be of significance to the results.
Conversely, and more importantly, the model assumes that the confined population is completely safe from infection, whereas in reality a vertical lockdown may not be feasible to implement as the elderly are not adequately distanced from the younger in their family and/or social circle, and infection cannot be avoided if the younger is exposed to COVID-19.
Finally, the analysis assumes that the data on fatalities is accurate. Underreported deaths should lead to an unknown source of error in the present study. Also, the MCMC produces error bars in the parameters that we did not take into account in the forward modeling.
Conclusion
Is this study we examine the strategy of vertical confinement as currently debated in Brazil. Since the fatality rate of COVID-19 is disproportionate among the elderly, the federal government has suggested confining only the at-risk age groups. Such a strategy would limit the economic impact of the pandemic while at the same time minimizing the risk of collapse of the health care system and number of fatalities. State governments have by and large resisted the idea, and tried in turn to enforce horizontal lockdowns, as also suggested by the majority of epidemiologists.
At the current rate of advance of the pandemics, Brazil should face collapse of the health care system by May 15, with 300 000 ICU beds needed (10 times more than the current capacity), and 106 fatalities. A decrease in the rate of confirmed cases is seen with respect to the rate of fatalities, which is indicative of the effect of the lockdown. A 60% lockdown reduces the number of deaths to 400 000 due to COVID-19, still not avoiding overwhelming of the health care system. An increase in lockdown to 70% is needed to avoid collapse.
An exit strategy of vertical confinement of individuals older than 60 year old by May 1st would see a second wave disproportionally affect the 50-60 age bin. The ICU cases in this age range alone would bring the health care system to collapse and result in over 100 000 deaths. If vertical confinement is contemplated, it should be of the population over 50 years old. However, this age range, 50-60, is also a part of the workforce, and thus defeats the purpose of a vertical confinement. Moreover, we emphasize that our model assumes an idealized lockdown where the confined are perfectly insulated from contamination, while in reality there would be several practical barriers to it as the confined elderly would depend on the young for most essential activities, and a perfect lockdown would not be achieved in a multi-generational household, especially in close quarters such as found in the low and even middle income neighborhoods common in Brazil. Our results therefore discourage vertical confinement as an exit strategy and point toward enforcing a higher degree of horizontal confinement than currently practiced. We urge Brazilian authorities to take action to prevent virus dissemination in the critical coming weeks.
Supporting information
Markov Chain Monte Carlo
To fit the best value to w and p, and to better constrain σ-1, γ-1, θ-1, ξ-1, we use the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) [21] to sample the parameter space around the solutions and evaluation of the parameter uncertainties. For the priors input, we use the values taken from the [18]. To search for the minimization of cumulative hospitalization Hc, we generated a cumulative error Cerr on the reported confirmed cases Cc.
As the JHU-CSSE reports on the confirmed cases are given daily with some fluctuations, we need to take this into account while weighing all solutions by adding a 1-day error matrix together with the confirmed cases (being conservative). In an ideal scenario, the cumulative number of hospitalization would be the same as the number of confirmed cases. In real life, not all confirmed cases are hospitalized so we do not expect to fit the Hc with Cc. Rather, we weigh the Cc array with the Hc array using:
is the weighed cumulative hospitalizations and n is the length of the data.
Following we get the residual between Cc and , and we used the negative binomial distribution to calculate each likelihood [22]:
Equations 21 to 26 are implemented in the likelihood function on the code. If in each run it returns a finite number, the algorithm parses the result, if not it returns a large number (1020) to discard as a bad fit.
We limit each parameter using a range cutoff in when feeding the probability function to restrict parameter space. That way, we do not run models with unrealistic physical parameters (such as e.g. symptomatic going to the hospital in –2 days), and also constrain the known range for the other parameters. The MCMC function feeds on 6 free parameters, 4 fixed parameters and 2 predetermined arrays as presented in Table 3.
Table notes: Priors taken from [18].
Data Availability
The code developed for the present work is public and made available at https://github.com/wlyra/covid19
Footnotes
2 The software is written in python 3.7, and is made public at https://github.com/wlyra/covid19