Quantification of plasma cell dynamics using mathematical modelling

Plasma cells (PCs) are the main antibody-producing cells in humans. They are long-lived so that specific antibodies against either pathogens or vaccines are produced for decades. PC longevity is attributed to specific areas within the bone marrow micro-environment, the so-called ‘niche’, providing the cells with required growth and survival factors. With antigen encounters, e.g. infection or vaccination, new PCs are generated and home to the bone marrow where they compete with resident PCs for the niche. We propose a parametrized mathematical model describing healthy PC dynamics in the bone marrow. The model accounts for competition for the niche between newly produced PCs owing to vaccination and resident PCs. Mathematical analysis and numerical simulations of the model allow explanation of the recovery of PC homoeostasis after a vaccine-induced perturbation, and the fraction of vaccine-specific PCs inside the niche. The model enables quantification of the niche-related dynamics of PCs, i.e. the duration of PC transition into the niche and the impact of different rates for PC transitions into and out of the niche on the observed cell dynamics. Ultimately, it provides a potential basis for further investigations in health and disease.

Days after TT boost ASCs secreting TT-IgG  Table S2. Dataset (D2): Dynamics of serum IgG to tetanus toxoid (TT). Serum anti-TT IgG is measured in µg per ml. Data refer to the same donor as in Table S1 and were obtained from [1].

1
A published dataset by Bernasconi et al. [1] is used to characterise the dynamics of vaccinespecific and vaccine-nonspecific plasma cells (PCs) after vaccination with tetanus toxoid (TT) for a follow-up of 145 days. Data for one human donor comprise two sets of immunoglobulin type G (IgG) measurements, i.e. (D1) the number of antibody-secreting cells (ASCs) secreting IgG specific to TT in the peripheral blood (see Table S1), and (D2) the dynamics of serum IgG specific to TT (see Table S2). Dataset (D1) is used as time series which describes TTspecific PBs created after TT boost, and dataset (D2) is used as time series for the number of TT-specific PCs by applying a transformation of the Ig data in order to match them with the observables of the vaccination model (M).
From estimates based on large cohorts of adult normal-weighted healthy individuals it can be deduced that about 95% of the Ig concentrations per litre of blood lie within the interval (7.65, 20.13) (in g/l) [2]. The average Ig concentration can be estimated to comprise approximately 13 g/l with fractions of common Ig types given by 10 g/l IgG, 2 g/l IgA and 1 g/l IgM [2,3]. Ig types with marginal incidence (IgD, IgE) are neglected. The fraction of an Ig type corresponds to the fraction of PCs secreting that particular Ig type [4]. Estimates for the average total number of PCs in the bone marrow are deduced from two different approaches: 1. Proposed estimates: The total number of PCs in the human bone marrow is estimated to be about 10 9 cells [5]. The bone marrow cellularity is estimated to comprise about 10 12 cells [6]. Assuming a fraction of PCs of 0.25 − 1% [5,7], the total number of PCs can be calculated to be in the range of 2.5 · 10 9 − 10 10 cells. A comparable result using that the number of PCs in the bone marrow comprises about 7.5 · 10 9 cells per litre of bone marrow can be obtained [8]. With a total bone marrow volume of 1.6 − 4.0 litres [9], this implies an estimate of PC counts in the range of 1.2 − 3 · 10 10 cells. Summarising reported values yields that the average total number of PCs in the human bone marrow lies within 2.5 · 10 9 − 3 · 10 10 cells.

2.
Estimates based on Ig synthesis rates of PCs: The total amount of IgG at homoeostasis is assumed to be 50 g. The half-life time of IgG molecules can be estimated to be 20 days [10][11][12]. IgG synthesis rate is in a range of 20 pg -35 pg per PC per day with mean of 26 pg per PC per day [13]. In order to relate the synthesis rates to the total number of PCs, we define G(t) as the amount of IgG molecules (in gramme) at time t (in days). The production and degradation of IgG is modelled as where C(t) is the total number of the IgG-producing PCs at time t, s is the synthesis rate, and d g is the degradation rate of IgG molecules (which is given by the reciprocal of the half-life time multiplied by ln (2)). At homoeostasis, it holds sC = d g G, where C and G are the total amount of IgG at equilibrium (which is known), and the total number of the IgG-producing PCs at equilibrium, respectively. Using this mean synthesis rate results in 26 · 20 · 50 · 10 12 PCs ≈ 6.7 · 10 10 PCs.
Taking into account that C reflects 10/13 of the total number of PCs [2,3], the total number of PCs can be estimated to be comprised of 8.7 · 10 10 cells.
Combining both ways of estimation results in a span of 2.5 · 10 9 − 8.7 · 10 10 cells. Taking the mean order of magnitude motivates the assumption that the average total number of PCs in the human bone marrow comprises 10 10 cells. Consequently, a measurement of 1 g/l Ig in the serum is equivalent to 1/13 · 10 10 PCs in the bone marrow.

