Data

The persistence data were time series of hSBA titers from adolescents vaccinated with two doses (priming vaccination) of 4CMenB vaccine given 1 to 6 months apart in six phase 2 and phase 3 clinical studies (Table 4). Only participants that satisfied the following criteria were included: they received both the priming doses, and had at least two datapoints for at least one antigen, with the last occurring at least two years after priming. The number of data points varied across countries, from a minimum of 136 in the USA (17 participants)16,17, up to 1,929 in Chile (131 participants vaccinated with three different 4CMenB schedules13,14,15). The longest follow-up was for participants recruited in Chile, where data were collected up to 7.9 years after complete vaccination; the shortest follow-up was just over two years after a 0–2-months schedule for participants from Poland and the USA16,17. Participants from Australia (38 participants12,15) and Canada (106 participants12,15) received two doses given one month apart and were followed for more than four years.

Table 4 Summary of the data to which all the models were fitted

Post-booster data (one month after a booster dose, given at any time to subjects that had already received the two priming doses) was available for 98% (n = 306) of the participants with post-priming persistence data. Participants from Poland and the USA did not receive the same vaccine given at priming (4CMenB), but an investigational vaccine that incorporates 4CMenB in combination with a multivalent meningococcal MenACWY-CRM conjugate vaccine16,17.

hSBA is the correlate of protection of meningococcal disease and is used to evaluate the efficacy of meningococcal vaccines10. Bactericidal titers were obtained for four MenB indicator strains, each exhibiting one, and only one, of the four vaccine antigens, thus indicating immune response to each vaccine component. Strain 44/76 was used as indicator for fHbp, strain 5/99 as indicator for NadA, strain M10713 (M07-0241084 for studies in USA and Poland) for NHBA, strain NZ98/254 for PorA. Titers were calculated as the reciprocal of a dilution (e.g., a dilution of 1:4 corresponds to titer 4). The lowest detectable titer was 2, values under the limit of detection were unknown and indicated as “<2”.

Only adolescents that received two doses of 4CMenB 1 to 6 months apart and with at least one hSBA result from a blood draw taken at least two years after the second dose were included in the analysis. Data from blood draws taken after booster doses were excluded. For each hSBA titer, time was calculated as number of days intercurred between second dose and blood draw.

All MenB typing data used here were derived from a distinct study where a panel of 442 strains (representative of MenB disease in the USA)20 were collected and tested using the MATS. This assay was developed to predict if a MenB strain is covered or not by 4CMenB, and by which specific antigens (coverage by one antigen being sufficient to predict bacterial killing)20.

Mathematical models for the evolution of immune response over time

Several mathematical models were developed and tested to predict four outcomes i.e., the four time-series of hSBA titers against 4CMenB vaccine antigenic indicator strains, measured in sera from adolescents that received two vaccine doses. Candidate models were characterized by different levels of heterogeneity of immune response and different combinations of two models reproducing the natural evolution of antibody titers after vaccination.

The two antibody evolution models were the power law model and the exponential model, which have been widely used to reproduce the decline of antibodies elicited by vaccines against different viral and bacterial diseases, including human papillomavirus, hepatitis, pertussis, and meningococcal disease caused by different serogroups21,22,29,30,31. In the power law model, for each participant i the bactericidal titer T evolves over time t after the last dose as:

$$_\left(t\right)=__}\right)}^_},}\,\,t\ge _.$$

At time t = tmin, T equals A. We set tmin = 15 days, so that A represents the bactericidal activity 15 days after vaccination, when its effect should be nearly at peak. Parameter b is the power law exponent that determines the evolution of T over time. Titer T(t) declines over time t for negative values of b. After log-transforming, the model becomes linear with respect to the logarithm of time logt:

$$log\, _\left(t\right)=_+_\left(-log\, _\right),$$

(1)

where \(_=\log _\) is the intercept of the corresponding logt-linear model.

In the exponential model, for each participant i the bactericidal titer T evolves over time t after the last dose as:

$$_\left(t\right)=_^_\frac_}_}}.$$

