Characterising antibody kinetics from multiple influenza infection and vaccination events in ferrets

The strength and breadth of an individual’s antibody repertoire is an important predictor of their response to influenza infection or vaccination. Although progress has been made in understanding qualitatively how repeated exposures shape the antibody mediated immune response, quantitative understanding remains limited. We developed a set of mathematical models describing short-term antibody kinetics following influenza infection or vaccination and fit them to haemagglutination inhibition (HI) titres from 5 groups of ferrets which were exposed to different combinations of trivalent inactivated influenza vaccine (TIV with or without adjuvant), A/H3N2 priming inoculation and post-vaccination A/H1N1 inoculation. We fit models with various immunological mechanisms that have been empirically observed but have not previously been included in mathematical models of antibody landscapes, including: titre ceiling effects, antigenic seniority and exposure-type specific cross reactivity. Based on the parameter estimates of the best supported models, we describe a number of key immunological features. We found quantifiable differences in the degree of homologous and cross-reactive antibody boosting elicited by different exposure types. Infection and adjuvanted vaccination generally resulted in strong, broadly reactive responses whereas unadjuvanted vaccination resulted in a weak, narrow response. We found that the order of exposure mattered: priming with A/H3N2 improved subsequent vaccine response, and the second dose of adjuvanted vaccination resulted in substantially greater antibody boosting than the first. Either antigenic seniority or a titre ceiling effect were included in the two best fitting models, suggesting a role for a mechanism describing diminishing antibody boosting with repeated exposures. Although there was considerable uncertainty in our estimates of antibody waning parameters, our results suggest that both short and long term waning were present and would be identifiable with a larger set of experiments. These results highlight the potential use of repeat exposure animal models in revealing short-term, strain-specific immune dynamics of influenza.

