The unexpected dynamics of COVID-19 in Manaus, Brazil: Herd immunity versus interventions

In December 2020, a vicious second wave of COVID-19 erupted in Manaus, Brazil, resulting in >100 residents dying per day, and triggering yet another collapse in the country’s rundown healthcare system¹. Awkwardly, this contradicts the recent Science paper of Buss et al. (2020)² which suggested that any second wave would be highly unlikely. The authors estimated that the first wave infected ∼76% of the city’s population by October 2020 (see also³). With an Attack Rate (AR)∼76%, it appeared that herd immunity had been reached, and no further epidemic should be expected, assuming reinfection is relatively rare. Currently there are only speculations as to why a major second wave then appeared (Sabino et al, 2021)⁴. For example, measurements of AR via seroprevalence testing could possibly be inaccurate given the biases of the blood donor data⁵ used by Buss et al. (2020). It is therefore useful to explore alternative and complementary epidemiological modelling analyses of publicly available mortality datasets to help understand why the first wave of the COVID-19 outbreak suddenly “crashed,” and why a secondary wave later unexpectedly appeared in Manaus in December 2020. Here we discuss and expand on our results presented in Ref 6. Elucidating the underlying processes that took place is of central importance for planning future interventions appropriately, not only for the benefit of Manaus which is currently experiencing the full fury of SARS-CoV-2, but also for numerous other potentially high-risk locations worldwide.

We first mention that it is not uncommon for potentially large-scale epidemic outbreaks to die out well before predictions of conventional epidemiological models, and well before herd-immunity can develop⁷. For example, the Ebola outbreak in Sierra Leone Guinea 2014⁸ began with an initially susceptible population of 5,364,000 residents. While a conventional SEIR-type model predicted 5,157,000 Ebola cases⁶, the epidemic resulted in “only” 14,124 infected individuals. We documented a similarly perplexing account of a catastrophic epidemic outbreak⁹ in WW2. Brauer (2019) ⁶ showed that the early demise of an epidemic can occur in the presence of: i) heterogeneous transmission and super-spreader-like events^7,10; ii) human behavioural responses which reduce the size of epidemics by reducing the reproductive number R₀(t) over time^7,9,11,12, creating a “herd-protection” effect based on behaviour rather than immunity. Brauer (2019)⁷ modelled behavioural interference by assuming that R₀(t) (which is a measure of transmission) reduces exponentially over time at rate e^-κt (See also¹²). Other studies of cholera and influenza use similar parametric approaches¹². Here we generalise this principle by fitting a “flexible” R₀(t) to the data without pre-imposing any specific functional form ^9,11.

Figure 1 plots the COVID-19 dynamics in Brazil for eight selected states (data from¹⁴). Each panel displays the number of COVID-confirmed deaths per day for the given state over the epidemic. In theory, mortality data should reflect trends in the disease dynamics more reliably than other surveillance data (e.g., confirmed cases), although under-reporting and disease confirmation issues remain sources of bias. We repeated the plots for all Brazilian states (see Appendix), and it became immediately obvious that the epidemic curve of the first wave was roughly synchronized over all of Brazil in the period May-July. All 27 regions peaked in this 2-month period, some earlier than others, and then declined within several months. The differences in timing between states presumably relate to spatial spread and differences in local NPI policies in terms of initiation and intensity. From this perspective, it would be unusual to argue that the epidemic crashed naturally in Amazonas after reaching herd immunity, while in other locations the outbreak turned around due to the involvement of non-pharmaceutical interventions (NPIs; social-distancing, facemasks, hand-sanitizing, work stoppages, lockdowns, travel bans, restrictions on gatherings, etc.).

• Download figure
• Open in new tab

Figure 1.

Fitting an SEIRD model to reported daily COVID-19 confirmed mortality data in eight Brazilian states¹³ with flexible transmission. We set n_m=7. We used the plug-and-play inference framework which is built on the iterated filtering method. More details on the methodology can be found in references^11,15.