S2 Simplified vaccination model
In this section, we consider the following system of ordinary differential equations (ODEs) for times t ≥ 0, referred to as simplified vaccination model, If not stated differently, without loss of generality, we consider initial conditions denoted as respectively. Note that system (SM1) corresponds to system (M) with g ≡ 0.

Reduced model
Since and the dynamics of the total numbers of PCs outside and inside the niche, and y(t) := y 0 (t) + y v (t), respectively, are described by the reduced model Initial conditions relate as x 0 = x 0 0 + x 0 v and y 0 = y 0 0 + y 0 v . This allows investigating x and y in the phase plane, and characterising the niching dynamics with respect to the function z given different initial conditions. The reduced model possesses exactly one equilibrium as it can be seen by setting the right-hand side of system (SM2) equal to zero, and solving for x and y. A straightforward calculation yields denoted as healthy equilibrium. Observe that E h corresponds to PC homoeostasis, which is characterised by n more PCs inside than outside the niche.

Vaccine-induced perturbation
Characteristics of the inflow kinetics of PCs due to vaccination (described by g(t) in the vaccination model (M)) vary within and among individuals [5,14]. Lacking data describing those kinetics imply the necessity for approximating the time-dependent inflow dynamics using a simpler formulation. A vaccine-induced perturbation of PC homoeostasis due to primary immunisation can be implemented into the simplified vaccination model (SM1) framework using a time-discrete formulation, which can be interpreted as initial condition.
Assuming that the system is at equilibrium prior to vaccination, and the perturbation occurs at time T , such perturbation may be characterised as follows: , meaning that the system is at healthy equilib- v > 0 corresponds to the amount of vaccine-specific PCs perturbing homoeostasis. Such formalism enables modelling vaccine boosts, i.e. further administrations of a vaccine following earlier doses. Prior to the boost, the simplified vaccination model (SM1) is at homoeostasis, which is characterised by a respective homoeostatic arrangement of vaccinespecific and vaccine-nonspecific PCs located on the manifold E (see main article). In particular, to investigate the number of vaccinations of antigens different from a previously used vaccine needed to reduce the newly established immunity characteristic to 50%, an iterative procedure is accomplished as follows: Starting at a healthy equilibrium E 0 ∈ E with compo- Here, vaccine-nonspecific PCs are generated by antigens different from the previously administered vaccine. This leads to a new healthy equilibrium E j ∈ E, and a further perturbation is carried forward. In accordance with the obtained results (see main article), we choose N = 3 · 10 7 and the iterative procedure is terminated as soon as

S3 Mathematical results
The presented models belong to the class of piecewise-smooth dynamical systems with switching manifold given by the set z(t) = 0 [15]. Existence and uniqueness of solutions for given initial conditions can be guaranteed due to Lipschitz continuity of the vector fields. The latter can be shown by straightforward calculations [16,17]. Likewise, elementary calculations applying the method of invariant sets yields non-negativity of solutions for biologically plausible parameter settings [17,18].

Stability of plasma cell homoeostasis
The healthy equilibrium E h of the reduced model and the manifold E of the (simplified) vaccination model are located on the switching manifold. Thus, linearised stability theory as applied for smooth dynamical systems [19] may lead to incorrect conclusions [15,20]. Instead, providing a common Lyapunov function, i.e. a function that is Lyapunov for each of the vector fields defining the systems dynamics in each of the phase space regions separated by the switching manifold guarantees global asymptotic stability of the homoeostatic equilibria [21].
(a) The healthy equilibrium E h of the reduced model (SM2) is globally asymptotically stable.
(b) The manifold of non-isolated equilibria E of the simplified vaccination model (SM1) is globally asymptotically stable, i.e. solutions asymptotically approach E.
is a strict Lyapunov function of system (SM2) for all (x, y) = E h . The proof is left to the reader. Thus, E h is globally asymptotically stable, which shows (a).
It holds that , where x(t) and y(t) are solutions of system (SM2). Due to (a), it follows that and in particular lim t→∞ z(t) = 0. It remains to show that lim t→∞ x v (t) = 0. Then, every Let ε > 0 arbitrary but fixed. Since lim t→∞ z(t) = 0, it follows that for 0 <ε < εd 2 max{b,c} there exists aT > 0 such that |z(t)| <ε for all t ≥T . Thus, for all t ≥T , it is Next, consider the following ordinary differential equation with initial condition, for all t ≥T . Using the comparison principle [16] implies x v (t) < w(t) < ε 2 + ε 2 = ε. Observe that the first inequality is true for all t ≥T , whereas the second holds for all t ≥ T , where T >T is chosen such that the second term in w(t) is smaller than ε 2 . This can be achieved since the latter is positive and exponentially decreasing in time. Due to non-negativity of x v (t) for all t ≥ 0, it follows that lim t→∞ x v (t) = 0. This proves (b).

Biological consequences
Biologically, Theorem S1 implies that any non-malignant perturbation of PC homoeostasis eventually results in PC homoeostasis again with the same number of PCs outside and inside the niche. In particular, this holds true for a vaccination-induced perturbation, where the composition of the recovered homoeostasis depends on the number of vaccine-specific PCs arriving in the bone marrow.