Introduction Natural infection with influenza stimulates a complex and multifaceted immune response to neutralise and clear the infection [1]. The adaptive immune response is of particular interest for seasonal epidemic and pandemic preparedness, as responses from previous exposures provide some long-term protection against reinfection and disease via antibody and T-cell mediated immunity [2,3]. Focusing on the adaptive immune response is also advantageous because it can be (i) induced in advance of an epidemic through vaccination and (ii) measured and compared against correlates of protection to improve public health forecasting [4][5][6][7].
Influenza is an antigenically variable virus and undergoes continual antigenic drift, whereby mutations in immunodominant epitopes are selected by immunological pressure, allowing influenza lineages to escape population herd immunity [8][9][10]. This results in the continual loss of long-term immunity as antibodies effective against past strains fail to neutralise novel variants [11]. The current strategy for combating antigenic drift is to regularly update the seasonal influenza vaccine to better represent circulating strains, resulting in a competition between virus and vaccine formulation. Vaccines are more effective in some years than others due to factors such as: the antigenic match between the vaccine and circulating strain, and prior exposure histories of the population [12,13]. Consequently, there has been a recent push towards a universal influenza vaccination strategy, either through new vaccines or improved strength and breadth of immunity using existing technologies [14,15].
Whilst there is some cross-reactivity and cross-protection within influenza A virus subtypes and within influenza B virus lineages, humans experience numerous infections over their lives [16][17][18]. Each successive influenza exposure, which may be vaccination or infection, can strengthen the available repertoire of T and B cells which target epitopes on circulating and previously encountered strains [19,20]. In the humoural response, this occurs by boosting antibodies produced by pre-existing long-lived plasma cells and activated memory B cells (MBCs), and through generating a novel B cell response targeting unrecognised epitopes [21][22][23]. Given that individuals experience repeated infections and vaccinations from antigenically varied influenza viruses, interpreting the composition of an observed antibody response is confounded by the complex interaction of an individual's immunity from prior infections with the infecting virus [24][25][26][27].
A large body of experimental and observational work exists describing the contribution of infection histories to observed influenza susceptibility profiles, antibody landscapes and vaccination responses, often under the terms "original antigenic sin" or "antigenic seniority" [1,[28][29][30][31][32][33][34]. Furthermore, next generation assays to characterise antibody diversity and B cell identity have provided a detailed understanding of immune dynamics, including short term immune kinetics, duration of the humoural response, and immunodominance of different antigenic sites [35][36][37][38][39][40]. However, few studies have integrated these mechanisms into quantitative frameworks which can be used to explain and predict serological data from human populations, which often rely on simpler and less finely resolved assays [18,22,41,42].
Animal models, in particular ferrets, have been used to generate much of our understanding of influenza immunology due to opportunity for intensive observation and control [36,[43][44][45][46][47][48]. Here, we exploit the experimental flexibility and transparency of a ferret model to find evidence for and quantify multiple immunological mechanisms that may be important in characterising antibody landscapes generated from complex exposure histories, yet are observable using only routine antibody assays. Quantifying short term mechanisms in a ferret system might reveal patterns that could be used to improve the predictability and interpretation of human antibody landscapes following exposure [49].
We developed a mathematical model of antibody boosting and biphasic waning to describe antibody kinetics using previously published antibody titre data from a group of ferrets with varied but known exposure histories [43]. Previous models of antibody kinetics have often focused on the response to a single immunogen following one exposure [50][51][52]. Here, we take into account previously described immunological phenomena to describe cross-reactive antibody titres arising from varied exposure histories. These phenomena included exposure typespecific homologous and cross-reactive antibody kinetics, the role of priming on subsequent vaccination, titre-dependent antibody boosting (or titre ceiling effects) and reduced antibody boosting with each subsequent exposure (antigenic seniority) [32,[53][54][55][56][57][58][59]. By fitting models with various combinations of these mechanisms to haemagglutination inhibition (HI) titre data from ferrets, we sought to identify immunological mechanisms that are important in describing observed antibody profiles arising from multiple exposures. Parameter estimates from these model fits allowed us to quantify the impact of prior infection and adjuvant inclusion on antibody levels following vaccination and to compare homologous and cross reactive boosting profiles of different exposure types.

Study data
Antibody titre data were obtained from a previously published ferret study [43]. The experimental protocol was originally designed to reflect different possible human infection and vaccination histories at the time of the 2009 pandemic. Here we present a secondary analysis of these data, with the intention of characterizing underlying immunological processes.
Briefly, five experimental groups each of three ferrets underwent different combinations of infection with seasonal influenza A and/or vaccination with Northern and Southern Hemisphere trivalent inactivated influenza vaccine (TIV), with or without Freund's incomplete adjuvant (IFA), over the course of 70 days (Table 1). Serum samples were collected at days 0, 21, 37, 49 and 70 from all ferrets (Fig 1). HI titres were used to determine antibody titres to each infection and TIV strain. Dilution plates with 12 wells were used, such that the highest possible recorded dilution was 1:40960, and the lowest detectable titre was 1:20. Undetectable titres were recorded as <1:20. All analyses here were carried out using log titres, defined as k ¼ log 2

Models of antibody kinetics
The mathematical model describes the kinetics of homologous and heterologous antibody titres following exposure. Fig 2 depicts the example of an individual becoming infected and later vaccinated, though the model may characterise any sequence of exposures. Conceptually similar mathematical models of boosting followed by biphasic waning have been used previously to describe antibody secreting cell (ASC) and antibody kinetics [50,52,60,61].
After an infection (start at time ξ 1 ), homologous antibody titres undergo boosting rising linearly (on the log scale) by μ 1 log units to a peak after time t p1 , ignoring any delay between exposure and the start of antibody production. Titres then quickly drop by a fixed proportion, d 1 , over t s1 days (in the timescale of a few days to weeks), representing the initial short-term waning phase as free antibodies and early short-lived ASCs begin to decay following clearance of the initial antigen dose [33]. Antibody waning then switches to a constant rate m 1 (log titre units lost per day) for the remainder of time (representing the population of persistent ASCs, lasting months to years) until subsequent vaccination (syringe at time ξ 2 ), when antibody dynamics become dominated by a new set of boosting and waning parameters. We did not include a third, steady state phase due to the short time frame of these experiments [62]. Antibodies effective against heterologous strains experience boosting and biphasic waning in proportion to the exposure strain, with the proportion dependent on the antigenic distance between the measured and exposure strains. We assumed that the lower bound of detection of the HI assay (a log titre of 0) was synonymous with a true absence of antibodies, such that model predicted titres could not wane below this level.
We then built on this base model to incorporate additional immunological mechanisms that are important in describing antibody boosting and waning. These included: biphasic or monophasic antibody waning; exposure-type specific or type non-specific cross-reactivity; antigenic seniority; the impact of priming infection on subsequent vaccine response; and titredependent boosting. We considered models with different numbers of exposure types to match the experimental design: either 3 (infection, TIV, TIV + adjuvant) or 6 (priming infection, secondary infection, initial TIV, secondary TIV, initial TIV + adjuvant, secondary TIV + adjuvant). The base boosting and waning model remains the same across model variants, but these mechanisms add complexity to the boosting parameter, μ, and link different exposures with common parameters. A full description of each of these mechanisms and their implementation is described in S1 Supporting Protocol.
We fit each of the 64 potential model variants in a Bayesian framework using parallel-tempering Markov chain Monte Carlo (PT-MCMC) to estimate the posterior medians and 95% credible intervals (CI) of all free model parameters. For each model, we ran 3 chains each for 5000000 iterations. Where the effective sample size (ESS) was <200 or the Gelman-Rubin diagnostic (R) was <1.1 for any estimated parameter (calculated using the coda R package [63]), we ran 5 chains each for 10000000 iterations and obtained upper 95% confidence intervals forR of <1.1 for all estimated parameters presented here. ESS andR estimates for all parameters are provided in S5 Table. We then performed a model comparison analysis using Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO) with the loo R package [64,65]. Briefly, the purpose of this analysis was to compare the expected log point-wise predictive density (ELPD) of different model fits to compare their out-of-sample prediction accuracy. Comparing ELPD estimates serves a similar purpose to comparing other information criteria, where a lower ELPD suggests greater predictive power penalised by model complexity. Results shown in the main text are from the most complex model (most free parameters) variant with δELPD<1 compared to the lowest ELPD. Parameter estimates from all model variants with a δELPD<20 are shown in S7-S14 Figs. Posterior parameter estimates are shown as medians and 95% CIs. Further details of the model fitting and comparison are described in S1 Supporting Protocol. All code and data are available as an R package at https://github.com/jameshay218/antibodyKinetics.

Antibody kinetics following a single exposure support biphasic waning
To validate our boosting and biphasic waning model for a single exposure, we fit the base model to HI titres against A/Panama/2007/1999 (H3N2) from group E alone (Fig 3). Ignoring the later exposure to A/Fukushima/141/06 (H1N1) at day 56, from which we do not expect any cross-subtypic antibody reactivity, these data in isolation reflect a typical antibody trajectory following exposure to a single immunogen and measurement of antibodies against it [34,50]. The models with biphasic waning (both with estimated long term waning rate, m, and fixed m = 0) were better supported than the models with monophasic waning or no waning After each exposure, antibody levels undergo linear boosting on a log scale followed by an initial, short waning phase and then a slower, long-term waning phase. This example demonstrates two exposures, initially with infection (star symbol) and subsequently with vaccine (syringe symbol), where antibody dynamics are governed by a set of parameters depending on the exposure type. Note that the y-axis is on a log scale and all observations are discrete and taken as the floor value. Model parameters are described in S1 Table. https://doi.org/10.1371/journal.pcbi.1007294.g002 (ELPD -20.6 (standard error (SE), 3.37) and -20.5 (SE, 3.77) compared to -23.4 (SE, 3.71) and -28.1 (SE, 3.48) respectively), although we note that these differences are small with respect to the standard error of the ELPD estimates. The biphasic waning models with estimated longterm waning m and fixed long-term waning m = 0 had a difference in ELPD of <1, suggesting that both models had similar predictive performance. Overall, these results suggest that the model with monophasic waning is justified over the version with no waning, and that the biphasic waning model is better still than the monophasic waning model. Posterior estimates for model parameters were: μ = 9.91 (median, 95% CI 7.08-12.7); d = 0.551 (median, 95% CI 0.183-0.695); t s = 19.5 days (median, 95% CI 6.39-27.4 days); m = 0.0414 (median, 95% CI 0.00405-0.103).

Variation in antibody kinetics driven by different exposure histories
Overall, ferrets that received more frequent and immunogenic exposures achieved the highest, most broadly reactive and long-lived antibody titres. The full data show substantial variation in observed antibody titres across the groups driven by different exposure types and combinations. Following two doses of unadjuvanted TIV, ferrets achieved only modest increases in titres against the vaccine strains ( Fig 4A), with 2 out of 3 ferrets failing to generate H3N2 titres that persisted past day 37. The addition of an adjuvant resulted in increased and persistent titres against the vaccine strains in all ferrets by day 49. Titres against A/Fukushima/ 141/2006 (H1N1), which is antigenically similar to A/Solomon Islands/3/2006 (H1N1), were also increased at this time point (Fig 4B). Similarly, priming infection resulted in higher and long-lived titres to the vaccine strains and A/Fukushima/141/2006 (H1N1) relative to ferrets in the unprimed, unadjuvanted TIV protocol ( Fig 4C). Observed titres at day 21 against A/Panama/2007/1999 (H3N2) were consistently high following priming infection in groups C-E, with one ferret in each of groups C and E also experiencing some boosting of antibodies against the other H3N2 strains. All ferrets were infected with A/Fukushima/141/2006 (H1N1) at day 56, leading to elevated titres to both H1N1 strains by day 70 in all ferrets.

Model comparison results
The top two model variants had ELPD estimates of -412.1 (SE, 21.4) and -412.6 (SE, 20.7) respectively. Both of these models included: a role for priming infection in increasing subsequent vaccine response; different boosting profiles between vaccination and infection; different boosting profiles with adjuvant versus without adjuvant; and biphasic antibody waning. The model with the lowest ELPD (model ID 21, S6 Table) had 30 free parameters and also included titre-dependent boosting, no antigenic seniority, and no exposure type-specific cross reactivity. The other model (model ID 62, S6 Table) had 33 free parameters and did not include titre-dependent boosting, but did include antigenic seniority and exposure type-specific cross reactivity. Fig 4 shows the latter (more complex) model variant fitted to the data. Parameter estimates for these two models are shown in S4 Table. The remainder of the results refer to the latter model with more free parameters.
Other model variants may provide latent titre predictions more in line with biological expectations, though we note that they are not as well justified based on the model comparison analysis. For example, the predicted latent titres at the time of secondary vaccination (day 42) were unexpectedly lower in group D than group C under the chosen model variant (Fig 4C &  4D). Oil-in-water adjuvants are hypothesised to increase recruitment of neutrophils, antibody presenting cells and antigen bearing B cells at draining lymph nodes, and we would therefore expect antibody titres following adjuvanted TIV to be higher than unadjuvanted TIV throughout this time frame [66,67]. These unexpected results are likely due to limitations of the flexible model structure, which finds the set of parameter estimates best supported by all of the data, potentially at the cost of some biological realism. For example, a model variant identical to the one used in Fig Table).
Overall, ELPD estimates ranged from -412.1 (SE, 21.4) in the highest ranked model to -543.6 (SE, 21.8) in the lowest ranked model (S6 Table). The simplest model with 8 free parameters was the third lowest ranked model (ELPD -539.0 (SE, 21.9)), whereas the most complex model with 35 free parameters was the the 7th highest ranked model (ELPD -417.0 (SE, 21.2)). We note that some of the simpler model variants may have similar predictive performance to the best fitting model and may therefore be more suitable in a general predictive application. For example, further constraining the antibody trajectories in Fig 4E from Table). Our aim was not to predict unseen data but rather to quantify immunological mechanisms.
As a crude measure of mechanism importance, we performed Pseudo-Bayesian model averaging (Pseudo-BMA+) to estimate the relative weights of each model variant and thereby weights of models with a particular mechanism relative to models without [68]. Although comparison of variable importance using information criteria must be interpreted with caution (for example, changing the sample size or experimental protocol may change the results), Pseudo-BMA+ serves as a rough estimate of which mechanisms are most important in explaining these data [69]. Variable weights were: 1.00 for the presence of priming; 0.999 for the presence of 6 exposure types; 0.836 for the presence of biphasic waning; 0.579 for the presence of titre dependent boosting; 0.572 for the presence of type specific cross reactivity; and 0.406 for the presence of antigenic seniority. The top two models had Pseudo-BMA+ weights of 0.331 and 0.303, with a drop off to the third model with a weight of 0.0977. The top two models included only titre-dependent boosting and antigenic seniority respectively, suggesting that inclusion of at least one of these mechanisms improved predictive performance. The consistency of parameter estimates across the best fitting model variants is demonstrated in S7-S14 Figs.