Exploring this more deeply, we fitted an SEIR-type epidemiological model (see Appendix) to COVID-19 mortality datasets and estimated R₀(t) as it changed, as shown in Figures 1&2 (blue). The model is described in the Appendix and in Refs.7,9&13. In all states, and in the cities of Manaus and Sao Paolo, we found evidence that R₀(t) fell substantially with time over the first wave as the epidemic proceeded, and it thus appeared to be a regional phenomenon. In line with Brauer (2019)⁷, the reduction in R₀(t) is likely to be due to behavioural reactions such as state-wide NPI’s, individual behavioural changes, and possibly also combined with seasonal and/or climatic factors. In addition, it was impossible to fit SEIR-type models to the mortality data with a constant R₀, as might be a reasonable assumption if interventions or seasonality had little impact, as has been suggested. Any such fit yielded outcomes that did not appear meaningful. For example, for Manaus, the best fit in Figure 2a predicted a large R₀(t) = 3.1 over the whole period through to January 2021, and provided a very poor fit to the data. However, with a flexible fitting of R₀(t), as in Figure 2b&c, an excellent fit to the data is obtained. Figure 2b fits COVID-19 confirmed mortality data obtained from Ref.16 while Fig.2c fits Severe Acute Respiratory Illness mortality data from Ref.21, both collections for Manaus. Moreover, the SEIR-type model, despite being famous for its ability to fit epidemics of all types, was unable to do so unless R₀(t) decreased rapidly in April, thereby allowing a low stationary endemic state to be maintained data from July to November, as had been the case. The model estimated an AR∼30% until October 2020 in Manaus, which would appear to be reasonable given the large ongoing second wave that followed. Mellan et al. (2020)²² predicted a smaller AR. In summary, the first wave of the Manaus epidemic appeared to crash due to a reduction in R₀(t), rather than because of a large attack rate and the attainment of herd immunity.

• Download figure
• Open in new tab

Figure 2.

Fitting an SEIRD model to reported daily Covid-19 confirmed deaths in Manaus¹⁶ with a flexible transmission rate R₀(t). We use the plug-and-play likelihood based inference framework¹⁵. (a) A constant transmission rate (constant R₀) fails to fit the major wave well (panel a; number of nodes in transmission rate spline n_m=2). (b) When the spline’s number of nodes increases to 7, this gives a flexible R₀(t), and the fitting performance is improved. Based on the 2^nd-order Akaike Information Criterion (AICc), n_m=7 yields the best performance. (c) Same as panel b but based on SARI Manaus mortality data from Ref.21. We assume the mean latent period at 2 days, mean infectious period at 3.5 days, thus the mean generation time 5.5 days (in line with previous studies on serial interval and generation time), mean time from loss of infectiousness to death at 12 days (thus from infection to death is about 17.5 days which is in line with observation²¹.) Estimated R₀ for Manaus 2.250 95% CI(1.787, 2.791). See Appendix for discussion of age-structure.

Buss et al. (2020)² showed that their estimate of a high AR∼76% implies that the infection fatality rate (IFR) must be extremely low for consistency (on average IFR∼0.018-0.26%) – but they did not mechanistically model the epidemic trajectory over time as done here. We were unable to fit the mortality dynamics sensibly with such a low IFR (see Figure S1). Our results indicate that intervention activities implemented in Manaus played an important role in the crash of the first wave seen in early May 2020, as also argued by Mellan et al. (2020)²². Although no national lockdown had officially been set early in the epidemic, states and cities adopted their own measures. Buss et al (2020)² document that physical distancing increased during March. By March 23, a state of emergency was declared in Amazonas and in fact all of Brazil. Soon after, Manaus introduced a range of NPI’s. There have been valid criticisms concerning how well these NPIs were enforced, their late implementation, and the levels of compliance; some sectors of society generally associated with the far-right, deliberately flouted the most basic rulings. In addition, those sectors most severely suffering from inherent structural inequalities had also showed a lower adherence to lockdowns initially, as many lacked financial security to survive for over a week without income. Yet quarantine was generally favoured in these sectors, and with increasing financial and social aid, adherence rose. Generally, in Manaus, the NPIs were implemented often with the strong support of the Mayor, who was often undermined by the Brazilian President Bolsarano. By the beginning of April, health systems were reaching dangerously close to their threshold of collapse. Nevertheless, by mid-May the hospitals were back on their feet, and the epidemic was declining rapidly (Figure 2), although only after many lives were lost. Similar NPI’s were implemented in nearly all other Brazilian states, roughly synchronously to combat the newly emerging countrywide epidemic, although they have been subject to criticism for poor organization, compliance, and disregard. Nonetheless, these results strongly suggest that the implementation of NPIs was sufficiently effective to provide “herd-protection”, despite the notable political and infrastructural challenges faced by the residence of Manaus.

Our work highlights the difficulties policy makers face when attempting to predict herd immunity, as well as the need for developing and using complementary mechanistic epidemiological modelling tools and datasets.

Supplementary Figures

• Download figure
• Open in new tab

Figure S1.

Reproduces Figure 2 with IFR <0.3%. The estimated R0(t) goes to unreasonably high values during the second wave.

• Download figure
• Open in new tab

Figure S2.

Fitting SEIRD model to daily SARI death in 26 Brazilian states and the federal district¹.

View this table:

• View inline
• View popup

Table S1.

27 major Brazilian cities: the infection attack rate (IAR) median is 28%, with a range (19% to 58.8%). Manaus 29.2% by 16 January 2021, as found from SEIR model fits. To SARI mortality data from Ref.1.

Appendix

Models and Results

