# Modeling the evolution of the neutralizing antibody response against SARS-CoV-2 variants after several administrations of Bnt162b2

Quentin Clairon, Mélanie Prague, Delphine Planas, Timothée Bruel, Laurent Hocqueloux, Thierry Prazuck, Olivier Schwartz, Rodolphe Thiébaut, Jérémie Guedj

<figure><img src="image_1.png" /></figure>

To cite this version:

Quentin Clairon, Mélanie Prague, Delphine Planas, Timothée Bruel, Laurent Hocqueloux, et al. Modeling the evolution of the neutralizing antibody response against SARS-CoV-2 variants after several administrations of Bnt162b2. PLoS Computational Biology, 2023. hal-03946556

HAL Id: hal-03946556

https://inria.hal.science/hal-03946556

Submitted on 19 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling the evolution of the neutralizing antibody response against SARS-CoV-2 variants after several administrations of Bnt162b2

Quentin Clairon<sup>1,2,3,✉</sup>, Mélanie Prague<sup>1,2,3,✉</sup>, Delphine Planas<sup>3,4,✉</sup>, Timothée Bruel<sup>3,4 ‡</sup>, Laurent Hocqueloux<sup>5‡</sup>, Thierry Prazuck<sup>5‡</sup>, Olivier Schwartz<sup>3,4 ‡</sup>, Rodolphe Thiébaut<sup>1,2,3 ‡*</sup>, Jérémie Guedj<sup>6</sup> <sup>✉</sup>

1 Université de Bordeaux, Inria Bordeaux Sud-Ouest

2 Inserm, Bordeaux Population Health Research Center, SISTM Team, UMR1219, F-33000 Bordeaux, France

3 Vaccine Research Institute, F-94000 Créteil, France

4 Virus and Immunity Unit, Institut Pasteur, Université de Paris Cité, CNRS UMR3569, Paris, France

5 Service des Maladies Infectieuses et Tropicales, Centre Hospitalier Régional, Orléans, France

6 Université Paris Cité, IAME, Inserm, F-75018 Paris, France

✉These authors contributed equally to this work.

‡These authors also contributed equally to this work.

✉Current Address: Bordeaux Population Health Research Center, Université de Bordeaux, 146 rue Léo Saignat, 33076 Bordeaux cedex, France.

* rodolphe.thiebaut@u-bordeaux.fr

## Abstract

Because SARS-CoV-2 constantly mutates to escape from the immune response, there is a reduction of neutralizing capacity of antibodies initially targeting the historical strain against emerging Variants of Concerns (VoC)s. That is why the measure of the protection conferred by vaccination cannot solely rely on the antibody levels, but also requires to measure their neutralization capacity. Here we used a mathematical model to follow the humoral response in 26 individuals that received up to three vaccination doses of Bnt162b2 vaccine, and for whom both anti-S IgG and neutralisation capacity was measured longitudinally against all main VoCs. Our model could identify two independent mechanisms that led to a marked increase in humoral response over the successive vaccination doses. In addition to the already known increase in IgG levels after each dose, we identified that the neutralization capacity was significantly increased after the third vaccine administration against all VoCs, despite large inter-individual variability. Consequently, the model projects that the mean duration of detectable neutralizing capacity against non-Omicron VoC is between 366 days (Beta variant, 95% Prediction Intervals PI [323; 366]) and 606 days (Alpha variant, 95% PI [555; 638]). Despite a very low protection after three doses, the mean duration of detectable neutralizing capacity against Omicron variants varies between 184 days (BA.5 variant, 95% PI [155; 215]) and 268 days (BA.1 variant, 95% PI [238; 302]). Our model shows the benefit of incorporating the neutralization capacity in the follow-up of patients to better inform on their level of protection against the different SARS-CoV-2 variants as well as their optimal timing of vaccine administration.
## Author summary

Developed vaccines against SARS-CoV-2 have been a turning point against the ongoing Covid-19 pandemic. When the Wuhan virus was dominant, they help to dramatically reduce the number of severe cases as well as infection and transmission rates. For mRNA vaccines, it was in great part explained by the high level of induced antibodies a few weeks/months after injection and linked to high neutralizing capacity, the ability to prevent viruses to enter and infect target cells. However, decreasing antibody concentration over time and apparition of variants escaping their neutralizing action dramatically reduced the initial vaccine efficacy. As a countermeasure, additional injections were used to re-establish significant antibody population and ensure a long-term neutralizing activity against emerging variants. To infer if this multi-dose strategy fulfils such task, we construct a model of the evolution of the induced antibodies and their neutralizing capacity against different variants. This model helps us to quantify the gain brought by each new injection on both antibody population and their neutralizing ability against all tested variants as well as the dramatic differences between them. We also predict the long-term evolution of neutralizing activity, years after last injection, and thus discuss the longevity of the induced protection by vaccine.

## Introduction

The discovery and the rapid availability of several vaccines against SARS-CoV-2 virus has been a turning point in the combat against Covid-19 [1]. Although their efficacy may vary to some extent, it is undisputable that large scale vaccination campaigns have dramatically reduced both the risk of severe diseases [2–4] and, to a lesser extent, the rates of transmission and disease acquisition [5–7], resulting in millions of saved lives [1, 8, 9].

However vaccine efficacy has been jeopardized by the apparition of various Variants of Concern (VoCs) that partially escape immune protection. A clear decrease in the neutralization capacity has been observed [10, 11] which has translated to a substantial reduction of efficacy against transmission and disease acquisition with Delta and Omicron variants, and, to a lesser extent, to a decrease of efficacy against severe COVID-19 disease [12, 13]. The concern caused by a potential loss of protection against VoC has been further enhanced by the natural waning immunity and the progressive reduction in antibody levels over time [14–16]. This has supported boosting strategies with one or two additional vaccine doses to maintain a high level of protection. However the optimal time to administer boosters, and how these times may vary for different VoC, remains unclear.

To characterize in detail the duration of protection against SARS-CoV-2, it is therefore essential to measure not only total anti-S IgG antibodies over time, as typically done in large observational studies, but also how this translates in terms of neutralization capacity. The latter requires intensive *in vitro* measurements, but it provides a much more accurate description of the level of protection present in the sera of COVID-19 vaccine recipients [17, 18]. Then, a detailed characterization of the immunological or virological factors modulating the duration of protection can be obtained by using mathematical models of immune marker dynamics [19].

Here we propose to use for the first time a mathematical model to analyze the joint kinetics of anti-S IgG antibodies and neutralization capacity after repeated vaccine injections against the main VoCs. For that purpose we relied on data from a cohort of Bnt162b2 vaccine recipients, in which both antibody kinetics and neutralizing activity were measured longitudinally [11, 20, 21]. We built on previous models of antibody kinetics [22, 23] to develop a novel approach to quantify the kinetics of neutralization
capacity, and we use this model to characterize the effects of repeated vaccine administrations on the neutralizing activity. We finally use the model to discuss the duration of protection conferred by Bnt162b2 against VoCs.