Comparison of homologous boosting by exposure type
The level of homologous boosting resulting from priming infection (Infection 1) and secondary infection (Infection 2) was similar, shown by similar estimates for μ from both infections ( Fig 5A). We inferred that antibody titres fell only marginally following the initial waning phase (μ(1 − d), Fig 5B). The antibody waning rate was not identifiable for secondary infection due to the lack of observations following this exposure. We found evidence for only low levels of homologous antibody boosting following both initial and secondary doses of unadjuvanted TIV (TIV 1 and TIV 2) that quickly waned to near undetectable levels during the initial waning phase.
The addition of an adjuvant appeared to have no significant impact on the homologous antibody response to the first vaccine dose, but did improve the response to a second dose of vaccine (TIV 1 compared to TIV 1 + adjuvant and TIV 2 compared TIV 2 + adjuvant, Fig 5B). Titres against A/Brisbane/10/2007 (H3N2) and A/Solomon Islands/3/2006 (H1N1) were similar following the first unadjuvanted vaccine dose and the first adjuvanted vaccine dose (TIV 1 compared to TIV 1 + adjuvant, Fig 4A & 4B). However, the second adjuvanted TIV dose appeared to elicit a significant persistent boost to the vaccine strains, which resulted in peak titres near the limit of detection of this assay (TIV 2 compared to TIV 2 + adjuvant, Fig 4A &  4B).

