ARTICLE **OPEN**

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

Check for updates

# Prediction of long-term humoral response induced by the two-dose heterologous Ad26.ZEBOV, MVA-BN-Filo vaccine against Ebola

Marie Alexandre<sup>1,2</sup>, Mélanie Prague<sup>1,2,8</sup>, Chelsea McLean<sup>3</sup>, Viki Bockstal<sup>3,6</sup>, Macaya Douoguih<sup>1,3</sup>, Rodolphe Thiébaut<sup>1,2,8✉</sup> and for the EBOVAC 1 and EBOVAC 2 Consortia*

The persistence of the long-term immune response induced by the heterologous Ad26.ZEBOV, MVA-BN-Filo two-dose vaccination regimen against Ebola has been investigated in several clinical trials. Longitudinal data on IgG-binding antibody concentrations were analyzed from 487 participants enrolled in six Phase I and Phase II clinical trials conducted by the EBOVAC1 and EBOVAC2 consortia. A model based on ordinary differential equations describing the dynamics of antibodies and short- and long-lived antibody-secreting cells (ASCs) was used to model the humoral response from 7 days after the second vaccination to a follow-up period of 2 years. Using a population-based approach, we first assessed the robustness of the model, which was originally estimated based on Phase I data, against all data. Then we assessed the longevity of the humoral response and identified factors that influence these dynamics. We estimated a half-life of the long-lived ASC of at least 15 years and found an influence of geographic region, sex, and age on the humoral response dynamics, with longer antibody persistence in Europeans and women and higher production of antibodies in younger participants.

npj Vaccines (2023)8:174; https://doi.org/10.1038/s41541-023-00767-y

## INTRODUCTION

The 2014–2016 Ebola virus disease (EBOV) outbreak in West Africa and the current SARS-CoV-2 pandemic have led to accelerated development of vaccines to control the spread of infection and reduce the severity of disease in infected individuals. As a result, effective vaccines were developed and became available quickly after the start of these two epidemics. In the case of Ebola, the recombinant replication-competent vesicular stomatitis viral vectored vaccine (Ervebo) was approved by the FDA in December 2019¹ and used during epidemics in a ring vaccination strategy. The two-dose heterologous strategy, combining immunizations with Ad26.ZEBOV (Zabdeno) and MVA-BN-Filo (Mvabea), were approved by the European Commission in July 2020² under exceptional circumstances for use in children and adults. An important question for those who have already been vaccinated, and for using the vaccines as a preventive strategy to control the occurrence of outbreaks, is the duration of protection conferred by vaccination.

In the context of rapid vaccine development, long-term follow-up in large populations of vaccinated persons, as with older vaccines, is not possible³,⁴. When data are sparse, mathematical modeling is helpful because it can provide estimates of the duration of response by using additional information from biological knowledge about the vaccine mechanism and biological parameters. It is also helpful in quantifying the effect of factors that influence the response to the vaccine. This type of work is performed by modeling the dynamics of one or several markers that could be considered good correlates of protection⁵. Vaccine efficacy and mechanisms of action, or optimal immunogenic vaccine doses, have been evaluated for various infectious diseases, such as influenza⁶⁻⁸, yellow fever⁹,¹⁰, Zika¹¹, tuberculosis¹², and

more recently SARS-CoV-2¹³. In the case of the Ad26.ZEBOV and MVA-BN-Filo vaccine strategy, the concentration of binding antibodies is considered a good correlate of protection based on work performed in non-human primates¹⁴. It was agreed with the FDA to be suitable for use in a Biological Licensing Application under the Animal Rule¹⁵ and it was the basis for marketing authorization in the EU.

In a previous publication¹⁶, we used a mathematical model for antibody-secreting cell (ASC) dynamics that distinguishes between short-lived and long-lived cells (SL and LL, respectively), and we estimated the model parameters using the data from the available first Phase I studies. We found that antibody production is maintained by the population of long-lived cells with an estimated half-life of at least 5 years. New data from three Phase II studies¹⁷⁻¹⁹ conducted in two international consortia (EBOVAC1 and EBOVAC2) provided an opportunity to validate the model and better characterize factors associated with the variation of the antibody response.

## RESULTS

### Descriptive analysis of the data

The baseline and demographic characteristics of the 487 participants included in the study are shown in Table 1. In all results hereafter, Benjamini and Hochberg correction²⁰ for multiple testing has been used (see the section “Methods” for more details). Comparable baseline characteristics in terms of age, body mass index (BMI), and weight are observed in European participants across the Phase I and II clinical studies (all p-values > 0.80). Similarly, no differences are observed in Africa across trials and sites in terms of weight, however, BMI appears

<footnote n="1">¹Department of Public Health, Bordeaux University, Inserm UMR 1219 Bordeaux Population Health Research Center, Inria SISTM, Bordeaux, France. ²Vaccine Research Institute, Créteil, France. ³Janssen Vaccines and Prevention, Leiden, the Netherlands ⁴Present address: ExeVir, Ghent, Belgium. ⁵These authors contributed equally: Mélanie Prague, Rodolphe Thiébaut. *A list of authors and their affiliations appears at the end of the paper. ✉email: rodolphe.thiebaut@u-bordeaux.fr</footnote>
<table><caption>**Table 1.** Demographic and baseline characteristics of participants.</caption><thead><tr><th rowspan="3"></th><th colspan="3">Phase I trials</th><th colspan="3">Phase II trials</th><th rowspan="2">Total</th></tr><tr><th>Europe</th><th colspan="2">East Africa</th><th>Europe</th><th>East Africa</th><th>West Africa</th></tr><tr><th>UK<br>EBL1001</th><th>Kenya<br>EBL1003</th><th>Tanz./Ug.<br>EBL1004</th><th>UK/France<br>EBL2001</th><th>Ken./Tanz.<br>EBL2002</th><th>BFA/IVC</th><th>West Africa<br>Sierra Leone<br>EBL3001</th></tr></thead><tbody><tr><td>Part., no.</td><td>14</td><td>15</td><td>15</td><td>71</td><td>79</td><td>58</td><td>235</td><td>487</td></tr><tr><td>Sex</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Men</td><td>4 (29%)</td><td>11 (73%)</td><td>10 (67%)</td><td>32 (45%)</td><td>45 (57%)</td><td>44 (75%)</td><td>203 (86%)</td><td>349 (72%)</td></tr><tr><td>Women</td><td>10 (71%)</td><td>4 (27%)</td><td>5 (33%)</td><td>39 (55%)</td><td>34 (43%)</td><td>14 (24%)</td><td>32 (14%)</td><td>138 (28%)</td></tr><tr><td>Age (yrs)</td><td>37.6 (9.3)</td><td>23.7 (2.8)</td><td>26.5 (6.8)</td><td>41.2 (14.7)</td><td>34.1 (13.5)</td><td>34.1 (10.8)</td><td>27.2 (10.0)</td><td>31.3 (12.4)</td></tr><tr><td>BMI (kg/m²)</td><td>26.1 (3.3)</td><td>22.5 (4.1)</td><td>22.9 (4.2)</td><td>25.4 (4.5)</td><td>23.8 (4.0)</td><td>23.0 (3.4)</td><td>21.9 (3.3)</td><td>23.0 (3.9)</td></tr><tr><td>Weight (kg)</td><td>73.7 (13.7)</td><td>63.3 (12.7)</td><td>63.5 (11.7)</td><td>74.7 (14.7)</td><td>63.9 (10.2)</td><td>67.2 (9.9)</td><td>62.4 (9.4)</td><td>65.4 (11.7)</td></tr></tbody></table>

Data are *n* (%) or mean (SD). Only healthy adults receiving Ad26.ZEBOV followed by MVA-BN-Filo 56 days later were selected within each of the 6 trials.

*Part.* participants, *no.* number, *yrs* years, *kg* kilograms, *m* meter, *UK* United Kingdom, *Tanz.* Tanzania, *Ug.* Uganda, *BFA* Burkina Faso, *IVC* Ivory Coast.

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Dynamics of Ebola GP-specific binding antibody concentrations, in log₁₀ scale (ELISA units/mL, EU/mL) for each clinical study from 7 days after the second vaccination. Each subplot represents the antibody dynamics measured in one of the three accredited laboratories: Battelle (left side), Focus (middle) and Q² Solutions (right side). Each color corresponds to a clinical study (red: EBL1001, dark blue: EBL1003, light blue: EBL1004, orange: EBL2001, turquoise: EBL2002, light green: EBL3001). Solid and dashed lines represent medians in European and African participants, respectively. Circles correspond to Phase I studies and triangles to Phase II studies. Error bars correspond to 25th–75th confidence intervals. The vertical dotted line represents the first year after the first vaccination (309 days after the 2nd vaccination). The horizontal dot-dashed lines represent the LLOQ values considered for each laboratory (36.60, 26.22, and 36.11 EU/mL, at Battelle, Focus, and Q² Solutions, respectively).</figcaption></figure>

significantly higher in East African participants (+6%, *p*-value = 0.007) than in West African ones. European participants were significantly older than Africans (41 vs. 29 years, *p*-value < 0.001) and consequently, participants in EBL2002 tended to be older (34 vs. 27 years, *p*-value < 0.001). BMI and weight (*p*-values < 0.001 in both cases) were significantly higher in European participants (+13% and +18%, respectively) than in African participants.

Figure 1 shows the dynamics of antibody concentrations (median and interquartile ranges) 7 days after the second vaccination for each study according to the assay used to quantify the binding antibodies. In addition, Table 2 summarizes antibody concentrations observed at predefined sampling time points. Only participants who had received both the first and second vaccinations were included in both the descriptive and the modeling analyses. Similar kinetics were observed in all studies, with the highest binding antibody concentrations observed at 21 days post-dose 2 (hereafter referred to as “peak”), followed by a biphasic decline up to 1 year after the first vaccination. Furthermore, the longer-term dynamics observed in EBL3001 suggest a durable immune response after the biphasic decline.

**Mechanistic model of the humoral response**