# Materials and methods

## Data

### Population Study

Data originate from a cohort of N=29 subjects who received up to three injections of Bnt162b2 (ClinicalTrials.gov:NCT04750720 and ClinicalTrials.gov:NCT05315583). In brief SARS-CoV-2 naive patients were recruited in Orléans, France between August 27, 2020 and May 24, 2022. Individuals were followed for up to 483 days after their first vaccine injections (see more details on the data in [11, 20, 21]). Two patients without longitudinal follow-up and 1 immunocompromised individual were not included in our analysis. In total, N = 26 individuals were analyzed (see Table 1). Briefly, all subjects received at least 2 doses, administered on average 27 days after the first injection. N = 22 subjects received a third injection, administered on average 269 days after the first injection. During the follow-up N = 12 had a positive PCR, and only data prior to infection were analyzed, leaving an average follow-up of 11 visits and a median follow-up time of 362 days.

<table><caption>**Table 1. Characteristics of the analysed population**</caption><thead><tr><th rowspan="2">Characteristics</th><th>Median</th><th colspan="2">Median Time of vaccination</th></tr><tr><th>[Min; Max] or n (%)</th><th>[Min; Max]</th><th></th></tr></thead><tbody><tr><td>Men</td><td>14 (54%)</td><td></td><td></td></tr><tr><td>Age</td><td>59 [33; 95]</td><td></td><td></td></tr><tr><td>Follow-up duration after first-dose (days)</td><td>368 [168; 483]</td><td></td><td></td></tr><tr><td>Number of follow-up visits</td><td>14 [2; 18]</td><td></td><td></td></tr><tr><td>Number of vaccination doses</td><td></td><td></td><td></td></tr><tr><td>1st</td><td>26 (100%)</td><td>-</td><td>-</td></tr><tr><td>2nd</td><td>26 (100%)</td><td>22 [17; 60]</td><td>-</td></tr><tr><td>3nd</td><td>23 (88%)</td><td>243 [175; 385]</td><td>221 [154; 361]</td></tr></tbody></table>

### Ethics statement

This study was approved by the Ethics Committee ILE DE FRANCE IV. The cohort was approved by the national external committee (CPP Ile-de-France- IV IRB No. 00003835). Study participants did not receive any compensation.

### Longitudinal markers of immune response

Two types of measurements were available at each visit: 1) anti-spike binding IgGs, measured in BAU/mL) neutralization titers of sera provided in ED50, which is the effective dilution required to neutralize 50% of an arbitrary viral load of reference (eg, the higher the ED50 the larger the protection level). Neutralization capacity was assessed against historical strain (D614G), Alpha, Beta, Delta, and Omicron variants (strains BA.1, BA.2, and BA.5).

In brief IgGs markedly increased after each dose, but rapidly declined over time, with a rate that did not substantially differ after the second or the third dose (Fig 1 - top).
In contrast the kinetics of neutralizing activity was much more heterogeneous, and was characterized by large differences against the different VoCs. Further the neutralizing activity was markedly increased after the third dose against all Omicron variants, albeit remaining at much lower levels than against the other VoC (Fig 1 - bottom).

<figure><img src="image_3.png" /><figcaption>**Fig 1.** Top: longitudinal evolution of the binding antibody concentration of anti-S IgG. Bottom: longitudinal evolution of the neutralizing activity against VoCs after the first (left), second (middle) and third (right) vaccination dose. Squares represent median values, and plain horizontal lines represent the minimal and maximal encountered values among subjects. The lower limit of detection (LOD) is equal to 6 BAU/mL for IgG and 30 for ED$_{50}$. Given the limited number of samples available, data were grouped, using a one week sliding window after the first dose, 20 days in the first 100 days following the second or third infection, and 50 days for the other data points.</figcaption></figure>

Mechanistic model of B cells response

Mechanistic model for antibody kinetics

We rely on a simplified and rescaled version of a previously published model in the context of vaccine against Ebola infection [23]. In brief, after each dose, cells transfected with Bnt162b2 generate antigen, noted V, which triggers the constitution of a memory
compartment, noted M, at a rate ρ. This memory compartment is a general one accounting for all cell populations able to differentiate into secreting cells upon antigen presence. These can be activated or memory B-cells either circulating or present in germinal centers. So, M can differentiate into secreting plasma cells, noted ˜S, at a rate μV. These cells then produce antibodies, noted Ab, at a rate θ. V, S, and Ab are degraded at rates δV, δS, and δAb, respectively, leading to the following ODE system:

$$
\begin{align}
\dot{V} &= -\delta_V V \\
\dot{M} &= \rho V - \mu V M \\
\dot{\tilde{S}} &= \mu V M - \delta_S \tilde{S} \\
\dot{A} b &= \theta \tilde{S} - \delta_{Ab} A b. \tag{1}
\end{align}
$$

Assuming that individuals are naive of infection, and noting t₁ the time of first injection, the initial conditions are given by: M(t₁) = ˜S(t₁) = Ab(t₁) = 0.

To model the effect of repeated doses, we consider that V is a function presenting discontinuity at time of first, second and third injection (t₁, t₂, t₃). By denoting k = 1, ..., 3 the dose number, on each interval [tₖ, tₖ+₁], solving previous ODE for V gives us Vₖ(t) = V₀e⁻⁵δV(t-tₖ) where V₀ is the initial antigen concentration, assumed equal from all doses.

Because this model is not identifiable when only Ab are measured, we derived a structurally identifiable approximated model described in Equation 2; see Appendix A in S1 Appendix for a description of this simplification. Briefly, it consists of rescaling the model for S = (μV₀M̄₁)⁻¹ ˜S and assuming that M can be replaced by its steady-state value M̄ₖ if equilibrium is reached quickly after each injection:

$$
\begin{align}
\dot{S} &= f_{\overline{M_k}} e^{-\delta_V (t-\mathbf{t}_k)} - \delta_S S \\
\dot{A} b &= \vartheta S - \delta_{Ab} A b \tag{2}
\end{align}
$$

where $f_{\overline{M_k}} = \overline{\frac{M_k}{M_1}}$ is the fold-change for steady-state memory compartment after $k^{th}$ injection compared to the first one (by definition $f_{\overline{M_1}} = 1$). Of note, [23] initial model contains two populations of secreting cells S and L, differing by their life expectancy. In our case, preliminary statistical analysis conclude that there was no statistical differences between model adjustments when accounting for S and L or S only (results no shown). This allows us to reduce the number of unknown parameters. This is crucial for parameters related to cell kinetics known to be very different for newly developed mRNA vaccines comparing to viral vector ones and for which no values have been previously inferred. Thus, the retained model (2) is complex enough to account for the effect of multiple injections on antibody concentration evolution while avoiding identifiability issues. We define $\eta^{ODE} = (f_{\overline{M_2}}, f_{\overline{M_3}}, \delta_V, \delta_S, \vartheta, \delta_{Ab})$ the vector of model parameters defining the dynamics of the system.

