Mechanistic modeling of SARS-CoV-2 immune memory, variants, and vaccines

Early waves of the SARS-CoV-2 pandemic were driven by importation events and subsequent policy responses. However, epidemic dynamics in 2021 are largely driven by the spread of more transmissible and/or immune-evading variants, which in turn are countered by vaccination programs. Here we describe updates to the methodology of Covasim (COVID-19 Agent-based Simulator) to account for immune trajectories over time, correlates of protection, co-circulation of different variants and the roll-out of multiple vaccines. We have extended recent work on neutralizing antibodies (NAbs) as a correlate of protection to account for protection against infection, symptomatic COVID-19, and severe disease using a joint estimation approach. We find that NAbs are strongly correlated with infection blocking and that natural infection provides stronger protection than vaccination for the same level of NAbs, though vaccines typically produce higher NAbs. We find only relatively weak correlations between NAbs and the probability of developing symptoms given a breakthrough infection, or the probability of severe disease given symptoms. A more refined understanding of breakthrough infections in individuals with natural and vaccine-derived immunity will have implications for timing of booster vaccines, the impact of emerging variants of concern on critical vaccination thresholds, and the need for ongoing non-pharmaceutical interventions.


27
Towards the end of 2020, several trends developed that challenged the structure of many existing 28 models of SARS-CoV-2 and COVID-19. First, variants of concern (VOC) emerged with phenotypic 29 differences from the ancestral wild-type virus. These differences conferred increased transmissi- In the absence of standardized assays to measure NAbs, normalization against a 125 convalescent serum standard has been suggested as a method for providing greater comparability 126 between results from different assays (Kristiansen et al., 2021). 127 Individuals are assumed to have no protective neutralizing antibodies against SARS-CoV-2 ini-128 tially (Gouma et al., 2021). Upon first infection, individuals draw an initial NAb level from a log-129 normal distribution (log 2 (NAb 0 ) ∼ (0, 2)), following the distribution used in Khoury et al. (2021). 130 Individuals' initial NAb levels are adjusted based on the symptoms they experience throughout 131 their infection, with asymptomatic infection assumed to reduce the average NAb level by 15% 132 and severe infection assumed to increase it by 50%, again following Khoury et al. (2021).  larly, vaccination is assumed to generate an initial NAb level drawn from a lognormal distribution 134 (log 2 (NAb 0 ) ∼ ( , 2)), where depends on the vaccine administered (see Table 2).

135
For individuals who have some prior immunity, both reinfection and vaccination are known 136 to boost the level of NAbs, although the exact degree of boosting is not yet well understood. We 137 provide a parameter that can be used to control the degree of boosting, and set it to a default value 138 of 150% for reinfections. Vaccination is assumed to boost by , where is a vaccine-specific 139 boosting factor (see Table 2).

140
NAbs grow linearly until they reach their peak after three weeks and then follow a two-part ex- parameters to explore different parameters or functional forms. 146 We keep track of the most recent NAb-conferring event for each individual, storing information 147 about the variant with which they were infected or the vaccine they received. In order to map NAb 148 level to protective efficacy, we use vaccine immunogenicity and efficacy trial data as well as data 149 on reinfection in health care workers (see Table 4). We re-normalized the average NAb for each of 150 the vaccine trial cohorts as well as the re-infection datasets to account for waning that may have 151 occurred between collection of NAb assays in the immunogenicity studies and vaccine efficacy trial 152 endpoints using the antibody kinetics model described above, yielding an average 50% reduction 153 in average NAb level compared to those reported in Khoury et al. (2021). 154 We are modeling three types of immunity: protection against infection, protection against 155 symptomatic disease, and protection from severe disease. Once we consider the impact of wan-156 ing immunity, we find that it's challenging to fit protection against infection with a single curve for was last infected with; cross-immunity for different variants is described in the following section. 173 We find that conditional on a breakthrough infection, reduction in symptoms is weakly corre- . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted June 1, 2021.  Different variants can also be associated with different relative probabilities of disease progres-191 sion. Since Covasim also includes age-linked disease progression and mortality probabilities (see 192   Table 2 in Kerr et al. (2020)), each variant is assigned a relative scale factor which is then applied 193 to these age-dependent probabilities.