To better identify the factors associated with the dynamic of the antibody response and to predict its duration, we used a model initially applied by Pasin et al.¹⁶ in Phase I trials evaluating the two-dose heterologous Ad26.ZEBOV, MVA-BN-Filo vaccine regi- men. In this mechanistic model, antibodies are assumed to be produced by plasma cells (antibody-secreting cells, ASCs) divided into two distinct sub-populations characterized by their lifespan: short-lived (SL) and long-lived (LL). For various infectious diseases, a rapid expansion of antigen-specific ASCs in blood peaking on day 7 post-infection or vaccination, followed by a fast depletion is observed²¹,²². Therefore, strictly decreasing dynamics were considered from 7 days after the second vaccination for the two compartments of plasma cells assuming no additional exposure to the antigen. A schematic diagram of the mathematical model used to describe the humoral response from 7 days after the second vaccination is displayed in Fig. 2. This simple model relied on three biological processes. LL and SL ASCs decay with time at rate *δ*<sub>L</sub> and *δ*<sub>S</sub>, respectively, and produce antibodies at rates *θ*<sub>L</sub> and
<table><caption>**Table 2.** Ebola glycoprotein-specific antibody concentrations (in log₁₀ ELISA units/mL) in each trial from 7 days after the second vaccination to study completion.</caption><thead><tr><th></th><th colspan="3">Phase I trials</th><th colspan="3">Phase II trials</th></tr><tr><th></th><th>EBL1001<br>(n = 14)<sup>a</sup></th><th>EBL1003<br>(n = 15)<sup>a</sup></th><th>EBL1004<br>(n = 15)<sup>a</sup></th><th>EBL2001<br>(n = 71)<sup>a</sup></th><th>EBL2002<br>(n = 137)<sup>b</sup></th><th>EBL3001<br>(n = 235)<sup>c</sup></th></tr></thead><tbody><tr><td colspan="7"><strong>Day 64 (7 days after the 2nd vaccination, MVA-BN-Filo)</strong></td></tr><tr><td>Number Part.</td><td>14</td><td>15</td><td>15</td><td></td><td></td><td></td></tr><tr><td>Positive Part.<sup>d</sup></td><td>14 (100%)</td><td>15 (100%)</td><td>15 (100%)</td><td></td><td></td><td></td></tr><tr><td>Missing data</td><td>0 (0%)</td><td>0 (0%)</td><td>0 (0%)</td><td></td><td></td><td></td></tr><tr><td>Mean [IQR]</td><td>3.19 [2.99; 3.48]</td><td>3.33 [2.90; 3.73]</td><td>3.09 [2.55; 3.58]</td><td></td><td></td><td></td></tr><tr><td colspan="7"><strong>Day 78 (21 days after the 2nd vaccination, MVA-BN-Filo)</strong></td></tr><tr><td>Number Part.</td><td>14</td><td>15</td><td>15</td><td>70</td><td>137</td><td>231</td></tr><tr><td>Positive Part.<sup>d</sup></td><td>14 (100%)</td><td>15 (100%)</td><td>15 (100%)</td><td>70 (100%)</td><td>137 (100%)</td><td>231 (100%)</td></tr><tr><td>Missing data</td><td>0 (0%)</td><td>0 (0%)</td><td>0 (0%)</td><td>1 (1%)</td><td>0 (0%)</td><td>4 (2%)</td></tr><tr><td>Mean [IQR]</td><td>3.88 [3.64; 4.10]</td><td>4.21 [3.96; 4.45]</td><td>4.03 [3.80; 4.31]</td><td>4.00 [3.79; 4.43]</td><td>3.88 [3.62; 4.16]</td><td>3.60 [3.34; 3.88]</td></tr><tr><td colspan="7"><strong>Day 156 (155 days after the 1st vaccination, Ad26.ZEBOV)</strong></td></tr><tr><td>Number Part.</td><td></td><td></td><td></td><td></td><td></td><td>42<sup>e</sup></td></tr><tr><td>Positive Part.<sup>d</sup></td><td></td><td></td><td></td><td></td><td></td><td>42 (100%)</td></tr><tr><td>Missing data</td><td></td><td></td><td></td><td></td><td></td><td>1 (2%)</td></tr><tr><td>Mean [IQR]</td><td></td><td></td><td></td><td></td><td></td><td>2.73 [2.54; 3.00]</td></tr><tr><td colspan="7"><strong>Day 180 (179 days after the 1st vaccination, Ad26.ZEBOV)</strong></td></tr><tr><td>Number Part.</td><td>12</td><td>15</td><td>15</td><td></td><td>23<sup>f</sup></td><td></td></tr><tr><td>Positive Part.<sup>d</sup></td><td>12 (100%)</td><td>15 (100%)</td><td>15 (100%)</td><td></td><td>23 (100%)</td><td></td></tr><tr><td>Missing data</td><td>2 (14%)</td><td>0 (0%)</td><td>0 (0%)</td><td></td><td></td><td></td></tr><tr><td>Mean [IQR]</td><td>3.47 [3.29; 3.60]</td><td>3.00 [2.67; 3.20]</td><td>2.97 [2.71; 3.25]</td><td></td><td>2.70 [2.47; 2.85]</td><td></td></tr><tr><td colspan="7"><strong>Day 240 (239 days after the 1st vaccination, Ad26.ZEBOV)</strong></td></tr><tr><td>Number Part.</td><td>13</td><td>15</td><td>15</td><td></td><td></td><td></td></tr><tr><td>Positive Part.<sup>d</sup></td><td>13 (100%)</td><td>15 (100%)</td><td>15 (100%)</td><td></td><td></td><td></td></tr><tr><td>Missing data</td><td>1 (7%)</td><td>0 (0%)</td><td>0 (0%)</td><td></td><td></td><td></td></tr><tr><td>Mean [IQR]</td><td>3.35 [3.20; 3.41]</td><td>2.66 [2.22; 2.93]</td><td>2.83 [2.56; 3.09]</td><td></td><td></td><td></td></tr><tr><td colspan="7"><strong>Day 360/365 (1 year after the 1st vaccination, Ad26.ZEBOV)</strong></td></tr><tr><td>Number Part.</td><td>12</td><td>15</td><td>15</td><td>51</td><td>134</td><td>207</td></tr><tr><td>Positive Part.<sup>d</sup></td><td>12 (100%)</td><td>15 (100%)</td><td>15 (100%)</td><td>51 (100%)</td><td>134 (100%)</td><td>205 (99%)</td></tr><tr><td>Missing data</td><td>2 (14%)</td><td>0 (0%)</td><td>0 (0%)</td><td>20 (28%)</td><td>3 (2%)</td><td>28 (12%)</td></tr><tr><td>Mean [IQR]</td><td>3.24 [3.09; 3.33]</td><td>2.61 [2.40; 2.96]</td><td>2.74 [2.46; 3.06]</td><td>3.07 [2.89; 3.28]</td><td>2.54 [2.26; 2.78]</td><td>2.44 [2.11; 2.68]</td></tr><tr><td colspan="7"><strong>Day 540 (539 days after the 1st vaccination, Ad26.ZEBOV)</strong></td></tr><tr><td>Number Part.</td><td></td><td></td><td></td><td></td><td></td><td>33<sup>e</sup></td></tr><tr><td>Positive Part.<sup>d</sup></td><td></td><td></td><td></td><td></td><td></td><td>33 (100%)</td></tr><tr><td>Missing data</td><td></td><td></td><td></td><td></td><td></td><td>10 (23%)</td></tr><tr><td>Mean [IQR]</td><td></td><td></td><td></td><td></td><td></td><td>2.43 [2.15; 2.68]</td></tr><tr><td colspan="7"><strong>Day 720 (2 years after the 1st vaccination, Ad26.ZEBOV)</strong></td></tr><tr><td>Number Part.</td><td></td><td></td><td></td><td></td><td></td><td>190</td></tr><tr><td>Positive Part.<sup>d</sup></td><td></td><td></td><td></td><td></td><td></td><td>184 (97%)</td></tr><tr><td>Missing data</td><td></td><td></td><td></td><td></td><td></td><td>45 (19%)</td></tr><tr><td>Mean [IQR]</td><td></td><td></td><td></td><td></td><td></td><td>2.45 [2.19; 2.69]</td></tr></tbody></table>

IQR: Interquartile range = 75% confidence intervals, Part.: Participants.

<sup>a</sup>Participants receiving the 2nd vaccination in the protocol-defined window of 57 ± 1 day.

<sup>b</sup>Participants receiving the 2nd vaccination in the protocol-defined window of 57 ± 3 days.

<sup>c</sup>Participants receiving the 2nd vaccination in the protocol-defined window of 57 ± 1 week.

<sup>d</sup>Refers to the number of participants with antibody concentration above the lower limit of quantification (LLOQ), expressed as n/N (%) where n is the number of participants with a concentration above the LLOQ (i.e., 36.11 EU/mL for Battelle lab, 26.22 EU/mL for Focus lab and 36.11 EU/mL for Q² Solutions lab) at that timepoint and N is the total number of participants with data at the first and the second vaccination and at that time point.

<sup>e</sup>Refers only to 43 participants enrolled in a substudy to receive a third dose (Ad26.ZEBOV) 2 years after the first vaccination.

<sup>f</sup>Refers to participants enrolled in EBL 2002 having an additional timepoint, initially scheduled for participants who do not receive a second vaccination because of a study pause.
<figure><img src="image_3.png" /><figcaption>**Fig. 2 Schematic diagram of the model describing the humoral immune response from 7 days after the 2nd vaccination.** S and L stand for short- and long-lived ASCs, respectively and Ab for antibodies. The parameters $\delta_S$, $\delta_L$, and $\delta_{Ab}$ are respectively the decay rates of SL ASCs, LL ASCs, and antibodies while $\theta_S$ and $\theta_L$ represent the production rates of antibodies by SL and LL ASCs.</figcaption></figure>

$\theta_S$, respectively. Finally, antibodies are assumed to decay over time at rate $\delta_{Ab}$. Since the baseline level of ASCs is unknown, the parameters $\phi_L = \theta_L L_0$ and $\phi_S = \theta_S S_0$ were defined, which represent the influx of LL and SL ASCs, respectively (see "Mathematical model of antibody kinetics" for more details).

## Quality of model prediction

Using parameter estimations obtained by Pasin et al.<sup>16</sup> on humoral response observed in Phase I trials, we evaluated the robustness of the model and its predictive abilities.

First, we looked at the capacity of the model to capture the dynamic of the antibodies during the early phase following vaccination, based on the previously estimated parameters, in a new population of participants. Fixing the antibody, short- and long-lived ASCs half-lives at 24 days, 3.0 days and 6.0 years, respectively, as well as the SL ASC influx parameter at 2755 ELISA units/mL/day and the LL ASC influx parameter at 16.6 and 70.7 ELISA units/mL/day for African and European participants (see the "Methods" subsection "Evaluation of the model quality of prediction" for more details), only random effects (i.e. individual deviation from population mean) for the 487 participants were evaluated using empirical Bayes estimates (EBEs).

When restricted from 7 days post-second vaccination to the peak of individual dynamics, the model predictions fit well with observed antibody concentrations (see Table 3). The overall percent of observations falling within the 95% individual prediction intervals (also referred to as the percent of coverage) was evaluated at 100%. The root mean squared error (RMSE) was consistent with a relatively small average difference between observations and predictions estimated at 0.012 $log_{10}$ ELISA units/ml (RMSE = 0.028 $log_{10}$ ELISA units/ml in Phase I studies accounting for two observations). The model provided similar goodness of predictions when observations extended to the first year following vaccination were used to update EBEs (Supplementary Table 1). The percent of coverage was evaluated at 100% and the RMSE at 0.045 $log_{10}$ ELISA units/ml (RMSE = 0.066 and 0.038 $log_{10}$ ELISA units/ml in Phase I and Phase II studies, respectively). These results confirmed the ability of the model,

estimated using only data from Phase I trials, to capture the antibody response in all additional participants included in Phase II trials.

Then, we looked at the ability of the model to predict antibody concentrations beyond the peak of the dynamics. Individual parameters assessed for the early part of the dynamics of humoral responses were then used to predict both short-term antibody responses between the peak and 1 year after the first vaccination, and long-term antibody responses between 1 year and 2 years. As described in Table 2, while participants from all trials were included in the analysis of short-term predictions, only participants from the EBL3001 clinical trial contributed to the analysis of long-term predictions because they were the only ones with a follow-up beyond 12 months. As shown in Fig. 3 and Table 3, the model demonstrated a high quality of short- and long-term predictions with a total of 98% of the observed antibody concentrations falling within the 95% individual prediction intervals. Nevertheless, the high width of 95% prediction intervals (2.297 and 2.479 $log_{10}$ ELISA units/mL for short- and long-term forecast, respectively) highlights a large uncertainty in individual model parameters and explains the high percentage of coverage. Plots of individual predictions for more participants having at least two observations from 7 days post-second vaccination are given in Supplementary Figs. 1–3. Similar work was done to evaluate the ability of the model to predict long-term antibody concentrations beyond 12 months when data from 7 days post-second vaccination to 1 year were used to estimate individual parameters (see Supplementary Figs. 4–6). These results highlight the benefit of using an additional short-term observation to improve long-term predictions (beyond 1 year after the second vaccination). The uncertainty of predictions was much lower (65% reduction in the size of prediction intervals) leading to a fair but smaller coverage (90% vs. 98%) and a significant improvement in the quality of the predictions (RMSE: 0.298 instead of 0.518; bias: -0.017 instead of 0.253).