Comparison of cross reactivity by exposure type
In models with type-specific cross-reactivity, we found differences in the width of cross reactivity elicited by the 6 exposure types shown in Fig 5. Secondary infection appeared to elicit a level of cross reactivity in line with that of the priming infection, whereas cross reactivity for both unadjuvanted and adjuvanted vaccination appeared to be narrower and only boosted antibodies that were effective against antigenically similar viruses (Fig 6). σ describes the degree by which antibody titre decreases as a function of antigenic distance, where higher values of σ suggest lower cross reactive breadth. When a single cross reactivity gradient was assumed for all exposure types (as in the highest ranked model), we estimated the cross reactivity gradient to be 2.33 (median; 95% CI: 1.74-3.01), suggesting narrower cross reactivity than would be expected given the definition for cross reactivity based on ferret antisera (an antigenic distance of 1 unit should see a reduction in antibody boosting of 1 log titre unit) [70]. Fig 6 demonstrates that homologous boosting (the y-intercept) was too small to elicit any measurable cross reactive boosting at these antigenic distances. The cross reactivity gradient parameter, σ, could therefore not be identified for the second dose of unadjuvanted TIV and first dose of adjuvanted TIV, and we were only able to recover the prior distribution for these parameters. These values were therefore excluded from Fig 5F.