Functional model for neutralizing activity

After modeling antibody concentration evolution in the previous section, we aim to model their neutralizing activity. This means in our case proposing a model describing the evolution of EDν₅₀ with respect to Ab. We consider the following linear model:

$ED_{50}^{\nu}(t) = F(\nu, t)Ab(t).$

The function F(ν, t) represents the relationship between the concentration of binding
antibodies in BAU/mL and its neutralization capacity against the VoC ν. It is
variant-specific and time-varying, let us first derive its expression for the strain D614G
before moving to any arbitrary VoCs. After t₁, we assume a proportional relationship
between Ab and neutralizing activity against D614G i.e $F(D614G, t) = \gamma$ (equivalently $ED_{50}^{D614G}(t) = \gamma Ab(t)$). After additional injections, we assume there is a neutralization gain quantified by the fold-change $f_2$ after $t_2$ and $f_3$ after $t_3$ i.e. $F(D614G, t) = \gamma f_2$ when $t \in [t_2; t_3]$ and $F(D614G, t) = \gamma f_3$ for $t \ge t_3$. Now, we account for VoCs specific neutralizing activity by modifying baseline value $\gamma$ by the fold-changes $f_\nu$ such that $F(\nu, t) = F(D614G, t)f_\nu = \gamma f_\nu$ when $t \in [t_1; t_2]$ and $F(\nu, t) = F(D614G, t)f_\nu = \gamma f_\nu f_2$ for $t \in [t_2; t_3]$. We assume that the relative gain brought by third injection can be also VoC-specific. That is why we introduce the fold-changes $g_\nu$ to quantify this gain i.e. $F(\nu, t) = F(D614G, t)f_\nu g_\nu = \gamma f_\nu f_3 g_\nu$ for $t \ge t_3$. This piece-wise constant function can be then expressed in a general form:

$$F(\nu, t) = \gamma f_\nu (1_{t<t_2} + f_2 1_{t \in [t_2; t_3]} + f_3 g_\nu 1_{t \ge t_3}).$$

The choice of this model is the result of exploration based on the minimization of an adjustment criteria. In particular, the current model only quantifies the effect of the repetition of injections on affinity enhancement. Other factors can play a role as the elapsed time since antigen presentation, for example to account for the progressive Memory B-cells repertoire expansion [24, 25]. An alternative neutralization model only considering the time factor has been developed. This supplemental analysis is described in Appendix B in S1 Appendix but lead to a less accurate model (in terms of AIC). A general model accounting for both factors, the number of injections and the elapsed time, has been also tested leading to non-significant improvements over the retained model and at the expense of identifiability problems (results not shown). We also investigate the possibility of a variant-specific fold-change after second injection. This was discarded due to practical identifiability issues.

We define $\eta = (\eta^{ODE}, \gamma, f_\nu, f_2, f_3, g_\nu)$ the vector of model parameters that have to be estimated from the observed data. Description of the model parameters can be found in Table 2.
<table><caption>Table 2. Model parameters and estimation. Fcn:fold-change in neutralisation</caption><thead><tr><th>Parameter</th><th>Description</th><th>Unit</th><th>Fixed Effect [IC95%]</th><th>Random Effect [IC95%]</th></tr></thead><tbody><tr><td>f<sub>M2</sub></td><td>Fold change for M equilibrium after second injection</td><td>dimensionless</td><td>4.5 [2.9; 7.0]</td><td>0.8 [0.5; 1.1]</td></tr><tr><td>f<sub>M3</sub></td><td>Fold change for M equilibrium after third injection</td><td>dimensionless</td><td>12.4 [11.0; 14.0]</td><td></td></tr><tr><td>θ</td><td>Initial acceleration for Ab production</td><td>[A].days<sup>-2</sup></td><td>18.7 [14.8; 23.8]</td><td>0.5 [0.3; 0.6]</td></tr><tr><td>γ</td><td>Proportion of neutralisation provided by first vaccination</td><td>[V].[A]<sup>-1</sup></td><td>0.3 [0.2; 0.4]</td><td>0.8 [0.6; 1.0]</td></tr><tr><td>f<sub>Alpha</sub></td><td>Fcn for variant Alpha</td><td>unitless</td><td>1.3 [1.0; 1.7]</td><td></td></tr><tr><td>f<sub>Beta</sub></td><td>Fcn for variant Beta</td><td>unitless</td><td>0.2 [0.1; 0.3]</td><td></td></tr><tr><td>f<sub>Delta</sub></td><td>Fcn for variant Delta</td><td>unitless</td><td>0.3 [0.2; 0.4]</td><td></td></tr><tr><td>f<sub>BA.1</sub></td><td>Fcn for variant BA.1</td><td>unitless</td><td>0.070 [0.003; 0.015]</td><td></td></tr><tr><td>f<sub>BA.2</sub></td><td>Fcn for variant BA.2</td><td>unitless</td><td>0.014 [0.007; 0.030]</td><td></td></tr><tr><td>f<sub>BA.5</sub></td><td>Fcn for variant BA.5</td><td>unitless</td><td>0.016 [0.010; 0.023]</td><td></td></tr><tr><td>f<sub>2</sub></td><td>Fcn for second vs. first injection</td><td>unitless</td><td>8.3 [5.6; 12.3]</td><td></td></tr><tr><td>f<sub>3</sub></td><td>Fcn for third vs. first injection in original strains D614G</td><td>unitless</td><td>16.6 [9.7; 28.4]</td><td></td></tr><tr><td>g<sub>Alpha</sub></td><td>Fcn for third vs. first injection in variant Alpha</td><td>unitless</td><td>7.2 [2.1; 24.8]</td><td></td></tr><tr><td>g<sub>Beta</sub></td><td>Fcn for third vs. first injection in variant Beta</td><td>unitless</td><td>2.8 [0.9; 8.9]</td><td></td></tr><tr><td>g<sub>Delta</sub></td><td>Fcn for third vs. first injection in variant Delta</td><td>unitless</td><td>1.3 [0.4; 3.3]</td><td></td></tr><tr><td>g<sub>BA.1</sub></td><td>Fcn for third vs. first injection in variant BA.1</td><td>unitless</td><td>13.7 [4.2; 45.8]</td><td></td></tr><tr><td>g<sub>BA.2</sub></td><td>Fcn for third vs. first injection in variant BA.2</td><td>unitless</td><td>5.3 [1.5; 22.0]</td><td></td></tr><tr><td>g<sub>BA.5</sub></td><td>Fcn for third vs. first injection in variant BA.5</td><td>unitless</td><td>1.6 [1.0; 5.3]</td><td></td></tr><tr><td>σ<sub>BAU</sub></td><td>Measurement error for Y<sup>BAU</sup> = log<sub>10</sub>(Ab) + e<sup>BAU</sup></td><td></td><td>0.23 [0.21; 0.25]</td><td></td></tr><tr><td>σ<sub>D614G</sub></td><td>Measurement error for Y<sup>D614G</sup> = log<sub>10</sub>(ED<sup>D614G</sup>) + e<sup>D614G</sup></td><td></td><td>0.48 [0.44; 0.52]</td><td></td></tr><tr><td>σ<sub>Alpha</sub></td><td>Measurement error for Y<sup>Alpha</sup> = log<sub>10</sub>(ED<sup>Alpha</sup>) + e<sup>Alpha</sup></td><td></td><td>0.59 [0.58; 0.60]</td><td></td></tr><tr><td>σ<sub>Beta</sub></td><td>Measurement error for Y<sup>Beta</sup> = log<sub>10</sub>(ED<sup>Beta</sup>) + e<sup>Beta</sup></td><td></td><td>0.45 [0.40; 0.50]</td><td></td></tr><tr><td>σ<sub>Delta</sub></td><td>Measurement error for Y<sup>Delta</sup> = log<sub>10</sub>(ED<sup>Delta</sup>) + e<sup>Delta</sup></td><td></td><td>0.44 [0.39; 0.49]</td><td></td></tr><tr><td>σ<sub>BA.1</sub></td><td>Measurement error for Y<sup>BA.1</sup> = log<sub>10</sub>(ED<sup>BA.1</sup>) + e<sup>BA.1</sup></td><td></td><td>0.44 [0.38; 0.50]</td><td></td></tr><tr><td>σ<sub>BA.2</sub></td><td>Measurement error for Y<sup>BA.2</sup> = log<sub>10</sub>(ED<sup>BA.2</sup>) + e<sup>BA.2</sup></td><td></td><td>0.45 [0.36; 0.54]</td><td></td></tr><tr><td>σ<sub>BA.5</sub></td><td>Measurement error for Y<sup>BA.5</sup> = log<sub>10</sub>(ED<sup>BA.5</sup>) + e<sup>BA.5</sup></td><td></td><td>0.36 [0.30; 0.42]</td><td></td></tr><tr><td>δ<sub>Ab</sub></td><td>Antibody degradation rate</td><td>days<sup>-1</sup></td><td>0.03 (fixed)</td><td></td></tr><tr><td>δ<sub>V</sub></td><td>Induced vaccine antigen declining rate</td><td>days<sup>-1</sup></td><td>2.7 (fixed)</td><td></td></tr><tr><td>δ<sub>S</sub></td><td>Death rate of S cells</td><td>days<sup>-1</sup></td><td>0.01 (fixed)</td><td></td></tr></tbody></table>
## Observation model