## Additional insight on longevity of the humoral immune response

The model performed well in forecasting short- and long-term humoral response. However, the increase in RMSE and the width of 95% prediction intervals, and the decrease in the percent coverage beyond 12 months (see Table 3 and Supplementary Table 1) motivated an update of the parameters using all available data from both the Phase I and Phase II trials. First, the model was modified to include laboratory effects (Battelle, Focus, or $Q^2$ Solutions) in the observation model. Among the four observation models tested (see the "Methods" subsection "Update and re-estimation of the model" for more details), none was able to outperform the model without adjustment for laboratories. However, since laboratory effects reflect an observable reality, the adjusted observation model was constrained and the model including a scaling factor between the observations and the compartment Ab in natural scale, providing the lowest corrected Bayesian information criteria (BICc), was chosen. We secondly focused on the half-life of LL ASCs, $log(2)/\delta_L$. The estimation of the lower bound of the loss rate of LL ASCs $\delta_L$ was performed with a profile likelihood. Thanks to the longer follow-up available, the previous estimation of 5 years for the lower bound of the half-life of LL ASCs was updated to 15 years (Fig. 4). The method used to achieve this estimation is described in Supplementary Methods. In other words, since long-lived antibody-secreting cells are non-proliferating cells<sup>23</sup>, half of these antibody-secreting cells, which are produced at 7 days after the second vaccination, should persist for at least 15 years. Given this result, further estimations were performed with the parameter $\delta_L$ as fixed at the value corresponding to a lifespan of 15 years.
<table><caption>**Table 3.** Evaluation of the robustness and the quality of prediction of the model developed by Pasin et al.¹⁶.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">All trials</th><th colspan="3">Phase I trials</th><th colspan="3">Phase II trials</th></tr><tr><th>EBL1001</th><th>EBL1003</th><th>EBL1004</th><th>EBL2001</th><th>EBL2002</th><th>EBL3001</th></tr></thead><tbody><tr><td>Time of peak<sup>a</sup></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Mean [95% CI]</td><td>15.0 [5.0; 21.0]</td><td>13.4 [11.6; 14.6]</td><td>14.4 [13.0; 18.2]</td><td>13.4 [12.0; 14.0]</td><td>14.3 [13.0; 16.3]</td><td>14.0 [12.4; 15.0]</td><td>16.0 [3.9; 22.2]</td></tr><tr><td>Number of timepoints<sup>b</sup></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Mean [Min-Max]</td><td>1.09 [1.0; 2.0]</td><td>2.0 [2.0; 2.0]</td><td>2.0 [2.0; 2.0]</td><td>2.0 [2.0; 2.0]</td><td>1.0 [1.0; 1.0]</td><td>1.0 [1.0; 1.0]</td><td>1.0 [1.0; 1.0]</td></tr><tr><td colspan="8"><em>Predictions from 7 days post-2nd vaccination to the peak</em></td></tr><tr><td>RMSE<sup>c</sup></td><td>0.012</td><td>0.034</td><td>0.023</td><td>0.026</td><td>0.005</td><td>0.006</td><td>0.006</td></tr><tr><td>Coverage (%)</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td></tr><tr><td>Bias<sup>c</sup></td><td>0.002</td><td>0.018</td><td>-0.003</td><td>0.006</td><td>0.001</td><td>0.001</td><td>0.000</td></tr><tr><td>95% PI width</td><td>0.410</td><td>0.547</td><td>0.502</td><td>0.510</td><td>0.389</td><td>0.389</td><td>0.387</td></tr><tr><td colspan="8"><em>Short-term forecast from the peak to 1 year</em></td></tr><tr><td>RMSE<sup>c</sup></td><td>0.471</td><td>0.218</td><td>0.539</td><td>0.361</td><td>0.442</td><td>0.460</td><td>0.515</td></tr><tr><td>Coverage (%)</td><td>98.1</td><td>100</td><td>93.3</td><td>100</td><td>100</td><td>98.7</td><td>97.6</td></tr><tr><td>Bias<sup>c</sup></td><td>0.251</td><td>0.014</td><td>0.239</td><td>0.077</td><td>0.293</td><td>0.224</td><td>0.328</td></tr><tr><td>95% PI width</td><td>2.297</td><td>2.097</td><td>2.124</td><td>2.184</td><td>2.502</td><td>2.358</td><td>2.298</td></tr><tr><td colspan="8"><em>Long-term forecast beyond 1 year</em></td></tr><tr><td>RMSE<sup>c</sup></td><td>0.518</td><td></td><td></td><td></td><td></td><td></td><td>0.518</td></tr><tr><td>Coverage (%)</td><td>97.8</td><td></td><td></td><td></td><td></td><td></td><td>97.8</td></tr><tr><td>Bias<sup>c</sup></td><td>0.253</td><td></td><td></td><td></td><td></td><td></td><td>0.253</td></tr><tr><td>95% PI width</td><td>2.479</td><td></td><td></td><td></td><td></td><td></td><td>2.479</td></tr></tbody></table>

The model was estimated on Phase I data and individual parameters were assessed, for each participant of Phase I and Phase II trials, using observation from 7 days post-second vaccination (day 64) to the peak of individual dynamics.

CI confidence interval, PI prediction interval, RMSE root mean squared error.

<sup>a</sup>Time delay in days (Mean [95% CI]) from 7 days post-vaccination (day 64) to the peak (first local maximum).

<sup>b</sup>Number of observations from 7 days post-vaccination to the peak.

<sup>c</sup>Criteria calculated on the median of individual predictions.

The application of three algorithms of covariate selection (SCM, COSSAC, and SAMBA; see the “Methods” subsection “Update and re-estimation of the model” for more details) enabled us to identify factors influencing the dynamics of the humoral response. Although methods of covariate search differ from one algorithm to another, the adjustment of the biological parameters of the model for baseline characteristics, including demographics, selected by the different methods were quite consistent. All procedures led to the selection of an effect of continent, sex, and age on antibody responses. The same covariates were identified in the model without adjustment for laboratories.

The best model estimated the effect of the continent on $\phi_L$ (see Table 4 presenting a summary of parameter estimations). The mean value of $\phi_L$ was estimated at 36.6 ELISA units/mL/days in Europe compared to 10.2 ELISA unit/mL/days in African participants. These results are in accordance with those previously obtained with Phase I trial data¹⁶.

By adding more information with Phase II trial data, we also identified sex as another significant covariate for explaining the inter-individual variability of the decay rate of antibodies. Indeed, we estimated that antibodies have a significantly higher half-life in women (p-value estimated by Wald test <0.001) with an increase of the decay rate of 41% (95% confidence interval (CI): [24%; 60%]) for men as compared to women.

We also found that older age was associated with a decrease in the influx of short-lived ASCs (parameter $\phi_S$). For example, a 31-year-old participant (ages are assumed to be centered, see Table 1) displayed a mean value of SL ASCs influx of 3057 ELISA units/mL/days (see Table 4). Each additional year from this mean age induces a division of the resulting influx of SL ASCs by 7% (95% CI: [5%; 9%]). Therefore, for a participant 10 years older, its influx of SL ASCs will then be divided by 49% (see Table 4

footnote), corresponding to a decrease of the peak of its dynamics of 0.23 [0.22; 0.25] $log_{10}$ EU/mL.

Once the optimal structure was identified, we estimated the value of the parameters of the model as shown in Table 4, providing the model parameters estimated by Pasin et al.¹⁶ on Phase I data as well as the model parameters obtained on combined Phase I and II data. Figure 5 displays the dynamics estimated by the model, highlighting the goodness-of-fit of the data (the reader can refer to Supplementary Figs. 7 and 8 for additional results about model estimation and its goodness-of-fit).

Compared with the estimates we obtained from the Phase I data, the new estimates show a decrease in the magnitude of the $\phi_L$ parameter. This decrease is likely due to the significant increase in the half-life of the LL ASC from 6 to 15 years. Nevertheless, the mean $\phi_L$ remained four times higher in Europe than in Africa, similar to the approximation previously obtained using Phase I data¹⁶, where the mean $\phi_L$ was 4.3 times higher in Europe than in Africa. For the dynamics of SL ASCs, the parameter estimates remained quite stable, between the newly estimated model and the earlier estimates. As noted above, the information gained from longer follow-up allowed an update of the lower bound of the LL ASC decay rate. Similarly, the use of 443 additional participants improved the precision of the model parameter estimates. Indeed, the confidence intervals in the new estimates have become narrower, for each parameter $\phi_L$, $\phi_S$, or $\delta_S$, and are mostly included within the confidence intervals of the old estimates. Moreover, a comparison of the model estimates showed a slight increase in inter-individual variability for parameters $\phi_L$ and $\delta_{Ab}$ in the new model compared with the old one. The latter may be due to the use of additional data collected in a more heterogeneous population than in Phase I studies. However, adjustment of the parameter $\phi_S$ for the age of the participants reduced the
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Individual antibody concentrations predicted by the model, estimated on Phase I data, for a random sample of participants from the six clinical studies. Each subplot represents the individual antibody dynamics (in log₁₀ ELISA units/mL) from 7 days after the 2nd vaccination. For each participant, the vertical dashed line represents the time limit (individual peak of dynamics) between the predictions (on the left) and the forecasts (short-term in blue and long-term in orange). Plain dots correspond to observations used to evaluate individual parameters while circles are observations not used in parameter estimation. Shaded areas correspond to 95% individual prediction intervals (accounting for the uncertainty on the individual parameter estimation and the measurement error) and the solid lines correspond to the prediction of the model.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Profile likelihood on parameter $\delta_L$. The y-axis corresponds to the non-penalized log-likelihood computed by importance sampling for several values of LL ASCs half-life which needs to be maximized. The blue dotted vertical line represents the lower bound of the LL ASCs half-life estimated by profile likelihood by Pasin et al.¹⁶ on Phase I data. The red dashed vertical line represents the newly estimated lower bound using both Phase I and II data.</figcaption></figure>

