We formulated different systems of nonlinear ODEs to model the observed antibody dynamics. We assumed that the incubation period is the same for all individuals and we take the value proposed in the literature. We assumed antibody levels to change over time due to a proliferation function (f ₁) and a decay function (f ₂). In all models, we assumed that proliferation depended on the number of antibody secreting cells (ASC, e.g., plasma cells), and decay on the number of antibodies (AB). For ease of presentation, we suppress the time dependence of the number of AB and the number of ASC from the notations below. Hence, we can express the progression of the number of antibodies in the following general differential equation:

Next, we describe the number of ASC using a second differential equation. We assume that proliferation of ASC occurs according to a function g ₁, during a time period [0, h] after which no new ASC will be generated. ASC decay is assumed to occur at all-time points according to a decay function g ₂. This leads to the second differential equation:

where I_t≤h represents an indicator function that takes value one if t ≤ h, and is equal to zero otherwise.

In all our proposed models, antibody decay is assumed to be proportional to the number of antibodies and can therefore be written as: f ₂ (AB) = u_AB × AB, with u_AB denoting a time-invariant decay rate (i.e., implying exponential decay). Similarly, the antibody production is assumed to be proportional to the number of ASC, and thus f ₁ (ASC) = p_AB × ASC, where p_AB is the constant antibody production rate by ASC. Equations (1) and (2) can be combined as follows:

where AB ₀ = AB(0) and ASC ₀ = ASC(0) are the antibody and plasma cell counts at time 0 (days), respectively. We can further divide the population of ASC into two sub-populations: one with a short lifespan and the other with a long lifespan (11, 12). Short-living antibody secreting cells (SASC) can be interpreted as B-cells (whether or not including plasmablasts). Long-living antibody secreting cells (LASC) can be seen as plasma cells, which can reside in the bone marrow with a lifespan of many years (13, 14).

Models (4), (5) and (6) treat the population of ASC as a single population without discriminating between short- and long-living plasma cells, while models (7), (8) and (9) look at the aforementioned sub-populations separately. by the choice of the proliferation function g ₁(ASC).

In model (4), we assume that ASC expansion occurs with the constant function g ₁(ASC) = p_ASC during time period [0, h]. In model (5), expansion happens proportional to the number of ASC, g ₁(ASC) = p_ASC × ASC. Model (6) assumes no proliferation g ₁(ASC) = 0 and thus the number of ASC monotonously decreases over time starting from the initial level ASC ₀. We can interpret this as proliferation that can be neglected over a short time period [0, h].

Similarly, we defined models (7), (8) and (9), explicitly distinguishing between LASC and SASC. First of all, in order to limit the number of model parameters, we assume that the number of LASC remains constant over time, i.e., dLASC / dt = 0. Moreover, we assume no SASC at time 0, i.e., SASC(0) = 0 and decay of SASC following g ₂(SASC) = u_SASC × SASC.

Models (7) to (9) rely on the use of different proliferation functions for SASC. More specifically, in model (7) a constant expansion g ₁(SASC) = p_SASC is assumed, in model (8) a proportional expansion g ₁(SASC) = p_SASC × SASC is considered and finally in model (9) no expansion g ₁(SASC) = 0 is presumed.

We can now formulate nonlinear mixed models based on the dynamic models described in the previous subsection. The nonlinear mixed effects model implemented is defined as:

where y _ij ∈ℝ is the jth observation in the ith subject, Ψ_i is the parameter vector of the structural model f for individual i. The residual error model is defined by the function g of the structural model f and an additional vector of parameters ξ. The residual errors (ε_ij ) are standard Gaussian random variables (mean 0 and standard deviation 1). In this case, it is clear that f(t _ij ,ψ _i ) and g(f(t _ij ,ψ _i ),ξ) are the conditional mean and standard deviation of y_ij , E(y _ij |ψ _i )=f(t _ij ,ψ _i ) and sd(y _ij |ψ _i )=g(f(t _ij ,ψ _i ),xi) .