The structural model used to describe the log-transformed concentration of binding antibodies in BAU/mL for the $i^{th}$ individual ($i = 1, ..., N$) at the $j^{th}$ time point ($j = 1, ..., n_i$) is:

$$Y_{ij}^{BAU} = \log_{10}(Ab(\eta_i, t_{ij})) + e_{ij}^{BAU},$$

where $e_{ij}^{BAU}$ is the residual additive error which follow a normal distribution of mean zero and constant standard deviation $\sigma_{BAU}$. The vector $\eta_i$ is the specific value for individual $i$ of vector $\eta$.

We also consider a log-transformation of $ED_{50}^{\nu}$ raw measurements for the variant in the list {D614G, Alpha, Beta, Delta, BA.1, BA.2, BA.5}. For the $i^{th}$ individual at the $j^{th}$ time point, we have:

$$Y_{ij}^{\nu} = \log_{10}(ED_{50}^{\nu}(\eta_i, t_{ij})) + e_{ij}^{\nu},$$

where $e_{ij}^{\nu}$ is the residual additive error for variant $\nu$ which follow a normal distribution of mean zero and constant standard deviation $\sigma_{\nu}$.

## Statistical model for parameters over time and injections

**Fixed parameters.** Because not all parameters can be estimated when only concentration of binding antibodies and antibody neutralizing activity are measured, some parameters had to be fixed based on the current literature. Here, we just fix one parameter, the clearance of antibodies found to be stable across studies and infections [26] with a half life of 21 days. Further, the model predictions were found largely insensitive to the choice of the degradation rate of V and S. Using a profiled likelihood approach [27], we fixed their half-life to 0.25 and 51 days, respectively.

**Inter-individual variability.** In the vector $\eta$, some parameters have to be individual-specific to account for inter-individual variability. It is the case for $\psi_i = (\vartheta, f_{\overline{M}_2}, \gamma)$. We suppose it follows a log-normal distribution such that:

$$\psi_i = \psi_0 \exp(u_i),$$

where $\psi_0$ is the fixed effect and average mean value in the population. The vector $u_i$ is individual random effects, which follow a normal distribution of mean zero and standard deviation $\Omega$, and account for heterogeneity across individual. We assume that other parameters in vector $\eta$ except error measurements are also estimated in log-transformation and are common to all individuals in the population. Altogether, the vector of parameters to estimate is given by $\theta = (\eta, \Omega, \sigma_{BAU}, \sigma_{\nu})$.

**Estimation procedure** Parameters were estimated (and named $\hat{\theta}$ in the following) with the SAEM algorithm implemented in MONOLIX software version 2022R1 [28] allowing to handle left censored data [29]. Likelihood was estimated using the importance sampling method and standard error were obtained by asymptotic approximation and inversion of the Fisher Information Matrix. Graphical and statistical analyses were performed using R version 3.4.3.

## Simulation of long-term humoral response

Next, we used the model to predict the long-term evolution of Ab and $ED_{50}^{\nu}$ over time. To account for uncertainty in our predictions, we used a Monte-Carlo sampling method, where $K = 1000$ replicates of parameters values $\theta^{(k)}$ were sampled in the posterior distribution of the parameter estimates to derive 95% prediction intervals (PI) of the predicted trajectories.
Finally we used these predictions to calculate the time to reach a given threshold value. To take into account between-subjects variability, we added a second layer to our Monte-Carlo sampling method and we sampled $N = 100$ replicates in the population parameter distribution. We used these predictions to derive the probability of having a concentration of binding antibodies greater than given thresholds, in particular higher than 264 BAU/mL, which corresponds to the standard threshold of protection defined by Feng et al. [30] and adopted by WHO. The level of neutralizing activity has been identified as a correlate of protection for vaccine efficacy against the historical strain [31, 32]. However, to date, no threshold for $ED_{50}^v$ value ensuring protection has been isolated for D614G, let alone for the new VoCs. So, for a range of threshold values, we calculated the probability that the neutralizing activity against each VoC would be higher than these values over time, especially if this activity was still detectable at a given time. In this way, we can compare the longevity of neutralizing activity between VoCs even in the absence of a clear threshold of protection for each of them.