As for the power law model, T = A when t = 15 days after vaccination. Parameter c is the evolution rate of T, rescaled by tmin. After log-transforming, the model becomes linear with respect to time:

$$_\left(t\right)=_+_\frac_}_}.$$

(2)

For this model, \(_=\log _\) is the intercept of the corresponding t-linear model.

In the different models tested, some or all of the parameters a, b and c were fixed for all the participants, or varied between individuals in a hierarchical way (also called mixed effects model). In hierarchical models, the parameters \(\theta =a,b,c\) were assumed to be normally distributed with mean \(_\) and variance \(_\):

$$\theta \sim N\left(_;_\right).$$

where \(\theta =a,b,c\). Models are no more hierarchical if dispersion is removed by setting \(_\) to zero or, equivalently, by identically assigning \(\theta =\,_\) to all the participants in a stratum. Moreover, model changes in performance were tested when implementing stratification by study country or not.

Post-booster models

All the available post-booster hSBA titers were measured one month after booster doses. The post-booster model was developed as a continuation of the post-priming model, thus conditioned on the inferred post-priming posteriors, to account for correlations between titers after booster and post-priming response and persistence of protection. It required the introduction of one additional parameter.

Through model comparison, two different hypotheses were tested when analyzing data collected one month after booster. The first was that immune response depends on time intercurred between priming vaccination and booster dose. This was obtained by adding the additional parameter \(_}}\) to Eqs. (1) and (2), that is switched on after booster:

$$_\left(t=t^\right)=_\left(t\right)+_}}.$$

Where \(t^\) is one month after booster.

The second hypothesis is that immune response one month after booster is not dependent on the time that intercurred between priming and booster. This was modeled through a k parameter that is not summed to Eqs. (1) and (2), but substitutes them at one month after booster:

$$_\left(t=t^ \right)=_$$

Moreover, three different levels of model hierarchy and stratifications for post-booster variance parameters h and k were evaluated, for a total of six models.

Model fitting and evaluation

Titers were all logarithmically transformed before fitting. The models were Bayesian, with non-informative priors and left-censored normal likelihoods. The likelihoods were left-censored to allow evaluation even when titers were reported to be under the limit of detection, thus avoiding any a-priori imputation of titer data for undetected bactericidal activity that can introduce bias when fitting longitudinal data32. Bayesian full posterior distributions of model parameters were inferred with a Markov Chain Monte Carlo (MCMC) algorithm from the Python PyMC3 package33.

Models were compared through their expected log pointwise predictive densities calculated using the WAIC. A higher WAIC score indicates a model with better predictive accuracy34.

Simulation of immunogenicity and persistence curves

Posterior distributions inferred from the best-fitting post-priming and post-booster models were used to generate immunogenicity curves after the two primary doses (in absence of booster) and after booster dose. For post-priming simulation, antigens log titers were simulated to evolve as in Eqs. (1) or (2), depending on model best fit antigen data. Post-booster log titers were simulated assuming to evolve in the same way as post-priming, but starting from a one-month-post-booster intercept. In practice, since the best model was the one in which booster immunogenicity did not depend on time between priming and booster, logT followed Eqs. (1) or (2) (depending on which was antigen’s best model after priming) with \(a\) substituded by \(k\).

For both post-priming in absence of booster and post-booster, the procedure consisted of computationally simulating 1000 times the vaccination of 5000 hypothetical participants from each of the five considered countries, and their bactericidal titers against each antigen at different times. Synthetic data were first generated by country, then merged into “overall” data.

The simulated hSBA data were then used to calculate the proportion of participants with bactericidal activity against the four antigens at each time point, by country, and overall. Participants with an hSBA titer equal or superior to four were considered protected10. Overall proportions represent the expected value of the bactericidal activity, when combining evidence from different countries.

From antigens persistence to vaccine persistence using MenB typing data

