An ODE-based mixed modelling approach for B- and T-cell dynamics induced by Varicella-Zoster Virus vaccines in adults shows higher T-cell proliferation with Shingrix compared to Varilrix

Clinical trials covering the immunogenicity of a vaccine aim to study the longitudinal dynamics of certain immune cells after vaccination. The corresponding immunogenicity datasets are mainly analyzed by the use of statistical (mixed effects) models. This paper proposes the use of mathematical ordinary differential equation (ODE) models, combined with a mixed effects approach. ODE models are capable of translating underlying immunological post vaccination processes into mathematical formulas thereby enabling a testable data analysis. Mixed models include both population-averaged parameters (fixed effects) and individual-specific parameters (random effects) for dealing with inter-and intra-individual variability, respectively. This paper models B-cell and T-cell datasets of a phase I/II, open-label, randomized, parallel-group study in which the immunogenicity of a new Herpes Zoster vaccine (Shingrix) is compared with the original Varicella Zoster Virus vaccine (Varilrix). Since few significant correlations were assessed between the B-cell datasets and T-cell datasets, each dataset was modeled separately. By following a general approach to both the formulation of several different models and the procedure of selecting the most suitable model, we were able propose a mathematical ODE mixed-effects model for each dataset. As such, the use of ODE-based mixed effects models offers a suitable framework for handling longitudinal vaccine immunogenicity data. Moreover, it is possible to test differences in immunological processes between the two vaccines. We found that the Shingrix vaccination schedule led to a more pronounced proliferation of T-cells, without a difference in T-cell decay rate compared to the Varilrix vaccination schedule. Author summary Upon vaccination, B-cells and T-cells are activated to induce an immune response against the vaccine antigen at hand. In this paper, we study and compare the longitudinal dynamics of the specific immune response based on a vaccine trial in which the immunogenicity of a new Herpes Zoster vaccine (Shingrix) is compared with the original Varicella Zoster Virus vaccine (Varilrix). We combine the use of ordinary differential equations (ODEs), i.e. mathematical models which are used to describe the dynamics of the immune response, with advanced regression analyses enabling us to infer the model parameters describing these dynamics. The resulting ODE-based mixed effects models enable describing the immune response dynamics allowing for both inter-and intra-individual variability; comparing the dynamics induced by the two vaccines and studying the B-and T-cell interactions. We found a more pronounced proliferation of T-cells for the Shingrix vaccination schedule as compared to the Varilrix vaccination schedule. The proposed methodology offers a suitable framework for better understanding the immunogenicity of vaccines.

T-cells have a direct cytotoxic function and can target host cells that are infected by a 10 pathogen. 11 The vaccine-induced B-cells and T-cells are hypothesized to be capable of preventing 12 or minimizing the morbidity related to the infectious disease against which the vaccine 13 is targeted. Vaccine immunogenicity trials aim to study the longitudinal dynamics of 14 the specific immune response following vaccination. These trials can range from several 15 months to several decades. The quantitative analysis of longitudinal immune response 16 data has evolved from between-group and time point comparisons to statistical 17 regression analyses [1,2,3]. Current state-of-the art statistical analyses of longitudinal 18 data consist of a mixed effects model approach in which a separation is made between 19 population-averaged parameters (so called fixed effects) and individual-specific 20 parameters (so called random effects). More recently, Andraud et al. [4] and Le et al. [5] published the first papers in which the mixed effects modeling approach was 22 combined with the use of ordinary differential equations (ODE), thereby more closely 23 resembling immune response dynamics post vaccination. Whereas [4] focused on the 24 long term dynamics following vaccination, [5] focused on the short term dynamics 25 following vaccination. 26 ODE-based mathematical models are capable of translating the underlying 27 immunological/biological theory into a testable data analysis. Moreover, the 28 combination with mixed effects modeling offers a methodology capable of dealing with 29 inter-and intra-individual variability. As such, the use of ODE-based mixed effects 30 models offers a suitable framework for handling longitudinal vaccine immunogenicity 31 data. 32 In this paper, we set out to use ODE-based mixed effects models to study B-cell and 33 T-cell dynamics following varicella-zoster virus (VZV) vaccinations in VZV-immune adults. In particular, this framework will allow us to disentangle the immunogenic 35 differences between two different VZV-specific vaccines. 36 In the Materials and methods section, we first present the immunogenicity data from 37 two VZV vaccine studies consisting of B-cells and CD4+ T-cells of participants at 38 different time points. We then describe the differential equations, the ODE and the 39 ODE-based mixed effects models used to describe the immune response dynamics 40 within each individual as well as the associated model selection procedures. In the 41 Results section we apply the above methods and select a suitable model for each 42 dataset. Next, we compare the results of the two VZV vaccines, using a group-related 43 effect on a chosen parameter. Correlations in and between the datasets are also 44 explored. Finally, we review our findings in the Discussion section. adults, two groups of young adults (18-30 years) were vaccinated with two vaccine doses 52 two months apart. The first group (GROUP 1; sample size: n 1 = 10) received one dose 53 of HZ/su and one dose of OKA concomitantly at month 0 and month 2 (i.e. four doses 54 in total), whereas the second group (GROUP 2; n 2 = 10) received a dose of HZ/su both 55 times (i.e. two doses in total).