# Results

## Mechanistic model for humoral response

We first aimed to investigate whether there is a proportional relationship between the evolution of concentration of binding antibodies and its neutralization capacity. Fig 2 displays the observed relationship from data between antibody concentration and $ED_{50}^v$ for each VoC after each injection. First, we notice that these ratios are different for the variants. Then, using a Student's t test, we compared the evolution of these ratios with respect to the previous vaccination. In most cases, the ratios improved significantly, indicating an intrinsic gain in neutralisation that cannot be explained by the variation in antibody concentration alone, justifying the need to quantify this phenomenon precisely.

Estimation of model parameters can be found in Table 2. This estimation indicates that multiple injections both increase antibody concentration and intrinsic affinity per constant antibody unit. Regarding antibody concentration, estimation of mechanistic parameters indicates a significant increase in the size of the memory compartment. It increased by $f_{M_2} = 4.5$ (95% Confidence interval CI [2.9; 7.0]) after the second injection and by $f_{M_3} = 12.4$ (95% CI [11.0; 14.0]) after the third injection compared to the first one.

Regarding neutralization per constant antibody concentration unit, we found that there are two main influencing factors: the repetition of the injections and the VoC. Regarding repeated injections effect for the original strain D614G, the second dose increases neutralization by a factor $f_2 = 8.3$ (95% CI [5.6; 12.3]) and the third one by $f_3 = 17$ (95% CI [9.7; 28.4]) compared to the first injection. Now regarding the neutralization capacities for emerging VoCs, they are significantly decreased compared to the original strain, with the exception of Alpha, where there is no significant change in neutralization compared to D614G. It ranges from a reduction of 70% (95% CI [60%; 80%]) for Delta to a dramatic reduction of 98.4% (95% CI [97.7%; 99%]) for BA.5. Still, we find that the sequential injection strategy confers a gain in long-term neutralizing capacities for all VoCs. The second injection increases neutralization against all VoCs by the same factor $f_2 = 8.3$ (same as D614G). The third injection increases neutralization in a VoC-specific manner, given by $f_3 g_v$. It ranges from an increase in fold change of 21 (95% CI [4; 94]) for Delta to 230 (95% CI [40; 1200]) for BA.1 times higher for the third injection than for the first injection. For comparison with D614G, the neutralization is $f_3/f_2 = 2.0$ (95% CI [1.7; 2.3]) times higher for the third injection compared to second injection. Transitively, the fold change is $f_3/f_2 g_{Delta} = 2.6$ (95% CI [0.7; 7.6]) for Delta to $f_3/f_2 g_{BA.1} = 27.4$ (95% CI [7.2; 105.2]) for BA.1 times higher for
<figure><img src="image_4.png" /><figcaption>**Fig 2.** Evolution of the predicted ratio $ED_{50}/BAU$ for each VoC after successive vaccine doses. Each circle represents a ratio $ED_{50}/BAU$ computed when both measurements for $ED_{50}$ and BAU where available for a given patient at a given observation time. Most of patients contribute several times due to the repeated measurements made over time after each dose. All predictions below the limit of detection for $ED_{50}$ were removed to avoid overoptimistic $ED_{50}/BAU$ ratio when replacing $ED_{50}$ values by detection threshold. This explains why very few values are available for Beta and Delta and none for Omicron strains for one dose case. Comparison between vaccine dose were done using a Student test, "**" indicate p-values lower than 0.01, "*" lower than 0.1 and ">0.1" higher than 0.1.</figcaption></figure>

the third injection than for the second injection.

Examples of fitted trajectories are given for four randomly selected patients in Fig 3. We observe a very good adequation with most of the observations lying in the 95% prediction intervals.

These trajectories are exemplified in Fig 4, that shows the mean markers trajectories
<figure><img src="image_5.png" /><figcaption>**Fig 3.** Individual fits for four representative individuals. The solid line is the subject-specific prediction and the shaded area is the 95% prediction interval. The plain dots and crosses represent the observed and censored data, respectively.</figcaption></figure>

for an average individual (i.e random effects $u_i$ set to 0). As expected, the level of the response is higher after a repeated number of injections for both binding antibody concentration and neutralization for all variants. Interestingly, the neutralization curves for BA.1, BA.2, and BA.5 are significantly lower than for the other variants, with no overlap in prediction intervals. The first and second doses elicit a neutralization response for Omicron (BA.1, BA.2, and BA.5) that remains below the detection limit in most individuals (which is consistent with the observed data) but is dramatically enhanced by the third injection. Regarding the $ED_{50}/BAU$ ratio (See S1 Fig), we find that for the same concentration of binding antibodies, neutralization is significantly increased after each new injections for all variants and is significantly different for Alpha, D614G, {Delta, Beta} and {BA.1, BA.2, BA.5} variants.

## Long-term predictions

As already shown in Fig 4, we can use the estimated models to predict the long-term trajectories of markers corresponding to the mean parameter values as well as 95%
<figure><img src="image_6.png" /><figcaption>**Fig 4.** Left: Predicted evolution of binding antibody concentration. The horizontal line corresponds to the value of 264 BAU/ml considered as a threshold against symptomatic infection. Right: Predicted kinetics of $ED_{50}$). The horizontal line corresponds to the LOD. In all panels, the shaded area is the 95% prediction interval.</figcaption></figure>

prediction intervals. It allows to derive an estimation of the time needed to reach a certain threshold after a three injections vaccination scheme with first vaccination at time $t_1 = 0$, second injection at time $t_2 = 27$ days and third injection at time $t_3 = 269$ days, corresponding to the mean observed time of injection in our cohort. Binding antibodies concentration is below 264 BAU/mL 158 (95/% PI [143; 174]) days after third vaccination. Neutralisation reaches undetectable levels between 184 days (95/% PI [155; 215]) for BA.5 to 606 (95/% PI [555; 638]) for Alpha after the third dose.

Fig 5 (top) displays the probability of having antibody concentration higher than the protection threshold established by [30] of 264 BAU/mL each days after the last injection in the counterfactual scenario where subjects only received one, two or three doses. The same is done for neutralizing activity again the VoCs (bottom, left: one, middle: two or right: three). It is possible to see the drastic effect of repeated injections on the levels reached by both binding antibodies concentration and neutralization for all variants. Strikingly, the full response duration is similar in length for the binding antibodies concentration after two or three injections. However, whereas response higher than 264 BAU/mL is reached in 100% (95% PI [99%; 100%]) of the population after three injections, it is only reached in 75% (95% PI [60%; 90%]) of the population after two doses, and not even reached in 2% (95% PI [0%; 3%]) of the population after the first injection. Table 3 provides the time needed for a proportion of a vaccinated population to return under a certain threshold. It explores multiple thresholds (100 BAU/mL, 264 BAU/mL, and 1000 BAU/mL for antibodies concentrations; and undetectability, 100 and 1000 for neutralization) that could be investigated when and if a clear level of correlate of protection is found. For all markers, there is a systematically and significantly higher duration of humoral activity after three compared to two injections. After three injections, duration of neutralization against Omicron variants (BA.1,BA.2 and BA.5) is significantly lower than for other variants for all thresholds.
<figure><img src="image_7.png" /><figcaption>**Fig 5.** Top: Predicted probability of having predicted antibody concentration (anti-S IgG) greater than 264 BAU/mL. Bottom: Probability of having detectable neutralizing activity against VoCs after the first (left), second (middle) or third (right) vaccination dose. Simulations were performed assuming that the second and third vaccination doses occurred at day 27 and 269, respectively.</figcaption></figure>

