Abstract
The ongoing COVID-19 epidemics poses a particular challenge to low and middle income countries, making some of them consider the strategy of “vertical confinement”. In this strategy, contact is reduced only to specific groups (like age groups) that are at increased risk of severe disease following SARS-CoV-2 infection. We aim to assess the feasibility of this scenario as an exit strategy for the current lockdown in terms of its ability to keep the number of cases under the health care system capacity. We developed a modified SEIR model, including confinement, asymptomatic transmission, quarantine and hospitalization. The population is subdivided into 9 age groups, resulting in a system of 72 coupled nonlinear differential equations. The rate of transmission is dynamic and derived from the observed delayed fatality rate; the parameters of the epidemics are derived with a Markov chain Monte Carlo algorithm. We used Brazil as an example of middle income country, but the results are easily generalizable to other countries considering a similar strategy. We find that starting from 60% horizontal confinement, an exit strategy on May 1st of confinement of individuals older than 60 years old and full release of the younger population results in 400 000 hospitalizations, 50 000 ICU cases, and 120 000 deaths in the 50-60 years old age group alone. The health care system avoids collapse if the 50-60 years old are also confined, but our model assumes an idealized lockdown where the confined are perfectly insulated from contamination, so our numbers are a conservative lower bound. Our results discourage confinement by age as an 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 increased reaching a total of 1.5 million confirmed cases and 88 thousand deaths in early April 2020. Recent disease 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 developed countries and challenging their health system. In Europe, Italy has been particularly affected.
With 140 thousand cases, the Italian national health system has been struggling to effectively respond to the exponentially increasing flow of patients in need of 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 likely 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 threat 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 on the health care system. Mathematical models are a powerful tool that proved important 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 containing simplifying assumptions, SIR models have been of great help on stopping epidemics in the past by e.g. informing effective vaccination protocols [16].
Here we develop a SIR type compartmental model for COVID-19 including both symptomatic and asymptomatic, quarantined, and hospitalized while taking into consideration differences by age groups. We also analysed 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.
In this work we focus on the case of Brazil, where the pandemics counts 16 000 confirmed 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 SARS-CoV-2 case in Brazil was confirmed in São Paulo on February 26th and the first official COVID-19 death was reported on March 19th. Shortly after, a lockdown was enacted first in Rio de Janeiro on March 22nd, then on other regional urban centers. There is no reliable measument of the percentage of the population that is currently in confinement; however, the number is estimated to be around 56% according to satellite data 1.
Given the 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. In the next section we will present the model, followed by validation. We then apply the model to the specific SARS-CoV-2 scenario in Brazil. Finally, we test the effect of both general and vertical confinement on the epidemic curve.
The model
We used a modified version of a 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 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 we consider that some flow rates between compartments 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 = ∑iXi and compartments are such that S + C + E + A + I + Q + H + R = N, with N =∑ i Ni 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. The exposed (E) become infectious at the rate σ; a fraction p of them becoming symptomatic (I) and a fraction (1 − p) becoming 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 and a fraction (1 − w) becoming 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, ai and bi are allowed to vary by age (e.g., 80% of the 40’s age group population are confined). The flow chart between compartments is shown in Fig. 1.
Schematic flow chart between compartments.
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 on how we validated the model and how to determine the 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, as well as 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 ζi) are shown in Table 2.
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 of 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 infection, we start the integration τ days before the first reported COVID-19 death, i.e., t = 0 means tr = τ. The initial conditions are fully specified when the initial number of exposed individuals is defined. 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.
Fatality rate as a function of median age. The fatality rate by age bins is taken from [18] and the population pyramids from UN data.
Left: The time series of fatalities for a number of countries. Right: the time derivative of the curve of fatalities is converted into ℛ(t) according to Eq. (19).
The timescales σ, γ, θ, and ξ, as well as the fractions p and w, are found by Markov chain Monte Carlo (MCMC) fitting, with the priors 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 up to mid-June. Parameters determined by the MCMC modeling are shown in Fig. 5, being 
, and 
. ℛ(t) at present is hovering around 2.
a) Evolution of the compartments of exposed (E), asymptomatic (A), symptomatic (I) and hospitalized (H). b) The same curve of H, and the fraction of hospitalizations needing ICU (U). The green dashed curve shows the total number of ICU beds in the country. At the current rate, the epidemics should peak in mid-May and collapse of the health care system should happen around May 1st. c) Cumulative number of reported cases and fatalities as orange and black dots, respectively. The number of hospitalizations closely matches the number of confirmed cases.
Posterior probabilities for the epidemic parameters determined by the MCMC modeling.
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 the number of symptomatic is predicted to end at July 1st. The peak of symptomatics is predicted for May 17th with 20 million symptomatics. Consequently, there is a predicted rise in the number of hospitalized, reaching saturation on May 3rd and peaking on 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 on June 1st is 106.
Fig. 4 contrasts the predicted cumulative numbers of infected persons (orange line), hospitalized persons (blue line), and deaths (black line). The figure also shows the cumulative number of confirmed cases (yellow dots) and actual 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. 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-by-state variation from a maximum of 64.7% to a minimum of 53.7%). At this low level of distancing, capacity should be reached in less than 50 days, which is in agreement with the dynamical ℛ(t) model in Fig. 4.
a) Number of hospitalizations. b) ICU cases. c) Fatalities. The three curves are shown as function of the degree of horizontal confinement. To not overwhelm the health care system capacity (≈3 × 104) ICU beds, the level of social distancing should be over 70%. Brazil is managing 56%.
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 rows of plots show the same analysis but confining 60% of the population up 50, 60, and 70 years old, respectively. As seen in the cyan line, the number of hospitalized rises from 30 to 60 years old and falls for 70 years old onwards. That is because even though 70+ are more likely to be hospitalized, the number of 30-60 is much higher in the population.
Upper: Number of hospitalizations in horizontal confinement with different proportions of the population, broken down by age. 2nd row: 60% of the population under 40 is confined, the population older than 40 is confined to a higher degree, at 90% (solid blue line) and 99% (dashed blue line). The 3rd, 4th, and 5th rows of plots show the same analysis but confining 60% of the population up to 50, 60, and 70 years old, respectively.
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 the 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, hundreds of thousands of deaths are seen in the 60-70, 70-80, and 80+ age bins. The number drops 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 a much lower death rate for this age segment.
Middle: Same as Fig. 7, but for fraction of hospitalized that need ICU.
Same as Fig. 7, but for the number of fatalities.
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 the general 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 the 40-50 age group.
Vertical confinement as exit strategy after a horizontal lockdown of 60% was held, from March 22nd to May 1st. Upper plots: S and C compartments, normalized by the number of individuals in the respective age bin. 2nd row: number of hospitalizations, 3rd row: number of ICU cases, lower plots: cumulative number of fatalities. The figure shows full release for the population under 50 (dashed line) and under 60 (solid line). The population past this age is kept at 90% confinement.
Limitations
As in any setting, the outbreak response strategy plays a crucial role on the quality of the outputs the models can give. Since the identification of the first case, the response strategy in Brazil 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 on April 1st. This is more than the number of expected cases in the country while we write this article, considering the estimated sub-notification of cases [20] and inaction on controlling the infection. It is possible that the actual number is lower, although it is also important to notice that Brazil has not done a real lockdown so far.
The model also ignores mobility, in the sense that it does not consider travel to and from the country. Given that right now we are at the stage of community transmission, 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 are 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
In this study we examine the strategy of vertical confinement as currently debated in Brazil. Since the fatality rate of COVID-19 appears to be higher among the elderly population, we studied how confinement by age groups (particularly 60 years old and beyond) affects the demand for hospital beds and intensive care beds.
Our model suggests that at the current rate of advance of the pandemics, Brazil should face collapse of the health care system by May 15th, 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 overload of the health care system. An increase in lockdown to 70% is needed to avoid the number of cases from overcoming the number of available critical care beds.
An exit strategy of confinement of individuals older than 60 years 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. Confinement by age group should consider the population over 50 years old as well. However, the age range 50-60 is also a part of the workforce, and thus defeats the purpose of a confinement by age. 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 those found in the low and even middle income neighborhoods common in Brazil. Our results therefore discourage confinement by age as the only exit strategy. 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 [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 (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.
Input parameters and arrays for MCMC
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
References