The assumption that the distribution of any observation y_ij is symmetrical around its predicted value is a very strong one. In order to achieve that, we may want to transform the data to make it more symmetric around its (transformed) predicted value. An extension of the statistical model to include a Log-normal distribution was therefore proposed as transformation since all observations are

In mixed models, both fixed and random effects are included. Consider a parameter P_pop (e.g., the decay rate of antibodies) which is constant across different individuals. Such a parameter can be interpreted as a population(-averaged) parameter, representing the average value across all individuals, while the use of random effects leads to an individual-specific interpretation of parameters and a description of possible variation between individuals. An individual-specific model parameter P_i can be written as P _i =u _i ×P _pop , where P_pop is the so-called population parameter and u_i is an individual-specific random effect. In this paper, we will assume that u_i follows a log-normal distribution with a unit mean (to ensure identifiability and the aforementioned interpretation for P_pop ) and variance ω ².

Categorical variables, such as sex or the use of antivirals, can be taken into account by adding a dummy variable with an additional parameter β_j describing how the parameter of group j deviates from the reference group.

For example, if we consider sex as a categorical variable assuming that the population values of antibodies are different for male and female, we implemented the following model: log(A _i )=log(A _pop )+β _A 1 _{sex i =F} +η _A,i , where 1 _{sex i =F} if individual i is a female and 0 otherwise. Then, A_pop is the population antibody for males while A _pop e ^{β _A} the population antibody for females. This allows investigating which particular parameter of the structural model (e.g., proliferation of antibodies, decay of antibody secreting cells,…) leads to the differences observed between different groups (e.g., male vs. female individuals). Dependency can be introduced between individual parameters by assuming that the random effects η_i are not independent.

The model parameters were estimated using the Monolix software ^©Lixoft (2021 version) (15). The estimation of the population parameters was achieved by a two-step algorithm with 10⁶ + 10⁵ iterations to assess convergence. This algorithm consists of a built in stochastic approximation of the standard expectation maximization algorithm (SAEM) with simulated annealing, combined with a Markov Chain Monte Carlo (MCMC) procedure which replaces the simulation step of the SAEM algorithm (16). Next, importance sampling was used to determine the log-likelihood per model, in which a fixed t-distribution is assumed with 5 degrees of freedom. The models were assessed using the Akaike Information Criterion (AIC) (17).

For the antibody data set, the following procedure was used for comparing and selecting the most suitable biologically plausible model to describe the data.

In order to compare the different models and select the most suitable and biologically plausible model to describe the antibody data, we deployed the following procedure. First of all, for each model, the model parameters were estimated using the Monolix software. In case of SAEM-MCMC convergence, we used the AIC (17) to compare the different models resulting in a list of six models to be compared in the next step. Models with poor SAEM convergence, likely because of abundant model complexity, were discarded. Subsequently, the model with the lowest AIC value on the list was selected as the first candidate model. We performed a non-parametric bootstrap with 1000 bootstrap samples to investigate model stability. If the bootstrap did not converge well, we excluded the model from the list of candidate models. Since a sequential approach based on the candidate models with the lowest AIC values was used, the need to perform bootstraps for all candidate models was avoided, in order to decrease the number of computations. It was found that for a bootstrap, 63%-77% of the samples had proper SAEM convergence. For this reason, the criterion for good bootstrap convergence was defined as having at least 65% of bootstrap samples with proper SAEM convergence. Finally, a sensitivity analysis of the (converging) candidate model’s bootstrap results was performed to determine whether the presence or absence of specific profiles of individuals in the bootstrap samples influenced model convergence. For example, if a single participant’s profile appeared more frequently in a non-converging dataset, a new bootstrap was run, excluding the specific participant’s profile. Again, if the bootstrap convergence was poor, the candidate model was rejected. If convergence remained robust enough, the candidate model was chosen as the final model. (for more details see Appendixes C and D ).