Discussion

We proposed here a modeling framework to characterize the kinetics of antibodies to successive doses of Bnt162b2 vaccine. The originality of our approach is that we relied on both the kinetics of anti-S IgG binding antibodies and their neutralization against the major VoCs that have emerged since 2021. Our model quantifies the benefit of successive injections and can be used to predict the optimal duration against each VoC. The model shows a large and significant benefit of the third dose against all VoCs in terms of neutralizing activity. However, both the maximum level achieved and the rate of decline could vary greatly between VoCs. Accordingly, the mean duration of detectable neutralizing activity after the third dose of vaccine was 20, 12.5, 9, 8.5 and 6 months for Alpha, Delta, and Omicron BA.1, BA.2 and BA.5 respectively. Our results also highlight the wide variability in patient response, with at least 5% of patients with undetectable neutralizing activity against Omicron BA.5 only 3 months after the third
<table><caption>Table 3. Predicted distribution for the duration of anti-S IgG and neutralization activity above different threshold levels [95% prediction interval].</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Population quantiles</th><th colspan="3">Time to anti-S IgG</th></tr><tr><th>&lt;100 BAU/mL</th><th>&lt;264 BAU/mL</th><th>&lt;1000 BAU/mL</th></tr></thead><tbody><tr><td rowspan="3">IgG</td><td>95%</td><td>181 [157; 197]</td><td>103 [86; 126]</td><td>0 [0; 0]</td></tr><tr><td>50%</td><td>238 [220; 255]</td><td>162 [149; 183]</td><td>60 [0; 80]</td></tr><tr><td>5%</td><td>300 [276; 315]</td><td>233 [228; 243]</td><td>120 [104; 145]</td></tr><tr><th rowspan="2">Variant</th><th rowspan="2">Population quantiles</th><th colspan="3">Time to ED50</th></tr><tr><th>Undetectable</th><th>&lt;100</th><th>&lt;1000</th></tr><tr><td rowspan="3">D614G</td><td>95%</td><td>360 [304; 375]</td><td>254 [213; 286]</td><td>100 [0; 117]</td></tr><tr><td>50%</td><td>441 [413; 481]</td><td>351 [326; 392]</td><td>183 [152; 223]</td></tr><tr><td>5%</td><td>541 [517; 592]</td><td>454 [429; 508]</td><td>283 [257; 338]</td></tr><tr><td rowspan="3">Alpha</td><td>95%</td><td>522 [449; 551]</td><td>416 [353; 463]</td><td>264 [190; 292]</td></tr><tr><td>50%</td><td>603 [557; 658]</td><td>512 [467; 566]</td><td>345 [298; 399]</td></tr><tr><td>5%</td><td>703 [666; 768]</td><td>615 [572; 679]</td><td>445 [401; 510]</td></tr><tr><td rowspan="3">Beta</td><td>95%</td><td>282 [212; 308]</td><td>176 [126; 218]</td><td>0 [0; 0]</td></tr><tr><td>50%</td><td>363 [327; 411]</td><td>272 [231; 323]</td><td>103 [60; 154]</td></tr><tr><td>5%</td><td>463 [433; 526]</td><td>376 [338; 432]</td><td>205 [174; 264]</td></tr><tr><td rowspan="3">Delta</td><td>95%</td><td>287 [230; 303]</td><td>181 [140; 211]</td><td>0 [0; 0]</td></tr><tr><td>50%</td><td>368 [340; 409]</td><td>277 [250; 317]</td><td>108 [80; 145]</td></tr><tr><td>5%</td><td>468 [446; 519]</td><td>381 [356; 431]</td><td>210 [187; 258]</td></tr><tr><td rowspan="3">BA.1</td><td>95%</td><td>185 [128; 200]</td><td>73 [0; 110]</td><td>0 [0; 0]</td></tr><tr><td>50%</td><td>266 [236; 304]</td><td>175 [148; 216]</td><td>0 [0; 0]</td></tr><tr><td>5%</td><td>365 [343; 418]</td><td>278 [254; 329]</td><td>106 [81; 161]</td></tr><tr><td rowspan="3">BA.2</td><td>95%</td><td>169 [107; 190]</td><td>50 [0; 94]</td><td>0 [0; 0]</td></tr><tr><td>50%</td><td>251 [218; 294]</td><td>160 [129; 200]</td><td>0 [0; 0]</td></tr><tr><td>5%</td><td>350 [323; 404]</td><td>263 [236; 313]</td><td>89 [49; 149]</td></tr><tr><td rowspan="3">BA.5</td><td>95%</td><td>98 [0; 111]</td><td>0 [0; 0]</td><td>0 [0; 0]</td></tr><tr><td>50%</td><td>181 [154; 218]</td><td>87 [54; 126]</td><td>0 [0; 0]</td></tr><tr><td>5%</td><td>281 [258; 332]</td><td>194 [169; 240]</td><td>0 [0; 63]</td></tr></tbody></table>

injection.

These results were obtained based on a number of hypotheses, which we summarize below. First, the model of antibody concentration dynamics remains simplified, with the memory compartment simply represented by a piecewise constant function over successive doses. In addition, the model assumes only one type of secreting cell population and thus overlooks the complexity of the B-cell response mechanism. A piecewise constant function was also used to model the effects of successive doses on neutralization activity and did not account for other potential phenomena, such as time since last injection. Finally, the model assumed that the second dose would result in a similar change in protection for all variants. In the future, application of such approaches to larger populations of individuals may allow some of these hypotheses to be relaxed without compromising the identifiability of the parameters.

We restrict our analysis to vaccination induced humoral response and discarded patient data after breakthrough infection. To this date, 11 out of 26 followed subjects were infected with Omicron and this proportion is likely to increase as it is the case in the global vaccinated population. Thus, it is of great interest to model the hybrid protection induced by vaccination followed by natural infection. Still, it requires to deeply modify our model to integrate two different antibody populations, one coming
from vaccination and targeting the historical strain and the other one targeting the Omicron variant. That is why this analysis is left to future works.

