Spreading of COVID-19 in Brazil: Impacts and uncertainties in social distancing strategies

Appendix A.1. Model Description

Mathematical models can help the investigation of different scenarios, devising alternative strategies to control virus spreading. The SARS coronavirus epidemic in 2002-2003 generated several publications using diversified modeling techniques (see [5] for a review). Some of these publications employed mathematical models to provide some understanding of the propagation pattern and ultimately to infer control measures for the disease [22, 23, 24]. Similar mathematical models were employed for other diseases, for example: (I) the 1918 influenza pandemic using data from cities of the USA [25, 26], England and Wales [27]; (ii) the 2009 influenza pandemic was modeled in local regions in Japan [28]; and (iii) the 2014-2016 Ebola epidemic was modeled for the West African countries with reported cases of the disease [29, 30]. Overall, the modeling efforts for representing virus spreading share common challenges that have been pointed out in the literature recently [31, 32], including those related to the emergence of “novel infectious agents” [33]. Issues like virus genetic diversity, within-host variability, social, demographic, and spatial heterogeneity, among others, can have a significant impact on model predictions and control programs. In any case, model granularity and ability to capture the desired phenomena are linked to available data for model calibration, which may hinder the development and practical use of more complex models.

Several mathematical models have already been published to better understand the disease spread and the effectiveness of different mitigation programs. Two variations of the Susceptible-Infectious-Recovered (SIR) model are used in [7] to describe the early evolution of COVID-19 in Brazil (BR), using data from February 25, 2020 to March 30, 2020. Short and long term forecasts are performed, which highlight the need of keeping social distancing policies for longer periods to actually flatten the epidemic curve, and the necessity of testing to have a better picture of the amplitude of the peak of symptomatic infected individuals. In [34], an eight compartment model with nine age stratification groups was developed to investigate the effect of social distancing policies, specifically those related to specific (older) age groups. Starting from a 60% lockdown scenario and many predefined parameterizations based on available references, their analysis indicates the ineffectiveness of older age group confinement policy as a way of controlling the confinement exit. A similar model approach is developed in [35] to investigate the transmission and isolation measures in São Paulo state. They found that both young and older people must be isolated to reduce hospitalization. The two latter studies explicitly included quarantined individuals, either due to confinement of susceptible individuals or because of a positive testing isolation measure. More recently, the work [36] develops a model that accounts for either reported and unreported cases. Besides analyzing the impact of different public health mitigation measures, the work suggests the use of a mixture approach for estimating model parameters. At the early onset of the disease, confirmed reported cases are assumed to follow an exponential law to ease getting initial conditions (ICs). Later on, a Bayesian approach is used for model calibration. Since the unreported infected cases are major issue to control outbreaks, the approach for identifying them seems promising and has already been successfully applied to the Hong Kong seasonal influenza epidemic in New York City in 1968–1969 for a SIR model [37].

Another class of compartmental models widely used to describe the spread of virus diseases worldwide is SEIR, that consider the total population (N) is split into sub-groups of Susceptible (S), Exposed (E), Infected (I), and Recovered or Removed (R) individuals, such that S+E+I+R = N. The compartments S and I contain all individuals likely to be infected and those already infected, respectively. E individuals are usually considered to have low viral load and are not likely to be infectious, while R individuals can or cannot be infected again depending on the virus disease. This class of dynamical model may be used to ultimately identify the processes that interplay in pandemic waves and has been largely used recently to assess the dynamics of COVID-19 in many countries [38, 39, 40, 41, 42, 43], including BR [44, 7, 34, 35, 36, 45, 46, 13].

Inspired in the generalized SEIR model proposed in [38], Jia et al. [47] employed the SEAIRPQ model with seven compartments for modeling the Chinese outbreak of COVID-19, including populations of positively diagnosed (P) individuals who are under medical treatment and of quarantined (Q) individuals. A schematic description of the SEAIRPQ is shown in Fig. A.6, with its mathematical formulation given by the following system of ordinary differential equations:

Model parameters are non-negative constants and most of them are described in Table A.2. The three remaining are the fraction of infection rates between A and I, θ, the time required for an individual to be quarantined, 1/p, and the quarantine period, 1/λ. An interesting feature of the SEAIRPQ model is the explicit inclusion of a compartment for Q individuals. Longer periods of isolation yields a decrease in S population, ultimately decreasing the infected population. S individuals may be exposed to contact with A and I individuals. For the Chinese considered scenarios, it was assumed that A individuals are less likely to spread the disease, so that θ << 1. The SEAIRPQ model assumes that both A and I individuals are eventually positively diagnosed at rates ε_a and ε_I, respectively. They can also recover without medical treatment at rates γ_a and γ_I, respectively.