56
After vaccine safety was confirmed, three groups of older adults (50-70 years) were 57 vaccinated two months apart, one group (GROUP 3; n 3 = 45) received twice a single 58 dose of HZ/su, the second (GROUP 4; n 4 = 45) twice a single dose of OKA and the last 59 (GROUP 5; n 5 = 45) twice two concomitant doses of HZ/su and OKA. So, all in all, 155 60 participants were divided over these 5 groups. The properties of each group are 61 summarized in Table 1.

62
Safety and immunogenicity were assessed for all groups up to 12 months 63 post-vaccination in the original study. In order to obtain long-term immunogenicity 64 data on the newly proposed HZ/su vaccine, 23 individuals from the groups solely 65 receiving HZ/su (i.e. GROUPS 2 and 3) were assessed up to 42 months post-vaccination 66 in the extension studies: BIO109671 and BIO109674. We refer to [6] for a more in 67 depth description of the design and results of these studies. are plotted in Fig. 1 for the Varilrix-specific B-cells. considering the second vaccination at month 2.

83
We observe an increase in B-cells, or antibody-secreting cells (ASC), after vaccination 84 at time t = 0 months. At time t = 2 months, the subjects were re-vaccinated, but no 85 data were collected at that time point. Fig 1 shows only the time points for which data 86 were available ( t = 0, 1 and 12 months). Since it is reasonable to expect a (higher) peak 87 in the data after the second dose at t = 2 months, we will assume a time period decreases. The data plots of gE-specific ASC show a similar pattern (see Fig S1).   Given that T-cell data were collected at more time points than B-cell data, we now 103 observe a second peak in the group-specific data plots, as is expected given the vaccine 104 administrations at month 2. Therefore we will use two time periods [0, h 1 ] and ]2, h 2 ] 105 (with 0 < h 1 < 2 < h 2 ) during which the level of T-cells first increases and then 106 decreases. In case only one peak is observed, we will assume h 1 = 2 < h 2 .  Mathematical methods 109 We used systems of (nonlinear) ODEs to model the B-cell and T-cell dynamics. We 110 applied a systematic approach to fit and compare several models in order to obtain the 111 models that best describe the available data, while providing sufficient biological in S1 Appendix. In the following subsections we provide the basic rationale of these 114 ODEs for both B-cell and T-cell dynamics, respectively. 115 B-cell dynamics models 116 We describe the dynamics of the antibody secreting B-cells (ASC), using the following 117 ODE: where ASC 0 = ASC(0) denotes the initial number of antibody secreting cells at time=0 119 (months) and f 1 and f 2 are smooth functions of the number of ASC at time t (months), 120 describing the change in the number of ASC due to activation and decay of ASC. We  Essentially, the former have a much shorter life span than the latter. Plasma cells 134 remain in the bone marrow for a long time period [4,7]. Models B3 and B4 incorporate 135 this distinction by including different equations in the ODE system for short living ASC 136 (SASC) and long living ASC (LASC). The dynamics of the SASC in models B3 and B4 137 are similar to those in model B1, but at time 0 months, no SASC are present in models 138 B3 and B4. LASC however, are present in models B3 and B4 at time 0 months.

139
To distinguish between models B3 and B4, the dynamics of LASC are considered.  The design of the T-cell models follows a similar procedure as that of the B-cell models. 148 The following ODE describes the basic dynamics of the stimulus-specific T-cell 149 population:

178
The dynamics of the short living T-cells are similarly described as in models T1 and 179 T2: a constant number of short living T-cells will be activated after each vaccination 180 (not necessarily the same number), while the decay of short living T-cells occurs at all 181 times at a constant decay rate. For models T3 and T4, the assumption is made that the 182 total number of long living T-cells remains constant over time. If we add the distinction 183 between models with equal and different functions f 1 and f 2 , we arrive at models T3 184 and T4, respectively. Adding proliferation rates of long living T-cells after each 185 vaccination, models T3 and T4 are extended to yield models T5 and T6. In order to 186 restrict the total number of parameters, the long living T cell proliferation rates are 187 assumed to be equal. Table 3 summarizes the parameters used in the functions of the 188 different models.  Although mathematical identifiability is guaranteed for the models presented in the 215 Mathematical methods section, the complexity of these models when combining them 216 with many random effects in view of the data limitations in terms of sampling times 217 and sample sizes resulted in non-convergence. Therefore, simplifying assumptions 218 needed to be made. One such simplifying assumption is presuming that the decay of B 219 or T-cells is identical for all individuals, implying that the random effect for that decay 220 parameter is omitted from the model.

221
For both the B-cell and T-cell data sets, the following procedure was used for 222 comparing and selecting the most suitable biologically plausible model to describe the 223 data. In a first step, a list of models was composed, consisting of models B1 to B4 for 224 the B-cell data and of models T1 to T6 for T-cell data, together with assumptions on 225 the parameters reflecting whether or not individual variation on these parameters is 226 present, i.e. whether or not random effects were included for the different parameters. 227 The model parameters were then estimated with the Monolix software.

228
Models with poor SAEM convergence, likely because of abundant model complexity, 229 were discarded. Next, the candidate models were compared using Akaike's Information 230 Criterion (AIC) and the model with lowest AIC value was selected as first candidate 231 model. Subsequently, a non-parametric bootstrap, using 1000 bootstrap re-samples, was 232 performed on the candidate model. Since a sequential approach based on the candidate 233 models with the lowest AIC values was used, the need to perform bootstraps for all 234 candidate models was avoided, in order to decrease the number of computations. It was 235 found that for a bootstrap, either 65%-77% of the samples had proper SAEM 236 convergence, or the proportion of bootstrap samples with proper convergence was less 237 than 15%. For this reason the criterion for good bootstrap convergence was defined as 238 having at least 65% of bootstrap samples with proper SAEM convergence. In case of 239 poor bootstrap convergence, the candidate model was rejected from the list of candidate 240 models.

241
Next, a sensitivity analysis on the bootstrap results of the (converging) candidate calculate the maximal information coefficient (MIC). This is a way to detect linear and 263 non-linear relations between variables, and can thereby be used to indicate whether a 264 linear relation is feasible by comparing it with the R-squared value [9]. In addition we 265 compute the Spearman correlation.

271
Model selection of B-cell datasets 272 We started by modeling the Varilrix-specific B-cell dataset, for which the model selection 273 procedure outlined in the Inference and model selection section was followed. The upper 274 part of Table 4 summarizes the differences between all models that were tested.    Consequently, we looked at models in which a distinction between SASC and LASC 286 was made. Model B3a assumed all parameters had random effects. In models B3b and 287 June 11, 2018 9/19

275
B3c, uASC and h, respectively, were set as fixed population parameter. Model B3b was 288 the only model that showed convergence, with an AIC value of 6541.

289
The last model examined was model B4, in which a proliferation rate for LASC was 290 added. Many assumptions on the parameters were made; the decay rate of ASC 291 (uASC), proliferation rate of LASC (aLASC) and activation period (h) were set as fixed 292 parameters in models B4b, B4c and B4d, respectively. Combinations of these fixed 293 parameters were considered as well in models B4f, B4g and B4h. Apart from this, we 294 also looked at different group specific parameters, not only the proliferation rate of 295 SASC (aSASC) was considered, but the proliferation rate of long living B-cells (aLASC) 296 as well. Model B4g was the only model that accounted for aSASC and aLASC and still 297 converged. In this model, both aLASC and uASC were set as fixed population 298 parameters, and aSASC was considered to be a group specific parameter. This model 299 had the lowest AIC value of 6524 among all aforementioned models, and was selected as 300 first candidate model.

