key: cord-0772571-8rdqv5vi
authors: Lyra, Wladimir; do Nascimento, José-Dias; Belkhiria, Jaber; de Almeida, Leandro; Chrispim, Pedro Paulo M.; de Andrade, Ion
title: COVID-19 pandemics modeling with modified determinist SEIR, social distancing, and age stratification. The effect of vertical confinement and release in Brazil
date: 2020-09-02
journal: PLoS One
DOI: 10.1371/journal.pone.0237627
sha: 04351603479d67407c7f340df5dec02a405a645e
doc_id: 772571
cord_uid: 8rdqv5vi

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 (e.g. 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. Sensitivity analysis shows the 95% confidence interval brackets a order of magnitude in cases or three weeks in time. 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.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

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 SEIR 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 SEIR 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 (April 9th, 2020). 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 measurement of the percentage of the population that is currently in confinement; however, the number is estimated to be around 56% according to satellite data.

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, and run a sensitivity analysis. Finally, we test the effect of both general and vertical confinement on the epidemic curve.

We used a modified version of a SEIR-type deterministic compartmental model to trace COVID-19 epidemic evolution in an isolated population of N individuals. 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, 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 � ∑ i X i and compartments are such that S + C + E + A + I + Q + H + R = N, with N � ∑ i N i being the total population; N i is the population in each age bin. The software is written in python 3.7, and is made public at https:// github.com/wlyra/covid19. 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 q i 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 latent period 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 RðtÞ via

Lock-down consists of having a fraction of the susceptible population removed from the epidemic dynamic by moving them from S i to C i at a rate ψ i . Similarly, lifting the lock-down is done by placing C i into S i at the rate ϕ i . We consider these functions to be Dirac deltas

where t lock and t lift are the time (in days) of lock-down and of lifting of the lock-down, respectively. To allow for partial demographic lock-downs, a i and b i 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. Other diagnostic quantities are the numbers U i of people in need of an intensive care unit (ICU) bed

where z i is the fraction of hospitalized patients that need critical care. Both z i and the hospitalization fraction q i are age-stratified.

For integration, we use a standard Runge-Kutta algorithm, with timesteps

In this section we present details on how we validated the model and how to determine the characteristic timescales and other parameters.

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 q i , and fraction of hospitalized that need critical care z 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 (https://www.populationpyramid.net), and split into the pre-defined age bins.

We derive Rð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 t r � t + τ for the retarded time. Summing over all age bins D � ∑ i D i we have the cumulative death rate on the LHS, which is an observable

and hμSi�∑ i μ i S i . We can then substitute Eq (9) and solve for Rð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 t r = τ. The initial conditions are fully specified when the initial number of exposed individuals is defined. This should be

where t 0 is the time of the first death and � m � N À 1 P i m i n i 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 (https://systems.jhu.edu/research/public-health/ncov). We plot in the left panel of Fig 2 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 RðtÞ, feeding this value into Eqs (1)-(8) to start the SEIR evolution. The populations I(t) and S(t) that enter in Eq (19) are then calculated to update RðtÞ. The resulting values are plotted in the right-hand-side of Fig 2. 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. , 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 10 6 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 PLOS ONE 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 10 6 . Fig 3c 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.

We perform a sensitivity analysis, shown in Fig 5, by varying the parameters of the models by -2, 0, and 2 standard deviations as given by the results of the MCMC analysis (Fig 4) . Given 7 parameters, we run 3 7 = 2187 simulations. The fiducial model, with zero standard in all parameters, is shown as the thick line; all other models are shown as thin lines. The 95% confidence interval brackets about an order of magnitude above or below the fiducial model, or about three weeks left or right of it.

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 × 10 4 ) 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 RðtÞ model in Fig 3. 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. 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.

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.

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 assumes, in Eq (10) that there is no difference in infectiveness between symptomatic and asymptomatic population. This is an assumption that should be updated as further knowledge of COVID-19 is unveiled. We stress also that by using the cummulative hospitalization data as a guide in the MCMC, we are incurring into the problem raised by [21] , of overestimating the confidence interval precision; fitting raw incidence data instead would enhance the statistical accuracy.

The model also ignores mobility, in the sense that it does not consider travel to and from the country. Given that Brazil is at the stage of community transmission (as of April 9th, 2020), 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.

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 https://doi.org/10.1371/journal.pone.0237627.g010 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 10 6 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. The 95% confidence interval spans two orders of magnitude in cases or a month and a half in time.

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.

Eqs 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 (10 20 ) 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 .

Wladimir Lyra, José-Dias do Nascimento, Jr., Jaber Belkhiria, Leandro de Almeida, Pedro Paulo M. Chrispim, Ion de Andrade.

COVID-19) Situation Report-75. World Health Organization

Characteristics of and Important Lessons From the Coronavirus Disease 2019 (COVID-19) Outbreak in China: Summary of a Report of 72314 Cases From the Chinese Center for Disease Control and Prevention

COVID-19: towards controlling of a pandemic

The Ebola epidemic in West Africa: Challenges, opportunities, and policy priority areas

The Next Epidemic-Lessons from Ebola

Major issues and challenges of influenza pandemic preparedness in developing countries. Emerging infectious diseases

COVID-19 and Italy: what next? The Lancet

COVID-19) Situation Report-79. World Health Organization

COVID-19) Situation Report-72. World Health Organization

Mathematical modeling of the West Africa Ebola epidemic

Modeling the spatiotemporal transmission of Ebola disease and optimal control: a regional approach

A Modified Sir Model to Study on Physical Behaviour among Smallpox Infective Population in Bangladesh

Modeling control strategies for influenza A H1N1 epidemics: SIR models

Simulating SARS: Small-World Epidemiological Modeling and Public Health Policy Assessments

A Contribution to the Mathematical Theory of Epidemics

A simple vaccination model with multiple endemic states

Nacional do Estabelecimentos de Saú de

Impact of non-pharmaceutical interventions (NPIs) to reduce COVID19 mortality and healthcare demand

Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The Lancet Infectious Diseases

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2)

Avoidable errors in the modeling of outbreaks of emerging pathogens, with special reference to Ebola. Proceedings Biological sciences / The Royal Society

Bayesian Data Analysis

We thank the referee, Qianying Lin, for his constructive comments that helped improve the manuscript. 

To fit the best value to w and p, and to better constrain σ −1 , γ −1 , θ −1 , ξ −1 , we use the affineinvariant ensemble sampler for Markov chain Monte Carlo (MCMC) [22] 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 H c , we generated a cumulative error C err on the reported confirmed cases C c .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 H c with C c . Rather, we weigh the C c array with the H c array using:H 0 c is the weighed cumulative hospitalizations and n is the length of the data. Following we get the residual between C c and H 0 c , and we used the negative binomial distribution to calculate each likelihood [23] :

Conceptualization: Wladimir Lyra, Jaber Belkhiria, Pedro Paulo M. Chrispim.

Formal analysis: Wladimir Lyra, Jaber Belkhiria, Leandro de Almeida.Investigation: Wladimir Lyra, José-Dias do Nascimento, Jr., Leandro de Almeida, Pedro Paulo M. Chrispim.