194
The Covasim model comes pre-populated with parameter values capturing available evidence 195 on the relative transmissibility and severity of several known variants (see Table 1). Likewise, the 196 model captures available evidence on the cross-immunity that infection with wild-type variants 197 confer against infection with other known variants (see Table 1 for parameter values and   was not yet sufficient data to justify tracking individuals' full immune history within the model, so 207 we take a more parsimonious modeling approach based on the most recent infection.

Variant
Relative Relative Relative Source transmissibility severity immunity Relative immunity refers to the reduction in neutralization with each variant relative to wild-type variants; the reduction in protective efficacy against infection or disease is then calculated using the mechanisms described in the Immunity section and depicted in Figure 2. We also incorporate other estimates of cross-immunity effects where available, e.g. Cele et al. (2021) estimate that infection with B.1.351 is associated with a 2.3-fold reduction in neutralization with wild-type variants.

209
The immunity model can account for the effect of single or multi-dose vaccines with variant-specific 210 protection against infection, symptoms, and severe disease in addition to reduced transmissibility 211 of breakthrough infections. As described in the Immunity section above, for individuals with exist-212 6 of 19 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted June 1, 2021. ; https://doi.org/10.1101/2021.05.31.21258018 doi: medRxiv preprint ing immunity, vaccination confers a boost to their current NAb level by a vaccine-specific scaling factor , while individuals without pre-existing immunity are assigned a peak NAb level drawn 214 from a log-normal NAb distribution. For each vaccine, NAbs increase linearly from the day of first 215 dose until the vaccine's peak efficacy, at which point they begin to decay. Second doses of the 216 vaccine boost current NAbs by the same vaccine-specific scaling factor .

217
Variant-specific protection is implemented as an additional scaling factor reducing an agent's 218 NAb level that is input to the mapping shown in Figure 2F. The model provides vaccine parameters 219 for the Pfizer, Moderna, AstraZeneca, Janssen, and Novavax vaccines (see Table 2) along with the 220 capability to add custom vaccine profiles. In addition, the model will be updated to include new 221 default vaccines as these are released.

222
The default first or single dose allocation scheme employed in Covasim is a constant per-day 223 probability configured to a specified day or set of days during the simulation.  Table 3. Description of epidemiological profiles for case study of the population have already been exposed to wild-type infection, B.1.351 gains a stronger hold 249 due to its immune-evading advantage. For this example, we again simulate a population of 100,000 agents under the "controlled trans-253 mission" setting described above. Specifically, we assume that 40% of the population have some . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted June 1, 2021. ; https://doi.org/10.1101/2021.05.31.21258018 doi: medRxiv preprint that AstraZeneca is being distributed to the population, and we calculate the impact of lifting these  The presence of B.1.351 makes a considerable difference to the vaccination threshold that 263 needs to be attained before NPIs could be lifted without resulting in a resurgence of infections.

264
In this scenario, our estimates show that even high levels of vaccination coverage would not com-265 pletely eliminate transmission or disease if B.1.351 were present in the population (see Fig 4).  To illustrate an approach for working with trees in Covasim, we consider three scenarios:   We can build transmission trees directly from Covasim using a built-in analyzer. We can also 280 build a tree by directly processing the transmission log. Depending on the number of nodes, the 281 size of the transmission tree, and the visualization tools available, it may be useful to focus on 282 a sampled tree instead of the full tree. The most basic sampling strategy is to randomly sample 283 transmission events or nodes from the full tree. We could also use more targeted sampling strate-284 gies. One of these strategies is to randomly sample events on specific dates. This could be useful 285 for comparisons with sampled diagnoses and sequencing data. To illustrate this strategy, assume 286 that we sample a number of events every 10 days starting on day 30. We sample 3 events on each 287 of those days from day 30 to 80. Then we increase the sampled events to 5 until day 130. After 288 that we keep sampling 10 events on each of the remaining sampling days. This sampling strategy 289 9 of 19 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted June 1, 2021. ; https://doi.org/10.1101/2021.05.31.21258018 doi: medRxiv preprint Figure 5. Epidemic curves for the scenarios in case study 3. The figures show how a scenario with no vaccination (i.e., Scenario 1) leads to a higher number of infections and deaths than one with vaccination campaigns (Scenarios 2 and 3). The effect of multiple variants is also shown by the increased number of cases in Scenarios 1 and 3, which include both wild and P.1 variants. results in the trees shown in Fig 6. One may also want to sample based on other criteria, for exam-290 ple, from individuals that would render a PCR test positive on a given date, from the death events, 291 etc.

292
It is worth noting that, when using the number of days between events as the branch length, 293 we should expect the sampled tree to look topologically similar to the corresponding phylogenetic 294 tree. We should be careful, however, on the way we join sub-trees that characterize different virus 295 variants, as the assumptions introduced in this process could lead to very different results about 296 common ancestors for nodes in the tree. In some cases, we may have some additional information 297 about some variants arising from mutations of some other variants. In this case, it would be more 298 accurate to join the sub-trees of a descendent sub-tree to an internal node from the parent sub-299 tree. For example, if we believe that P.1 resulted from a mutation of an existing wild infection 300 around day 60, we could select one of the infections occurring on day 60 as the parent of the seed 301 P.1 infections. In this case, the tree would look as is shown in the tree on the right of  . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted June 1, 2021. ; Figure 6. Sampled transmission trees for the scenarios in case study 3. We can randomly sample a given number of events in particular dates to match, for example, actual samples collected from genome sequencing (see for example, Hadfield et al. (2018)). Branch lengths are proportional to number of days between infection events. A visual inspection of the trees allow us to identify patterns in the transmission of competing variants. For example, these trees show only infections with the wild variant early in the epidemic, with P.1 infections starting to show more prominently after 60 days. The tree on the left is built assuming that the wild and P.1 variants have a common ancestor on day 0. The tree on the right simulates a mutation giving rise to a new variant from an existing transmission chain (in this example, we assume that P.1 started from a mutation on one of the infections that occurred around day 60). our scenarios (see Fig 7). Note that these statistics allow us to infer (from our scenarios) that the 309 P.1 variant is more transmissible than the wild variant, showing that the tree topology reflects the 310 (assumed) dynamics. We can also verify that vaccination campaigns can help reduce transmissibil-311 ity, but they have to be timely (in the figure, the mild vaccination campaign that begins after heavy 312 transmission with the wild type variant is ongoing has little effect on the epidemic; however, it does 313 have an important effect on reducing transmission of P.1).

315
As efficacy results from COVID-19 vaccine trials began to appear in late 2020, there was new optimism that the worst of the pandemic may have passed. However, this optimism was quickly Relating NAbs to protective immunity 584 We broke vaccine efficacy into conditional parts to match the stages of the infection process. Vac-585 cine efficacy against infection, VE inf , is the first stage that modulates the probability of infection 586 given exposure. For people who get infected, symptomaticity is modulated by the conditional vac-587 cine efficacy given infection, VE symp|inf . This is defined from the more commonly reported marginal 588 efficacy against symptoms, VE symp in (Eq 1) as: Conditional on being infected and having symptoms, the (conditional) efficacy against severity, 590 VE sev|inf,symp is defined similarly from the commonly reported marginal efficacy against severity In order to estimate the relationship between neutralizing antibodies and correlates of pro-593 tection, we jointly fit four regressions and used empirical bootstrapping with 5,000 samples to 594 estimate the 90% confidence interval: 595 logit(VE natinf, ) = nat inf + inf * (NAb ) + natinf, logit(VE inf, ) = inf + inf * (NAb ) + inf, logit(VE symp, ) = logit(1 − (1 − invlogit( symp|inf + symp|inf * (NAb ))) * (1 − invlogit( inf + inf * (NAb )) + symp, logit(VE sev, ) = logit(1 − (1 − invlogit( inf + inf * (NAb ))) * (1 − invlogit( symp|inf + symp|inf * (NAb )) * (1 − invlogit( sev|symp + sev|symp * (NAb )) + sev, where natinf, , inf, , symp, , and sev, represent the model error terms VE natinf, , VE inf, , VE symp, , and 596 VE sev, respectively represent the reported vaccine efficacy against infection, symptomatic infection, 597 and severe disease from trial , and represents the average level of neutralizing antibodies 598 across participants in trial . This formulation assumes both the efficacy against infection and the 599 conditional efficacy against symptoms given infection are logit-log, and thus the marginal efficacy 600 against symptoms is the more complex function above. The conditional efficacy against severity 601 was modeled similarly.

602
In the absence of individual-level data, we leveraged the cohort estimates from various vac-603 cine trials and convalescent data sets that report both immunogenicity as well as an associated 604 reduction in the risk of infection and symptoms (Table 4). In order to compare the immunogenicity 605 data with the efficacy endpoints, we accounted for any waning that may have occurred across the 606 timescales reported. We find that accounting for waning would reduce the relative NAbs by nearly