Magnitude and duration of waning phases
Our model provided support for the presence of an initial short-term, rapid waning phase followed by a secondary long-term, sustained waning phase. For all vaccine doses, we estimated that the majority of the antibody boost waned within two weeks of reaching the peak (upper 95% CI 17.3, 5.15, 12.5 and 2.16 days for TIV 1, TIV 2, TIV 1 + adjuvant and TIV 2 + adjuvant respectively, Fig 5A & 5B). Conversely for priming infection, we estimated that the antibody titre was maintained at near peak levels with an estimated initial waning phase duration of 18.9 days (median; 95% CI 11.0-29.3) and a 21.6% (median; 95% CI: 2.02-48.8%) drop in log titre relative to the peak. We estimated similar long-term waning rates for second unadjuvanted TIV, second adjuvanted TIV and priming infection (Fig 5E).  We could not produce constrained estimates for the waning phases that take place following infection with A/Fukushima/141/2006 (H1N1) at day 56, given that only one subsequent observation was made at day 70. Although the 95% CI does not exclude biphasic waning rates consistent with the other exposures, any single trajectory in this range that passes through the single observation is similarly likely given these data.

Impact of priming
Prior to receiving non-adjuvanted TIV, experimental group C was infected with H3N2 Panama/2007/1999 at day 0, which represented a host being primed by natural infection prior to vaccination. Our model allowed us to identify additional homologous and cross-reactive antibody boosting that resulted from priming improving the subsequent vaccine response, as comparable experimental groups were given the same vaccination schedule with or without priming infection. The model suggested that priming infection added a substantial boost (7.28 log units (median, 95% CI 5.85-8.92)) to antibodies against the A/H1N1 and A/H3N2 vaccine strains at the time of vaccination in addition to that provided by the vaccine itself ( Fig 5A). Given our assumption that a log titre of 0 represents the true absence of antibodies, we cannot be certain that the higher titres observed at day 37 are due to a single large boost from primed vaccination rather than an antibody boost below the limit of detection from the priming infection followed by a small subsequent TIV boost. However, previous antibody kinetics results showing higher vaccine-induced antibody boosting following priming from the same detectable starting titre suggest that the former explanation is likely [71].
We estimated the cross reactivity of this additional boost to be broad with a gradient of 0.882 (median; 95% CI: 0.531-1.49), suggesting that priming increases the cross-reactive breadth of the vaccine response. It should be noted that whilst additional priming-induced vaccine boosting is well supported by the model fit, the model overestimates the antibody titre to A/Fukushima/141/2006 (H1N1) at day 37 elicited by initial dose of TIV following priming by H3N2 infection (Fig 4C). This may be a result of subtype-specific interactions that are not captured by our model.