301
The bootstrap that subsequently was performed did not converge for model B4g, 302 and it was therefore in the end rejected. Model B3b was selected as next-candidate   Table 5. Parameter estimates and corresponding 95% confidence intervals (CI) of final model B1a. No distinction between SASC and LASC is presumed. ASC 0 = ASC(0) denotes the initial number of ASC. The proliferation of ASC is constant in time period (0, h), at rate aASC and assumed to be group-specific with effects β i (i = 1, 2, 3, 5). Decay of ASC happens at rate uSASC.  Model T4 assumed different activation rates of T-cells after each vaccination. When 339 all parameters had random effects, and both a 1 T and a 2 T were group specific, SAEM 340 convergence was not reached (model T4a). When only a group specific effect on a 2 T 341 was assumed (model T4b), convergence was achieved resulting in an AIC value of 342 11,615. Subsequently, uT was set as fixed parameter in models T4c and T4d, again with 343 group specific effects on both activation rates (T4c) and on a 2 T only (T4d), 344 respectively. Both models showed SAEM convergence with an AIC value of 11,646 and 345 a lower AIC value of 11,626, respectively.

346
When assuming LT activation according to a constant proliferation rate (equal after 347 each vaccination in order to limit the number of parameters to be estimated), models 348 T5 and T6 were reached. In models T5, the activation rates of ST were presumed equal 349 after each vaccination. Together with the assumption that all parameters were random, 350 and aST was a group specific parameter, this leaded to model T5a, where an AIC value 351 of 11,631 was found. We note that setting aLT as a group specific parameter was tested 352 as well, but none of these models (including the following) showed convergence and thus 353 were omitted from Table 7 T4b  no  T6f  no  T4d  no  T6d  no  T5a  no  T5b  no  T3a  66 %  no  T3a' 67 % no Above: Overview of the different models used to model the Varilrix-specific T-cell data. First column: model identifier. Second column: parameters selected as fixed population parameter. Third column: parameter which is chosen to be group-specific.  considered. When uT was fixed, and both activation rates were group specific, 363 convergence was not reached (model T6c). With a 2 ST being group specific (model 364 T6d), convergence was obtained with an AIC value of 11,630. In case of setting aLT as 365 a fixed parameter, similar results were found; assuming both activation rates to be 366 group specific did not lead to convergence, but assuming only a 2 ST was group specific, 367 did, with a slightly lower AIC value equal to 11,624. The last scenario assumed both 368 aLT and uT were fixed population parameters, though no convergence was obtained.

369
Since model T4b was the model with lowest AIC (11,615), it was selected as first 370 candidate model. However, model T4b was subsequently rejected as a 1000 sample 371 bootstrap failed to converge. Likewise models T6f, T4d, T6d, T5a and T5b did not have 372 proper bootstrap convergence. Next, model T3a was selected as candidate model and 373 showed bootstrap convergence; 66% of the bootstrap samples reached SAEM 374 convergence. As before, a search for frequently deviant profiles in the converging and 375 non-converging bootstrap datasets was performed, but no such profile was identified.

376
Taking into account that the assumption that a 1 ST = a 2 ST = aST might not be a 377 realistic assumption in a model that described a real life cellular process, a difference in 378 proliferation rates after each vaccination was inserted, assuming a 2 ST was proportional 379 to a 1 ST : a 2 ST = k × a 1 ST . Adding this parameter to the pool of parameters to be 380 estimated in the SAEM algorithm, did not yield a converging model. In order to limit 381 the number of parameters to be estimated by the SAEM-algorithm, different values of k 382 ware set fixed and for each subsequent model, AIC values were compared. A model with 383 k = 0.15 showed the lowest AIC value (11,623), which we named model T3a', and a  Table 8.   Table 9. Model T1a" 396 does not differentiate between ST and LT and assumes a constant activation of T-cells 397 in time periods (0, h 1 ) and (2, h 2 ) with a 2 ST = 0.66 × a 1 ST . All parameters were 398 assumed to have random effects with the activation rate being a group-specific 399 parameter. Moreover, individuals 89 and 149 were left out of this dataset since they had 400 too much influence on model convergence. Table 9. gE-specific T-cell results.   Table 10. In view of the sample size of groups 1 and 2 (and age), we were 419 mainly interested in β 3 and β 5 . A β i higher than zero indicates a higher activation rate 420 of cells in the groups receiving the HZ/su vaccine, compared to the activation rate in 421 the reference group which received the OKA vaccine. In the case of Varilrix-and 422 gE-specific T-cells, both groups 3 and 5 showed a significant higher activation rate 423 (p < 0.05). The activation rate of gE-specific B-cells also was significantly higher 424 compared to the reference group. Varilrix-specific B-cells also seemed to have a higher 425 proliferation rate, though in the case of groups receiving solely the HZ/su vaccine, not 426 significantly so (p > 0.05).

438
To further examine this, we calculated the Spearman correlation and the maximal 439 information coefficient (MIC) as a way to assess and measure (non)linear relationships 440 between datasets.