One of the main advantages of the model is its flexibility to easily incorporate information on new VoCs and to use the strength of information obtained on other viral variants to update the model as data become available. In fact, we continuously updated the model to include successive Omicron variants. Interestingly, despite the small number of samples available, a high degree of precision was achieved for all variants. For example, although patients had on average only two data points with detectable neutralization against the BA.2 variant, this was sufficient to achieve a good precision of the model parameters (Table 2). Moreover, by using all available data (eg, by analyzing anti-S IgG and neutralizing activity of all patients simultaneously), the model reveals some signals of kinetics that were not visible when analyzing the individual markers separately. For example, we identified different slopes of antibody decline that directly affect the prediction of protection duration. Using the same data set and a simpler single-slope model, the time to undetectable neutralizing levels after the third dose of vaccine for D416G was estimated to be 11.5 months [21], which is much shorter than our estimate of 15 months (derived from Fig 4). In their approach, Planas et al. [21] chose to adjust the anti-S IgG and $ED_{50}^\nu$ decline separately for the different VoCs without considering causal relationships between them. On the contrary, our model assumes an influence of antibody concentration on the development of $ED_{50}^\nu$. Similar results are shown for the Delta variant (11.5 vs. 12.5 months, respectively), demonstrating the importance of a model-based approach to predict neutralizing activity in the long term. Predictions for the Omicron strains were similar for BA.1 and BA.2, respectively 9 and 8.5 months, and show a reduction for BA.5 strains with 6 months. Interestingly, our results also suggest a longer duration of detectable neutralizing activity than what has been directly extrapolated from other observational studies [33–35], although this difference may also be due in part to the different experimental procedure used to measure neutralization.

The large and VoC-dependent variability in neutralization values argues for the use of individualized approaches to identify patients most at risk. Although it should be acknowledged that such an approach is hampered by the lack of an established neutralizing activity threshold as correlates of protection, its level was found to be associated with the risk of breakthrough infection. In a cohort of elderly nursing home residents, none of those the individuals with ED50 above 2136 had Omicron BA.1 breakthrough infection [36]. A model-based study found that a threshold of 1000 dramatically reduced peak viral load, suggesting that such a threshold may be a good indicator of protection against infection (unpublished results). Interestingly, our results show that neutralizing levels for all Omicron variants remain largely below this value (Fig 5), consistent with the current understanding that BNT162b2 is poorly effective against disease acquisition in the Omicron era [36, 37]. Fortunately, the vaccine has shown high efficacy against severe disease to date [38, 39].

To date, the use of a fourth dose of vaccine to increase efficacy in France has been limited to high-risk patients who were not represented in this cohort. Nevertheless, we used the model to predict the neutralization levels that could be achieved after a fourth vaccine dose. Under the conservative assumption, yet consistent with available observational study [40], that this injection does not increase affinity or maturation parameters, our model predicts a similar duration of detectable neutralization as after the third dose, ranging from 183 to 268 days for the Omicron variants. Assuming that the fourth dose allows a similar increase in maturation and affinity as after the third dose, the model predicts that the duration of detectable neutralization could be much longer, ranging from 596 to 680 days for Omicron variants (see this supplementary analysis in Appendix C in S1 Appendix).
## Supporting information

<figure><img src="image_8.png" /><figcaption>Fig 6. Evolution of $ED_{50}^\nu$:BAU ratio.</figcaption></figure>

<caption>S1 Appendix. Appendixes for "Modeling the evolution of the neutralizing antibody response against SARS-CoV-2 variants after several administrations of Bnt162b2"</caption>

## Acknowledgments

We thank Isabelle Staropoli, Florence Guivel-Benhassine, Françoise Porrot and all the members of the Virus and Immunity Unit for discussion and help, as well as Fabienne Peira, Vanessa Legros, Barbara De Dieuleveult, Aurelie Theillay, Sandra Pallay and Daniela Pires Roteia (CHR Orléans) for their help with the cohorts. Part of the experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr). We thank Simulations Plus, Lixoft division for the free academic use of the MonolixSuite.

## References

1. Watson OJ, Barnsley G, Toor J, Hogan AB, Winskill P, Ghani AC. Global impact of the first year of COVID-19 vaccination: a mathematical modelling study. The Lancet Infectious Diseases. 2022;22(9):1293-1302.

2. Thomas SJ, Moreira Jr ED, Kitchin N, Absalon J, Gurtman A, Lockhart S, et al. Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine through 6 months. New England Journal of Medicine. 2021;385(19):1761-1773.

3. Polack FP, Thomas SJ, Kitchin N, Absalon J, Gurtman A, Lockhart S, et al. Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine. New England journal of medicine. 2020;.
4. Baden LR, El Sahly HM, Essink B, Kotloff K, Frey S, Novak R, et al. Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine. New England journal of medicine. 2020;

5. Harris RJ, Hall JA, Zaidi A, Andrews NJ, Dunbar JK, Dabrera G. Effect of vaccination on household transmission of SARS-CoV-2 in England. New England Journal of Medicine. 2021;385(8):759-760.

6. Pritchard E, Matthews PC, Stoesser N, Eyre DW, Gethings O, Vihta KD, et al. Impact of vaccination on new SARS-CoV-2 infections in the United Kingdom. Nature medicine. 2021;27(8):1370-1378.

7. Eyre DW, Taylor D, Purver M, Chapman D, Fowler T, Pouwels KB, et al. Effect of Covid-19 vaccination on transmission of alpha and delta variants. New England Journal of Medicine. 2022;386(8):744-756.

8. Gupta S, Cantor J, Simon KI, Bento AI, Wing C, Whaley CM. Vaccinations Against COVID-19 May Have Averted Up To 140,000 Deaths In The United States: Study examines role of COVID-19 vaccines and deaths averted in the United States. Health Affairs. 2021;40(9):1465-1472.

9. Meslé MM, Brown J, Mook P, Hagan J, Pastore R, Bundle N, et al. Estimated number of deaths directly averted in people 60 years and older as a result of COVID-19 vaccination in the WHO European Region, December 2020 to November 2021. Eurosurveillance. 2021;26(47):2101021.

10. Cele S, Jackson L, Khoury DS, Khan K, Moyo-Gwete T, Tegally H, et al. Omicron extensively but incompletely escapes Pfizer BNT162b2 neutralization. Nature. 2022;602(7898):654-656.

11. Planas D, Veyer D, Baidaliuk A, Staropoli I, Guivel-Benhassine F, Rajah MM, et al. Reduced sensitivity of SARS-CoV-2 variant Delta to antibody neutralization. Nature. 2021;596(7871):276-280.

12. Bernal JL, Andrews N, Gower C, Gallagher E, Simmons R, Thelwall S, et al. Effectiveness of Covid-19 vaccines against the B. 1.617. 2 (Delta) variant. New England Journal of Medicine. 2021;

