Using vaccine Immunostimulation/Immunodynamic modelling methods to inform vaccine dose decision-making

Unlike drug dose optimisation, mathematical modelling has not been applied to vaccine dose finding. We applied a novel Immunostimulation/Immunodynamic mathematical modelling framework to translate multi-dose TB vaccine immune responses from mice, to predict most immunogenic dose in humans. Data were previously collected on IFN-γ secreting CD4+ T cells over time for novel TB vaccines H56 and H1 adjuvanted with IC31 in mice (1 dose groups (0.1–1.5 and 15 μg H56 + IC31), 45 mice) and humans (1 dose (50 μg H56/H1 + IC31), 18 humans). A two-compartment mathematical model, describing the dynamics of the post-vaccination IFN-γ T cell response, was fitted to mouse and human data, separately, using nonlinear mixed effects methods. We used these fitted models and a vaccine dose allometric scaling assumption, to predict the most immunogenic human dose. Based on the changes in model parameters by mouse H56 + IC31 dose and by varying the H56 dose allometric scaling factor between mouse and humans, we established that, at a late time point (224 days) doses of 0.8–8 μg H56 + IC31 in humans may be the most immunogenic. A 0.8–8 μg of H-series TB vaccines in humans, may be as, or more, immunogenic, as larger doses. The Immunostimulation/Immunodynamic mathematical modelling framework is a novel, and potentially revolutionary tool, to predict most immunogenic vaccine doses, and accelerate vaccine development.


Supplementary material for: Using vaccine Immunostimulation/Immunodynamic modelling methods to inform vaccine dose decision-making Additional Results
Analysis 1: Fitting the IS/ID model to the mouse data stratified by dose group and the human data As described, the LRT was used to establish the mouse covariate model in analysis 1i. The selected model for analysis 1i was the one which satisfied the LRT against the pooled model (see supplementary methods for the model fit to the pooled mouse data, table S6), had all estimated model parameters RSE <30%, and had the lowest -2LL.  Table S1 shows that the best covariate model is when dose group was indexed on model parameter bTEM with all model parameter standard deviations fixed to 0.5 (highlighted)(allowing the standard deviations to be estimated led to RSE of one or more parameters >30%).

Diagnostic Plots
The VPC plot, prediction distribution and observed versus predicted response plots can be found in Figures S1-S3.
The VPC shows that for each dosing group (low, middle and high), the model predicts the data well ( Figure S1), although with less data per group the VPC is not as well defined. This is due    Table S1, estimated parameters in Table 1 Table S1, estimated parameters in Table 1). The blue points represent the data. The bands represent the 25 th to 75 th percentiles of the theoretical predictions using the estimated population parameters and associated variation for analysis 1i ( Table 1). The black line shows the median total cell response prediction. Note, Y-axis not on the same scale.
6 Figure S3. Mouse observed data versus model predicted IFN-γ responses stratified by dose group (dose group indexed on parameter βTEM, see Table S1, estimated parameters in Table 1).
The results of the model fit to the human data can be found in Table 1. As one of the parameters RSE =30% the effects of estimating model parameter standard deviations (to account for BSV) was not tested as it was clear there was not enough data to estimate further parameters.