<table><caption>**Table 4.** Model parameters estimated on Phase I participants by Pasin et al.¹⁶ and the new estimates obtained on pooled Phase I and Phase II data and considering adjustment for laboratory effects.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Meaning</th><th colspan="2">Phase I data</th><th colspan="2">Phase I & II data</th></tr><tr><th>Mean</th><th>95% CI</th><th>Mean</th><th>95% CI</th></tr></thead><tbody><tr><td colspan="6"><strong>Fixed effects</strong></td></tr><tr><td>δ<sub>Ab</sub></td><td>Antibody decay rate (day<sup>-1</sup>)</td><td>0.029</td><td>[0.027; 0.033]</td><td></td><td></td></tr><tr><td>Women</td><td></td><td></td><td></td><td>0.0251</td><td>[0.0223; 0.0283]</td></tr><tr><td>Men</td><td></td><td></td><td></td><td>0.0353</td><td>[0.0296; 0.0421]</td></tr><tr><td>log(2)/δ<sub>Ab</sub></td><td>Antibody half-life (days)</td><td>24</td><td>[22; 26]</td><td></td><td></td></tr><tr><td>Women</td><td></td><td></td><td></td><td>27.6</td><td>[24.5; 31.1]</td></tr><tr><td>Men</td><td></td><td></td><td></td><td>19.6</td><td>[16.4; 23.4]</td></tr><tr><td>δ<sub>S</sub></td><td>SL ASCs decay rate (day<sup>-1</sup>)</td><td>0.231</td><td>[0.15; 0.36]</td><td>0.333</td><td>[0.326; 0.340]</td></tr><tr><td>log(2)/δ<sub>S</sub></td><td>SL ASCs half-life (days)</td><td>3.0</td><td>[1.9; 4.7]</td><td>2.08</td><td>[2.04; 2.13]</td></tr><tr><td>δ<sub>L</sub></td><td>LL ASCs decay rate (year<sup>-1</sup>)</td><td>3.16 × 10<sup>-4</sup></td><td>[1.46; 7.03] × 10<sup>-4</sup></td><td>1.25 × 10<sup>-4</sup></td><td></td></tr><tr><td>log(2)/δ<sub>L</sub></td><td>LL ASCs half-life (years)</td><td>6.0</td><td>[2.7; 13]</td><td>15.0</td><td></td></tr><tr><td>φ<sub>S</sub></td><td>SL ASCs influx (EU/mL/day)</td><td>2755</td><td>[1852; 4100]</td><td></td><td></td></tr><tr><td>Mean Age (31.3 years)</td><td></td><td></td><td></td><td>3057</td><td>[2418; 3865]</td></tr><tr><td>FC ΔAge = + 1 year<sup>a</sup></td><td></td><td></td><td></td><td>0.934</td><td>[0.915; 0.954]</td></tr><tr><td>φ<sub>L</sub></td><td>LL ASCs influx (EU/mL/day)</td><td></td><td></td><td></td><td></td></tr><tr><td>African part.</td><td></td><td>16.6</td><td>[13.7; 20.1]</td><td>10.2</td><td>[9.01; 11.4]</td></tr><tr><td>Eur. part.</td><td></td><td>70.7</td><td>[54.0; 92.7]</td><td>36.6</td><td>[27.3; 49.2]</td></tr><tr><td>α</td><td>Scaling factor—lab effects</td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>focus</sub></td><td></td><td></td><td></td><td>1.04</td><td>[0.93; 1.16]</td></tr><tr><td>a<sub>Q2sol</sub></td><td></td><td></td><td></td><td>1.00</td><td>[0.98; 1.02]</td></tr><tr><td colspan="6"><strong>Random effects</strong></td></tr><tr><td>ω<sub>φS</sub></td><td>Sd of RE on φ<sub>S</sub></td><td>0.92</td><td>[0.83; 1.01]</td><td>0.84</td><td>[0.56; 1.13]</td></tr><tr><td>ω<sub>φL</sub></td><td>Sd of RE on φ<sub>L</sub></td><td>0.85</td><td>[0.78; 0.92]</td><td>0.88</td><td>[0.81; 0.96]</td></tr><tr><td>ω<sub>δAb</sub></td><td>Sd of RE on δ<sub>Ab</sub></td><td>0.30</td><td>[0.24; 0.36]</td><td>0.35</td><td>[0.29; 0.41]</td></tr><tr><td colspan="6"><strong>Error model</strong></td></tr><tr><td>σ<sub>Ab</sub></td><td>Sd of error model</td><td>0.10</td><td>[0.10; 0.10]</td><td>0.107</td><td>[0.101; 0.112]</td></tr></tbody></table>

CI confidence interval, EU ELISA units, Eur. European, FC fold change, LL ASCs long-lived antibody-secreting cells, Part. participants, RE random effects, SL ASCs short-lived antibody-secreting cells, Sd standard deviation.

aRepresents the multiplicative factor to apply to the value of $\phi_S$ obtained for the mean age, for an increase in participant age of 1 year: $\phi_S(\text{Mean Age} + 1\text{year}) = \phi_S(\text{Mean Age}) \times \text{FC}(\Delta\text{Age} = +1)$. Therefore, the percentage of decrease of $\phi_S$ for a participant X years older than the mean age is given by $100 \times (1 - \text{FC}(\Delta\text{Age} = +1))^X$.

unexplained inter-individual variability for the same parameter by 24%. Finally, it can be noticed the absence of effects of laboratory adjustment, with the two scaling factors $a_{focus}$ and $a_{Q2sol}$ estimated as non-significantly different from 1.00, at 1.04 and 1.00, respectively.

In order to evaluate the impact of fixing the parameter $\delta_L$, we performed a model averaging analysis²⁴,²⁵, which integrates model uncertainty in the value of $\delta_L$ in the calculation of the parameter confidence intervals. Results shown in Supplementary Tables 2 and 3 indicate very stable estimation.
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Individual antibody concentrations estimated by the model for a random sample of participants from the six clinical studies. Each subplot represents the individual antibody dynamics (in log₁₀ ELISA units/mL) from 7 days after the 2nd vaccination to 2 years. Colored circles correspond to observations used to estimate the model. The thick solid lines correspond to the individual dynamics and the 95% individual confidence intervals (accounting for the uncertainty of the estimation of the individual parameters only) are delimited by the shaded areas.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Evaluation of the ability of the model to predict unseen data using Monte-Carlo cross-validation. The predictive quality was assessed by the evaluation of two criteria: the RMSE (left side) and the percent coverage (right side). The x-axis corresponds to the percent of participants randomly selected for the training dataset and the y-axis to the value of criteria calculated on the testing dataset. One hundred replicates were performed for each train-test split percentage. Solid lines display the values of criteria and dashed lines, the 95% confidence intervals. The horizontal red dotted line on the right side displays the threshold of 95%.</figcaption></figure>

We examined the ability of our model to predict the response for new participants by performing Monte Carlo cross-validation (MCCV) and using RMSE and percent coverage as quality criteria for prediction. (see the “Methods” subsection “Update and re-estimation of the model” for more details). The results of this analysis are summarized in Fig. 6, where these two criteria are displayed as functions of the percent of participants used in the training dataset. Finally, despite the wide range of percentages tested for the split of train and test, the quality of the model prediction was very stable. The mean RMSE gradually decreased from 0.0870 to 0.0828 log₁₀ ELISA units/mL until it reached the value of 0.0843 log₁₀ ELISA units/mL when 100% of the data were used to estimate the models. The mean percent coverage remained higher than 95% even when only 20% of participants were used to estimate the model. Consequently, the model showed reasonably good quality in predicting the humoral immune response from 7 days after the second vaccination to two years after the first vaccination.

DISCUSSION

With this modeling work, we evaluated the quality of long-term predictions of the mechanistic model developed by Pasin et al.<sup>16</sup> which considers two populations of ASCs. We demonstrated with new data and a longer follow-up from phase 2 trials that the model, developed with a small number of participants from phase 1 trials, provides long-term predictions of the antibody response with high validity. Re-estimation of the model with a longer-term follow-up allowed us to update the value of the lower limit of the LL ASC half-life and showed that the longevity of LL plasma cells is much longer than previously estimated.

One advantage of this modeling work is the ability to distinguish the half-life of the antibodies from that of ASCs. The estimated values obtained for the half-life of the antibodies (δ<sub>Ab</sub>) were in the range of those reported in the literature of 20–50 days<sup>26–32</sup>. References on ASC half-lives are more difficult to figure out because these cells are circulating in various compartments (lymphoid, bone marrow) and the survival of the cells may vary greatly according to the infectious agents<sup>33</sup>. Hence, estimates of short-lived ASC vary from days<sup>34</sup> to several weeks<sup>3,35</sup>. The half-life of long-lived ASCs is even more variable and reflects the variability of the antibody dynamics according to infectious agents and type of vaccines<sup>36</sup>. The updated estimate for the long-lived ASCs calculated in the present study is particularly

encouraging with a half-life of at least 15 years. This is comparable to the half-life estimated for long-lived ASCs induced by the Hepatitis A virus vaccine³. Long-term humoral immunity is maintained through the memory B cells and the long-lived plasma cells³⁷. These latter cells, residing preferentially in the bone marrow, produce antibodies in the absence of antigenic stimulation. Interestingly, it has been recently demonstrated that these cells are intrinsically long-lived and can maintain the secretion of antibodies without replenishment of the pool of memory B cells³⁸. Specific niches in the bone marrow promote cell survival through various factors³⁶.

The model was able to capture inter-individual variation of the antibody dynamics. A part of this variability was associated with the geographic region, age, and sex. The cause of the influence of geographic regions on the humoral response to vaccines is still unknown. Concomitant malaria infection is suspected to play a role in compromising the immune response<sup>39–41</sup>. Nevertheless, these results should be interpreted cautiously, as simple nonspecific cross-reactivity within the assays used could also be responsible for this association<sup>40</sup>. Chronic parasitic infections, such as schistosomiasis<sup>42</sup>, have also been suspected to play a role in dampening immune response to vaccination<sup>43</sup>. Both hypotheses are currently explored in the EBOVAC consortia. The less pronounced decline of antibodies observed in women is consistent with several reports from other vaccines, including SARS-CoV-2<sup>44,45</sup>, Ebola<sup>46</sup>, and Flu<sup>47</sup>. The biological mechanism behind the lower clearance of binding antibodies in women remains an area of research and may differ from vaccine to vaccine<sup>47</sup>. It could also be a limitation of the modeling work in identifying precisely enough which compartment is influenced by sex differences in the absence of more immunological measurements. The influence of age on the response to the vaccine is known<sup>48</sup>, but the characterization of its effect through the production of short-lived antibody-secreting cells is novel and requires further confirmation. Pritz et al.<sup>49</sup> have noted an age-related decline in the number of plasma cells in human bone marrow. Here also, the identifiability of this effect could have been compromised because of a restricted range of participant ages because none of them were older than 65 years.

In conclusion, the dynamical model constructed from early Phase I data has demonstrated its predictive capacity, with longer follow-up and updated estimates giving promising results for the duration of the immune response. Nevertheless, the simplicity of this model significantly limits its use to fully describe the immune
response induced by a multi-dose vaccine strategy. Indeed, the model focuses exclusively on the decrease in antibody concentrations following vaccination, without modeling the establishment of the immune response triggered by each injection. Thereafter, more complex versions of this model have been proposed to model the establishment, reactivation, and persistence of the humoral immune responses induced by vaccination⁵⁰,⁵¹. In particular, these models integrate the stimulation of the immune system triggered by vaccine antigens and the role of immunological memory resulting from interactions with memory B cells and plasma cells.

The main findings obtained in this modeling work are not necessarily specific to the Ad26.ZEBOV, MVA-BN-Filo vaccine strategy. Some parameters, such as the long-lived ASCs half-life, could be found with other vaccine regimens. Nevertheless, the value and the interpretation of the parameters are dependent on the model structure and assumptions. For instance, in the modeling work conducted by Clairon et al.⁵¹ to model antibody dynamics induced by mRNA vaccine strategies against SARS-CoV-2, the absence of long-term persistence of antibody responses required an adaptation of the model for monophasic trajectories. In our work, the biphasic nature of the model fitted antibody trajectories perfectly, resulting in the estimation of the long-lived ASCs half-life.

The sporadic nature of Ebola outbreaks makes the identification of correlates of protection difficult⁵². Currently, neither a universal Ebola immune correlate of protection nor a binding antibody concentration threshold which would ensure a high probability of protection against EBOV, has been identified. Nevertheless, studies performed in non-human primates¹⁴ identified vaccine-induced binding antibody concentrations as the immune parameter most highly correlated with survival after the EBOV challenge. Due to the shorter disease course and full lethality of the Ebola disease model in non-human primates, the extrapolation of these results from animals to humans remains difficult. Consequently, no protective threshold of post-vaccination binding antibody concentration was derived for the Ad26.ZEBOV, MVA-BN-Filo vaccine.

# METHODS

## Ethics statement