Although SEAIRPQ model was able to describe quite well the outbreak in China, the Brazilian reality may prevent its use in Brazilian localities. BR testing capacity is still very limited, one of the lowest in the world. In April 2020, BR tested around 296 per million people. According to [3], this number increased to 4,104 in May 2020, although it can be significantly lower in some Brazilian localities due to huge demographic and social heterogeneities. Thus, available data on testing is very scarce: the present Brazilian policy is to test only people with severe symptoms. Such evidence has driven the development of the model proposed in this paper.

Here we consider a modification of model (A.1) taking into account protective quarantine measures as an implicit variable. The proposed model, called SEIARPD-Q, has also seven compartments since it includes the compartment of dead (D) individuals to explicitly keep track of them. The notation used indicates that Q individuals are not separated in a specific compartment. Instead, quarantine mechanisms are implicitly considered, with individuals in isolation included in R compartment. Likewise SEIARPQ model, the infected class is modeled as two separate compartments encompassing I and A individuals. R compartment includes individuals who have recovered from the disease as well as those who have been removed to this compartment as a result of isolation measures.

For completeness, we repeat here the mathematical formulation of the SEAIRPD-Q model, which is given by the following set of differential equations:

The model is schematically shown in Fig. A.7. All model parameters meanings are described in Table A.2. Data sources are informed in Table A.3 and Table 1 (Materials and Methods). Information regarding the numerical differential equation solver employed can be seen in SI-SAppendix A.2.

Looking at the mathematical description of model (A.2), it is important to note that the transformation of S into E individuals occurs through interaction with individuals from the infected groups (I and A) according to β and μ rates, respectively. The remaining mechanisms are linear, depending only on the corresponding transfer rates. As already mentioned, an important mechanism introduced in the model describes the application of social distancing measures: individuals from compartments S, E, A, and I are removed to compartment R at a rate ω. In this way, it is possible to exclude R individuals from the infection dynamics. A similar approach was adopted in [11]. We emphasize that we have chosen an implicit quarantine/social distancing approach, that is, we do not model a specific compartment of Q individuals, unlike other studies [47, 38]. The advantage of this implicit representation is to obtain a simpler model, with less computational cost and fewer parameters.

The previous assumptions reflect our current knowledge of the COVID-19 in BR. We remark that the adaptive immunity to SARS-CoV-2 of people who recovered from the disease has not been established yet as well as how long the immunity could eventually last. Actually, a follow-up study on 262 recovered patients in China revealed that 15% were tested positively for SARS-Cov-2. They were all younger than 14-years old and had mild and moderate conditions at the time of the primary infection. During the follow-up, they did not show clinical symptoms and it is speculated if either more sensitive test is required or indicates a sign of relapse rather than reinfection. Another study performed in previously infected rhesus macaques found that reinfection of SARS-Cov-2 could not occur [48]. Overall, the current understanding, based on what is known on coronaviruses in general, is that antibodies produced by infected patients may give them immunity to the specific virus for short term or even years. For example, previously infected SARS patients can become susceptible to reinfection only after around three years [49].

Appendix A.2. Ordinary Differential Equations Solving Method

The mathematical formulation of the SEAIRPD-Q model described by the system of seven ordinary differential equations (A.2) is solved using the LSODA method [50] from SciPy [51]. This package provides an interface for such routine from ODEPACK [52].

Appendix A.3. Evaluation of the Reproduction Number

Following the ideas introduced in [53], we apply the Next Generation Matrix (NGM) method to obtain for the SEAIRPD-Q model. Firstly, we use a reduction process to arrive at a linearized infection model, under the disease-free steady state. The identification of the infected variables (E, A, and I) allows defining the vector x^T = [E, I, A], where the superscript T stands for the transposition of a matrix. Thus, we consider that S, R, P, and D individuals are not able to transmit the disease. Denoting by T and Σ the transmission and transition matrices, respectively, the three-dimensional linearized infection sub-system is given by: where

The reproduction number is the dominant eigenvalue of the next generation matrix K = −TΣ⁻¹ [53], which yields:

The corresponding effective reproduction number is:

Appendix A.4. Bayesian Calibration

In this work we adjust our model parameters to make results consistent with available observations (i.e., C and D for both BR and RJ scenarios). With Bayesian calibration, we are capable of determining the most likely uncertainties for input parameters as well as their MAP values described by the a posteriori probability distribution π_post. Considering that our aim is to adjust a set of parameters θ given y observations, Bayes’ theorem states: where π_prior is the a priori probability distribution that represents the initial (a priori) knowledge over θ, π_like is the likelihood function, and π_evid is the evidence (information) encompassed in the data (see [54] for a more detailed description). We assume that all data follows a Gaussian distribution centered on simulated of observable quantities values with variances and , which yields the following form for the likelihood function:

The standard deviations for the cumulative confirmed (σ_C) and dead cases (σ_D) are tuned as hyperparameters within the calibration process. Regarding π_prior, a list of prior distributions and MAP estimates of marginal distributions can be seen in Table 1 (Materials and Methods). The frequency histograms for the calibrated parameters of both BR and RJ scenarios are shown in Figs. A.8 and A.12, respectively. The posterior distributions are obtained by a tempered adaptive Bayesian updating procedure proposed in [16], the Cascading Adaptive Metropolis in Parallel, used here through PyMC3 [15] implementation in Python language (version 3.7).

Despite the fact that we have mentioned which observed data are used as observations for the calibration procedure, a proper definition of such quantities is absent. The estimate of C is calculated from the model as: which is the cumulative value of inflow individuals in the P compartment, since only diagnosed individuals are recorded on Brazilian confirmed case databases due to the lack of large-scale testing in the overall population. Using Eq. A.10, we can estimate the amount C at t days after day t₀. Analogously, in order to compute the D from SEAIRPD-Q, we calculate: which is the integration of Ḋ (see system (A.2) for reference) over the time range of interest. C and D time-stamps are gathered from data sources specified in the Materials and Methods section. Such quantities are used as observed data in the calibration procedure.

Appendix A.5. Sensitivity Analysis

Due to paucity data and possible model identifiability issues, not every parameter and unknown initial conditions (ICs) are estimated using the Bayesian calibration procedure. Based on previous studies, we assume values for every parameter (with the exception of β, μ, ω, d_I, and d_P), and IC (with the exception of S(0), P(0), D(0), and R(0)). Relevant information on model settings is provided in Tables 1 (Materials and Methods), A.3 and A.2. In order to quantify how these choices affect the outcome of the dynamics, we perform a sensitivity analysis. Specifically, we are interested in the sensibility of the quantities:

As previously stated in the main text, we employ the Elementary Effects method [20] to carry out this analysis. In contrast to local sensitivity methods, which evaluate how small changes around a singular point affect a QoI [55], the Elementary Effects method is a global sensitivity method. Considering the θ ∈ ℝ^k as the set of analysed parameters, we construct a sample of m initial points θ⁽¹⁾, ···, θ^(m) from a k-dimensional p-level grid. Afterwards, a metric defined as an elementary effect EEi, i = 1,…,k, is computed for every parameter: where is the magnitude of change in the ith parameter in the k-dimensional p-level grid. Values of p as an even number and are usual suggested choices to ensure an equal sampling probability in the parametric space, and are the default option in the SALib [21] library used. After computing the elementary effects, the process ends by calculating the global sensitivity coefficients for each parameter. Our work dedicates the analysis only to the first order sensitivity coefficient:

Parameters with high scores are the most influential to the QoI. In turn, the first sensitivity index values allow ranking the parameters with respect to their order of importance. Sensitivity analysis is performed for both scenarios (BR and RJ) considering all parameters and ICs E(0), A(0), and I(0). The settings p = 4 and m = 40 are used for all experiments. A 50% variation around the fixed and MAP values in Tables 1 (Materials and Methods) and A.3 is considered to construct parameter ranges. The QoIs are observed in days 2, 25, 50, 75, 100, 125, 150, 175, 200, 225, and 250 in order to evaluate the scores s_i along all simulation period. For the sake of presentation, the scores are normalized (sum to 1).

The sensitivity analysis results for BR are depicted in Fig. A.10. The sensitivity indexes for the cumulative C and D (QoI₂(t)) are heavily influenced by the initial conditions in the early stages of the disease propagation, although their influence decreases as time evolves. In contrast, they are not influential for along all the simulation time. The influence of the parameter related to social distancing (ω) on both QoIs is notable. More importantly, its influence increases with the evolution of time, which emphasizes the need for care when establishing social distancing relaxation measures. On the other hand, the proportion of I individuals is also one of the most influential parameters on both QoIs, although it is more relevant considering C and D. This points out the need of having a more widespread testing policy. Regarding the RJ scenario, Fig. A.14 reveals very similar scores. One noteworthy distinction is the rank of influence of E(0) and I(0) associated with QoI₂(t). The former is more influential in the BR scenario (Fig. A.10b), while the latter has a higher score for RJ (Fig. A.14b).

Appendix A.6. BR: Additional Results