Diagnostic Plots
The VPC plot, model prediction distribution plot and the observed versus predicted (for the population and individual participants) for the pooled human model can be found in Figures The VPC plot shows that the simulated model cover the data well except for an underestimation of the median response at the latest time point ( Figure S4). This may be due to the fact that there is less data at the latest time point (these are only the H1 responses, not the pooled H56 and H1 responses). Again, due to the small sample size, this VPC plot does not summarise all responses, either observed data (green line) and model simulations (blue and orange regions) for all times points which is why the green line, blue and orange regions do not reflect the shape of the model prediction in Figure 2D of the main paper. Similarly, this is not a reflection of the fitting of the model, but an artefact of the default settings for the VPC plot in Monolix, where model predictions for small sample sizes are misrepresented.
However, the expected profile from the IS/ID model can be seen better in the model prediction distribution plot, which suggest the percentiles of the data are adequately covered ( Figure S5) despite widely variable responses over time in the human data set. Figure 2D in the main paper shows how the model predictions follow the trend of this variable data.
However, similar to the mouse pooled model, as all parameter standard deviations are fixed at 0.5, this may be underestimating the responses in some cases (although the observed versus predicted individual responses suggests the model is a good fit ( Figure S6)). The individual plots for each human participant can be found in Figure S7 and S8.   Table 1). The black line shows the median total cell response prediction  Analysis 2: Use fitted mathematical models in analysis 1, and a vaccine dose allometric scaling assumption, to predict the human immune response dynamics and predict the most immunogenic dose in humans In analysis 2, the estimated model parameters identified for the dose groups in mice and for the one dose in humans (analysis 1) were used to predict the IFN-γ response in humans for a range of doses. The steps to achieve this are outlined in the methods.
We found in analysis 1i, that the dose-dependent parameter was bTEM ( To achieve this step, we estimated three values for bTEM from the mouse data stratified by dose (for the low, middle and high dose group) by fitting the IS/ID model to the empirical mouse data using NLMEM (analysis 1i, Table 1). For step one we assume the average dose value for the low dose group (average of 0.1, 0.5, 1). We also included the bTEM value for the control mouse data (mice who received no H56) which we assumed was low as the IFN-γ response for the zero group was flat ( Figure S11). We verified this assumption by fitting the model to the zero dose data keeping all parameters except bTEM fixed to the estimated population bTEM values in Table 1  For the remaining steps, see Table S2.

As the current (antigen) dose allometric scaling factor between mouse and humans
for the H-series vaccines is assumed to be ten [1][2][3]  Step 1 Step 2 Step 3 Step 4 Step 5 For the remaining scaling factors, which for the H-series could potentially be between one and ten (personal communication, T Evans), we chose doses that would produce approximate values from 1 to 10. For example, for a scaling factor of 9, the dose 50/9=5.556 was predicted in step 1. For the remaining scaling factors (1 to ten), steps 2-5 were repeated. These calculations are not included here.
The human predicted dose response curves are in Figure 3 for Scaling factor 1, 5 and 10. For the dose response curves for the remaining scaling factors, see Figure S10.

Model Assumptions
Key model assumptions from the IS/ID model are outlined in Table S3.

Baseline responses were fixed at the median value
In this model, the initial values for the Transitional Effector Memory cells (TEM0) were not estimated. This is due to the fact that all mice IFN-γ responses at baseline were based on measurements from one unvaccinated mouse and therefore were all zero. As all human participants in the clinical trials were previously BCG vaccinated and no other human covariates were considered that could impact on a baseline response, the baseline responses were fixed to the median value. This also aided in avoiding over parameterisation compared to the small sample size of the human data.

Central Memory (CM) cells do not die
The central memory cell population is assumed to be maintained be a constant turnover, so we assumed the death rate could be omitted from the both the human and mouse model 4 . Although there is evidence to suggest CD4+ longterm memory cells turnover may diminish with time 5,6 , we assumed this does not affect the time frame of the model.

Introducing a death rate of memory cells would result
in a decline of the long-term responses.

transition, b CM
In the model, after revaccination, the CM cells replicated at a fixed rate for time t, which was estimated in the model fitting stage. Only after replication had occurred, the cells transitioned back to TEM cell type at a rate bCM, which was assumed to be fast. Although this may be a simplification of the host immune response dynamics, it was necessary to assume as we did not have information on bCM. We therefore considered the transition of CM cells to TEM cells as a result of revaccination to be a proliferation followed by a "burst" as opposed to a slower gradual transition (where proliferation and transition occur simultaneously). We believe this assumption is justified as the purpose of CM cells are to mount an immune (in our case, IFN-γ) response faster than a primary response as a result of re-exposure to the antigen (revaccination) 7 and a "burst" response is an effective method to represent this dynamic.

IFN-γ responses are not scaled to host body size
The ELISPOT assay readout is conventionally measured per million cells in all species and we considered the model to represent a systemic response regardless of host blood volume, it was not necessary to scale the ELISPOT readout to reflect body size. As our focus was on translating the change in dynamics due to change in dose between mouse and human, therefore this scaling the ELISPOT readout was not essential.

CD4+ T cell stimulation greatly simplified
The immune response to vaccination is a complex network of cells and cytokines behaving nonlinearly over time. In the Th1 response to Mtb. infection (or vaccination), innate and adaptive cells interact to optimise and maintain a protective response 8 . Very simply, cytokines secreted by innate cells after infection or vaccination, such as IL-12, work to stimulate adaptive cells to produce IFN-γ that both encourages innate cells to phagocytose bacteria and produce more IL-12 9,10 . As such, a feedback stimulation loop is established. In addition, to avoid an over-inflammatory response (which is harmful to the host) cytokines such as IL-10 are produced to regulate and dampen the immune response 11 . In the model, function δ is used to represent the delay of T cell initiation due to processes such as antigen processing and presentation and the decline of T cell responses due to depreciation of the required stimulation (creating a "n-shaped" curve). However, δ neglects the influence of stimulation amplification as a result of cytokine feedback loops, amongst other co-stimulation factors. As such, δ is a generalization of the complex networks required to protect against infection or vaccination and may not be as prolonged as required to generate a response to vaccination.
If data were available on IL-12 or other cytokines believed to be important to an immune response to BCG, It is possible that δ could be modelled as a parallel "innate response" compartmental model.
Incorporating such a model would provide insight into the innate cell mechanisms and thus strengthen the conclusions we draw on the T cell dynamics.

Transition and replication of transitional effector cells happens in Lymph node
before entering the blood

Data
Mouse IFN-γ ELISPOT data The methods and materials used to generate the mouse IFN-γ response data following H56+IC31 vaccination are outlined below. These methods are published in 12 .

Ethics Statement
All animal work was carried out in accordance with the Animals (Scientific Procedures) Act 1986 under a license granted by the UK Home Office (PPL 70/8043), and approved by the LSHTM Animal Welfare and Ethics Review Body.

Animals
Female CB6F1 mice were acquired from Charles River UK at 6-8 weeks of age. Animals were housed in specific pathogen-free individually vented cages, were fed ad libitum, and were allowed to acclimatize for at least 5 days before the start of any experimental procedure.

Vaccination
The  Figure S11 shows the ELISPOT results using the 24 hour incubation time for each dosing group.
Each coloured dot represents the responses of one mouse, the lines indicate the median responses. Human IFN-γ ELISPOT data Table S4 summarizes the two H-series trials from which the human ELISPOT data was taken. Figure S12 shows the individual IFN-γ responses (measured using ELISPOT assay) over time for both trials and the pooled median response across both trials.

H1+IC31 phase I trial 19 : ClinicalTrials.gov no NCT00929396
The ELISPOT methods for the H1+IC31 clinical trial are outlined in 19 . We summarise the methods below.
Frozen cells were pre-stimulated for 16-18 hours, followed by 24 hours in the ELISPOT plate.
1x10 6 thawed cells/well were stimulated in 24 well plates with H1 antigens (Ag85B and ESAT-6 proteins) as well as PPD, separate peptide pools and positive and negative controls (see 19 ).
All samples were assayed in triplicate. Incubation was done overnight in a fully humidified incubator at 37 •C, 5% CO2. Subsequently, cells were resuspended and divided over 3 wells BCG vaccination status was determined by tuberculin-skin-test (TST), whereby a reaction range 6-15mm or any documented value between 6 and 15mm on medical file in the past, indicated the participant was BCG vaccinated. To determine LTBI status, a QuantiFERON®-TB Gold In Tube test and a 6-day lymphocyte stimulation test (as described in 20 ) in addition to chest X-rays, were conducted at screening. HIV status was determined by reviewing recorded medical history and conducting standard blood tests.

Mathematical vaccine Immunostimulation/Immunodynamic (IS/ID) Model
The conceptual mechanisms of the IS/ID model for the IFN-γ secreting CD4 T cells as a result of H56+IC31 vaccination can be found in Figure 1. The mathematical model used in Monolix to estimate the parameters represents these exact mechanism, however we separated the compartments into TEM1 and CM1 and TEM2 and CM2 corresponding to TEM and CM after primary vaccination and TEM and CM after revaccination. Here, once primary vaccination occurs, cells are recruited at rate ! "#$%&' into TEM1 and these cells then either die (at rate ( )*+ ) or transition to CM1 at rate -)*+ . At time of revaccination two process are initiated simultaneously. The first is recruitment into TEM2 at rate ! "#$%&.%/011 (this is the same recruitment, but adjusted for revaccination time; parameter values are the same as in the case of primary vaccination). The second is the replication of CM1 cells. After t days, the CM1 cells stop replicating and transition to TEM2 at rate -2+ which is fixed to an arbitrarily high value so that this process is instantaneous. TEM2 cells then either die (at rate ( )*+ ) or transition to CM2 at rate -)*+ . The output of the model is the sum of all TEM and CM compartments over time (see main methods for justification of this).
In summary, this model essentially separates the dynamics of the cells corresponding to which vaccination time they follow, but retains the overall conceptual mechanisms outlined in  Figure 1 was adopted. In the conceptual schema in Figure 1, the only way to eliminate this from happening was to apply a condition on -2+ to be 0 once 99% of CM cells had left the CM compartment after replication. This would have added complexity to the Monolix code we felt was unnecessary. We felt the representation of the model dynamics in Figure 1 was easier for a reader to understand conceptually than the model of the separated compartments.
As such, the equations for the model are as follows: A schema of the model is below. Where dotted arrows correspond to rates activated at time of (or following in the case of -2+ ) revaccination.

Statistical (NLMEM) model
In To conduct NLMEM to estimate the IS/ID model parameters for the H56 IFN-γ response data for mice and humans, the following is required: 1. Statistical models to account for individual responses.
To account for the intra and inter-variation, the NLMEM framework statistical model requires two components:

c. Observed versus predicted data plots
In these plots, the observed data is plotted against the predicted value on a population and individual level. Plots where the observed and predicted values are similar should show a diagonal distribution along a line of unity (i.e. if the prediction is the same as observed value that point will lie on the diagonal line extending from the origin of the plot).

Monolix code
Alongside the data, which is outlined in figures S11 and S12, the statistical model outlined in the methods of the main paper, the monolix code used to fit the IS/ID model to the data can be found below:

Justification for fixing the random effects due to between subject variability (BSV)
Given the sample size of the data and what we aimed to achieve (using the model to describe the IFN-γ dynamics by dose), we prioritised estimating the "average" response for as many subpopulations (dose groupings) as possible (the population typical model parameter values for each dose group), over estimating variation in response for those subpopulations (the BSV of the model parameters).
As such we implemented the following tasks until power to achieve well estimated parameters was lost due to sample size:

Testing the structural model for parameter d and fit of the IS/ID model to the pooled mouse data
Three different mathematical forms were used to represent parameter d:  As outlined in the main paper: the Gaussian equation is as follows (Figure 1, Table S5):  Table S6.    The model predicted IFN-γ responses for the parameter set in table S6 are plotted in Figure   S13. The visual predictive check (VPC) plot showed the model predictions represented the median pooled data well ( Figure S14).  Table S6) is plotted against the median data (blue line). The orange and purple dashed lines are the model prediction (total cells) of the 75th and 25th percentiles of the data, a result of the variation in the estimated parameters (standard deviation fixed to 0.5 for all parameters (Table S6)).
The VPC plot shows that the simulated model predictions cover the data well and there are little red areas (red areas indicate the simulated model predictions did not adequately cover the observed data) ( Figure S14). The red areas in the early response stages may be due to variable responses at this stage. The red area for the 25 th percentile prediction indicates the model is under predicting the data. This could be due to the 0.5 value constraint placed on the standard deviation of the parameters which limit the degree to which the predictions can