It is possible to introduce dependency between individual parameters by assuming that the random effects η_i are not independent. This can be done by introducing linear correlation between the random effects. To test for this type of correlation, we performed Pearson’s correlation tests. In our scenario, we found a significant correlation between η _{AB 0} (the initial value of antibodies) and η _{μ AB} (the decay rate of antibodies). This indicates that the distributions of these two random effects are not independent and that the correlation must be included in the model and estimated.

Once the correlation is included in the model, the random effects for AB ₀ and μ_AB are drawn from a joint distribution rather than two independent distributions. In our scenario, including the correlation in the model resulted in an AIC of 4232.25, compared to 4243.84 for the original model. No other significant correlations were found.

In order to understand the drivers of inter-individual variability, building a covariate model is a critical effort to explain antibody population dynamics. Many runs are usually required to find a solid covariate model (18). Several ways to automate this operation have been proposed in the past. Stepwise Covariate Modeling (SCM) is the most often utilized. We apply a new stepwise strategy based on statistical tests between individual parameters taken from their conditional distribution and covariates validated by Lixoft and published elsewhere (18).

The Conditional Sampling usage for Stepwise Approach based on Correlation testing (COSSAC) uses the information in the current model to determine which parameter-covariate relationship to test next. This technique drastically decreases the number of covariate models examined while keeping the models that improve the log-likelihood on the search path. The Pearson correlation test for continuous covariates and ANOVA for categorical covariates can be used to calculate p-values. The p-values are used to sort all of the random effect-covariate correlations, regardless of whether or not they are included in the model. Depending on the results of the correlation tests, COSSAC iterations switch between forward and backward selection. Group-specific effects on chosen parameters make it possible to investigate the influence of covariates and correlation in the selected model 9. More specifically, we investigated whether participants’ sex, age, viral load, and use of antivirals affect model parameters, as well as whether there is any correlation between random effects.

We acquired the best fitting model after testing and running all possible parameter-covariate relationships, resulting in an AIC value of 4214.14 improving the AIC value of the original model 9 by 29.7 points (see Table 2 for the estimates of the model parameters). Since there is a significant correlation between AB ₀ and u_AB in this model, where AB ₀ is the initial value of the antibodies and u_AB is the decay rate allowing the correlation to be included in the model and approximated.

This model, on the other hand, yielded viral load and age as significant covariates. More particularly, an increasing viral load was found associated with a much higher proliferation rate p_ABS (see Figure 6 ). The decay rate of SASC, u_SASC , was shown to decline considerably with increasing age.

Viral load refers to the amount of virus present in a person’s body and is typically measured in terms of the concentration of virus in a sample of blood or other bodily fluid (19). High viral load can indicate a more severe infection or a greater risk of transmission to others. Age is another important factor that can influence immune function and the body’s response to infection or vaccination (20, 21). As people get older, their immune system may become less efficient at mounting a response to foreign substances, such as viruses or vaccines.

Specifically, we found that increasing viral load was associated with a higher production rate of antibodies by SASC (p_ABS ). This finding is supported by the study of (22), which found that increased levels of gp350-specific neutralizing activity were directly correlated with higher peripheral blood Epstein-Barr virus DNA levels during acute infectious mononucleosis. This suggests that individuals with higher viral loads may have a more pronounced immune response to varicella-zoster virus, resulting in the production of more antibodies. This has important implications for understanding the immune response to varicella-zoster virus and could potentially inform vaccine development efforts. It is worth noting, however, that while higher viral loads may lead to increased antigen exposure and potentially more need for antibody production, it is also possible that individuals with higher viral loads may have less effective antibodies and therefore require a stronger immune response to control the virus. Further research is needed to confirm and better understand these findings and the underlying mechanisms involved.

We found that sex and antiviral load did not significantly affect the results of the model. This may be due to the small sample size or the short duration of the study. Further research with larger sample sizes and longer follow-up periods would be needed to more fully examine these questions.

Additionally, we found that the decay rate of certain antibody-secreting cells was shown to decrease considerably with increasing age. This finding is supported by a previous study that demonstrated a positive correlation between VZV IgG and ageing (22). Further research is needed to confirm this finding and to better understand the potential underlying mechanisms involved.