Limited evidence for antigenic seniority and titre dependent boosting
Despite the relatively short duration of these experiments, we found some evidence for a trend of decreasing antibody response with increasing number of prior exposures and/or higher preexposure titres. In the best fitting model with antigenic seniority, we estimated τ to be 0.213 (median; 95% CI: 0.134-0.300), suggesting that antibody boosting decreased substantially with increasing number exposures after taking into account exposure type and priming. τ measures the proportion of the full boost that is lost with each successive exposure experiences relative to the first (ie. boosting decreases linearly as a function of increasing prior exposures). A higher value of τ therefore indicates more boosting suppression with an increasing number of prior exposures. Based on these estimates, the amount of antibody boosting would be reduced by over 50% following 4 exposures.
In the best fitting model with titre-dependent boosting, we estimated the titre-dependence gradient γ to be 0.0898 (median; 95% CI: 0.0788-0.102) applying to all titres below 10.9 (median; 95% CI: 8.62-11.95). γ gives the proportion of full boost that is lost per unit increase in log titre at the time of exposure (with no suppression from a starting log titre of 0). The full posterior estimate for the titre-dependent boosting mechanism is shown in S5 Fig. We note that the inferred titre-dependent boost relationship may be different if the limit of detection of the HI assay was lower. The top two model variants incorporated one of antigenic seniority or titre dependent boosting, suggesting that either one significantly improves model fit relative to the model variants with neither. The two mechanisms are correlated in these experiments, and antigenic seniority was not well identified with estimates for τ that did not exclude 0 for models with both titre-dependent boosting and antigenic seniority. However, all of the top models with antigenic seniority but no titre-dependent boosting give constrained estimates for τ away from 0 (S14 Fig).

