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 × 1010 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 × 108 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.53 and Winslow et al.54 for the UK (study registered at ClinicalTrials.gov, NCT02313077, and labeled EBL1001 here), Mutua et al.55 for Kenya (study registered at ClinicalTrials.gov, NCT02376426, and labeled EBL1003 here), and Anywaine et al.56 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.17 for a detailed description of results in the European trial and to Barry et al.18 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.19 (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 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 Q2 Solutions Laboratory (San Juan Capistrano, CA, USA; formerly Focus Diagnostics; hereafter referred to as Q2 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.16, we mainly focused our analysis on immunogenicity measurements assessed after the second vaccination.

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).

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 correction20 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 antibodies3. 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/vaccination21,22, 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.16, the model can be written as a single equation (1).

$$\left\\frac}}}t}=_}}}}^_}}t}+_}}^_}}t}-_}}}\\ }(t=0)=}}_=}}_\end\right.$$

(1)

with δS, δL and δAb representing the average decay rates of SL ASCs, LL ASCs and antibodies, respectively. The parameters ϕS and ϕL are, respectively, the influx of SL and LL ASCs defined as ϕS = θSS0 and ϕL = θLL0, where S0 = S(t = 0) = SD64 and L0 = L(t = 0) = LD64 are the initial conditions at 7 days after the second vaccination and θS and θL are their respective antibody production rates. The initial antibody concentration Ab0 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 Ti 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.

$$\left\\frac}}}t}=_}}}}^_(t+_)}+_}}}}^_}}(t+_)}-_}}}\\ }(t=0)=}}}_=}}_\end\right.$$

(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 (AbD64+), given: time = time observation−64−Ti. Therefore, for a participant with a first measurement at day 64 and observations at days , the lag-time is null (Ti = 0), rescaled time points of observations are given by and AbD64+ is equal to measurement at day 64, AbD64. For a participant with a first measurement at day 78 and observations at days , the lag-time Ti = 78−64 = 14, rescaled time points of observations are then given by and AbD64+ is equal to measurement at day 78, AbD78. We estimated the five following biological parameters Ψ = (ϕS, δS, ϕL, δL, δ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 \(}}_^\) for the participant i can be defined as follows, for k = .

$$\log \left(}}_^\right)=\log (}}_)+__^+_^$$

(3)

where Ψ0 is the fixed effect, Zk and βk are, respectively, the vectors of explanatory variables and regression coefficients related to the biological parameter Ψk, and \(_^\) is the individual random effect assumed to be normally distributed with the variance \(_^\). Random effects were assumed to be independent from each other. Based on results obtained in the previous work16, we assumed random effects on the influx parameters, ϕL and ϕS, and on the decay rate of antibodies δ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 log10 value of the antibody concentrations, with a variance \(_}}^\). The antibody concentration for patient i at the jth time is given by

$$Y(_)=_\left[}(}}^,_)\right]+_\quad _ \sim }}}(0,_^)$$

(4)

Model estimation

Mathematical and practical identifiability has been assessed in previous work16. Thus the parameter δL was estimated by profile likelihood57 which consists of defining a grid of values for the parameter, sequentially setting the parameter δL at one of those different values, and estimating the model by maximizing the log-likelihood, given that the value of δL. The resulting profile shows the maximum possible log-likelihood for each value of δL and has its maximum at the maximum likelihood estimate \(}_}}\). 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 (Ψ0), the regression coefficients (β), the standard deviation of random effects (ω) and the standard deviation of the error model (σAb). Model estimation was performed by the Monolix ®software versions 2019R1 and 2019R2. This software uses the Stochastic Approximation Expectation-Maximization (SAEM) algorithm58,59 to estimate the population parameters with likelihood computed by importance sampling60 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) procedure61 using the Metropolis–Hasting algorithm62 to compute the conditional distribution and the Nelder–Mead Simplex algorithm63 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.16, 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 (Ψk,0 and βk, ∀ k ∈ ), distribution of random effects (ωk, ∀ k ∈ ), as well as the standard deviation of the error model (σAb) as fixed to previously obtained values, we evaluated individual parameters for the 487 participants, via the variables \(_^\), using the empirical Bayes estimates (EBEs) approach implemented in Monolix. As shown in Table 4, we fixed the decay rate of antibodies (δAb), SL ASCs (δS) and LL ASCs (δL) at values corresponding to half-lives of 24 days, 3 days and 6 years, respectively. The parameter ϕS was fixed at 2755 ELISA units/mL/day while ϕ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 ϕS, ϕS and δAb were chosen as \(__}}}=0.92,__}}}=0.85\) and \(__}}}=0.30\). The parameter σ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 ϕL, as shown in the following equation:

$$\log (_}}^)=\log (_},0})+__},}}}\times }}_}}+__}}}^$$

(5)

For each individual, the 95% prediction interval64 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 data65 was performed to identify distinct trajectories of the dynamics of the humoral response. Using the kml R package66, 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 Q2 Solutions laboratories, we assumed in the observation model an adjustment for laboratory effects, as shown on the following equation:

$$\beginY(_)\,=\,_\left[\alpha \times }(}}^,_)\right]+_\quad _ \sim }}}(0,_}}^)\\ \qquad\alpha \,=\,\left\1\,\qquad}}\,\,i\in \,\,}}\,\\ _}}\,\,}}\,\,i\in \,\,}}\,\\ _}}\,\,}}\,\,i\in \,}}}^\,\,}}\,\end\right.\end$$

(6)

with α representing the proportional scaling factor (in natural scale) between the three laboratories, considering Battelle as the reference, and where the two parameters αfocus and αQ2sol are estimated with the five other biological parameters (Ψ). 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 lifetime23,67,68,69, therefore only an approximation of the lower bound of the confidence interval of their half-life (\(\log (2)/_}}\)) 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 δL. Considering the statistical model found by Pasin et al.16, 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 sampling60, 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 estimation70, 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 δS and ϕ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 δAb as well as on the ASCs influx, ϕL and ϕS.

The statistical model was updated by performing a covariate selection. We applied the classic stepwise covariate modeling (SCM) algorithm71,72 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)72, and (2) the stochastic approximation for model building algorithm (SAMBA)73. 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 δ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 αfocus and α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\\log (_}}^)=\log (_},0})+__}}}\times }}_}}+__}}}^\\ \log (_}}^)=\log (_},0})+__}}}\times \left(}}_-\overline}}\right)+__}}}^\\ \log (_}}^)=\log (_},0})+__}}}\times }}_}}+__}}}^\end\right.$$

(7)

where Agei and \(\overline}}\) are, respectively, the age of the participant i and the average age of the participants and with \(__}}}^ \sim }}}(0,__}}}^),__}}}^ \sim }}}(0,__}}}^)\) and \(__}}}^ \sim }}}(0,__}}}^)\). 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-validation74. 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 , 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.