The Phase I UK trial protocol and study documents were approved by the UK National Research Ethics Service. The Phase I Kenya trial protocol and study documents were reviewed and approved by the local Ethics Committee and the Kenyan regulatory authority. The Phase I Uganda/Tanzania trial protocol and study documents were reviewed and approved by the Tanzanian Medical Research Coordinating Committee of the National Institute for Medical Research, the Tanzania Food and Drugs Authority, the Uganda Virus Research Institute Research and Ethics Committee, the Uganda National Council for Science and Technology, the Uganda National Drug Regulatory Authority, and the Ethics Committee of the London School of Hygiene and Tropical Medicine. The Phase II UK/France trial protocol and study documents were approved by the French National Ethics Committee (CPP Ile de France III; 3287), the French Medicine Agency (150646A-61), the UK Medicines and Healthcare Products Regulatory Agency (MHRA), and the UK National Research Ethics Service (South Central, Oxford; A 15/SC/0211). The Phase II Kenya/Uganda/Burkina Faso/Ivory Coast trial protocol and study documents were approved by local and national independent Ethics Committees and Institutional Review Boards. The Phase II Sierra Leone trial protocol and study documents were approved by The study was approved by the Sierra Leone Ethics and Scientific Review Committee, the Pharmacy Board of Sierra Leone, and the London School of Hygiene & Tropical Medicine ethics committee.

These trials were conducted in accordance with the principles of good clinical practice and the Declaration of Helsinki, and all participants gave formal, written consent before undergoing any trial-related procedure.

## Immunogenicity measurements

We considered data from six studies aiming at evaluating the safety, tolerability, and immunogenicity of two-dose vaccine regimens with Ad26.ZEBOV and MVA-BN-Filo. Ad26.ZEBOV is a monovalent, recombinant, E1/E3-deleted, replication-defective, adenovirus type 26 vector vaccine encoding Ebola virus Mayinga variant GP, produced in PER.C6 human cells and injected as a single dose of 5 × 10¹⁰ viral particles. MVA-BN-Filo is a recombinant, replication-defective, modified vaccinia Ankara vector vaccine encoding Mayinga variant GP, Sudan virus Gulu variant GP, Marburg virus Musoke variant GP, and Tai Forest nucleoprotein. This multivalent vaccine was produced in chicken fibroblasts and injected at a dose of 1 × 10⁸ Infectious Units (Inf. U). Three of the six studies are randomized, observer-blinded, placebo-controlled Phase I trials on healthy volunteers aged 18–50 years. These studies were performed in four countries: the United Kingdom (UK), Kenya, Tanzania and Uganda. Results of the trials were previously described by Milligan et al.⁵³ and Winslow et al.⁵⁴ for the UK (study registered at ClinicalTrials.gov, NCT02313077, and labeled EBL1001 here), Mutua et al.⁵⁵ for Kenya (study registered at ClinicalTrials.gov, NCT02376426, and labeled EBL1003 here), and Anywaine et al.⁵⁶ for Tanzania/Uganda (study registered at ClinicalTrials.gov, NCT02376400, and labeled EBL1004 here). In addition, we considered data from two randomized, observer-blinded placebo-controlled, parallel-group Phase II trials on healthy volunteers aged 18–65 or 75 years. These studies were performed in six countries: the UK, France, Kenya, Uganda, Burkina Faso and Ivory Coast. We refer to Pollard et al.¹⁷ for a detailed description of results in the European trial and to Barry et al.¹⁸ for the African trial (two studies registered at ClinicalTrials.gov, NCT02416453 and NCT02564523, and labeled EBL2001 and EBL2002 here respectively for the European and African studies). The last study is a combined open-label, non-randomized stage 1, and a randomized, observer-blinded, placebo-controlled stage 2 Phase II trial on healthy adults. This study conducted in Sierra Leone also aimed to evaluate the long-term immunogenicity and the humoral immune memory induced by the vaccine regimen. Results of this trial were described by Ishola et al.¹⁹ (study registered at ClinicalTrials.gov, NCT02509494, and labeled EBL3001 here).

In Phase I trials, participants were equally randomized into four vaccination regimens: two with MVA-BN-Filo as the first vaccination on day 1, followed by Ad26.ZEBOV on day 29 or 57, and two with Ad26.ZEBOV was the prime vaccine on day 1, followed by MVA-BN-Filo on day 29 or 57. Within each regimen, participants received either an active vaccine or placebo in a 5:1 ratio. In the study EBL2001, participants in Cohorts I–III were equally randomized into three parallel groups in which they received Ad26.ZEBOV was the first vaccine on day 1, followed by MVA-BN-Filo on day 29, 57, or 85. This first cohort was excluded from the analysis as participants were enrolled to provide data only on safety and the timing of anti-Ebola virus GP ASCs responses. Within each group, participants received active vaccines or placebo in a 14:1 or 10:3 ratio in cohorts II and III respectively. In the study EBL2002, healthy adults (Cohort I) were equally randomized into the same three parallel groups with an active vaccine: placebo ratio of 5:1. Adults HIV-infected patients (Cohort IIa) and healthy children (Cohorts IIb and III) were not included in the analysis. Finally, in the study EBL3001, participants received either Ad26.ZEBOV as first vaccination on day 1 followed by MVA-BN-Filo on day 57, or MenACWY vaccine on day 1 and placebo on day 57, with a ratio of 1:0 and 3:1 in stage 1 and 2 respectively. In
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Design of EBOVAC 1 (EBL 1001, 1003, 1004, and 3001) and EBOVAC 2 (EBL 2001 and 2002) trials for participants receiving Ad26, MVA D57 as vaccine regimen. Immunogenicity measurements provide the concentration of IgG-binding antibodies against Ebola, as measured by ELISA (ELISA units/mL).</figcaption></figure>

this work, only participants receiving Ad26.ZEBOV as the first vaccination on day 1 and MVA-BN-Filo as the second vaccination in the protocol-defined window of 57 ± X days (Ad26/MVA D57; with X = 1 for Phase I trials and EBL2001, 3 for EBL2002 and 7 for EBL3001) were included. Based on these criteria, a total of 487 participants over all studies were included (among the 725 participants enrolled to receive Ad26/MVA D57, a total of 238 participants were excluded for not receiving their second dose (n = 108) or outside the protocol-defined window (n = 130)), 44 of whom where in Phase I studies, 71 in EBL2001, 137 in EBL2002 and 235 in EBL3001. In addition, the 168 participants receiving a placebo as a vaccine strategy were excluded.

Participants were followed up to 1 year after the first vaccination in all the studies, with longitudinal immunogenicity measurements performed on blood samples. As shown in Fig. 7, for the vaccine regimen of interest, immunogenicity samples were collected in all participants immediately before the administration of the first vaccination (Ad26.ZEBOV) on day 1, before the second vaccination (MVA-BN-Filo) on day 57, then 21 days after the second dose at day 78 and 1 year after the first dose (at day 360 or 365 according to the trial). In Phase I trials, additional samples were taken at days 7, 29, 64, 180, and 240, while immunological assays were done on blood samples taken at day 180 and day 156 in EBL2002 and EBL3001 respectively. Participants enrolled in EBL3001 were additionally followed up to 2 years after the first vaccination, with blood samples collected every 6 months after the first year. We analyzed total IgG Ebola virus GP-specific binding antibody concentrations measured by an Ebola virus GP (Kikwit strain) Filovirus Animal Non-Clinical Group (FANG) ELISA assay. The FANG ELISA assays were performed at three different

accredited laboratories: (a) at Battelle Biomedical Research Center (Columbus, OH, USA; hereafter referred to as Battelle) for the studies EBL1001 and EBL1004, (b) at Focus Diagnostics (San Juan Capistrano, CA, USA; hereafter referred to as Focus) for the study EBL1003, and (c) at Q² Solutions Laboratory (San Juan Capistrano, CA, USA; formerly Focus Diagnostics; hereafter referred to as Q² Solutions) for the studies EBL2001, EBL2002 and EBL3001. Particular attention has been paid in this work to account for a possible systematic difference in measurements induced by the distinct ELISA assays and thus between studies. Being interested in the longevity of the long-term immunity induced by the two-dose heterologous vaccine, similarly to Pasin et al.¹⁶, we mainly focused our analysis on immunogenicity measurements assessed after the second vaccination.

## Statistical analysis