-2*LogLikelihood sensitivity analysis on the pooled mouse data
Likelihood sensitivity analysis was conducted by varying one parameter whilst holding the remaining parameters fixed to the values in Table S6 (the fit of the model to the pooled mouse data). The sensitivity range tested was [parameter value/10, parameter value*10] (except for parameter b, which was bound by day 15 to ensure the peak of the recruitment parameter, delta, occurred before revaccination). For the varied parameter, the standard deviation was fixed at 0.5 as a result the lowest -2LL values (those when the parameter is at its best value; those recorded in table S6) are different across the parameters in the sensitivity analysis. The results are plotted and discussed in the figure   S15 and table S7 below.
In all cases, the slope of the -2LL sensitivity curves are steeper for the values of the parameters smaller than that of the optimal value (the parameter value where the -2LL is lowest), suggesting that model parameters were fix too low, the description of the data would be likely be worse than if fixed too high. The likelihood is most sensitive to parameter c in the given parameter range explored (the difference between the maximum and minimum -2LL value (optimal parameter value) is largest for parameter value both below and above the optimal value). Parameter c is the variance in the Gaussian equation that describes the TEM recruitment parameter delta. The effect of increasing c on the predicted IFN-γ response over time is a an earlier and more sustained TEM cell recruitment which acts to increase the overall magnitude of response. As such, it is possible that the model is very sensitive to this parameter, so very low and very high values of this parameter will affect the likelihood against all data points (for all times points).
a The likelihood is sensitive to parameter a in the given range. Parameter a scales the shape of the recruitment rate parameter, delta. An increase in a increases the magnitude of the TEM cell population and therefore the response over all, although unlike the parameter c, a does not regulate the timings of the recruitment only the absolute number of cells. Although similar to parameter c, for the given parameter range in the sensitivity analysis, there is less of a likelihood penalty at larger values of a then for the equivalent change in c. b We assumed the parameter b could only take values less than 15, as it was an assumption in the model that the TEM recruitment would be delayed after time of vaccination, but occur before revaccination. Potentially, very early (small) values for this parameter suggest that the response clears too quickly (for all other parameters fixed), which is the likelihood is noticeably worse for smaller values than for larger values (up to day 15).

Tau
The parameter tau represents how long the CM replicate for before they transition back to effector type. For the model output, this dictates how high the height of the peak after vaccination. For all other parameters fixed, this means, too high values of tau will lead to an overall response after revaccination that is too high, the converse is true for too small values of tau. The model is least sensitive to the value of tau. This is potentially as it is the only parameter in the model that only effects the response after revaccination and therefore, only the fit of the model to three time points.

Btem
The result of changing Btem -the transition rate of TEM to CM cells -is the evident in the later time points, before revaccination and in late time points after revaccination. A higher value of Btem suggests more TEM cells are transitioning to CM type, which then contributes to the overall output of the model (as CM cells do not die). This is most apparent for the very late time points of the model as the model output will be 100% CM cells (given a fixed death rate of TEM cells). A lower value of Btem will mean TEM cells may die out before they can transition, lowering the magnitude of response once the initial recruitment of TEM cells has waivered. This is important, as if the CM cells are minimal just before revaccination time, with a fixed Tau, the response may not be boosted high enough to fit the data after revaccination. This may be the reason why the -2LL sensitivity slope increases fast for smaller values of BTem and to higher -2LL values.