SARS-CoV-2 correlates of protection from infection against variants of concern

Serum nAb titers against SARS-CoV-2 D614G and BA.1 variants (Lentiviral Pseudovirus Production and Neutralization Assay)

Virus production and pseudovirus neutralization assays were done as previously described ⁵⁰. Briefly, 293T/ACE2.MF cells modified to overexpress human ACE2 (kindly provided by M. Farzan (Scripps Research)) were cultured in DMEM (Gibco BRL Life Technologies) containing 10% heat-inactivated serum (FBS) and 3 μg ml⁻¹ puromycin at 37 °C, 5% CO₂. Cell monolayers were disrupted at confluency by treatment with 0.25% trypsin in 1 mM EDTA (Gibco BRL Life Technologies). The SARS-CoV-2, Wuhan-1 spike, cloned into pCDNA3.1 was mutated using the QuikChange Lightning Site-Directed Mutagenesis kit (Agilent Technologies) and NEBuilder HiFi DNA Assembly Master Mix (NEB) to include D614G (wild-type) or lineage defining mutations for Delta (T19R, 156-157del, R158G, L452R, T478K, D614G, P681R and D950N) and), Omicron BA.1 (A67V, 69-70del, T95I, G142D, 143-145del, 211del, L212I, 214EPE, G339D, S371L, S373P, S375F, K417N, N440K, G446S, S477N, T478K, E484A, Q493R, G496S, Q498R, N501Y, Y505H, T547K, D614G, H655Y, N679K, P681H, N764K, D796Y, N856K, Q954H, N969K, L981F), Omicron BA.2 (T19I, L24S, 25-27del, G142D, V213G, G339D, S371F, S373P, S375F, T376A, D405N, R408S, K417N, N440K, S477N, T478K, E484A, Q493R, Q498R, N501Y, Y505H, D614G, H655Y, N679K, P681H, N764K, D796Y, Q954H, N969K). Pseudoviruses were produced by co-transfection in 293T/17 cells with a lentiviral backbone (HIV-1 pNL4.luc encoding the firefly luciferase gene) and either of the full-length SARS-CoV-2 spike plasmids with PEIMAX (Polysciences). Culture supernatants were clarified of cells by a 0.45μM filter and stored at −70 °C. Plasma samples were heat-inactivated and clarified by centrifugation. Pseudovirus and serially diluted plasma/sera were incubated for 1 h at 37 °C, 5% CO₂. Cells were added at 1 × 10⁴ cells per well after 72 h of incubation at 37 °C, 5% CO₂, luminescence was measured using PerkinElmer Life Sciences Model Victor X luminometer. Neutralization was measured as described by a reduction in luciferase gene expression after single-round infection of 293T/ACE2.MF cells with spike-pseudotyped viruses. Titers were calculated as the reciprocal plasma dilution (ID₅₀) causing 50% reduction of relative light units.

Noting that we measured neutralization titer using a lentiviral-backboned pseudovirus neutralization assay. A systematic review of Omicron neutralization data showed that pseudovirus neutralization assays tend to report higher neutralizing titers compared to live-virus assays. The titer drops from wild type to Omicron also tend to be less pronounced for pseudovirus platforms, suggesting the pseudovirus assay may underestimate Omicron’s capability to escape neutralization

sars-CoV-2 spike enzyme-linked immunosorbent assay (ELISA)

For ELISA, Hexapro SARS-CoV-2 full spike protein with the D614G substitution was expressed in Human Embryonic Kidney (HEK) 293F suspension cells by transfecting the cells with the respective expression plasmid. After incubating for 6 days at 37°C, proteins were first purified using a nickel resin followed by size exclusion chromatography. Relevant fractions were collected and frozen at -80°C until use. Two µg/mL of D614G spike protein was used to coat 96-well, high-binding plates (Corning) and incubated overnight at 4°C. The plates were incubated in a blocking buffer consisting of 1x PBS, 5% skimmed milk powder, 0.05% Tween 20. Plasma samples were diluted to 1:100 starting dilution in a blocking buffer and added to the plates. IgG secondary antibody (Merck) was diluted to 1:3000 in blocking buffer and added to the plates followed by TMB substrate (Thermofisher Scientific). Upon stopping the reaction with 1 M H2SO4, optical density (OD) was measured at 450 nm. The monoclonal antibodies (mAbs) CR3022 and Palivizumab were used as the positive and negative controls respectively.