Since there are four 4CMenB antigens which can be co-expressed, there are 16 theoretically possible types of MenB bacteria, relative to the vaccine-preventability of the disease they may cause. Four types correspond to strains that harbor one and only one antigen (like the indicator strains used for the hSBA). Six possible types include strains predicted to be killed by two antigens e.g., type fHbp-NHBA. Another four types are covered by three antigens (e.g., fHbp-NHBA-PorA), and one by all the four 4CMenB antigens. Finally, the sixteenth type simply includes all the MenB strains not covered by the vaccine, thus those MenB bacteria expressing no antigen corresponding to 4CMenB.

Bacterial killing (binary outcome) for the 15 covered types were derived at each time point for each simulated participant, in each country, for each posterior trajectory, assuming that a titer over threshold for one antigen is sufficient to elicit killing. For example, considering only two antigens for simplicity: if one year after vaccination a participant has an fHbp titer of 3 (<4, therefore not enough for killing) and an NHBA titer of 5 (enough for killing), the participant will be protected from NHBA and fHbp-NHBA types, but not from fHbp types, and of course not from “no-antigen” types. If after some time NHBA also falls below the protective threshold, the participant is no longer protected from any strain. This criterion was systematically applied to all simulations and derived immunogenicity curves (proportion of protected participants) for each of the 15 types. The sixteenth type has no curve, as there is no protection.

The prevalence of MenB types assessed in the USA panel of 442 strains was then used to weight such curves by the relative abundance of each strain type (with or without including the non-covered types, representing 8.8% of the panel)20. The sum of the weighted curves represents the final 4CMenB persistence, at each time point after vaccination. This 4CMenB-typing-based weighted average combines the previously calculated antigen persistence curves. The procedure was repeated, including and excluding non-covered types. In the first case, the results show protection against any MenB strain, while in the second case, they represent protection strictly from 4CMenB-preventable MenB strains. The same approach was applied for both post-priming and post-booster persistence.

Comparison with a compartmental model

The predicted proportion of protected subjects over time based on an assumed antibody titer evolution model (either exponential or power law), was compared to a compartmental model. In scenarios where longitudinal data on antibody titers are not available, the persistence of protection is usually assumed, such as by estimating the time at which the proportion of immunized subjects is halved. A commonly used model to predict the proportion of immune subjects over time is a SIR-like (Susceptible–Infected–Recovered) compartmental model that includes vaccination. This model can be described by the following differential equation (for simplicity, only the vaccinated compartment V is reported):

where the variation in the number of immunized subjects V over time is proportional to V itself, multiplied by a rate of waning immunity w. The solution of this equation is an exponential decay for the protected population. This is distinct from the exponential titer evolution model (Eq. 2), which does not generally correspond to an exponentially declining protected population.

To allow a direct comparison between models, all data were computationally generated for the same model population of 10,000 vaccinated individuals, with two constraints: i) following vaccination, 100% of the population is initially protected, and ii) after an arbitrary time of 500 time units (which could represent days, weeks, etc.), the proportion of the protected population declines to 50%. Specific details regarding the three models are provided below.

The population was composed of 10,000 immunized subjects, with 50% no longer protected after 500 time units.

For the two evolution models, 10,000 antibody titers were generated at time tmin = 1 using Python’s random generator for the lognormal distribution (parameters: mean = 5.03 and sigma = 1). The protective threshold was set to log(4). Given the randomly generated titers, at time tmin 9998 subjects were over threshold (99.98%). Titers for t > tmin were then generated using Eqs. 1 and 2, with the titers at tmin serving as the intercept and computationally solving for the slope parameter at which the proportion of titers over a protective threshold of log(4) was exactly 50% at time t = 500 units (Table S1). The solutions found for the rate and exponent parameters were, respectively for exponential and power law evolution models, c = 0.007311 and b = 0.5882.

For the compartmental model, two compartments were used to represent the number of unprotected subjects U and the number of protected subjects V, with U + V = 10,000 at any given time. The initial number of protected subjects was set to V(tmin) = 9,998, to match the two evolution models. The waning parameter w was determined such that the proportion of protected individuals at time 500 would decrease to 50% (i.e., V(t = 500) = 5000). The resulting waning parameter was w = 0.001388.