Fig. A.9a shows the prediction of P, I, A, D, and C in BR. A total of 75.7 thousand (95% CI: 67.7–85.8) active cases are expected on the simulation day 81 (95% CI: 78–83), with D and C expected to be around 34.7 thousand (95% CI: 32.4–37.3) and 364.2 thousand (95% CI: 325.7–414.8), respectively. The posterior distribution of the peak position is displayed in Fig. A.9b, with the vertical dashed lines corresponding to those displayed in Fig. A.9a. Fig. A.9c depicts the temporal changes of and displays two vertical lines identifying the credible interval (95% CI: 69–73) of the time above which . The same lines are depicted in Fig. A.9d that shows the uncertainty associated with that time value.

Fig. A.11 provides the model forecasts for the considered hypothetical scenarios of social distancing relaxation. Fig. A.11a shows what to expect in case of a sudden release from quarantine after the simulation day 100 (t_1/2 = 0.1 days as half-life decay). In this case, C reaches 434.8 thousand (95% CI: 372.1–526.8) and D 41.4 thousand (95% CI: 37.8–46.0) at the end of the simulation (day 262), which corresponds to an increase of approximately 19% in both values when compared to the original scenario with ω fixed at its MAP estimate. For a gradual release after the simulation day 100 (t_1/2 = 15 days as half-life decay), the numbers of C and D at the end of simulation are 378.0 thousand (95% CI: 334.7–436.0) and 36.0 thousand (95% CI: 33.6–38.9), respectively, as shown in Fig. A.11b. Figs. A.11c and A.11d compare the variability of C and D at the end of simulation in the form of box plots for all scenarios considered. Outliers appear as individual points and the samples medians are depicted in red.

Appendix A.7. RJ: Additional Results

Fig. A.13a shows the prediction of P, I, A, D, and C in RJ. A total of 7.3 thousand (95% CI: 6.38.7) active cases are expected on the simulation day 80 (95% CI: 76–85), with D and C expected to be around 7.0 thousand (95% CI: 5.8–8.7) and 45.0 thousand (95% CI: 37.7–55.5), respectively. The posterior distribution of the peak position is displayed in Fig. A.13b, with the vertical dashed lines corresponding to those displayed in Fig. A.13a. Fig. A.13c depicts the temporal changes of and displays two vertical lines identifying the credible interval (95% CI: 67–76) of the time above which . The same lines are depicted in Fig. A.13d that shows the uncertainty associated with that time value.

Fig. A.15 provides the model forecasts for the considered hypothetical scenarios of social distancing relaxation. Fig. A.15a shows what to expect in case of a sudden release from quarantine after the simulation day 100 (t_1/2 = 0.1 days as half-life decay). In this case, C reaches 57.1 thousand (95% CI: 44.4–80.1) and D 8.8 thousand (95% CI: 6.8–12.3) at the end of the simulation (day 252), which corresponds to an increase of approximately 27% in both values when compared to the original scenario with u fixed at its MAP estimate. For a gradual release after the simulation day 100 (t_1/2 = 15 days as half-life decay), the numbers of C and D at the end of simulation are 47.8 thousand (95% CI: 39.3–61.1) and 7.4 thousand (95% CI: 6.0–9.5), respectively, as shown in Fig. A.15b. Figs. A.15c and A.15d compare the variability of C and D at the end of simulation in the form of box plots for all scenarios considered. Outliers appear as individual points and the samples medians are depicted in red.

• Download figure
• Open in new tab

Figure A.6:

A schematic description of the SEAIRPQ model, modified from [47]

• Download figure
• Open in new tab

Figure A.7:

A schematic description of the SEIARPD-Q model.

• Download figure
• Open in new tab

Figure A.8:

Frequency histograms for the calibrated parameters (BR).

• Download figure
• Open in new tab

Figure A.9:

Dynamics of the COVID-19 in BR modeled with the available data (Materials and Methods).

• Download figure
• Open in new tab

Figure A.10:

Temporal changes of the first order sensitivity index of model factors for the BR scenario (Materials and Methods).

• Download figure
• Open in new tab

Figure A.11:

Model forecasts in BR considering sudden and gradual releases from quarantine (t₀ = 100 days and t_1/2 = 0.1 and 15 days). The original (baseline) scenario of the main text is also presented to ease comparison.

• Download figure
• Open in new tab

Figure A.12:

Frequency histograms for the calibrated parameters (RJ).

• Download figure
• Open in new tab

Figure A.13:

Dynamics of the COVID-19 in RJ modeled with the available data (Materials and Methods).

• Download figure
• Open in new tab

Figure A.14:

Temporal changes of the first order sensitivity index of model factors for the RJ scenario (Materials and Methods).

• Download figure
• Open in new tab

Figure A.15:

Model forecasts in RJ considering sudden and gradual releases from quarantine (t₀ = 100 days and t_1/2 = 0.1 and 15 days). The original (baseline) scenario of the main text is also presented to ease comparison.