A preliminary descriptive analysis was performed on the baseline and demographic characteristics of the 487 participants to describe and summarize the basic features of the data. Statistical differences among groups of participants were evaluated using classic t-tests (Welch's *t*-test in case of unequal variance, identified by a *F*-test, and Student *t*-test otherwise) implemented in *R*, and *p*-values were adjusted for test multiplicity with Benjamini and Hochberg correction²⁰ using the built-in *R* function *p-adjust*. Because of the difference in antibody concentrations measured by the distinct laboratories, comparisons of immunogenicity data between trials were not possible.

Finally, Spearman correlations between antibody concentrations measured 21 days after the second vaccination and
longer-term humoral responses were evaluated by integrating adjustment for test multiplicity on p-values.

## Mathematical model of antibody kinetics

To analyze the humoral immune response induced by the two-dose heterologous vaccine regimen Ad26.ZEBOV, MVA-BN-Filo against Ebola virus and evaluate the long-term immunogenicity, we used a mechanistic model divided into three parts. First of all, a mathematical model based on ordinary differential equations is defined to describe the dynamics of plasma cells and antibodies³. As shown in Fig. 2, antibodies are assumed to be produced by two plasma cell populations differentiated by their lifespan: short- and long-lived antibody-secreting cells (ASCs). Consequently, the ordinary differential equation (ODE) system contains three compartments: the short-lived ASCs (labeled S), the long-lived ASCs (labeled L) and the antibodies (Ab). Based on the hypothesis that antibody-secreting cells peaked at day 7 post-infection/ vaccination²¹,²², time was rescaled to consider only the antibody dynamics from 7 days after the second vaccination (day 64) during which plasma cells only decreased over time (t = time observation – 64). As demonstrated by Pasin et al.¹⁶, the model can be written as a single equation (1).

$$ \begin{cases} \frac{d\text{Ab}}{dt} = \phi_S e^{-\delta_S t} + \phi_L e^{-\delta_L t} - \delta_{\text{Ab}} \text{Ab} \\ \text{Ab}(t=0) = \text{Ab}_0 = \text{Ab}_{D64} \end{cases} \quad (1) $$

with $\delta_S$, $\delta_L$ and $\delta_{\text{Ab}}$ representing the average decay rates of SL ASCs, LL ASCs and antibodies, respectively. The parameters $\phi_S$ and $\phi_L$ are, respectively, the influx of SL and LL ASCs defined as $\phi_S = \theta_S S_0$ and $\phi_L = \theta_L L_0$, where $S_0 = S(t=0) = S_{D64}$ and $L_0 = L(t=0) = L_{D64}$ are the initial conditions at 7 days after the second vaccination and $\theta_S$ and $\theta_L$ are their respective antibody production rates. The initial antibody concentration $\text{Ab}_0$ is defined by the individual measure of antibody concentration at 7 days after the second vaccination. Keeping in mind that the antibody concentration can be unobserved at day 64 for some participants (see Fig. 7) while the decrease of the dynamics of ASCs is still assumed to start 7 days post-vaccination, an individual lag-time $T_i$ was introduced in Eq. (1). This lag-time represents the individual time interval between day 64 and the first observation following this specific time. The equation can then be written as follows.

$$ \begin{cases} \frac{d\text{Ab}}{dt} = \phi_S e^{-\delta_S(t+T_i)} + \phi_L e^{-\delta_L(t+T_i)} - \delta_{\text{Ab}} \text{Ab} \\ \text{Ab}(t=0) = \tilde{\text{Ab}}_0 = \text{Ab}_{D64+} \end{cases} \quad (2) $$

Based on this equation, time was rescaled for each individual such that the initial condition (Ab(t=0)) coincided with the first observation following day 64 ($\text{Ab}_{D64+}$), given: time = time observation – 64 – $T_i$. Therefore, for a participant with a first measurement at day 64 and observations at days {64, 78, 180, 240, 360}, the lag-time is null ($T_i=0$), rescaled time points of observations are given by {0, 14, 116, 179, 296} and $\text{Ab}_{D64+}$ is equal to measurement at day 64, $\text{Ab}_{D64}$. For a participant with a first measurement at day 78 and observations at days {78, 180, 240, 360}, the lag-time $T_i = 78 – 64 = 14$, rescaled time points of observations are then given by {0, 102, 162, 282} and $\text{Ab}_{D64+}$ is equal to measurement at day 78, $\text{Ab}_{D78}$. We estimated the five following biological parameters $\Psi = (\phi_S, \delta_S, \phi_L, \delta_L, \delta_{\text{Ab}})$. To account for inter-individual variability, we used a statistical model on which the five model parameters are assumed to be log-transformed, to ensure their positivity. Each parameter is then described by a mixed-effects model which depends on covariates. Each individual parameter $\Psi'_k$ for the participant i can be defined as follows, for $k = \{1, \dots, 5\}$.

$$ \log(\Psi'_k) = \log(\Psi_{k,0}) + \beta_k Z'_k + u'_k \quad (3) $$

where $\Psi_0$ is the fixed effect, $Z_k$ and $\beta_k$ are, respectively, the vectors of explanatory variables and regression coefficients related to the

biological parameter $\Psi_{k,i}$ and $u'_k$ is the individual random effect assumed to be normally distributed with the variance $\omega_k^2$. Random effects were assumed to be independent from each other. Based on results obtained in the previous work¹⁶, we assumed random effects on the influx parameters, $\phi_L$ and $\phi_S$, and on the decay rate of antibodies $\delta_{\text{Ab}}$.

For the observation model, we modeled the observed IgG binding antibody concentrations against the Kikwit glycoprotein from the six studies by the antibody ODE-compartment. We assumed an additive error model normally distributed on the $\log_{10}$ value of the antibody concentrations, with a variance $\sigma_{\text{Ab}}^2$. The antibody concentration for patient i at the jth time is given by

$$ Y(t_{ij}) = \log_{10}[\text{Ab}(\Psi^i, t_{ij})] + \varepsilon_{ij} \quad \varepsilon_{ij} \sim \mathcal{N}(0, \sigma_{\text{Ab}}^2) \quad (4) $$

## Model estimation

Mathematical and practical identifiability has been assessed in previous work¹⁶. Thus the parameter $\delta_L$ was estimated by profile likelihood⁵⁷ which consists of defining a grid of values for the parameter, sequentially setting the parameter $\delta_L$ at one of those different values, and estimating the model by maximizing the log-likelihood, given that the value of $\delta_L$. The resulting profile shows the maximum possible log-likelihood for each value of $\delta_L$ and has its maximum at the maximum likelihood estimate $\delta_L$. Other parameters were estimated by a population approach in which the model estimation relies on the estimation of the vector of population parameters including the fixed effects ($\Psi_0$), the regression coefficients ($\beta$), the standard deviation of random effects ($\omega$) and the standard deviation of the error model ($\sigma_{\text{Ab}}$). Model estimation was performed by the Monolix® software versions 2019R1 and 2019R2. This software uses the Stochastic Approximation Expectation-Maximization (SAEM) algorithm⁵⁸,⁵⁹ to estimate the population parameters with likelihood computed by importance sampling⁶⁰ and the Fisher information matrix calculated by stochastic approximation. Once population parameters are re-estimated, individual parameters are computed as empirical Bayes estimates (EBEs) representing the most likely values of the individual parameters, given individual data and population parameters. EBEs are calculated as the mode of the conditional parameter distribution by Markov-Chain Monte-Carlo (MCMC) procedure⁶¹ using the Metropolis-Hasting algorithm⁶² to compute the conditional distribution and the Nelder-Mead Simplex algorithm⁶³ to maximize it.

## Evaluation of the model quality of prediction

The mechanistic model described by Eqs. (2)–(4), initially estimated on Phase I data by Pasin et al.¹⁶, was validated on data from the six trials according to its quality of prediction. To this end, a two-step approach was applied: first, the robustness of the model was assessed by evaluating its ability to predict antibody dynamics from 7 days post-second vaccination to the peak of the dynamics (i.e., the first local maximum) for all Phase I and Phase II participants. Then, the ability of the model to forecast short-term (i.e., from the peak to 1 year after the first vaccination) and long-term antibody concentration (i.e., beyond 1 year following the first vaccination) was evaluated. Because validation of the mechanistic model estimated on Phase I data is sought here, no modification of the observation model defined in Eq. (4) was considered here to account for possible laboratory-induced effects in the measurement of antibody concentrations.

To investigate the robustness of the model initially estimated on Phase I data, only data restricted to the first year following the first vaccination were used to stay in the scope of applicability of the model (see Table 2 for a detailed description of the number of observations available at each time point). Consequently, for each
participant, the peak of its dynamics was sought during the first year (see Table 3 for a detailed description of individual times of peak in each trial). Assuming fixed effects and regression coefficients of the population parameters ($\Psi_{k,0}$ and $\beta_{k,}$, $\forall k \in \{1, \dots, 5\}$), distribution of random effects ($\omega_{k,}$, $\forall k \in \{1, \dots, 5\}$), as well as the standard deviation of the error model ($\sigma_{Ab}$) as fixed to previously obtained values, we evaluated individual parameters for the 487 participants, via the variables $u_{k}^{j}$, using the empirical Bayes estimates (EBEs) approach implemented in Monolix. As shown in Table 4, we fixed the decay rate of antibodies ($\delta_{Ab}$), SL ASCs ($\delta_{S}$) and LL ASCs ($\delta_{L}$) at values corresponding to half-lives of 24 days, 3 days and 6 years, respectively. The parameter $\phi_{S}$ was fixed at 2755 ELISA units/mL/day while $\phi_{L}$ was fixed at 16.6 ELISA units/mL/day for African participants and 70.7 ELISA units/mL/day for Europeans. East and West African participants were assumed to share the same value of LL ASCs influx. Finally, standard deviations of the inter-individual variability on the three parameters $\phi_{S}$, $\phi_{S}$ and $\delta_{Ab}$ were chosen as $\omega_{\phi_{S}} = 0.92$, $\omega_{\phi_{L}} = 0.85$ and $\omega_{\delta_{Ab}} = 0.30$. The parameter $\sigma_{Ab}$ was fixed at 0.10 (see ref. 16 or Table 4). To stay consistent with the model built on Phase I data, we included an adjustment for geographic region in the statistical model (binary variable equal to 0 in Africa and 1 in Europe) on $\phi_{L}$, as shown in the following equation:

$$ \log(\phi_{L}^{j}) = \log(\phi_{L,0}) + \beta_{\phi_{L,Eur}} \times 1_{i \in Eur} + u_{\phi_{L}}^{j} \quad (5) $$

For each individual, the 95% prediction interval⁶⁴ of the antibody dynamics was calculated and the percent coverage, defined as the percent of observations falling within the prediction interval, was assessed. Through these results, we highlighted the ability of the model to predict the very first antibody concentration measurements from 7 days post-second vaccination. Once these predictions were validated, individual parameters estimated on the early phase of the follow-up were used in the second step to quantify both the short- and long-term forecast skills of the model. To this end, we used the model to make individual predictions of antibody concentration between the peak and 2 years after the first vaccination. Predictions were then compared to observations and the percent of observations falling within the 95% individual prediction intervals was quantified.

Thereafter, the two-step approach was also applied for evaluating: first, the ability of the model to predict antibody dynamics from 7 days post-second vaccination to 1 year after the first vaccination (instead of the peak), and second, its ability to forecast antibody concentration beyond 1 year. This additional analysis was performed to identify whether the estimation of individual parameters on a longer follow-up can improve long-term predictions.

K-means clustering for longitudinal data⁶⁵ was performed to identify distinct trajectories of the dynamics of the humoral response. Using the kml R package⁶⁶, trajectories of antibody concentration from 7 days after the second vaccination to 2 years after the first vaccination were sequentially clustered into two and more clusters. Thereafter, we evaluated the percent coverage and the RMSE to investigate potential differences in prediction abilities according to underlying trajectories for each resulting partition.

## Update and re-estimation of the model

Once the quality of prediction of the mechanistic model was evaluated, an update of the model was performed in order to improve biological knowledge about the longevity of the long-term immune response induced by the two-dose heterologous vaccine regimen, Ad26.ZEBOV, MVA-BN-Filo. The low number of participants included in the three Phase I trials (177 participants, of whom only 44 received the Ad26/MVA D57 vaccine regimen) as well as the short-term follow-up of their immune response up to 1 year after the first vaccination tended to limit the precision of the estimation of the model parameters in the work conducted only

on Phase I trials. Despite the validation of the model according to its quality of prediction on additional data coming from the three Phase II trials (EBL 2001, 2002 and 3001), a re-estimation of the model using antibody dynamics from the 487 participants was performed to enhance and reinforce our understanding of the underlying biological processes leading to the long-term immunity following vaccination against Ebola. To account for the difference in measurements induced by the three distinct ELISA assays performed at Battelle, Focus and Q² Solutions laboratories, we assumed in the observation model an adjustment for laboratory effects, as shown on the following equation:

$$ Y(t_{ij}) = \log_{10}[a \times \text{Ab}(\Psi^{j}, t_{ij})] + \varepsilon_{ij} \quad \varepsilon_{ij} \sim \mathcal{N}(0, \sigma_{Ab}^{2}) \\ \alpha = \begin{cases} 1 & \text{if } i \in \text{Battelle} \\ a_{\text{focus}} & \text{if } i \in \text{Focus} \\ a_{\text{Q2sol}} & \text{if } i \in \text{Q}^{2} \text{ Solutions} \end{cases} \quad (6) $$

with $\alpha$ representing the proportional scaling factor (in natural scale) between the three laboratories, considering Battelle as the reference, and where the two parameters $a_{focus}$ and $a_{Q2sol}$ are estimated with the five other biological parameters ($\Psi$). To ensure their positivity, both of them are assumed to be log-transformed. Further investigations with other link-functions between $Y$ and Ab were conducted to model laboratory effects, such as a proportional relationship in the log10 scale, or with more complex functions like nonlinear sigmoid functions applied either in the natural or log10 scale. The function leading to the best model (i.e., lowest BICc value) before any covariate adjustment was kept. We also tested whether the accuracy of the three assays differed using different measurement error models. However, this modeling did not improve the fit and was therefore not retained (result not shown).

Participants from Phase I clinical studies, being monitored only during the first year following the first vaccination, provided information only on the early phase of the humoral response. In particular, the lack of information on long-term immunity made the estimation of the decay rate of the long-lived ASCs difficult. Long-lived ASCs are persistent plasma cells with a lifespan ranging from several months to the end of an individual's lifetime²³,⁶⁷⁻⁶⁹, therefore only an approximation of the lower bound of the confidence interval of their half-life ($\log(2)/\delta_{L}$) was possible. Using additional data from Phase II studies and, in particular, the humoral response measurements beyond 1 year, we performed a profile likelihood to identify whether enough information was available to precisely estimate the parameter $\delta_{L}$. Considering the statistical model found by Pasin et al.¹⁶, the model was estimated for multiple values of LL ASCs half-life ranging from 1 to 40 years. The profile likelihood was then drawn by maximizing the log-likelihood, computed by importance sampling⁶⁰, for each of those related models.

As a first estimation, a sequential Bayesian estimation was envisaged, that is using information provided by Phase I studies only through informative prior distribution for parameters. Maximum a posteriori (MAP) estimates, corresponding to a penalized maximum likelihood estimation⁷⁰, should then be obtained using humoral responses from only the 443 Phase II participants. However, the difference in sampling between Phase I and II studies, in particular the absence of data from 7 to 21 days after the second vaccination for Phase II participants (see Table 2), made estimation of the model difficult. The lack of information at the early stage of the dynamics induced practical identifiability issues for the parameters $\delta_{S}$ and $\phi_{S}$. To tackle this difficulty, all data were used to update the model. Random effects found on Phase I trials were kept, considering inter-individual variability on the parameter $\delta_{Ab}$ as well as on the ASCs influx, $\phi_{L}$ and $\phi_{S}$.

The statistical model was updated by performing a covariate selection. We applied the classic stepwise covariate modeling
(SCM) algorithm<sup>71,72</sup> which is a stepwise procedure with a forward selection followed by a backward elimination. In the forward selection, each parameter-covariate relationship is tested in turn and the relationship improving the model criteria (a corrected version of the Bayesian information criterion, BICc) the most is kept (the lower the better). Then the addition of a second covariate is tested. In the backward elimination, the removal of each parameter-covariate relationship selected in the first step is tested in an univariate manner. To verify the robustness of the results, two other algorithms of covariate selection in non-linear mixed effects models were performed, using BICc as model selection criteria: (1) the conditional sampling use for a stepwise approach based on correlation tests (COSSAC)<sup>72</sup>, and (2) the stochastic approximation for model building algorithm (SAMBA)<sup>73</sup>. The three algorithms were independently applied on an initial model without any covariates and tested the addition of the seven following potential covariates: Sex (=0 for women and =1 for men), Age, Weight, BMI, Continent (=0 for Africa and =1 for Europe), Region (=0 for East Africa, =1 for West Africa and =2 for Europe) and EBL3001 (=1 for participants from EBL3001 and 0 otherwise). Covariates such as Age, BMI, and Weight were centered around the mean value of the studied population (see Table 1). The parameter $\delta_L$, facing some identifiability issues due to the lack of measurements beyond two years, was removed from the covariate selection procedure. Based on their definition, the parameters $\alpha_{focus}$ and $\alpha_{Q2sol}$ were also excluded from this selection. The statistical significance of selected covariates was then evaluated using a Wald test. EBL3001 was the only study that had a follow-up beyond 1 year after the first vaccination and which was conducted in a single country (Sierra Leone). Therefore, the robustness of the results was analyzed to verify the short-term relevance of the selected covariates. To this end, the same procedure was performed on the model already adjusted for the selected covariates but considering only data up to 1 year after the first vaccination. At the end of the covariate selection procedure, an optimal model was obtained with the following statistical model (see section "Results" subsection "Additional insight on longevity of the humoral immune response" for more details).

$$ 
\left\{ 
\begin{aligned}
\log(\phi'_{L}) &= \log(\phi_{L,0}) + \beta_{\phi_{L,Eur}} \times \mathbb{1}_{i \in \text{Eur}} + u'_{\phi_L} \\
\log(\phi'_{S}) &= \log(\phi_{S,0}) + \beta_{\phi_{S,Age}} \times (\text{Age}_i - \bar{\text{Age}}) + u'_{\phi_S} \\
\log(\delta'_{Ab}) &= \log(\delta_{Ab,0}) + \beta_{\delta_{Ab,Men}} \times \mathbb{1}_{i \in \text{Men}} + u'_{\delta_{Ab}}
\end{aligned}
\right. 
\quad (7) 
$$

where $\text{Age}_i$ and $\bar{\text{Age}}$ are, respectively, the age of the participant $i$ and the average age of the participants and with $u'_{\phi_L} \sim \mathcal{N}(0, \omega_{\phi_L}^2)$, $u'_{\phi_S} \sim \mathcal{N}(0, \omega_{\phi_S}^2)$ and $u'_{\delta_{Ab}} \sim \mathcal{N}(0, \omega_{\delta_{Ab}}^2)$. Once the optimal model selected, its goodness of fit was checked and the robustness of the convergence of the estimation was assessed by using the convergence assessment tool implemented in Monolix which evaluated the robustness of the SAEM algorithm for numerous initial conditions.

The predictive quality of the newly estimated model was assessed by performing a Monte-Carlo cross-validation<sup>74</sup>. Participants from the overall dataset were randomly split into a training and a testing dataset, given a particular train-test split percentage. We ensure that the same ratio of participants in each trial was maintained within each of the two sub-datasets. Once the model was fitted on training data, EBEs resulting from this model were evaluated on test data, followed by the prediction of the individual antibody dynamics. Two criteria were then calculated on the testing dataset to estimate how accurately the predictive model performs: the percent coverage (the higher the better) and the RMSE (the lower the better). For each of the seven train-test split percentages {20%, 30%, 40%, 50%, 60%, 70%, 80%}, the procedure was replicated 1000 times.

## Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## DATA AVAILABILITY

Janssen has an agreement with the Yale Open Data Access (YODA) Project to serve as the independent review panel for the evaluation of requests for clinical study reports and participant-level data from investigators and physicians for scientific research that will advance medical knowledge and public health. Data will be made available following publication and approval by YODA of any formal requests with a defined analysis plan. For more information on this process or to make a request, please visit the Yoda Project site at http://yoda.yale.edu. The data sharing policy of Janssen Pharmaceutical Companies of Johnson & Johnson is available at https://www.janssen.com/clinical-trials/transparency.

## CODE AVAILABILITY

The original code (mlxtran models and R codes) developed in this work is available and free-of-cost on github (Inria SSTM Team) at the following link: https://github.com/s sistm/ModelingEbola.git.

Received: 23 March 2023; Accepted: 13 October 2023; Published online: 08 November 2023

## REFERENCES

1. FDA. First FDA-approved Vaccine for the Prevention of Ebola Virus Disease, Marking a Critical Milestone in Public Health Preparedness and Response (FDA, accessed 1 March 2019); https://www.fda.gov/news-events/press-announcements/first-fda-approved-vaccine-prevention-ebola-virus-disease-marking-critical-milestone-public-health (2019).

2. EC. Vaccine against Ebola: Commission Grants New Market Authorizations (EC, accessed 1 March 2022); https://ec.europa.eu/commission/presscorner/detail/en/IP_20_1248 (2021).

3. Andraud, M. et al. Living on three time scales: the dynamics of plasma cell and antibody populations illustrated for hepatitis a virus. PLoS Comput. Biol. 8, e1002418 (2012).

4. Monath, T. P. Yellow fever vaccine. Expert Rev. Vaccines 4, 553-574 (2005).

5. Plotkin, S. A. & Gilbert, P. B. Nomenclature for immune correlates of protection after vaccination. Clin. Infect. Dis. 54, 1615-1617 (2012).

6. Pawelek, K. A. et al. Modeling within-host dynamics of influenza virus infection including immune responses. PLoS Comput. Biol. 8, e1002588 (2012).

7. Lee, H. Y. et al. Simulation and prediction of the adaptive immune response to influenza a virus infection. J. Virol. 83, 7151-7165 (2009).

8. Xie, X.-T. et al. Within-host model of respiratory virus shedding and antibody response to h9n2 avian influenza virus vaccination and infection in chickens. Infect. Dis. Model. 6, 490-502 (2021).

9. Bonin, C. R. B., Fernandes, G. C., dos Santos, R. W. & Lobosco, M. A simplified mathematical-computational model of the immune response to the yellow fever vaccine. In 2017 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), 1425-1432 (IEEE, 2017).

10. Bonin, C. R., Fernandes, G. C., dos Santos, R. W. & Lobosco, M. A qualitatively validated mathematical-computational model of the immune response to the yellow fever vaccine. BMC Immunol. 19, 1-17 (2018).

11. Best, K. & Perelson, A. S. Mathematical modeling of within-host Zika virus dynamics. Immunol. Rev. 285, 81-96 (2018).

12. Rhodes, S. J. et al. Using vaccine immunostimulation/immunodynamic modelling methods to inform vaccine dose decision-making. Npj Vaccines 3, 36 (2018).

13. Perelson, A. S. & Ke, R. Mechanistic modeling of sars-cov-2 and other infectious diseases and the effects of therapeutics. Clin. Pharmacol. Ther. 109, 829-840 (2021).

14. Roozendaal, R. et al. Nonhuman primate to human immunobridging to infer the protective effect of an ebola virus vaccine candidate. NPJ Vaccines 5, 1-11 (2020).

15. FDA. Product Development Under the Animal Rule (FDA, accessed 1 March 2022); https://www.fda.gov/media/88625/download (2015).

16. Pasin, C. et al. Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation. J. Virol. 93, e00579-19 (2019).

17. Pollard, A. J. et al. Safety and immunogenicity of a two-dose heterologous ad26. ZEBOV and MVA-BN-Filo ebola vaccine regimen in adults in Europe (ebovac2): a randomised, observer-blind, participant-blind, placebo-controlled, phase 2 trial. Lancet Infect. Dis. 21, 493-506 (2021).
18. Barry, H. et al. Safety and immunogenicity of 2-dose heterologous ad26. ZEBOV, MVA-BN-Filo ebola vaccination in healthy and HIV-infected adults: a randomised, placebo-controlled phase ii clinical trial in Africa. *PLoS Med.* **18**, e1003813 (2021).

19. Ishola, D. et al. Safety and long-term immunogenicity of the two-dose heterologous ad26. ZEBOV and MVA-BN-Filo ebola vaccine regimen in adults in Sierra Leone: a combined open-label, non-randomised stage 1, and a randomised, double-blind, controlled stage 2 trial. *Lancet Infect. Dis.* **22**, 97–109 (2022).

20. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. *J. R. Stat. Soc.: Ser. B (Methodological)* **57**, 289–300 (1995).

21. Odendahl, M. et al. Generation of migratory antigen-specific plasma blasts and mobilization of resident plasma cells in a secondary immune response. *Blood* **105**, 1614–1621 (2005).

22. Huang, K.-Y. A. et al. Virus-specific antibody secreting cell, memory b-cell, and sero-antibody responses in the human influenza challenge model. *J. Infect. Dis.* **209**, 1354–1361 (2014).

23. Khodadadi, L., Cheng, Q., Radbruch, A. & Hiepe, F. The maintenance of memory plasma cells. *Front. Immunol.* **10**, 721 (2019).

24. Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding aic and bic in model selection. *Sociol. Methods Res.* **33**, 261–304 (2004).

25. Gonçalves, A., Mentré, F., Lemenuel-Diot, A. & Guedj, J. Model averaging in viral dynamic models. *AAPS J.* **22**, 1–11 (2020).

26. O'Dempsey, T. et al. Meningococcal antibody titres in infants of women immunised with meningococcal polysaccharide vaccine during pregnancy. *Arch. Dis. Childhood-Fetal Neonatal Ed.* **74**, F43–F46 (1996).

27. Zeitlin, L., Cone, R. A., Moench, T. R. & Whaley, K. J. Preventing infectious disease with passive immunization. *Microbes Infect.* **2**, 701–708 (2000).

28. Rosenthal, K. S. & Tan, M. J. *Rapid Review Microbiology and Immunology* E-Book 3rd edn (Elsevier Health Sciences, 2010).

29. Brinkhof, M. W., Mayorga, O., Bock, J., Heininger, U. & Herzog, C. Kinetics of maternally acquired anti-hepatitis a antibodies: prediction of waning based on maternal or cord blood antibody levels. *Vaccine* **31**, 1490–1495 (2013).

30. Vilajeliu, A. et al. Pertussis vaccination during pregnancy: antibody persistence in infants. *Vaccine* **34**, 3719–3722 (2016).

31. Voysey, M., Pollard, A. J., Sadarangani, M. & Fanshawe, T. R. Prevalence and decay of maternal pneumococcal and meningococcal antibodies: a meta-analysis of type-specific decay rates. *Vaccine* **35**, 5850–5857 (2017).

32. Adaken, C. et al. Ebola virus antibody decay-stimulation in a high proportion of survivors. *Nature* **590**, 468–472 (2021).

33. Amanna, I. J., Carlson, N. E. & Slifka, M. K. Duration of humoral immunity to common viral and vaccine antigens. *N. Engl. J. Med.* **357**, 1903–1915 (2007).

34. Ho, F., Lortan, J. E., MacLennan, I. C. & Khan, M. Distinct short-lived and long-lived antibody-producing cell populations. *Eur. J. Immunol.* **16**, 1297–1301 (1986).

35. Nguyen, D. C., Joyner, C. J., Sanz, I. & Lee, F. E.-H. Factors affecting early antibody secreting cell maturation into long-lived plasma cells. *Front. Immunol.* **10**, 2138 (2019).

36. Bhattacharya, D. Instructing durable humoral immunity for COVID-19 and other vaccinable diseases. *Immunity* **55**, 945–964 (2022).

37. Akkaya, M., Kwak, K. & Pierce, S. K. B cell memory: building two walls of protection against pathogens. *Nat. Rev. Immunol.* **20**, 229–238 (2020).

38. Langley, W. A. et al. Persistence of virus-specific antibody after depletion of memory b cells. *J. Virol.* **96**, e00026–22 (2022).

39. McArdle, A. J., Turkova, A. & Cunningham, A. J. When do co-infections matter? *Curr. Opin. Infect. Dis.* **31**, 209 (2018).

40. Mahon, B. E. et al. Baseline asymptomatic malaria infection and immunogenicity of recombinant vesicular stomatitis virus-zaire ebola virus envelope glycoprotein: Vaccine: the Sierra Leone trial to introduce a vaccine against ebola (strive). *J. Infect. Dis.* **224**, 1907–1915 (2021).

41. Williamson, W. & Greenwood, B. Impairment of the immune response to vaccination after acute malaria. *The Lancet* **311**, 1328–1329 (1978).

42. Dzhivhuho, G. A. et al. Chronic schistosomiasis suppresses HIV-specific responses to dna-mva and mva-gp140 env vaccine regimens despite antihelminthic treatment and increases helminth-associated pathology in a mouse model. *PLoS Pathog.* **14**, e1007182 (2018).

43. Borkow, G. & Bentwich, Z. Chronic parasite infections cause immune changes that could affect successful vaccination. *Trends Parasitol.* **24**, 243–245 (2008).

44. Levin, E. G. et al. Waning immune humoral response to bnt162b2 covid-19 vaccine over 6 months. *N. Engl. J. Med.* **385**, e84 (2021).

45. Grzelak, L. et al. Sex differences in the evolution of neutralizing antibodies to severe acute respiratory syndrome coronavirus 2. *J. Infect. Dis.* **224**, 983–988 (2021).

46. Simon, J. K. et al. Immunogenicity of rVSV6G-ZEBOV-GP ebola vaccine (ervebo®) in African clinical trial participants by age, sex, and baseline GP-ELISA titer: a post hoc analysis of three phase 2/3 trials. *Vaccine* **40**, 6599–6606 (2022).

47. Furman, D. et al. Systems analysis of sex differences reveals an immunosuppressive role for testosterone in the response to influenza vaccination. *Proc. Natl Acad. Sci. USA* **111**, 869–874 (2014).

48. Siegrist, C.-A. & Aspinall, R. B-cell responses to vaccination at the extremes of age. *Nat. Rev. Immunol.* **9**, 185–194 (2009).

49. Pritz, T. et al. Plasma cell numbers decrease in bone marrow of old patients. *Eur. J. Immunol.* **45**, 738–746 (2015).

50. Balelli, I. et al. A model for establishment, maintenance and reactivation of the immune response after vaccination against ebola virus. *J. Theor. Biol.* **495**, 110254 (2020).

51. Clairon, Q. et al. Modeling the kinetics of the neutralizing antibody response against SARS-CoV-2 variants after several administrations of Bnt162b2. *PLoS Comput. Biol.* **19**, e1011282 (2023).

52. Medaglini, D., Santoro, F. & Siegrist, C.-A. Correlates of vaccine-induced protective immunity against ebola virus disease. in *Seminars in Immunology*, Vol. 39 (ed. Ottenhoff, T. H. M.) 65–72 (Elsevier, 2018).

53. Milligan, I. D. et al. Safety and immunogenicity of novel adenovirus type 26-and modified vaccinia Ankara-vectored Ebola vaccines: a randomized clinical trial. *JAMA* **315**, 1610–1623 (2016).

54. Winslow, R. L. et al. Immune responses to novel adenovirus type 26 and modified vaccinia virus Ankara-vectored Ebola vaccines at 1 year. *JAMA* **317**, 1075–1077 (2017).

55. Mutua, G. et al. Safety and immunogenicity of a 2-dose heterologous vaccine regimen with ad26. ZEBOV and MVA-BN-Filo ebola vaccines: 12-month data from a phase 1 randomized clinical trial in Nairobi, Kenya. *J. Infect. Dis.* **220**, 57–67 (2019).

56. Anywaine, Z. et al. Safety and immunogenicity of a 2-dose heterologous vaccination regimen with ad26. ZEBOV and MVA-BN-Filo ebola vaccines: 12-month data from a phase 1 randomized clinical trial in Uganda and Tanzania. *J. Infect. Dis.* **220**, 46–56 (2019).

57. Cole, S. R., Chu, H. & Greenland, S. Maximum likelihood, profile likelihood, and penalized likelihood: a primer. *Am. J. Epidemiol.* **179**, 252–260 (2014).

58. Kuhn, E. & Lavielle, M. Maximum likelihood estimation in nonlinear mixed effects models. *Comput. Stat. Data Anal.* **49**, 1020–1038 (2005).

59. Delyon, B., Lavielle, M. & Moulines, E. Convergence of a stochastic approximation version of the em algorithm. *Ann. Stat.* **27**, 94–128 (1999).

60. Pinheiro, J. C. & Bates, D. M. Approximations to the log-likelihood function in the nonlinear mixed-effects model. *J. Comput. Graph. Stat.* **4**, 12–35 (1995).

61. Kuhn, E. & Lavielle, M. Coupling a stochastic approximation version of em with an MCMC procedure. *ESAIM: Probab. Stat.* **8**, 115–131 (2004).

62. Chib, S. & Greenberg, E. Understanding the metropolis-hastings algorithm. *Am. Stat.* **49**, 327–335 (1995).

63. Lagarias, J. C., Reeds, J. A., Wright, M. H. & Wright, P. E. Convergence properties of the Nelder-Mead simplex method in low dimensions. *SIAM J. Optim.* **9**, 112–147 (1998).

64. Kümmel, A., Bonate, P. L., Dingemanse, J. & Krause, A. Confidence and prediction intervals for pharmacometric models. *CPT: Pharmacomet. Syst. Pharmacol.* **7**, 360–373 (2018).

65. Genolini, C. & Falissard, B. Kml: k-means for longitudinal data. *Comput. Stat.* **25**, 317–328 (2010).

66. Genolini, C., Alacoque, X., Sentenac, M. & Arnaud, C. kml and kml3d: R packages to cluster longitudinal data. *J. Stat. Softw.* **65**, 1–34 (2015).

67. Slifka, M. K., Antia, R., Whitmire, J. K. & Ahmed, R. Humoral immunity due to long-lived plasma cells. *Immunity* **8**, 363–372 (1998).

68. Radbruch, A. et al. Competence and competition: the challenge of becoming a long-lived plasma cell. *Nat. Rev. Immunol.* **6**, 741–750 (2006).

69. Hammarlund, E. et al. Plasma cell survival in the absence of b cell memory. *Nat. Commun.* **8**, 1–11 (2017).

70. Drylewicz, J., Commenges, D. & Thiebaut, R. Maximum a posteriori estimation in dynamical models of primary HIV infection. *Stat. Commun. Infect. Dis.* **4** (2012).

71. Jonsson, E. N. & Karlsson, M. O. Automated covariate model building within nonmem. *Pharm. Res.* **15**, 1463–1468 (1998).

72. Ayral, G., Si Abdallah, J.-F., Magnard, C. & Chauvin, J. A novel method based on unbiased correlations tests for covariate selection in nonlinear mixed effects models: the Cossac approach. *CPT: Pharmacomet. Syst. Pharmacol.* **10**, 318–329 (2021).

73. Prague, M. & Lavielle, M. Samba: a novel method for fast automatic model building in nonlinear mixed-effects models. *CPT: Pharmacomet. Syst. Pharmacol.* **11**, 161–172 (2022).

74. Picard, R. R. & Cook, R. D. Cross-validation of regression models. *J. Am. Stat. Assoc.* **79**, 575–583 (1984).

## ACKNOWLEDGEMENTS

We thank the members of the EBOVAC1, EBOVAC2 and EBOVAC3 consortia. We also thank Thierry Van Effelterre (Janssen Pharmaceutica N.V., Beerse, Belgium) and Laura Solforosi (Janssen Vaccines & Prevention B.V., Leiden, Netherlands) for contribution in an earlier stage of the modeling project. We thank Anna Dari (Janssen Research and
Development, Beerse, Belgium) for extensive review. We thank Lixoft SAS for their support. Numerical computations were in part 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). This project has received funding from the Innovative Medicines Initiative 2 Joint Undertaking under grant agreement EBOVAC1 (grant nr. 115854), EBOVAC2 (grant no. 115861) and EBOVAC3 (grant no. 800176). This Joint Undertaking receives support from the European Union's Horizon 2020 research and innovation program and EFPIA. The funder of the study had no role in the study design, data collection, data analysis, data interpretation, or writing the report. The corresponding author had full access to all the data in the study and had final responsibility for the decision to submit for publication. This work was also supported by the Investissements d'Avenir program managed by the ANR under reference ANR-10-LABX-77. The funder had no role in the study design, data collection, data analysis, data interpretation, or writing the report.

## AUTHOR CONTRIBUTIONS

R.T. and M.D. coordinated the project and acquired funding, R.T. and M.P. supervised the work and formulated research goals and aims, R.T., M.P. and M.A. designed the methodological framework, M.A. analyzed the data, carried out the implementation of the computer code and performed formal analysis, V.B. performed data curation, M.A., M.P. and R.T. wrote the original manuscript. All authors contributed to the writing of the manuscript, read and approved the final manuscript.

## COMPETING INTERESTS

C.M. and M.D. declare no competing interests but are full-time employees of Janssen Vaccines and Prevention, B.V., and may hold shares of Johnson & Johnson. V.B. declares no competing interests but was a full-time employee of Janssen Vaccines and Prevention at the time of the study and may hold shares of Johnson & Johnson. The other authors declare no competing interests.

## ADDITIONAL INFORMATION

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41541-023-00767-y.

**Correspondence** and requests for materials should be addressed to Rodolphe Thiébaut.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023

## FOR THE EBOVAC 1 AND EBOVAC 2 CONSORTIA

Thierry Van Effelterre<sup>4,7</sup>, Laura Solforosi<sup>3</sup>, Viki Bockstal<sup>3</sup>, Macaya Douoguih<sup>1,3</sup>, Chelsea McLean<sup>3</sup> and Anna Dari<sup>5</sup>

<footnote n="4"><sup>4</sup>Janssen Pharmaceutica N.V., Beerse, Belgium. <sup>5</sup>Janssen Research and Development, Beerse, Belgium. <sup>7</sup>Janssen Pharmaceutica N.V., Global Commercial Strategy Organization, Beerse, Belgium.</footnote>