Discussion
In this study, we used a mathematical model of antibody kinetics to describe boosting and waning following influenza vaccination or infection in a group of well characterised ferrets. We fit various subsets of the model with different immunological mechanisms and found that the two best supported models both included: type-specific antibody boosting; type-specific biphasic waning; 6 distinct exposure types; and a role for priming in increasing subsequent vaccine response. Antigenic seniority, antigenic distance-mediated cross reactivity specific to each exposure type and titre-dependent boosting were also included amongst these top models, suggesting that these mechanisms may be important in accurately describing observed antibody titres following multiple exposures. We found quantitative differences in the level of homologous and cross-reactive antibody boosting between vaccination, infection and adjuvanted vaccination in this ferret model. A single TIV dose with or without adjuvant elicited negligible levels of homologous and cross reactive boosting. A second dose of TIV with adjuvant resulted in significant, broadly reactive antibody boosting, whereas a second dose of TIV without adjuvant did not elicit significant antibody boosting. The profile of boosting for primary infection was consistent across experimental groups, and similar in magnitude to secondary infection. Furthermore, we found that priming infection induced a significantly broader and stronger boosting profile following subsequent vaccination.
Our work has a number of limitations. The model predicted latent antibody titres were broadly in line with expected immune dynamics, though there were some exceptions. The model variant presented in the main text was selected based on the optimal balance between fewest parameters and accuracy of predicted antibody trajectories with respect to the full set of observed titres. Other model variants or designs may provide results more in line with biological expectations, but would require estimation of more parameters or collection of additional data. Given the relative sparsity of samples across time, some aspects of the model were poorly identified or included spurious features for some subsets of the data. In particular, sampling around the biphasic waning period of the vaccinations and following the final exposure event was limited, resulting in poor identifiability for some of the waning and timing parameters. We therefore restricted our reported results to estimates that were consistent across the best supported model variants (S7-S14 Figs). Experiments of a similar design with fewer exposures and more frequent sampling would power the model to elucidate these waning phases further and look for differences in response longevity by exposure type.
Our experimental timeline was much shorter than the typical human exposure timescale of months to years, with a minimum gap of 14 days between the two vaccine doses and 28 days between infection and vaccination [18]. No further increase in antibody titre was detected from 14 days post TIV 1 in group A or post A/Panama/2007/1999 (H3N2) infection in group E, suggesting that serum antibody titres consistently peaked within 14 days and therefore before each subsequent exposure. Furthermore, germinal centre (GC) structures and GCderived ASCs had likely developed within this time, as it has been shown in other small mammals (mice) that GC B cells are present within 14 days post infection [39, 72,73].
However, given that GC responses peak after 4 weeks and persist for months in mice, there may not have been sufficient time for the MBC population, which is a significant contributor to post-vaccination ASCs in humans [23], to fully develop. Recruitment of naïve B cells may also be impacted by both the presence of pre-existing GC structures from primary infection and the inclusion of IFA [66,67,74]. The inferred antibody kinetics of the two TIV doses might therefore be different if they were further apart. It is possible that the time to peak response parameter would also differ depending on the relative contributions of memory and de novo ASCs; however, models that allowed the time to peak parameter to vary were less well supported based on ELPD, resulted in reduced identifiability of some parameters, and did not change estimates for identifiable parameters (S1 Supporting Protocol).
We captured the impact of immune memory on subsequent antibody responses in two ways: antigenic seniority and titre-dependent boosting [26]. Titre-dependent boosting was a function of only homologous antibody titres, whereas antigenic seniority was assumed to be a function of all previous exposures regardless of exposure strain, subtype and administration route. If antigenic seniority is primarily a result of HA reactive antibodies, then it would only act within a subtype in contrast to our assumed mechanism [32]. However, it is not clear that immune memory effects do not extend between influenza subtypes [75,76]. For example, immune imprinting to a particular subtype has been proposed as an explanation for age-specific mortality in the 2009 H1N1 pandemic and to explain the age distribution of avian H5N1 and H7N9 cases, though these observations do not necessarily suggest any cross-subtype or cross-group impact on subsequent antibody responses [30,77]. It has been proposed that cross-subtype effects might act through cross-reactive memory T-cell responses that act to deplete heterosubtypic antigen load, which may in turn lead to lower antibody boosting [24].
HI assays measure only the aggregated activity of polyclonal antibodies targeting multiple epitopes on the haemagglutinin head, and we were therefore unable to investigate epitope-or B cell-specific contributions to serum antibody titres [78]. Understanding the immunological mechanisms of imprinting effects, short-term kinetics due to GC overlap, and the relative contributions of MBC derived and de novo antibody boosting would require data either on epitope-specific antibodies or single-cell assays [23,39,79]. Epitope masking, wherein preexisting antibodies targeting recognised, but poorly conserved epitopes (ie. the epitopes that generate cross-reactivity in the HI assay between eg. A/Panama/2007/1999 (H3N2) and A/ Wisconsin/67/2005 (H3N2) antibodies) sterically mask access to previously unseen and conserved, immunosubdominant epitopes (eg. at the receptor binding site), has been proposed as a mechanism for the recall bias of subsequent responses [19,31,[80][81][82]. A model that captures the contribution of MBCs and naïve B cells to antibodies targeting a variety of epitopes may explain how immune imprinting contributes to observed antibody titres, and may include additional insights such as decreased cross-reactive breadth and magnitude with each repeated exposure.
Our data included only trivalent vaccination with and without IFA, and estimates of any boosting parameters are therefore conditional on the presence of three antigens in a single vaccination. It would be interesting to compare the inferred homologous and cross-reactive boost of different vaccination strategies (eg. a three antigen TIV compared to a monovalent vaccine, or comparison by inoculum dose), and for adjuvants more relevant to human vaccination such as MF59 [83,84]. There may also be underlying heterogeneities in antibody response between and within influenza subtypes as well as between vaccine types [85,86]. For example, Live Attenuated Influenza Vaccines (LAIV), as well as newer DNA vaccines may provide different antibody kinetic profiles and may elicit broader antibody responses, or provide different priming effects [58,87].
We found evidence for biphasic waning following both primary infection and secondary vaccination. There was some evidence that the magnitude and duration of waning differed between exposure types: TIV 1 and 2 and adjuvanted TIV 1 waned very quickly, whereas Infection 2 and TIV 2 + adjuvant were more persistent. Heterogeneity in antibody waning rates between individuals and vaccine types have been shown for other pathogens [60,88]. Although studies of influenza antibody response duration have been carried out in humans, quantifying waning rates independent of subsequent exposures that cause repeated boosting is difficult [2,61,[89][90][91]. Our model fit to the single exposure ferrets provides an estimation of the waning rate of homologous antibodies in the absence of further exposure, but the cut off of 70 days limits the applicability of this waning rate to a timescale more relevant to humans. Extrapolating our estimated waning rate following primary infection would suggest that antibody titres would wane to non-detectable within a few months, whereas antibody responses against many viruses are known to persist for decades [88,92]. Longer term studies investigating the longevity of the antibody response in the absence of repeated exposure would be useful to quantify a long-term, steady state antibody waning rate [52]. Further mechanisms such as differential waning rates between cross-reactive and homologous antibodies are likely to be important, but were not identifiable here [21,59]. Although animal models are potentially useful, identification of these mechanisms in human populations is likely possible given long-term, frequent sampling of human sera combined with robust statistical methods [18,22].
Our results have implications for comparing different vaccination strategies. Achieving high HI titres against currently circulating strains is a key endpoint in influenza vaccine trials due to its correlation with clinical protection [4,[93][94][95]. However, there are a number of obstacles to achieving these high titres in some populations including antigenic interactions, age specific responses and antibody waning [11,[96][97][98][99]. One approach to improving vaccine effectiveness may therefore to elicit a broader antibody response to compensate for potential strain mismatch [100]. Adding adjuvants such as MF59 and AS03 has been shown to induce higher antibody titres that have greater cross-reactive properties [55,56,101,102]. Quantitative comparisons of cross reactivity profiles, as we have provided here, could be a useful tool in comparing the effectiveness of different adjuvants, which would provide a measurable benefit to trade-off against safety and immunogenicity concerns [103,104].
In addition to modelling boosting suppression due to prior immunity, we considered potential enhancement via priming infection. "Prime-boosting" has been described previously as a strategy to induce broadly reactive immune responses that may be rapidly boosted in advance of exposure to an antigenically novel virus [53,54,71]. Models that included a priming mechanism were ranked systematically higher in our model comparison analysis than those that did not, suggesting that this phenomena is important in explaining titres arising from repeated exposures. We found that vaccine responses to A/H1N1 strains were higher and more broadly reactive in A/H3N2 primed ferrets compared to unprimed ferrets, though our model did not account for subtype specific interactions and subsequently overestimated post vaccination A/H1N1 titres in primed ferrets. Although the phylogenetic relationship between the priming and subsequent boosting strain is likely to be important, heterosubtypic protection has been shown previously in animal models, potentially via cytotoxic T lymphocyte responses [45,47].
Our results suggest that mathematical models of antibody kinetics that explicitly consider immunological mechanisms and exposure-type specific parameters would be valuable for the prediction of antibody landscapes in human populations. Human cohort studies tracking infants from birth as they experience their first few influenza exposures are also now underway [105]. Combining these studies with single-cell immune profiling and mathematical models of multiple exposure kinetics will help to elucidate the role of these immunological mechanisms in building human antibody profiles. Direct inference from long-term observational data in humans may be difficult, but experimental models, such as the ferret system described here, provide an excellent alternative data source for the inference of short-term immunological mechanisms that may map onto models recovered using human sera [18, 21, 41, 42].  Table. Description of models with δELPD <20. Table is  Probability of observing a particular log titre given an underlying true, latent titre. Note that the true titre is a continuous value, whereas observations are discrete. Furthermore, truncation of the distribution at the upper and lower limit of the assay results in an asymmetrical distribution when the true value is at either of these limits. True values outside of these limits will be observed as a value within the assay limits. (TIF) Estimates are stratified by exposure type and ordered in order of increasing ELPD. Estimates are coloured according to whether or not cross reactivity was assumed to be a universal parameter or type-specific. Dashed horizontal lines represent uniform prior range. Model codes on x-axis relate to the first letter of each mechanism as described in S2 Table. (TIF) Estimates are coloured according to whether or not titre-dependent boosting was included. Dashed horizontal lines represent uniform prior range. Model codes on x-axis relate to the first letter of each mechanism as described in S2 Table. (TIF) S10 Fig. Summary of posterior distribution estimates for long-term waning rate, m from models with δELPD <20. Points show posterior median; line ranges show 95% credible intervals. Estimates are stratified by exposure type and ordered in order of increasing ELPD. Estimates are coloured according to whether or not waning was assumed to be biphasic or monophasic. Dashed horizontal lines represent uniform prior range. Model codes on x-axis relate to the first letter of each mechanism as described in S2 Table. (TIF) S11 Fig. Summary of posterior distribution estimates for cross reactivity gradient, σ from models with δELPD <20. Points show posterior median; line ranges show 95% credible intervals. Estimates are stratified by exposure type and ordered in order of increasing WAIC. Estimates are coloured according to whether or not cross reactivity was assumed to be a universal parameter or type-specific. Plots are truncated from above at 10 for clarity, but upper prior bound was 100. Red dashed line shows the fixed value of σ = 1 for priming infection. Blue dashed line shows value above which a homologous boost of μ = 5 would give an observed boost of 0 against a strain with an antigenic distance of 1. Model codes on x-axis relate to the first letter of each mechanism as described in S2 Table. (TIF)