Selection of blood donors for estimation of seroprevalence
Each of the eight cities had a monthly quota of 1,000 kits for testing selected donation samples in this study. In order to select more representative samples, we selected blood samples so that the spatial distribution of residential location of selected donors match the spatial distribution of population density in each municipality. More specifically, each city was divided into sub-municipal administrative zones, and the original quota (1,000 kits) was divided into sub-quotas following the populational distribution of the city administrative zones. Starting from the second week of each month, we selected consecutive blood donors based on the geolocation of their residential postcode to fill the sub-quotas. In this way, donations with missing or wrong postal code were considered ineligible for selection.
In Manaus, however, donor postcodes were not reliably collected, so that the number of missing and wrong values make this strategy unfeasible. So, samples were selected consecutively with no postal code restrictions. We also developed a study management system to operationalize this sampling strategy, whereby blood donor postcodes and epidemiological data were automatically extracted and selected. After that, the selected donation sample IDs were released for the research assistant to be separated for testing.
In Brazil, blood donation samples are usually saved for six months, so when serological test kits were made available in July 2020, we could retrospectively select and test frozen samples from February to July. After this period samples were selected and tested in real time. Antibody tests results were not made available to the blood donors themselves.
Blood donors are a convenience sample, and thus may not be representative of the wider population in terms of their risk of SARS-CoV-2 exposure. Supplementary Figs. 17–19 shows a comparison between recorded blood donor demographics and the last available Brazilian census conducted in 2010. Donors differ systematically in age, sex and self-reported skin colour compared to the population, but the income per capita is similar. To account for the differences in the age-sex structure of blood donors, we divide donors in age-sex groups and estimate the prevalence of each age-sex group separately. Then, we calculate the seroprevalence of the population as a weighted sum of the seroprevalences of each age-sex group.
SARS-CoV-2 serology assays
We applied chemiluminescent microparticle immunoassays (CIMA, AdviseDx, Abbott) that detect IgG antibodies against the SARS-CoV-2 nucleocapsid (N) because it was the only automated commercially available kit in Brazil when the study started (July 2020) and we used this kit throughout the study until March in all eight cities except Recife, where we used the kit until February 2021. Supplementary Figs. 20–22 contain the number of tests disaggregated by month, age, sex and the monthly signal-to-cutoff distribution. We also decided to validate the results observed in Manaus, as this represents a unique sentinel population, by retesting all samples in November 2020 using the chemiluminescent microparticle immunoassay (CIMA, AdviseDx, Abbott) that detect IgG antibodies against the SARS-CoV-2 spike (S) protein (see Supplementary Information for the validation analyses).
To determine the test sensitivity, we considered a cohort of 208 non-hospitalised symptomatic SARS-CoV-2 PCR-positive convalescent plasma donors tested within 60 days after symptom onset (Supplementary Table 1). These donors had symptomatic COVID-19 with PCR-confirmed SARS-CoV-2 infection and were recruited to provide convalescent plasma. We found a sensitivity of 90.6% for the anti-N assay using a threshold of 0.49 signal-to-cutoff and 94.0% for the anti-S assay. Specificity for the anti-N assay was 97.5%, with 801 negative results in 821 pre-pandemic blood donation samples5. Sensitivity and specificity for other assay thresholds are shown in Supplementary Table 1. The anti-S assay has a specificity of > 99% 13,14, and we assume 100% in this study. Although the sensitivity of both assays declines through time due to waning of the detected antibodies below the positivity threshold, the anti-S IgG antibodies wane more slowly13,14. As previously highlighted, the convalescent plasma donors cohort overestimates sensitivity for our use case due to spectrum bias. This is because convalescent donors had moderate-to-severe SARS-CoV-2 infection, and thus differ from the whole blood donor population (used to estimate seroprevalence), who are more likely to have had asymptomatic or mild disease.
We subsequently estimated the distribution of time-to-seroreversion, and thus the sensitivity decreasing through time, for the anti-N assay. We first calculated this in the convalescent donors, in whom the date of symptom onset is known, and whose blood samples were collected longitudinally during convalescence. As such the time-to-seroreversion distribution was computed after accounting for right censoring. However, again due to spectrum bias, the extrapolation of antibody waning from convalescent donors to whole blood donors is unlikely to be valid. As such, we obtained a second cohort of repeat blood donors in Manaus that provided multiple donations during the 2020–2021 period. These donors are expected to have the same antibody dynamics as the seroprevalence cohort, as they are drawn from the same population and have predominantly mild or asymptomatic infections. However, in this group the time of infection is unknown, as infection is inferred by serostatus alone. The procedure to manage this problem is described below.
Methods used to estimate the time-dependent sensitivity
Let \({se}_{0}\) be the sensitivity measured shortly after symptomatic infection (that is, the probability of an infected individual seroconverting to an S/C above the threshold), and \({p}^{+}\left[m\right]\) the probability of a donor remaining positive \(m\) weeks after seroconversion (given that the donor seroconverted). Then, the sensitivity of the test \(m\) weeks after seroconversion is then \(s{e}_{0}\times {p}^{+}\left[m\right]\) for a given donor. In this section we describe the procedure used to determine \({p}^{+}\left[m\right]\) from repeat blood donor data, for which time of infection and time of seroreversion are unknown. The seroreversion correction model described in the next section uses the estimate of \({p}^{+}\left[m\right]\) to calculate the seroprevalence accounting for seroreversion.
The criteria to select repeat blood donors were: 1. at least one positive test, indicating SARS-CoV-2 infection, 2. at least one subsequent blood sample, in order to interpolate the date of seroreversion, and 3. falling S/C between these two samples, because one of the samples used to define the interpolation curve may have occurred before the peak S/C, hence the half-life and the date of seroreversion cannot be estimated. Therefore, all selected donors had at least one positive sample and at least one subsequent sample (positive or negative) with smaller S/C.
To calculate \({p}^{+}\left[m\right]\), we first estimate the date of seroreversion for each repeat blood donor using an exponential interpolation (a linear interpolation in the log scale). We choose an exponential interpolation because an exponential decay is frequently used to model antibody dynamics16. When seroreversion is interval-censored, i.e., a donor that has a positive test subsequently becomes negative, we interpolate an exponential curve that passes through the last positive sample and the first negative sample. Otherwise, when seroreversion is not interval-censored, then it is right-censored (a donor remains positive on their last sample), in which case we extrapolate an exponential line through the last two positive samples and project this forward. As such, the estimated instant of seroreversion for blood donor \(i\) (denoted as \({{t}_{i}}^{-}\)) is the point where the interpolation curve crosses the threshold for a positive test. The interpolation procedure is illustrated in Supplementary Fig. 23. The proposed method may overestimate \({{t}_{i}}^{-}\) if an S/C used to define the interpolation curve was sampled shortly after seroconversion before the peak S/C was reached, since in this case the S/C curve does not behave as an exponential, leading to an overestimated half-life. To partially overcome this problem, we discard donors that do not serorevert within 106 weeks (2 years) after seroconversion.
After estimating \({{t}_{i}}^{-}\), for each blood donor \(i,\) we compute the probability distribution of the date of seroconversion for that donor, \({p}_{i}\). For this, we identify the earliest and latest possible date of seroconversion \({t}_{min}\) (the date of the last negative result before seroconversion or March 1st 2020 if the donor has no positive results before seroconversion) and \({t}_{max}\) (the date of the first positive result). The relative probability of seroconversion within this window depends on the incidence of seroconversions due to SARS-CoV-2 infection for the cohort of repeat donors, denoted \({u}_{repeat}\left[n\right]\). To estimate this quantity, we calculate the histogram of the date of first positive donation for repeat blood donors and then apply a 30-day moving average. As a sensitivity analysis, we also calculate \({u}_{repeat}\left[n\right]\) by computing the histogram of the date of onset of SARI deaths observed in Manaus, and applying to it a 7-day window moving average, yielding similar results (Supplementary Figs. 3 and 7).
The distribution of the date of seroconversion is obtained by truncating the incidence curve of repeat blood donors \({u}_{repeat}\left[n\right]\) in the interval \([{t}_{min},{t}_{max}]\) and renormalizing the distribution. We then generate 1,000 samples of the instant of seroconversion \({t}_{i}^{+}\sim {p}_{i}\left[n\right]\) and compute the 1000 sample delays between seroconversion and seroreversion \(\varDelta {t}_{i}={{t}_{i}}^{-}-{{t}_{i}}^{+}\).
The probability of the delay between seroconversion and seroreversion being \(n\ge 1\) days (denoted as \({p}_{day}^{-}\left[n\right])\) is calculated with the empirical histogram of the \(1000\times {N}_{donors}\) samples of \(\varDelta {t}_{i}\), \(i=1,\cdots , {N}_{donors}\). The distribution \({p}_{day}^{-}\left[n\right]\) is then binned into weeks by taking the average of \(({p}_{day}^{-}\left[7n\right], {p}_{day}^{-}\left[7n+1\right],\cdots ,{p}_{day}^{-}\left[7n+6\right])\) for \(n\ge 1\). The resulting distribution, denoted as \({p}_{week}^{-}\left[n\right]\), represents the probability of seroreversion exactly \(n\) weeks after seroconversion.
Finally, the probability of a donor remaining positive \(n\) weeks after seroconversion
\({p}^{+}\left[n\right]\) (i.e., the probability of a donor seroconverting after week \(n\)) obtained through
$${p}^{+}\left[n\right]=1-{\sum }_{k=1}^{n}{p}_{week}^{-}\left[k\right].$$
The presented method is summarised in Algorithm 1 in Supplementary Information.
Our proposed seroreversion correction model
Here, we present a Bayesian model that draws posterior samples from the incidence over time corrected by sensitivity, specificity and seroreversion using as input the estimated curve for \({p}^{+}\left[n\right]\), the number of weekly positive tests and total number of tests. Even though the main output of the model is the incidence at week \(n\) for age-sex group \(a\) (denoted as \(u[n,a]\)), the seroprevalence at week \(n\) for group \(a\) can be calculated from \(u[n,a]\) as \(\rho \left[n,a\right]={\sum }_{k=1}^{n}u[k,a]\). For simplicity, the proposed model ignores the delay between infection and seroconversion, as it should have small impact on the estimate of \(u[n,a]\). To define the age-sex groups, age was discretized in the intervals 16–24, 25–34, 35–44, 45–54 and 55–70.
Assuming that the sensitivity \(s{e}_{0}\) and specificity \(sp\) of the assay are independent of the age-sex group, the probability of a random person from age-sex group \(a\) being tested positive at week \(n\), denoted as \(\theta [n,a]\), is
$$\theta \left[n,a\right]=s{e}_{0}{\sum }_{k=1}^{n}{p}^{+}[n-k]u[k,a]+\left(1-sp\right)\left(1-{\sum }_{k=1}^{n}u[k,a]\right).$$
The derivation of the expression above is presented in Supplementary Information. The left term \(s{e}_{0}{\sum }_{k=1}^{n}{p}^{+}[n-k]u[k,a]\) represents true positives (previously infected donors that are still seropositive), while the right term \(\left(1-sp\right)\left(1-{\sum }_{k=1}^{n}u[k,a]\right)\) represents false positives (uninfected donors that test positive).
Let us denote as \({T}^{+}[n,a]\) and \(T\left[n,a\right]\) respectively the number of positive tests and the total number of tests for week \(n\) and age-sex group \(a\). Given \(\theta [n,a]\), the probability distribution of \({T}^{+}[n,a]\) is
$${T}^{+}\left[n,a\right]|\theta \left[n,a\right]\sim Binomial\left(T\left[n,a\right],\theta \left[n,a\right]\right).$$
We use a Bayesian framework to draw posterior samples from \(u\) assuming a non-informative prior, but limiting the final seroprevalence in the interval \([0,b]\), where \(b\) is a fixed input of the algorithm that can be 1 or 2 depending on whether reinfections are allowed, and we use \(b=2\) in this work. Instead of defining a prior distribution for \(u[n,a]\) directly, we decompose it into \(u\left[n,a\right]={\rho }_{max}\left[a\right]{u}_{norm}\left[n,a\right],\) where \({\rho }_{max}\left[a\right]\sim Uniform(0,b)\) sets the upper bound of the final prevalence to \(b\) and \({u}_{norm}\left[:,a\right]\sim Dirichlet(\text{1,1},\dots ,1)\) is the normalised incidence which sums to 1. This decomposition is equivalent to assuming a uniformly distributed prior for \(u\left[:,a\right]\) in the simplex \(0\le {\sum }_{n=1}^{N}u\left[n,a\right]\le b\) with \(u\left[n,a\right]\ge 0 \forall n\).
After drawing posterior samples from \(u[n,a]\), we calculate the seroprevalence at week \(n\) for age-sex group \(a\)as \(\rho \left[n,a\right]= {\sum }_{k=1}^{n}u[k,a]\)and then compute the age-sex weighted seroprevalence \(\rho \left[n\right]\), given by
$$\rho \left[n\right]={\sum }_{a=1}^{M}w\left[a\right]\rho \left[n,a\right], w\left[a\right]=\frac{\text{p}\text{o}\text{p}\left[a\right]}{{\sum }_{k=1}^{M}\text{p}\text{o}\text{p}[k]},$$
where \(\text{p}\text{o}\text{p}\left[a\right]\) is the population for the age-sex group \(a\) in the corresponding city.
The presented Bayesian model is summarised in Algorithm 2 in Supplementary Information. The posterior samples are drawn using a Monte-Carlo Markov Chain (MCMC) algorithm with 100,000 iterations.
The incidence returned by the model was validated through posterior predictive checks by randomly selecting 1,000 samples from \(u\left[n,a\right]\)and drawing samples from \({T}^{+}\left[n,a\right]|\theta \left[n,a\right]\sim Binomial\left(T\left[n,a\right],\theta \left[n,a\right]\right)\). The resulting crude seroprevalence is then compared with the measured crude seroprevalence (Supplementary Fig. 24).
The proposed Bayesian seroreversion correction model can be seen as an improvement on that presented in Buss et al5. The model in Buss et al. assumed a parametric form for time to seroreversion and derived the parameters by assuming an increasing cumulative incidence in the repeated cross-sectional samples of blood donors in Manaus. Here we derived the distribution directly from repeat blood donors without assuming any parametric form. Also, Buss et al. applied the seroreversion correction method to the measured seroprevalence corrected for sensitivity, specificity and reweighted by age and sex, while here we estimate the seroprevalence in each age group separately, allowing the identification of non-homogeneous incidence in different age groups.
Despite these differences, the results presented here (Table 1, Fig. 2D) are compatible with the seroprevalence estimates of 28.8% and 76.0% respectively for São Paulo and Manaus in Buss et al. The proposed seroreversion method also differs from other methods in the literature16,28 in that we use the incidence curve to estimate the time-dependent sensitivity instead of the deaths or confirmed cases curve, producing a seroprevalence that does not depend on case reporting and that can be reliably inferred in epidemics where the IFR changes with time, as was the case in Manaus.
Estimating the infection fatality rate for December 2020
We estimate the infection fatality rate using total deaths due to severe acute respiratory infection (SARI) which include PCR- and clinically-confirmed SARS-CoV-2 infection as well as SARI deaths without a final diagnosis, excluding SARI deaths confirmedly caused by other aetiologies. This approach reduces under-reporting, particularly in 2020 when testing was not widely available. To estimate the IFR in 2020, we use the seroprevalence estimated by our model for December 16 2020, and select only SARI deaths with symptoms between March 1st and December 15, 2020. For the first wave of COVID-19 that occurred in the eight cities, we estimate the number of cases as the age-specific population size (https://demografiaufrn.net/laboratorios/lepp/) multiplied by the estimated seroprevalence in the corresponding age group. We propagate the uncertainty in the prevalence estimate through the calculation of IFR.
Let \(\rho \left[a\right]\) and \(pop\left[a\right]\) be the prevalence and the population estimated for age group \(a\). We assume a uniform distribution in the interval \(\left[\text{0,1}\right]\) as a non-informative prior for \(IFR\left[a\right]\), and the number of deaths \(D\left[a\right]\) observed for each age group \(a\) is Binomial-distributed with size \(\lfloor \rho \left[a\right]\times pop\left[a\right]\rfloor\) (the number of infections) and probability \(IFR\left[a\right]\). For each sample of \(\rho \left[a\right]\), we draw a sample of the posterior distribution of \(IFR\left[a\right]\), given by
$$Beta(1+D\left[a\right], 1+\lfloor \rho \left[a\right]\times pop\left[a\right]\rfloor -D[a\left]\right)$$
and compute the median, interquartile ranges and 95% confidence intervals of the IFR by retrieving the quantiles of the posterior distribution. The same analysis was repeated to compute the Infection-Hospitalisation Rate (IHR), by using the number of hospitalisations with SARI instead of the number of deaths.
To validate the obtained IFRs, we also estimate the IFRs using the measured prevalence corrected only by the sensitivity and specificity of the assay, without explicitly accounting for seroreversion. In this validation analysis, we use a small threshold of 0.1 signal-to-cutoff to avoid underestimating the prevalence due to seroreversion (see Supplementary Information).
Estimating the infection fatality rate for the Gamma VOC
We estimate the IFR and the attack rate separately for the second, Gamma-dominant, SARS-CoV-2 wave that occurred in Manaus. The Gamma variant was first detected in Manaus in November 2020, and its prevalence among PCR-positive patients grew rapidly to 87.0% on January 4 20216. For this reason, it is reasonable to assume that all infections in Manaus that occurred after December 15 are due to the Gamma VOC. The Gamma-dominated wave was characterised by a non-negligible proportion of reinfections6,21,29. It is estimated that 13.6–39.3% of the infections in the second wave of COVID-19 epidemic in Manaus were reinfections21, which is explained by the higher in-vitro reinfection potential of Gamma30 and partial immunity waning eight months after the first surge. Thus, to calculate the attack rate and IFR of the Gamma-dominated wave, reinfections must be considered.
However, estimating the incidence of reinfections among positive donors is not straightforward - as a positive result may be either primary infection or reinfection, and these cannot be distinguished using a single test result. For this reason, it was not possible to obtain a point estimate for the number of infections that happened in the second wave in Manaus. To overcome this problem, we calculate upper bounds for the attack rate of the Gamma-dominated wave in Manaus (that is, the incidence between December 2020 and March 2021), and conversely lower bounds for the IFR of the Gamma VOC.
We additionally computed the IFR obtained using the seroprevalence estimated by the model. It is worth noting that our seroreversion correction model only estimates the incidence among seronegative individuals, thus an S/C boosting due to reinfection is not detected by our method. As such, our model estimates the seroprevalence assuming there are no reinfections among positive individuals, underestimating the size of the second wave in Manaus.
To limit the maximum number of reinfections allowed among positive donors, we estimate the incidence using a threshold of 1.4 signal-to-cutoff (the upper threshold recommended by the manufacturer) instead of 0.49 signal-to-cutoff (the lower threshold recommended by the manufacturer), and correct for sensitivity based on 163 true positives and 30 false negatives in the plasma donors cohort. Since the specificity of the test using a threshold of 1.4 is 99.9%, and since it is not straightforward to take the specificity into account when reinfections are allowed, we do not correct for specificity in this analysis.
We compute the monthly number of positive tests \({T}^{+}\left[n\right]\) from December 2020 to March 2021 for a given age-sex group, as well as number of true positive (TP) and false negatives (FN) from convalescent plasma donors and the number of False Positives (FP) and true negatives (TN) from the pre-pandemic blood donors cohort in Manaus (Supplementary Table 1). Then, we draw posterior samples from the raw estimated monthly incidence and the seroprevalence in December and correct for the sensitivity and specificity of the assay. The lower bound for the attack rate of the Gamma VOC is the estimated incidence between December 2020 and March 2021, and the upper bound is the incidence in this period added by the seroprevalence measured in December. Hence, the upper bound is obtained assuming that all individuals that were seropositive in December get reinfected or that all previously seropositive individuals seroreverted in March 2021. The lower bound for the IFR is then calculated using the upper bound of the attack rate and the number of deaths with symptoms onset between December 16 and March 15. This procedure is repeated for each age-sex group independently, and is summarised in Algorithm 3 in Supplementary Information.
The IHR for the Gamma VOC was estimated using the same procedure, but using the number of hospitalisations by SARI instead the number of deaths.
Data availability
All serological data required to reproduce the analyses and the codes used for the main analyses are available at https://github.com/CADDE-CENTRE/seroprevalence_eight_cities.