13. Andrews N, Stowe J, Kirsebom F, Toffa S, Rickeard T, Gallagher E, et al. Covid-19 vaccine effectiveness against the Omicron (B. 1.1. 529) variant. New England Journal of Medicine. 2022;386(16):1532-1546.

14. Mizrahi B, Lotan R, Kalkstein N, Peretz A, Perez G, Ben-Tov A, et al. Correlation of SARS-CoV-2-breakthrough infections to time-from-vaccine. Nature communications. 2021;12(1):1-5.

15. Goldberg Y, Mandel M, Bar-On YM, Bodenheimer O, Freedman L, Haas EJ, et al. Waning immunity after the BNT162b2 vaccine in Israel. New England Journal of Medicine. 2021;385(24):e85.

16. Levin EG, Lustig Y, Cohen C, Fluss R, Indenbaum V, Amit S, et al. Waning immune humoral response to BNT162b2 Covid-19 vaccine over 6 months. New England Journal of Medicine. 2021;385(24):e84.

17. Cromer D, Juno JA, Khoury D, Reynaldi A, Wheatley AK, Kent SJ, et al. Prospects for durable immune control of SARS-CoV-2 and prevention of reinfection. Nature Reviews Immunology. 2021;21(6):395-404.
18. Zhu F, Althaus T, Tan CW, Costantini A, Chia WN, Chau NVV, et al. WHO international standard for SARS-CoV-2 antibodies to determine markers of protection. The Lancet Microbe. 2022;3(2):e81-e82.

19. Korosec CS, Farhang-Sardroodi S, Dick DW, Gholami S, Ghaemi MS, Moyles IR, et al. Long-term predictions of humoral immunity after two doses of BNT162b2 and mRNA-1273 vaccines based on dosage, age and sex. medRxiv. 2021;doi:10.1101/2021.10.13.21264957.

20. Planas D, Bruel T, Grzelak L, Guivel-Benhassine F, Staropoli I, Porot F, et al. Sensitivity of infectious SARS-CoV-2 B. 1.1. 7 and B. 1.351 variants to neutralizing antibodies. Nature medicine. 2021;27(5):917-924.

21. Planas D, Staropoli I, Porot F, Guivel-Benhassine F, Handala L, Prot M, et al. Duration of BA. 5 neutralization in sera and nasal swabs from SARS-CoV-2 vaccinated individuals, with or without Omicron breakthrough infection. Med. 2022;.

22. Pasin C, Balleli I, Effelterre TV, Bockstal V, Solforosi L, Prague M, et al. Dynamic of the humoral immune response to a prime-boost Ebola vaccine: quantification and source of variation. Journal of Virology. 2019;93(18):e00579-19.

23. Balelli I, Pasin C, Prague M, Crauste F, Thiébaut R. A model for establishment, maintenance and reactivation of the immune response after two-dose vaccination regimens against Ebola virus. Journal of Theoretical Biology. 2020; p. 110254.

24. Muecksch F, Wang Z, Cho A, Gaebler C, Tanfous TB, DaSilva J, et al. Increased potency and breadth of SARS-CoV-2 neutralizing antibodies after a third mRNA vaccine dose. bioRxiv. 2022;.

25. Wang K, Jia Z, Bao L, Wang L, Cao L, Chi H, et al. Memory B cell repertoire from triple vaccinees against diverse SARS-CoV-2 variants. Nature. 2022;603(7903):919-925.

26. Fatehi F, Bingham RJ, Dykeman EC, Stockley PG, Twarock R. Comparing antiviral strategies against COVID-19 via multiscale within-host modelling. Royal Society Open Science. 2021;8(8):210082.

27. Murphy SA, Van der Vaart AW. On profile likelihood. Journal of the American Statistical Association. 2000;95(450):449-465.

28. Lixoft. Monolix version 2020R1. Antony, France. 2020;

29. Samson A, Lavielle M, Mentré F. Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: Application to HIV dynamics model. Computational Statistics & Data Analysis. 2006;51(3):1562-1574.

30. Feng S, Phillips DJ, White T, Sayal H, Aley PK, Bibi S, et al. Correlates of protection against symptomatic and asymptomatic SARS-CoV-2 infection. Nature medicine. 2021;27(11):2032-2040.

31. Gilbert PB, Montefiori DC, McDermott AB, Fong Y, Benkeser D, Deng W, et al. Immune correlates analysis of the mRNA-1273 COVID-19 vaccine efficacy clinical trial. Science. 2022;375(6576):43-50.

32. Alexandre M, Marlin R, Prague M, Coleon S, Kahlaoui N, Cardinaud S, et al. Modelling the response to vaccine in non-human primates to define SARS-CoV-2 mechanistic correlates of protection. Elife. 2022;11:e75427.
33. Brockman MA, Mwimanzi F, Lapointe HR, Sang Y, Agafitei O, Cheung PK, et al. Reduced magnitude and durability of humoral immune responses to COVID-19 mRNA vaccines among older adults. The Journal of Infectious Diseases. 2022;225(7):1129-1140.

34. Mwimanzi F, Lapointe HR, Cheung PK, Sang Y, Yaseen F, Umviligihozo G, et al. Older adults mount less durable humoral responses to two doses of COVID-19 mRNA vaccine, but strong initial responses to a third dose. medRxiv. 2022;.

35. Mwimanzi FM, Lapointe HR, Cheung PK, Sang Y, Yaseen F, Kalikawe R, et al. Brief Report: Impact of age and SARS-CoV-2 breakthrough infection on humoral immune responses after three doses of COVID-19 mRNA vaccine. medRxiv. 2022;.

36. Bruel T, Pinaud L, Tondeur L, Planas D, Staropoli I, Porrot F, et al. Neutralising antibody responses to SARS-CoV-2 omicron among elderly nursing home residents following a booster dose of BNT162b2 vaccine: A community-based, prospective, longitudinal cohort study. EClinicalMedicine. 2022;51:101576.

37. Lapointe HR, Mwimanzi F, Cheung PK, Sang Y, Yaseen F, Kalikawe R, et al. Serial infection with SARS-CoV-2 Omicron BA. 1 and BA. 2 following three-dose COVID-19 vaccination. medRxiv. 2022;.

38. Hacisuleyman E, Hale C, Saito Y, Blachere NE, Bergh M, Conlon EG, et al. Vaccine breakthrough infections with SARS-CoV-2 variants. New England Journal of Medicine. 2021;384(23):2212-2218.

39. Kuhlmann C, Mayer CK, Claassen M, Maponga T, Burgers WA, Keeton R, et al. Breakthrough infections with SARS-CoV-2 omicron despite mRNA vaccine booster dose. The Lancet. 2022;399(10325):625-626.

40. Regev-Yochay G, Gonen T, Gilboa M, Mandelboim M, Indenbaum V, Amit S, et al. Efficacy of a fourth dose of COVID-19 mRNA vaccine against omicron. New England Journal of Medicine. 2022;386(14):1377-1380.