We fit a susceptible-exposed-infectious-recovered type model with a flexible time-varying transmission rate β (t) to the reported Severe Acute respiratory syndrome (SARI) mortality data ¹: where the compartments S, E, I, and R denote the conventional susceptible, exposed, infectious, and recovered individuals, respectively. The variable N is the total population size of the city or state, H denotes the hospitalized cases (or severe cases), and D denotes the total of SARI-deaths. Parameters β (t), σ, γ, κ denote the transmission rate, the infectiousness emergence rate, the infectiousness disappearance rate, and the removal rate (due to death or recovery) of hospitalized cases, respectively. Parameters θ and π denote the ratio of hospitalized cases out of all infected cases and the proportion of deaths out of hospitalized cases, respectively. Thus, the overall case fatality rate (or infection fatality rate) equals θ π. All parameters are constant except β (t) being time-varying. A similar model was used for the Amazonas and Manaus COVID confirmed mortality data.

We follow closely our previous work for fitting a “flexible” time varying transmission rate β(t) (see ^2–4). It requires defining β(t) = exp (cubic_spline) as an exponential cubic spline with n nodes evenly distributed over the study period. We set time step size as one day and integrated for one day and yielded the simulated daily deaths D_t. We defined the reported deaths as C_t, and then we have

Here, τ denotes the overdispersion, and accounts for the measurement noise due to surveillance and heterogeneity among individuals.

The parameter values of, σ, γ, κ are taken as 365/2, 365/3.5, and 365/12 per year, respectively. Thus the mean latent period and mean infectious period (the reciprocals of, σ and γ)) are 2 and 3.5 days, respectively. The generation time GI, the sum of mean latent period and mean infectious period ⁵, equals 5.5 days, which is in line with three key studies which provided estimated GI ^6,10. From the SARI data we noticed a delay between the first symptom onset and the death at 14 days, which justifies the choice of 12 days (the reciprocal of κ) from loss of infectiousness to death. Namely the onset of infectiousness is 2 days ahead of first symptom onset.

Although we do not explicitly model age-class structure here, the model nevertheless has direct relevance for two reasons:

1. The important Science paper of Earn et al ⁸ was designed especially to show that the changes in R₀(t) have much more impact on the epidemic trajectory than age-structure and should be the first consideration, while the effects of age-structure are in comparison far smaller.

2. Manaus can be qualitatively approximated as a two age-class system. Some 20% of the population are older than 45 years of age and contributes significantly to COVID-19 related deaths, while the remaining 80% of the population are under 45 years of age and have such low IFR that they barely contribute to the deaths. (In fact <7% of total COVID-19 deaths.) Yet the infection dynamics is far more homogeneous. (Buss et al. suppose all age classes have the same uniform prevalence.) It can be shown that the homogeneous system found from averaging two age classes (themselves homogenous in transmission) can explain the epidemic dynamics of the two age-class system (as shown in eg., Magpantay et al. 2019⁹).

We used the Euler-multinomial simulation approach to simulate our model with a time-step size 1 day. We use the popular iterated filtering method (note that at least 25 applications use this method, see a list at https://en.wikipedia.org/wiki/Iterated_filtering) to fit our model to the observed data and yield the maximum likelihood estimates of unknown parameters, including values of nodes in β (t) and π while holding θ=0.1. Namely, we assume 10% of cases are hospitalized (or severe cases), as widely reported. Namely in this work, we only fit death data rather than case data. The fitting is insensitive to the choice of θ, but sensitive to the product of π and of θ, which is equivalent to the infection fatality rate (IFR). We assume uniform prior for all parameters. We hypothesized that the infection fatality rate (IFR) lay in the interval (0.006, 0.01), and fit the daily SARI mortality data with a flexible transmission rate. The method finds the best fitting parameters including IFR.

Reproductive number R₀(t)

For each dataset we plotted , which is the blue line seen in the figures as a function of time (/date). The reproductive number R₀(t) should not be confused with the effective reproductive number R_e(t) = R₀(t)S (t) mentioned in our Note, which we do not plot. The plot of Ro (t) allows us to visualise how the reproductive number changes e.g., as a result of behaviour, or interventions.

For Amazonas and Manaus, we also arrived at estimates for the initial reproductive number in the first phase of the epidemic, using standard techniques⁷. This is shown in Figure S3 which is based on the COVID-19 mortality data from Amazonas government health department reports (Boletim Diario de casos COVID-19 no Amazonas. http://www.fvs.am.gov.br/publicacoes) ¹².

For Figure A1, we used a standard approach, , where M is the moment generating function of the distribution of the generation time which attains a Gamma distribution with a mean of 5.5 days and s.d. 3 days^6,10. We used the internal from March 30, 2020 (when there were 2 deaths) to April 12, 2020. Our results for estimates of R₀ were:

Amazonas 2.202 [1.712 2.782]

Manaus 2.250 [1.787 2.791]

• Download figure
• Open in new tab

Figure A1.

The dashed vertical lines indicate the window in which R₀ was estimated.