Mediation analyses and household transmission model fitted to observed infections in the cohort

Here we blend concepts from causal inference and infectious disease transmission models. The no interference assumption in causal inference stipulates that the outcome of an individual does not depend on the outcome of others, which is often violated in infectious disease dynamics ^36,52,53. This is because the spread of infectious diseases requires pathogens to be transmitted from one host to another. In other words, the infection outcome of one individual inherently depends on the infection outcome of others, and this is particularly pronounced in a household setting ³⁶. The “dependent happening” nature of infectious disease dynamics violates the no interference assumption. As a result, the traditional regression approach for causal inference analysis cannot be applied to infectious disease outcomes among individuals who can in theory transmit the disease from one to another. To overcome this, Halloran and Struchiner. ³⁶ introduced the probability of infection conditional on exposure to already infected individuals (transmission probability), as the causal parameter. Using this proposed framework, we can investigate how the presence/absence of pre-existing immunity along with the immunologic marker of interest could modulate probability of infection, after adjusting for levels of exposure to the infectious source(s). The corresponding causal inference framework requires modelling the transmission process explicitly. Under this framework, we conduct mediation analyses to investigate how nAb titers against variants at the start of a SARS-CoV-2 wave correlate with SARS-CoV-2 transmission risk, using the Delta and Omicron BA.1/2 waves as examples ^34,35. We focus on the Delta and Omicron subgroup participants who have had a single or no prior infection, and fit a chain-binoal model to their infection outcomes during the corresponding Delta/Omicron wave ⁵⁴. Specifically, we introduce the causal parameters:

• : the per-contact SARS-CoV-2 household transmission probability from infected individual i to individual j in household k.

• : the overall probability of acquiring SARS-CoV-2 infection from outside the household by individual j of household k (probability of infection from the community).

We use e_j to indicate individual j’s prior SARS-CoV-2 infection history, with e_j = 0 representing no prior infection reported before the start of Delta/Omicron BA.1/2 and e_j = 1 representing one prior infection by the start of Delta/Omicron BA.1/2 wave. A prior SARS-CoV-2 infection (e_j = 1) would induce immunologic responses, measured by a set of immune markers (i.e., candidate mediators) {m_j|e_j = 1} (e.g., nAb titers level). Then the household transmission probability can be expressed as a function of prior infection status e_j, immunologic mediators of SARS-CoV-2 transmission probability {m_j |e_j = 1} and additional adjustment terms {c_i, c_j, c_k}, representing a set of potential confounding factors of individual i, individual j, and household k (eg, age of the donor and/or recipient, comorbidities, household size, etc). Similarly, the community infection probability can be expressed as a function of individual j ‘s prior exposure history e_j, immunological markers {m_j |e_j = 1}, and additional adjustment terms {c_j}, representing a set of potential confounding factors of individual j (e.g., age or comorbidities).

The causal diagram of the mediation analysis framework is shown in Fig. 3. We fit a household transmission model to the imputed household transmission chains based on an Expectation-maximization (EM) algorithm (detailed in Section 4). Specifically, for the Delta/Omicron BA.1/2 wave, if we look into a specific household k of size N, there are a total of n individuals infected belonging to L distinct chains of transmission due to L independent introductions of SARS-CoV-2 into the household. The uninfected individuals are N − n. We denote the likelihood of any individual j of household k having the observed infection status over the Delta/Omicron BA.1/2 wave (i.e., either infected or not) in a particular realization of the model. There are a few scenarios to write down :

• Within a given transmission chain l ∈ L, the initial generation always has an individual j acquiring infection from the general community (outside the household k). Thus, the probability of individual j being infected is if j is the first individual to be infected in the chain.

• For infected individual j in the first generation of transmission chain l, i.e.,, this individual would have to escape infection risk from the general community but get infected by the infected household member of . Thus, the probability of individual j being infected can be written as .

• For infected individual j in transmission chain l with generation greater than 1, i.e., g_I > 1, this individual has escaped infection risk from the general community as well as infected individuals i two generations away but got infected by an infector i^′ of j’s previous generation on the same transmission chain l. Thus, the probability of individual j being infected can be written as .

• For uninfected individual j within household k, this individual has escaped infection risk from the general community as well as all the n infected individuals within the same household. Thus, the probability of individual j remaining uninfected can be written as .