441
The Spearman correlation between Varilrix-specific T 01 and B 01 was -0.0274 442 (p = 0.7914), which rejected the hypothesis that increases in Varilrix-specific T-cells 443 were associated to increases in Varilrix-specific B-cells. A non-significant MIC of 0.2161 444 (5% significance level) confirmed this result. 445 We also assessed MIC and Spearman correlations on the same datasets per subgroup, 446 reaching the same conclusion (Table 11).

MIC coefficients, Spearman correlation and corresponding p-values between increase in
Varilrix-specific B-cells (B01) and increase in Varilrix-specific T-cells (T 01), calculated for groups 3, 4 and 5. As the sample sizes of groups 1 and 2 were too small (n 1 = 4 and n 2 = 6), those groups were omitted from the analysis.
Potential associations between increases in gE-specific T-cells and B-cells were 448 investigated. The scatterplot of T 01 plotted against B 01 is found in Fig S3. 449 The Spearman correlation of 0.3833 (p = 1.0647e −04 ) suggested there was indeed an 450 association. The MIC score was calculated as 0.4107, which was significant (p < 0.001). 451 In case of a linear relation, the R-square is expected to be close to this MIC score. As 452 the R-square was equal to 0.05593, we could exclude a linear relationship.

453
As before, we also studied the relations between T 01 and B 01 per subgroup (see Table 454 12). The Spearman p-values showed that only in group 3 the Spearman correlation 455 could be considered significant, together with MIC, implying a nonlinear relationship. MIC coefficients, Spearman correlation and corresponding p-values between increase in gE-specific B-cells (B01) and increase in gE-specific T-cells (T 01), calculated for groups 3, 4 and 5. As the sample sizes of groups 1 and 2 were too small (n 1 = 4 and n 2 = 6), those groups were omitted from the analysis.

456
Correlations between the initial number of B-cells B(0), the initial number of T-cells 457 T(0), the number of B-cells at month 1 B(1) and the number of T-cells at month 1 T(1) 458 have also been investigated. 459 We started with the Varilrix-specific B-cell and T-cell data. Again, we only included 460 individuals for whom we had data points B(0), B(1), T(0) and T(1). Spearman individuals and then split by group. We observed that some correlations seemed to be 477 higher compared to the Varilrix-specific data, however, when examining the Spearman 478 matrices by group, the values between the different groups seemed to vary widely. For 479 this reason, we did not find decisive evidence to include the number of gE-specific  In this study we used a nonlinear mixed modeling approach using ordinary differential 486 equations to describe B-cell and T-cell dynamics in adults following 2-dose vaccination 487 against VZV by means of the novel subunit VZV gE vaccine (Shingrix, GSK) and 488 live-attenuated VZV (Varilrix, GSK). Our study was motivated by the difficulties in 489 attributing differences between vaccines and vaccination schedules to underlying 490 immunological processes when using "classical" statistical techniques that do not take 491 into account the underlying immunology. Recently, Andraud et al. [4] and Le et al. [5] 492 showed that nonlinear ODE mixed modeling in the setting of vaccinations was capable 493 of giving estimates on several biological parameters.  Importantly, this way of modeling allowed us to directly compare several vaccination 511 schedule specific parameters. We found that the Shingrix vaccination schedules led to a 512 more pronounced proliferation of T-cells, however without a difference in T-cell decay 513 rate between vaccination schedules. This shows the benefit of using of mathematical 514 mixed models instead of performing a statistical analysis of the datasets: in the latter 515 case it is possible to prove significant differences between vaccines, however, it is not 516 determinable whether that increase is the result of either a higher proliferation of cells, 517 a lower decay (mainly in the case of a restricted number of data points), or a longer 518 time period (0, h) in which cells are activated. 519 We note that the adjuvant used for the Shingrix vaccine has been reported to be a 520 very potent adjuvant [10] and our modeling approach thus confirms the increased 521 proliferation of T-cells for the Shingrix vaccine. 522 We also assessed whether a correlation existed between the B-cell and T-cell counts, 523 but we did not find a significant association between the two immune response types.

524
This confirms previous findings concerning the glycoprotein-E adjuvant, part of the 525 AS01 B Adjuvant System family, in which it has been shown that this family has been 526 reported to show the lowest correlations between B-cells and T-cells of all families [11]. 527 During our modeling analyses we encountered several limitations. First, we noted 528 that given the limited sample size only models with moderate complexity could be