Plasma cell dynamics after a perturbation of homoeostasis
Characterisation of the dynamics of z(t) is provided using the concept of invariant sets [18].
The below defined sets R 1 , R 2 (M 1 , M 2 ) correspond to the gray (dark gray) sets visualised in Figure 3 (see main article).

Theorem S2.
(a) Consider the sets R 1 and R 2 defined by where R i : R 2 → R, i = 1, 2, are given by Then R 1 and R 2 are positively invariant with respect to the flow of system (SM2), i.e.
solutions starting in R 1 or R 2 , respectively, remain therein for all times. In particular, the function z(t) switches its sign at most once.
(b) Consider the sets M 1 ⊂ R 1 and M 2 ⊂ R 2 defined by where M i : R 2 → R, i = 1, 2, are given by Then M 1 and M 2 are positively invariant with respect to the flow of system (SM2), i.e. solutions starting in M 1 or M 2 , respectively, remain therein for all times. In that case, the function z(t) is monotonically decreasing or increasing, respectively.
Proof. First, consider the set R 1 . Let (x, y) ∈ ∂R 1 \ E h . The boundary of R 1 splits up into two parts, i.e. (i) then y = f d + n, and thus (x, y) = E h . For the first part, a straightforward calculation shows that and for the second part, it is This implies positive invariance of the set R 1 . The proof of the assertion for R 2 is similar and left to the reader. A change in sign of z(t) is only possible for a solution starting in . Then, the corresponding trajectory (x(t), y(t)) eventually enters R 1 ∪R 2 as t → ∞. This is due to global asymptotic stability of E h ∈ R 1 ∪R 2 (see Theorem S1). If (x(t), y(t)) enters R 1 ∪ R 2 at a finite time T > 0, then (x(t), y(t)) enters Note that E h cannot be reached in finite time due to uniqueness of solutions. This implies a sign switch in z(t). It follows that (x(t), y(t)) ∈ R 1 ∪ R 2 for all times t ≥ T , and no further switch in the sign of z(t) occurs.
Next, consider the set (1) Then, positive invariance of M 1 follows by positive invariance of R 1 . Next, z(t) is shown to be monotonically decreasing. For that, let (x 0 , y 0 ) ∈ M 1 . It follows that M 1 (x(t), y(t)) ≤ 0 for all t ≥ 0. Therefore, it is

Biological consequences
Theorem S2 allows studying the mode of decay of the surplus of PCs (outside or inside the niche) after an initial perturbation of PC homoeostasis. In particular, the dynamics of PCs after vaccination can be described. After an initial flow of PCs into the niche due to a high number of PCs outside the niche, the concomitant surplus of PCs inside the niche leads to PCs being expelled from the niche until PC homoeostasis is recovered. Such dynamics comes along with a switch in the sign of the function z(t) (see Figure 3 in the main article).

Characterisation of the switching time
For a vaccine-induced perturbation of the healthy equilibrium at time T , the following result determines the switching timet. It is z(t) > 0 for T < t <t, and z(t) ≤ 0 for t >t.
Moreover, it provides an upper bound for the number of vaccine-specific PCs in the niche at the re-established healthy equilibrium.
Theorem S3. Consider system (SM1) with initial conditions given by a vaccine-induced perturbation of the healthy equilibrium at .
Proof. Since system (SM1) is a refinement of system (SM2), it follows that z(0) > 0. Observe used that z(t) ≥ 0 for 0 ≤ t ≤t. This enables the derivation of upper estimates for the solutions y 0 (t) and y v (t) of system (SM1), respectively. In the second step, evaluation of these estimates at the switching timet yields upper estimates for y 0 (t) and y v (t) for all t ≥ 0, respectively. This is due to the following reasoning: On the one hand, it is z(t) ≥ 0 for all 0 ≤ t ≤t implying y 0 (t) ≥ 0 and y v (t) ≥ 0 for all 0 ≤ t ≤t. On the other hand, it is z(t) ≤ 0 for all t ≥t, and consequently y 0 (t) ≤ 0 and y v (t) ≤ 0 for all t ≥t.
Step 1: Since Consequently, the following estimates hold true: where u(y 0 (t), y v (t)) : Due to the comparison principle for systems of ODEs [22], it follows that Solving this system gives Observe that the solutions are strictly positive for t > 0.
Step 2: The switching timet is calculated using explicit solution formulae for solutions x(t) and y(t) of system (SM2) for z(t) ≥ 0, where s = √ 4b 2 + d 2 with initial conditions x 0 = f d + x 0 v and y 0 = f d + n. They are derived via elementary calculation. This implies Note thatt > 0 is well-defined. In particular,t only depends on the parameters b and d. As a consequence, a switch in the sign of z(t) occurs at timet. Combining steps 1 and 2 yields  For assessing practical identifiability of the parameters, additional estimations were performed using different optimisation methods and varying initialisations within a biologically reasonable range (see documentation of Mathematica (Version 9, Wolfram Research)). As a result, estimates and the widths of the corresponding confidence intervals did not significantly differ, indicating that the parameter values could uniquely be identified on the basis of the data.