Then, within household k of size N, we can express the likelihood of transmission chain l as ; the likelihood of observing all infections within k can be expressed as ; the likelihood of observing N − n uninfected individuals can be expressed as . Putting these together, the likelihood of observing one realization of the imputed (details of the EM imputation method described in the next section) households’ transmission trees for Delta/Omicron wave can be expressed as:

Where the likelihood of a given household transmission chains configuration can be expressed as:

In the remainder of the section, we will consider a few versions of the transmission model with slightly different implementations for and .

Model 1: waning model for prior exposure with serologically ascertained Delta and Omicron wave infections

This is the transmission model presented in the main analysis of the manuscript (results of the model shown in Table 2. In this model, we consider that protection from prior infection unexplained by nAb titers wanes over time but is not dependent on the variant responsible for prior infection (i.e. prior D614G or Beta infections for the Delta wave analysis, and prior D614G, Beta, or Delta infections for the Omicron wave analysis). Additionally, in this model, both the Delta and Omicron wave infections were ascertained by serology based on approach describe in a prior session in Methods.

More specifically, for the Delta wave, and can be expressed as:

As described before, e_j indicates individual j’s prior SARS-CoV-2 infection history, with e_j = 0 representing uninfected individuals at the start of the Delta wave, e_j = 1 representing one prior infection, and ϵ representing the effect size of the immune protection by prior infection not mediated through anti-D614G nAbs (direct effect, Table 2). Δt is the elapsed time between prior infection and blood draw 5 (the blood draw taken prior to the Delta wave which we use in this model) and τ is the waning half-life of ϵ (direct effect, Table 2). represents the anti-D614G nAb titer at blood draw 5 and represents the effect size of in mediating infection probability against the Delta wave infection (mediator effect, Table 2) at blood draw 5. While represents the quantity of anti-D614G nAbs waned from peak level (measured as the highest anti-D614G nAb titer level among the first 5 blood draws) to that at BD5 and represents the effect size of in mediating transmission probability against the Delta wave infection (mediator effect, Table 2) at blood draw 5. Note that the term only exists when e_j = 1.

We further evaluate whether breakthrough infections have reduced infectiousness compared to primary infections and may in turn affect . We use e_i to indicate individual i’s (the donor) prior SARS-CoV-2 infection history (e_i = 0 means no infection, and e_i = 1 represents one prior infection at the start of Delta wave). Further, λ represents the effect size of prior infection (in i) in reducing the infectiousness of reinfections.

We also consider confounding factors for donor i and recipient j, where c_i and represent infector i’s confounding factor (i’s age, sex) and effect size, respectively; c_j and represent j’s confounding factor (j’s age/sex-specific susceptibility (biology), age/sex- and site-specific susceptibility (behavioral), HIV infection status) and effect size, respectively; c_k and represent household k’s confounding factor (household size) and effect size, respectively. Lastly, α_s and β_s are logits of the baseline risks for household and community exposures. All parameters’ effect sizes are measured in the log of odds ratios.

Similarly, for the Omicron BA.1/2 wave, and can be expressed as:

As described before, e_j indicates individual j’s prior SARS-CoV-2 infection history, with e_j = 0 representing individual j remained naïve to SARS-CoV-2 at the start of Omicron BA.1/2 wave while e_j = 1 representing individual j had one prior infection at the start of Omicron BA.1/2 wave and ϵ represents the effect size of the immune protection by prior infection not mediated through anti-D614G nAbs (direct effect, Table 2). Δt is the elapsed time between prior infection and blood draw 8 (the blood draw taken prior to the Omicron BA.1/2 wave) and τ is the waning half-life of ϵ (direct effect, Table 2). Here we consider that parameter τ is shared between the Delta and Omicron wave and will be jointly estimated (described in the next session in Methods). represents the anti-BA1 nAb titer at blood draw 8 and represents the effect size of in mediating transmission probability against the Omicron BA.1/2 wave infection (mediator effect, Table 2) at blood draw 8. While represents the difference in titer from anti-D614G nAb to anti-BA1 nAb at blood draw 8 and represents the effect size of in mediating transmission probability against the Omicron BA.1/2 wave infection (mediator effect, Table 2) at blood draw 8. Note that the term only exists when e_j = 1. All other parameters have the same definition of the Delta wave.

are estimated through maximizing the likelihood function L for each of the 100 bootstrapped realizations and bootstrap mean and confidence intervals are calculated for each of the parameters.

Model 2: Sensitivity analysis considering variant-specific prior exposure for the direct effects

A potential confounding factor in understanding the waning of protection through direct effects is the diversity of prior SARS-CoV-2 exposures, with the dominance of D614G variant in the first wave, Beta variant in the second wave, and Delta variant in the third wave (Fig. 1). The effectiveness of protection may vary depending on the specific variant of prior exposure that induced the immune response at play. We conducted a sensitivity analysis (Model 2) employing a variant-specific model for the direct effects, which accounted for distinct types of SARS-CoV-2 variants conferring prior immunity, instead of considering generic a waning model. Specifically, in Model 2, we considered a more complex version of Model 1, where protection from prior infection depends on the type of infecting variant (i.e. prior D614G or Beta infections for the Delta wave analysis, and prior D614G, Beta, or Delta infections for the Omicron wave analysis). We consider waning in neutralizing titers as in Model 1, but we eliminate waning in the effect of prior infection that is not captured by neutralizing titers. More specifically, for the Delta wave, and can be expressed as:

Here, indicates individual j, prior to the Delta wave, was infected with D614G (Beta) variant. If , individual j was naïve at the beginning of the Delta wave. and represent the effect size of immune protection by prior D614G and Beta infection not mediated through anti-D614G nAbs, respectively.

For the Omicron wave, and can be expressed as:

Here, indicates individual j, prior to the Omicron wave, was infected with D614G (Beta, Delta) variant. If , individual j was naïve at the beginning of the Omicron wave. and represent the effect size of the immune protection by prior D614G, Beta and Delta infection not mediated through anti-D614G nAbs, respectively.

Additionally, similarly to Model 1, both the Delta and Omicron wave infections were ascertained by serology for Model 2. All other settings of Model 2 were kept the same as Model 1. The results of the Model 2 are presented in Extended Data Table 2.

Our analysis revealed that for both the Delta and Omicron waves, more recent variants conferred stronger protection than earlier variants, albeit with overlapping confidence intervals (Extended Data Table 2). This temporal trend aligns with the expectations of the waning model. Both waning and variant-specific immunity may modulate the direct effects of prior immunity; however, our study lacked sufficient statistical power to jointly estimate the relative contributions of these two factors. Full estimates of this sensitivity analyses are presented in Extended Data Table 2.

Model 3: Sensitivity analysis with Delta wave infections ascertained by PCR and/or serology

For Model 1, both the Delta and Omicron wave infection outcomes were inferred using the kinetics of anti-nucleocapsid antibodies from longitudinal serologic sampling, as detailed in previously published studies of the PHIRST-C cohort ^30,32. This approach for inferring infections based on serology was calibrated against virological evidence of infection during the Delta wave, established through twice-weekly rRT-PCR tests regardless of symptom presentation. However, it should be noted that this calibration did not achieve perfect concordance; the serology approach demonstrated 93% sensitivity and 89% specificity when compared to infections identified by rRT-PCR tests ³⁰. To address the uncertainties arising from the imperfect concordance between the two approaches for ascertaining infections, we conducted a sensitivity analysis (Model 3) for the Delta wave, where we considered infections based on rRT-PCR positivity and/or anti-nucleocapsid antibody serology. We identified an additional 17 infections during the Delta wave through this more sensitive infection ascertainment approach, bringing the total number of Delta wave infections to 290. All other settings of Model 3 were kept the same as Model 1. The results of the Model 3 are presented in Extended Data Table 3.

Notably, estimates of the direct and indirect effects of the mediation analysis were comparable between this sensitivity analysis and the main analysis (compare Extended Data Table 3 to Table 1). These findings provide support for the utilization of anti-nucleocapsid serology to ascertain Omicron BA.1/2 wave infections in the studied cohorts, in a period where twice-weekly rRT-PCR testing was not available and confirms the robustness of our CoP analyses.

Model 4: D614G spike binding antibodies as mediators of protection

We conducted sensitivity analysis (Model 4) to explore the role of D614G spike binding antibodies (referred to as bAb hereafter), as potential correlates of protection for both Delta and Omicron infections. Employing an in-house enzyme-linked immunosorbent assay (ELISA), we quantified the level of D614G spike bAb based by measuring absorbance at 450nm at an optical density (OD) at peak levels and BD5 (DB8) for the Delta (Omicron) wave analysis (Extended Data Fig. 3). The reduction in binding antibody levels from peak (ΔbAb^W) was determined as the difference between OD values at peak and BD5 (BD8) for the Delta (Omicron) wave (Extended Data Fig. 3).

Model 4 builds on Model 2 but replaces nAb titers with D614G spiking binding ELISA readouts as mediators of protection, in order to compare the protection afforded by neutralizing vs binding antibodies. More specifically, for the Delta wave, and can be expressed as:

Here, represents the D614G spike binding antibodies ELISA readout at blood draw 5 and represents the effect size of in mediating transmission probability against the Delta wave infection at blood draw 5. Further, represents the drop from peak D614G spike binding antibodies readout prior to blood draw 5 (measured as the highest D614G spike binding Ab titer level among the first 5 blood draws) to that at blood draw 5 and represents the effect size of in mediating transmission probability against the Delta wave infection at blood draw 5.

For the Omicron wave, and can be expressed as:

Here, represents the D614G spike binding antibodies ELISA readout at blood draw 8 and represents the effect size of in mediating transmission probability against the Omicron wave infection at blood draw 8. While represents the drop from peak D614G spike binding antibodies readout prior to blood draw 8 (measured as the highest D614G spike binding Ab titer level among the first 8 blood draws) to that at blood draw 8 and represents the effect size of in mediating transmission probability against the Omicron wave infection at blood draw 8. All other settings of Model 4 were kept the same as Model 2. The results of the Model 4 are presented in Extended Data Table 4.

We found that binding antibody levels at BD5 (BD8) correlate with protection against Delta (Omicron) wave infections: the risk of infection decreased by 74% (95% CI 41% – 88%) and 40% (95% CI 33% – 54%) per unit increase in OD value for the Delta and Omicron wave analyses, respectively. Conversely, the decline in bAbs from peak levels to BD5/BD8 (ΔbAb^W) demonstrated no contribution to the overall protection, with risk reduction per 10-fold increase: - 2% (95%CI: -91% – 55%) for Delta wave infections and -2% (95%CI: -87% – 55%) for Omicron wave infections. These findings underscore the correspondence between waning of binding antibodies and a waning of protection. Furthermore, our estimations indicate that the proportion of protection conferred through D614G spike bAbs at BD5 is 35% (95%CI: 32% – 38%) against Delta wave infections, a figure comparable estimation based on anti-D614G nAbs (37%, 95%CI: 34% – 40%, Extended Data Table 4). Notably, D614G spike bAbs at BD8 accounted for 27% (95% CI: 25% – 29%) of protection against Omicron wave infection, representing a larger proportion compared to anti-BA.1 nAbs (11%, 95%CI: 9% – 12%, Extended Data Table 4).

Transmission chains imputation and parameters estimation based on an Expectation-maximization (EM) algorithm

Here we describe the process to fit the models described in Section 3 to the household infection data. The serologic data available for the Delta and Omicron only provides information on the total number of infections within the household between two blood draws collected before and after the SARS-CoV-2 wave. The data does not provide the details of the transmission chains within the household, the order of infections among infected individuals, nor the infection dates. To account for the uncertainties of the transmission tree structure within households given only the total number of infections, we enumerate and reconstruct all possible transmission chains among the infected individuals, where each infected individual may have been infected by members of their own household or the general community. Supplementary Fig. 1 illustrates all 16 possible configurations of transmission chains for a household with 3 infected individuals. We limited our analysis to households with no more than 6 infected individuals, as the possible configurations of transmission chains among 6 infected individuals already reaches 16,807. Enumeration of all possible transmission chain configurations would be computationally intractable for households with more than 6 infected individuals. Additionally, the probability of each possible transmission chain depends on the parameter estimates of the transmission model described in the previous session in Methods. To address the statistical uncertainties due to unresolved transmission chains (which would affect the statistical confidence of mediation analysis detailed in the prior section), we jointly fit the household transmission model and impute the topological structure of the transmission trees. We use an EM algorithm, as described below ⁵⁵.

To resolve who infected whom within the household in a probabilistic manner, we considered an EM algorithm that iteratively estimates the transmission model parameters through maximizing the likelihood function L as described in Equation (1) in the previous section then updates the imputed probability of each transmission tree configuration within each household based on the fitted transmission model. The process is as follows:

1. Initial imputation of the household transmission trees with equal sampling probability for all configurations: For each household, we randomly sample one transmission tree with equal probability among all transmission tree configurations that are compatible with the number of infections. We iterate through all households so that each household has a simulated transmission tree. We then repeat the imputation 1000 times to obtain 1000 realizations of each household’s transmission tree.

2. Maximization step: We consider the waning parameter τ a hyper-parameter (nonlinear term in equations (3-6), cannot be estimated by logistic regression). For a fixed value of τ, for each of the 1000 realizations of the simulated household transmission chains, we estimate transmission model parameters through maximizing the likelihood function L described in Equation (1). The maximization of the likelihood function is achieved through fitting a logistic regression of the infection/exposure outcomes for all participants using R package “brglm” (version 0.7.2). We then pool the estimates from the 1000 realizations using the “pool” function in the R package “mice” (version 3.16.0). The full likelihood of the combined Delta and Omicron waves fitting in this EM step m can be expressed as

3. Expectation step: for a fixed value of hyper-parameter τ, based on the pooled estimates of the transmission model parameters we calculate the likelihood all configurations of transmission chains within each household based on Equation (2). We use these configuration-specific likelihoods to resample transmission chains: For each household, we randomly sample one transmission tree among all transmission tree configurations with probability proportional to transmission tree likelihood prescribed in Equation (2), given the parameters estimated by the most recent maximization step. We iterate through all households so that each household is assigned one simulated transmission tree. We repeat the process 1000 times to obtain 1000 realizations of the household transmission trees.

4. For each of the fixed value of hyper-parameter τ over a plausible range (30 - 500 days), we iterate over the EM steps (2) and (3) until L_m(τ) converge to the maximum value of the EM algorithm. We scan through the values of τ from 30 to 500 days at 10 days step. The EM algorithm convergence curve is shown in Supplementary Fig. 2 for each of the τ values. The EM algorithm converges at step 50, irrespective of the value of τ. The marginal likelihood of the model at τ, L(τ) is estimated by taking the average of L_m(τ) for EM steps 50 through 100. Supplementary Fig. 3 shows the log of the likelihood L(τ) as a function of τ, based on a spline interpolation. The point estimate of τ is taken from the maximum of log(L(τ)) while the 95% confidence interval is estimated by finding τ values with log-likelihood value at the maximum minus 1.92 (Supplementary Fig. 3).

5. We then take the best estimate of hyper-parameter τ and repeat the EM algorithm till convergence to estimate transmission model parameters as show in Table 2. Same EM algorithm were applied to Model 2-4 for the sensitivity analysis as well.

The “treatment effect” by prior infection is estimated by simulating from the best-fit model. We first sample 1000 realizations of the imputed household transmission trees, with imputation probability proportional to the best estimates of transmission model using the EM algorithm and hyper-parameter τ. For each of the 1000 realizations, we focus on the subset of individuals who had one prior SARS-CoV-2 infection, denoted as S_j = {j| e_j = 1}). We use the fitted transmission model to predict the probability of infection (i.e. or , with) of these non-naïve subsets under three scenarios:

1. Scenario 1: the probability of infection estimated with predictors as reported in the data, denoted as .

2. Scenario 2: a counterfactual scenario where the probability of infection is estimated with predictor e_j = 0 (i.e. a counterfactual naïve individual) whereas all other covariates (confounders) are the same as observed, removing both direct and mediator effects. We denote the infection probability in this counterfactual scenario as .

3. Scenario 3: a counterfactual scenario where the probability of infection is estimated with predictor e_j = 1, but setting (or ), effectively removing the mediator effect of nAb on preventing transmission, but keeping the direct effect. We denote the infection probability in this counterfactual scenario as .

We then calculate the total protection conferred by prior infection as the population average of , based on bootstrap resampling with replacement (maintaining the same number of observations) of each of the 1000 realizations of the household transmission chains. Point estimates and 95% confidence intervals are based on the median and 95% quantiles of 1000 realizations’ estimates.

Similarly, we calculate the proportion of protection mediated by nAbs as the population average of . We use the same bootstrapping approach as for total protection.