Parameter Identification for a Model of Neonatal Fc Receptor-Mediated Recycling of Endogenous Immunoglobulin G in Humans

Salvage of endogenous immunoglobulin G (IgG) by the neonatal Fc receptor (FcRn) is implicated in many clinical areas, including therapeutic monoclonal antibody kinetics, patient monitoring in IgG multiple myeloma, and antibody-mediated transplant rejection. There is a clear clinical need for a fully parameterized model of FcRn-mediated recycling of endogenous IgG to allow for predictive modeling, with the potential for optimizing therapeutic regimens for better patient outcomes. In this paper we study a mechanism-based model incorporating nonlinear FcRn-IgG binding kinetics. The aim of this study is to determine whether parameter values can be estimated using the limited in vivo human data, available in the literature, from studies of the kinetics of radiolabeled IgG in humans. We derive mathematical descriptions of the experimental observations—timecourse data and fractional catabolic rate (FCR) data—based on the underlying physiological model. Structural identifiability analyses are performed to determine which, if any, of the parameters are unique with respect to the observations. Structurally identifiable parameters are then estimated from the data. It is found that parameter values estimated from timecourse data are not robust, suggesting that the model complexity is not supported by the available data. Based upon the structural identifiability analyses, a new expression for the FCR is derived. This expression is fitted to the FCR data to estimate unknown parameter values. Using these parameter estimates, the plasma IgG response is simulated under clinical conditions. Finally a suggestion is made for a reduced-order model based upon the newly derived expression for the FCR. The reduced-order model is used to predict the plasma IgG response, which is compared with the original four-compartment model, showing good agreement. This paper shows how techniques for compartmental model analysis—structural identifiability analysis, linearization, and reparameterization—can be used to ensure robust parameter identification.

1. Introduction

Immunoglobulin G (IgG) is the most abundant immunoglobulin (Ig) isotype in the circulation in humans, with a plasma concentration in healthy adults of 10–16 g l⁻¹ (1). Its high concentration is facilitated by the neonatal Fc receptor (FcRn), which binds IgG in intracellular endosomes and transports it to the plasma membrane to be returned to the circulation. A proportion of IgG molecules that are not bound by FcRn are degraded in lysosomes. In this way, FcRn continually protects a proportion of the circulating IgG from degradation. The recycling mechanism is saturable, such that at high plasma IgG concentrations a greater proportion of plasma IgG is degraded. Conversely, at depleted plasma IgG concentrations, a greater proportion is recycled and the half-life is extended beyond the normal 23 days (2).

Recent publications have drawn attention to the importance of FcRn-mediated recycling of endogenous IgG in the bone marrow cancer multiple myeloma. In multiple myeloma, clonal plasma cells secrete an excess of monoclonal Ig into the circulation. Patients undergoing therapy are primarily monitored by quantification of Ig in blood serum samples (3). Mills et al. (4) have suggested that FcRn-mediated recycling of IgG may result in different response rates between patients with IgG-producing multiple myeloma and patients with IgA-producing multiple myeloma. Yan et al. (5) have also suggested that FcRn-mediated recycling of endogenous IgG in patients with multiple myeloma may shorten the half-life of the therapeutic monoclonal antibody daratumumab. These studies highlight the need for a parameterized model of endogenous IgG kinetics for investigating these clinical scenarios.

Numerous mathematical models of IgG kinetics have been presented in the literature, mostly with the aim of describing the pharmacokinetics of therapeutic monoclonal antibodies (mAbs) that are also regulated by FcRn. Many of these models are therefore pharmacokinetic in nature: their parameter values are obtained from animal experiments and they may be physiologically-based, with up to around 10 organs explicitly represented in the model (6–14). Pharmacokinetic models developed for specific mAbs may not be generalizable to endogenous IgG if, for example, they include details such as binding of the mAb to its target. In addition, mAb disposition may be adequately described by linear models in many cases where the plasma concentration of therapeutic mAb is substantially smaller than the plasma concentration of endogenous IgG and the latter is constant (13, 14). However, the assumption of a constant plasma concentration of IgG is not always appropriate; for example, in multiple myeloma the plasma IgG concentration typically shows large changes during the course of therapy. Relative to a less complex model, the more complex model will usually provide a better fit to observed data. However, this alone does not imply that all the parameters in the complex model can be estimated consistently, nor does it imply that the underlying assumptions of the complex model are valid (15).

In this paper we study a mechanism-based model with a single plasma compartment, rather than separate plasma compartments for different organs, which is accessible to measurement in humans. The model, which has been previously shown by Kim et al. (16) and Hattersley (17), has in total four compartments, representing IgG in plasma, IgG in a peripheral compartment (representing less rapidly perfused tissues), unbound IgG in intracellular endosomes and IgG bound to FcRn receptors in intracellular endosomes. The IgG-FcRn interaction is represented by nonlinear receptor-ligand binding kinetics (6, 7, 9). With reliable parameter values for humans, it may be possible to use this model to predict the responses of plasma IgG under various clinical conditions.

The aim of this study is to determine whether the model parameter values can be obtained using the limited in vivo human data that are available in the literature. The data are from studies of the kinetics of administered small doses of radiolabeled IgG when the subject's endogenous IgG is in steady state. We consider two measured outputs: the timecourse of the proportion of an administered dose of radiolabeled IgG remaining in plasma and in the body; and the relationship between the fractional catabolic rate and the quantity of endogenous IgG in plasma. Structural identifiability analysis is performed with respect to these outputs and structurally identifiable parameters are estimated from the data.

2. Mathematical Models and Data Description

2.1. The Four Compartment Model

The model of IgG metabolism under study (16, 17) has four state variables, nine parameters, and an input function, I(t), representing the synthesis of IgG. The model equations are given by

$\begin{array}{ll}\begin{array}{c}{ẋ}_1(t)=-(k_{21}+k_{31})x_1(t)+k_{12}{x}_2(t)+k_{14}{x}_4(t)+I(t)\hfill \\ {ẋ}_2(t)=k_{21}{x}_1(t)-k_{12}{x}_2(t)\hfill \\ {ẋ}_3(t)=k_{31}{x}_1(t)-k_{03}{x}_3(t)-\frac{k_{\mathrm{\text{on}}}}{v_3}{x}_3(t)(R_{\mathrm{\text{tot}}}-x_4(t))+k_{\mathrm{\text{off}}}{x}_4(t)\hfill \\ {ẋ}_4(t)=\frac{k_{\mathrm{\text{on}}}}{v_3}{x}_3(t)(R_{\mathrm{\text{tot}}}-x_4(t))-(k_{14}+k_{\mathrm{\text{off}}})x_4(t),\hfill \end{array}\hfill & (1)\end{array}$

where x₁(t), x₂(t), x₃(t), and x₄(t) represent the quantities in μmol of IgG in plasma, IgG in a peripheral compartment, unbound IgG in endosomes and IgG bound to FcRn in endosomes, respectively. I(t) represents the rate of synthesis of IgG in µmol day⁻¹. The rate constants, k_ij, represent the rate of material flow from compartment j to compartment i, with the convention that 0 represents the environment outside the system. k_on and k_off are the receptor-ligand binding constants of IgG and FcRn. We denote the volumes of plasma, the peripheral compartment and the endosomes by v₁, v₂, and v₃, respectively. We assume a constant total (bound and unbound) quantity of FcRn, R_tot (6). This means that the quantity of unbound FcRn is represented by [R_tot − x₄(t)]. The state variables of the model and physiological interpretations of the parameters are summarized in Table 1. Note that all states and parameters can only take non-negative values. We refer to Figure 1 for a schematic of the model.

TABLE 1

Table 1. States and parameters of four-compartment model of IgG metabolism, with parameter values sourced in the literature.

FIGURE 1

Figure 1. Schematic of four-compartment model of IgG metabolism. Arrows represent material flow between compartments and hooked arrows represent nonlinear receptor-ligand binding. The fifth compartment shown (R_tot − x₄) represents unbound FcRn receptors and has been eliminated from the model equations. The arrow, labeled k₁₄, from the IgG-FcRn complex compartment (x₄) to the unbound FcRn receptor compartment (R_tot − x₄), represents internalization of FcRn receptors from the plasma membrane to the endosome, after releasing IgG.

When the production rate of IgG is constant, I(t) = I₀, the system has a stable equilibrium point given by

$\begin{array}{ll}\begin{array}{c}{\widehat{x}}_1=\frac{I_0(k_{03}{k}_{14}{v}_3+k_{03}{k}_{\mathrm{\text{off}}}{v}_3+k_{\mathrm{\text{on}}}{I}_0+k_{14}{k}_{\mathrm{\text{on}}}{R}_{\mathrm{\text{tot}}})}{k_{31}(k_{03}{v}_3(k_{14}+k_{\mathrm{\text{off}}})+k_{\mathrm{\text{on}}}{I}_0)}\hfill \\ {\widehat{x}}_2=\frac{k_{21}}{k_{12}}{\widehat{x}}_1\hfill \\ {\widehat{x}}_3=\frac{I_0}{k_{03}}\hfill \\ {\widehat{x}}_4=\frac{k_{\mathrm{\text{on}}}{I}_0{R}_{\mathrm{\text{tot}}}}{k_{03}{v}_3(k_{14}+k_{\mathrm{\text{off}}})+k_{\mathrm{\text{on}}}{I}_0}.\hfill \end{array}\hfill & (2)\end{array}$

A stability analysis for this equilibrium point is provided in the Supplementary Material.

2.2. In vivo Human Data From the Literature

The data available in the literature were obtained from tracer experiments. These studies entailed intravenous administration of a bolus dose of radiolabeled IgG (the tracer) and monitoring the proportion of the dose remaining in the blood and in the body over time. In this way the administered dose is distinguishable (by the experimenter) from the subject's own endogenous IgG. The quantity of administered tracer is small, so as not to perturb the steady state of the endogenous IgG. The purpose of tracer experiments is to enable observation of processes such as distribution and elimination undergone by the endogenous protein, whilst it is in steady state. The methods are described fully by Waldmann and Strober (21).

The data for an individual subject consist of the timecourse of the proportion of the injected dose of IgG remaining in plasma and the timecourse of the proportion of dose remaining in the body. In this paper we use the data from six such plots available in the literature. We refer to the individuals as subjects A–F. The timecourse data for subjects A–D are from Solomon et al. (18), for subject E from Waldmann and Terry (23), and for subject F from Waldmann and Strober (21). Several of the individuals have health conditions which may result in an increased or decreased plasma IgG concentration. Subjects A and C have IgG multiple myeloma and subject D has macroglobulinemia. Subjects B, E, and F are referred to as “normal” subjects. A spaghetti plot of the data is shown in Figure 2A. Subjects A and D show slower dynamics and subject C shows faster dynamics. The dynamics of IgG in these subjects is assumed to be described by the same model, as given by Equations (1), however they may have had altered production rates of IgG due to the diseases.

FIGURE 2

Figure 2. (A) Spaghetti plot of the proportion of administered IgG remaining in plasma (blue) and the body (red) in six subjects; data from Waldmann and Strober (21), Solomon et al. (18), and Waldmann and Terry (23). (B) Plasma concentration dependence of the fractional catabolic rate (FCR); redrawn from Waldmann and Strober (21), with permission from S. Karger AG, Basel.

Also available in the literature is a plot of the fractional catabolic rate (FCR) vs. the subject's plasma concentration of endogenous IgG, obtained from a group of individuals with a range of plasma IgG concentrations (21). The FCR is defined as the elimination rate of IgG as a fraction of the quantity of IgG in plasma. In practice the FCR is calculated from the rate at which the tracer dose leaves the body at time t divided by the proportion of tracer dose remaining in plasma at time t. The relationship between the FCR and the timecourse data is described further in section 2.6. A plot of the FCR vs. the plasma concentration of endogenous IgG for 41 individuals provided by Waldmann and Strober (21) is shown in Figure 2B. Each data point was derived from the timecourse data of an individual subject. All of the data described in this section were extracted from plots in the literature using the Digitizer tool in OriginPro 2016 (24).

2.3. Nonlinear Model of Coupled Tracer and Endogenous IgG Dynamics

The administered tracer and the endogenous IgG are assumed to be indistinguishable by the human body, that is they exhibit identical kinetic (input/output) behavior—a standard assumption in tracer studies (25). We therefore assume that the kinetics of both tracer and endogenous IgG are described by the model given by Equations (1). From Equations (1), letting x_i(t) = x_i,T(t) + x_i,E(t), where x_i,T(t) and x_i,E(t) denote the quantities in μmol in compartment i of radiolabeled and endogenous IgG, respectively, gives

$\begin{array}{l}{\dot{x}}_{1,\text{T}}(t)=-(k_{21}+k_{31})x_{1,\text{T}}(t)+k_{12}{x}_{2,\text{T}}(t)+k_{14}{x}_{4,\text{T}}(t)\\ {\dot{x}}_{2,\text{T}}(t)=k_{21}{x}_{1,\text{T}}(t)-k_{12}{x}_{2,\text{T}}(t)\\ {\dot{x}}_{3,\text{T}}(t)=k_{31}{x}_{1,\text{T}}(t)-k_{03}{x}_{3,\text{T}}(t)-\frac{k_{\text{on}}}{v_3}{x}_{3,\text{T}}(t)(R_{\text{tot}}-x_{4,\text{E}}(t)\\ -x_{4,\text{T}}(t))+k_{\text{off}}{x}_{4,\text{T}}(t)\\ {\dot{x}}_{4,\text{T}}(t)=\frac{k_{\text{on}}}{v_3}{x}_{3,\text{T}}(t)\left(R_{\text{tot}}-x_{4,\text{E}}(t)-x_{4,\text{T}}(t)\right)-(k_{14}+k_{\text{off}})x_{4,\text{T}}(t)\\ {\dot{x}}_{1,\text{E}}(t)=-(k_{21}+k_{31})x_{1,\text{E}}(t)+k_{12}{x}_{2,\text{E}}(t)+k_{14}{x}_{4,\text{E}}(t)+I_{\text{E}}\\ {\dot{x}}_{2,\text{E}}(t)=k_{21}{x}_{1,\text{E}}(t)-k_{12}{x}_{2,\text{E}}(t)\\ {\dot{x}}_{3,\text{E}}(t)=k_{31}{x}_{1,\text{E}}(t)-k_{03}{x}_{3,\text{E}}(t)-\frac{k_{\text{on}}}{v_3}{x}_{3,\text{E}}(t)(R_{\text{tot}}-x_{4,\text{E}}(t)n\\ -x_{4,\text{T}}(t))+k_{\text{off}}{x}_{4,\text{E}}(t)\\ {\dot{x}}_{4,\text{E}}(t)=\frac{k_{\text{on}}}{v_3}{x}_{3,\text{E}}(t)\left(R_{\text{tot}}-x_{4,\text{E}}(t)-x_{4,\text{T}}(t)\right)\\ -(k_{14}+k_{\text{off}})x_{4,\text{E}}(t).& (3)\end{array}$

I_E (μmol day⁻¹) represents the production rate of endogenous IgG, which is assumed constant. All other parameters are defined in Table 1.

The dose of tracer administered at time t = 0 days is treated as a non-zero initial condition for x_1,T(t). Tracer is administered to the plasma compartment only; therefore the initial conditions of the remaining tracer compartments are zero. The endogenous IgG is assumed to be in steady state throughout the experiment, such that the initial conditions of the endogenous IgG are given by the steady states in Equations (2), with I₀ = I_E. In summary, the initial conditions are given by

$\begin{array}{ll}\begin{array}{c}{x}_{1,\mathrm{\text{T}}}(0)=D\hfill \\ x_{2,\mathrm{\text{T}}}(0)=x_{3,\mathrm{\text{T}}}(0)=x_{4,\mathrm{\text{T}}}(0)=0\hfill \\ x_{1,\mathrm{\text{E}}}(0)={\widehat{x}}_1\hfill \\ x_{2,\mathrm{\text{E}}}(0)={\widehat{x}}_2\hfill \\ x_{3,\mathrm{\text{E}}}(0)={\widehat{x}}_3\hfill \\ x_{4,\mathrm{\text{E}}}(0)={\widehat{x}}_4,\hfill \end{array}\hfill & (4)\end{array}$

where ${\widehat{x}}_i$ is the steady state quantity of endogenous IgG in compartment i, given by Equations (2), and D (μmol) is the administered dose of tracer.

The experimenter observes the proportion of the dose remaining in plasma [denoted by y₁(t)] and in the body [denoted by y₂(t)] during the experiment. The observation functions are thus given by

$\begin{array}{ll}\begin{array}{c}{y}_1(t)=\frac{x_{1,\mathrm{\text{T}}}(t)}{D}\hfill \\ y_2(t)=\frac{x_{1,\mathrm{\text{T}}}(t)+x_{2,\mathrm{\text{T}}}(t)+x_{3,\mathrm{\text{T}}}(t)+x_{4,\mathrm{\text{T}}}(t)}{D}.\hfill \end{array}\hfill & (5)\end{array}$

2.4. Linearized Model of Tracer Dynamics

Provided that the administered dose of tracer is sufficiently small, the tracer kinetics can be approximated using the Taylor series expansion of the model state about the equilibrium point. In this way a linear model of the experiment, valid in a neighborhood of the equilibrium point, is derived. Our derivation is provided in the Supplementary Material. The derivation of a linearized model for tracer dynamics from a general compartmental model is provided by Anderson (26).

The linear equations describing the tracer kinetics are given by

$\begin{array}{ll}\begin{array}{c}{ẋ}_{1,\mathrm{\text{T}}}(t)=-(k_{21}+k_{31})x_{1,\mathrm{\text{T}}}(t)+k_{12}{x}_{2,\mathrm{\text{T}}}(t)+k_{14}{x}_{4,\mathrm{\text{T}}}(t)\hfill \\ {ẋ}_{2,\mathrm{\text{T}}}(t)=k_{21}{x}_{1,\mathrm{\text{T}}}(t)-k_{12}{x}_{2,\mathrm{\text{T}}}(t)\hfill \\ {ẋ}_{3,\mathrm{\text{T}}}(t)=k_{31}{x}_{1,\mathrm{\text{T}}}(t)-k_{03}{x}_{3,\mathrm{\text{T}}}(t)-k_{43}{x}_{3,\mathrm{\text{T}}}(t)+k_{34}{x}_{4,\mathrm{\text{T}}}(t)\hfill \\ {ẋ}_{4,\mathrm{\text{T}}}(t)=k_{43}{x}_{3,\mathrm{\text{T}}}(t)-(k_{14}+k_{34})x_{4,\mathrm{\text{T}}}(t)\hfill \end{array}\hfill & (6)\end{array}$

where x_1,T(t), x_2,T(t), x_3,T(t), and x_4,T(t) represent the quantities of radiolabeled IgG in the central compartment, in the peripheral compartment, unbound in intracellular endosomes, and bound to FcRn in intracellular endosomes, respectively. The new parameters k₃₄ and k₄₃ are given by

$\begin{array}{ll}\begin{array}{c}{k}_{34}=k_{\mathrm{\text{off}}}\hfill \\ k_{43}=\frac{k_{\mathrm{\text{on}}}(R_{\mathrm{\text{tot}}}-{\widehat{x}}_4)}{v_3}=\frac{k_{\mathrm{\text{on}}}{R}_{\mathrm{\text{tot}}}{k}_{03}(k_{14}+k_{\mathrm{\text{off}}})}{I_{\mathrm{\text{E}}}{k}_{\mathrm{\text{on}}}+k_{03}{v}_3(k_{14}+k_{\mathrm{\text{off}}})}.\hfill \end{array}\hfill & (7)\hfill \end{array}$

All other parameters are defined in Table 1. The initial conditions are given by the first two equations of Equations (4) and the observation functions are given by Equations (5).

2.5. Comparison of Nonlinear Model and Linearized Model for Large Tracer Doses

The linearization of the model of timecourse observations relies on the assumption of a sufficiently small dose of tracer, such that the endogenous IgG can be assumed to remain in steady state. A typical tracer dose is between 3 · 10⁻³ and 7 · 10⁻³ μmol (18). Simulations of the quantity of tracer in each compartment are shown in Figure 3. In Figure 3A, a dose of D = 1 µmol is assumed and in Figure 3B, a dose of D = 100 µmol is assumed. The value of 1 μmol was chosen to show that the linear model is a valid approximation of the nonlinear model, even when the dose is more than 100 times typical tracer doses. The extremely large value of 100 μmol was chosen specifically to show the dynamics of the linearized model when it is not a valid approximation of the nonlinear model. The parameter values in Table 1 are used. A normal IgG synthesis rate of I_E = 15 µmol day⁻¹ was used; however the linearized model was still valid for D = 1 µmol when comparatively very small values of I_E were used. We find that, for a dose of 1 μmol and the particular parameter values used, the linearized model is a valid approximation of the full nonlinear model over a 25-day simulated time course. When the dose is increased to 100 μmol, the assumption that the steady state is not perturbed by the administered dose no longer holds and the two models give different simulation results for the quantities of tracer.

FIGURE 3

Figure 3. Simulations of the quantities of tracer in each compartment after administration at t = 0 days, for a tracer dose of (A) 1 μmol and (B) 100 μmol. The nonlinear model (Equations 3–4) is represented by solid lines and the linearized model (Equations 6) by dashed lines. The linearized model is valid for the smaller dose but not for the larger dose. Note the different scales for x_1,T(t) and x_2,T(t), and x_3,T(t) and x_4,T(t), respectively.

2.6. Fractional Catabolic Rate

We recall that the FCR (μmol day⁻¹) is defined as the elimination rate of IgG as a fraction of the quantity of IgG in plasma and can be defined with respect to the tracer or with respect to the endogenous IgG. The FCR with respect to the tracer is therefore given by

$\begin{array}{ll}{\mathrm{\text{FCR}}}_{\mathrm{\text{T}}}(t)=\frac{k_{03}{x}_{3,\mathrm{\text{T}}}(t)}{x_{1,\mathrm{\text{T}}}(t)},\hfill & (8)\hfill \end{array}$

where x_3,T(t) and x_1,T(t) are given by the solution of Equations (6).

Whilst a single value of the FCR is measured for an individual subject (see Figure 2B), in actuality FCR_T(t) is not constant, as shown by the dependence on time in Equation (8). A simulation of FCR_T(t) during the experiment is shown in Figure 4. After around day 5, for the particular parameter values used, FCR_T(t) approaches a steady state value, which is denoted here by FCR_T,∞:

$\begin{array}{ll}{\mathrm{\text{FCR}}}_{\mathrm{\text{T}},\infty }={\displaystyle \underset{t\to \infty }{lim}}\frac{k_{03}{x}_{3,\mathrm{\text{T}}}(t)}{x_{1,\mathrm{\text{T}}}(t)}.\hfill & (9)\hfill \end{array}$

Solving Equations (6) gives

$\begin{array}{l}{x}_{i,\mathrm{\text{T}}}(t)=A_{i1}exp({\lambda }_{1}t)+A_{i2}exp({\lambda }_{2}t)+A_{i3}exp({\lambda }_{3}t)\hfill \\ +A_{i4}exp({\lambda }_{4}t),i=1,\dots ,4,\hfill & (10)\hfill \end{array}$

where A_ij and λ_j (j = 1, …, 4) are expressions in terms of the model parameters and |λ₁| > |λ₂| > |λ₃| > |λ₄|. After sufficient time, x_{i, T}(t) can be approximated by A_i4 exp(λ₄t); thus, FCR_T,∞ is given by

$\begin{array}{ll}{\mathrm{\text{FCR}}}_{\mathrm{\text{T}},\infty }=k_{03}\frac{A_{34}exp({\lambda }_{4}t)}{A_{14}exp({\lambda }_{4}t)}=k_{03}\frac{A_{34}}{A_{14}}.\hfill & (11)\hfill \end{array}$

The expressions for A₃₄ and A₁₄ in terms of the model parameters are extremely long. The Mathematica (27) code for generating the expressions for A_ij and λ_j is provided in the Supplementary Material.

FIGURE 4

Figure 4. Simulation of FCR_T(t) given by Equations (12) and (8), for the parameter values in Table 1 and dose D = 0.01 μmol.

Noting that there is only elimination from the system and no input for t > 0, FCR_T(t) is equal to the rate of change of radiolabeled IgG in all compartments, divided by the quantity of radiolabeled IgG in plasma:

$\begin{array}{ll}{\mathrm{\text{FCR}}}_{\mathrm{\text{T}}}(t)=\frac{-({ẋ}_{1,\mathrm{\text{T}}}(t)+{ẋ}_{2,\mathrm{\text{T}}}(t)+{ẋ}_{3,\mathrm{\text{T}}}(t)+{ẋ}_{4,\mathrm{\text{T}}}(t))}{x_{1,\mathrm{\text{T}}}(t)}=\frac{-{ẏ}_2(t)}{y_1(t)}.& (12)\end{array}$

From Equation (12), it can be seen that FCR_T(t) is equal to the slope of the observation y₂(t) divided by y₁(t), showing how FCR_T(t) can be obtained from the observations y₁(t) and y₂(t). In practice, the experimenter obtains a value for FCR_T(t_N), where t_N is a time toward the end of the experiment, such that FCR_T(t_N) can be assumed a close approximation of FCR_T,∞. Henceforth, the quantity obtained from experiments, FCR_T(t_N), is referred to simply as FCR_T.

It is also possible to derive an expression for the FCR with respect to the endogenous IgG, FCR_E. If the endogenous IgG is assumed to remain in steady state, then from the definition of the FCR,

$\begin{array}{ll}{\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}=\frac{k_{03}{\widehat{x}}_3}{{\widehat{x}}_1},& (13)\end{array}$

where ${\widehat{x}}_1$ and ${\widehat{x}}_3$ are the quantities of IgG in compartments 1 and 3 in steady state, given by Equations (2). Substituting the expression for ${\widehat{x}}_3$ from Equations (2) into Equation (13), eliminating I₀ in favor of ${\widehat{x}}_1$ using the first equation of Equations (2), and setting ${\widehat{x}}_1=x_{1,\mathrm{\text{E}}}$, gives the following expression for the FCR_E in terms of the quantity of IgG in plasma, x_1,E:

$\begin{array}{l}{\text{FCR}}_{\text{E}}=\frac{1}{2k_{\text{on}}{x}_{1,\text{E}}}(k_{31}{k}_{\text{on}}{x}_{1,\text{E}}-k_{14}{k}_{\text{on}}{R}_{\text{tot}}-k_{03}{k}_{14}{v}_3\\ -k_{03}{k}_{\text{off}}{v}_3+\{4k_{03}{k}_{31}(k_{14}+k_{\text{off}})k_{\text{on}}{x}_{1,\text{E}}{v}_3\\ +(-k_{31}{k}_{\text{on}}{x}_{1,\text{E}}+k_{14}{k}_{\text{on}}{R}_{\text{tot}}+k_{03}{k}_{14}{v}_3\\ +{{k_{03}{k}_{\text{off}}{v}_3)}^2\}}^{1/2}).& (14)\end{array}$

3. Results

3.1. Parameter Identification Using Tracer Timecourse Data

In this section we investigate whether it is possible to estimate unknown model parameter values by fitting the linear approximation described in section 2.4 to the timecourse data described in section 2.2. Firstly, a structural identifiability analysis is performed. Parameter values are then estimated from the data by fitting the linearized model described by Equations (6) to the data.

3.1.1. Structural Identifiability Analysis

Structural identifiability addresses the question of whether model parameters can be uniquely identified from the available observations, under the assumption of the availability of ideal (i.e. noise-free) and continuous observational data. Here we determine which of the model parameters are structurally uniquely identifiable from the observations y₁(t) and y₂(t), given by Equations (4–6). The unknown parameter vector is given by $\theta ={(k_{21},k_{31},k_{12},k_{14},k_{03},k_{43},k_{34})}^{\mathrm{\text{T}}}$.

The transfer function method is used (28). To apply this approach the system described by Equations (4–6) is re-written in vector-matrix notation as

$\begin{array}{ll}\begin{array}{c}{\dot{\text{x}}}_{\mathrm{\text{T}}}(t,\theta )=\text{A}(\theta ){\text{x}}_{\mathrm{\text{T}}}(t)+\text{B}(\theta )u(t)\hfill \\ {\text{x}}_{\mathrm{\text{T}}}(0,\theta )=0\hfill \\ \text{y}(t,\theta )=\text{C}(\theta ){\text{x}}_{\mathrm{\text{T}}}(t),\hfill \end{array}\hfill & (15)\end{array}$

where ${\text{x}}_{\mathrm{\text{T}}}(t,\theta )={(x_{1,\mathrm{\text{T}}}(t),x_{2,\mathrm{\text{T}}}(t),x_{3,\mathrm{\text{T}}}(t),x_{4,\mathrm{\text{T}}}(t))}^{\mathrm{\text{T}}}$ and $\text{y}(t,\theta )={(y_1(t),y_2(t))}^{\mathrm{\text{T}}}$ are column vectors representing the state vector and the observation vector, respectively, and u(t) represents the single input to the system, an impulse at time t = 0, given by u(t) = δ(t). A(θ) is a 4 × 4 matrix, B(θ) is a column vector and C(θ) is a 2 × 4 matrix. A(θ), B(θ), and C(θ) are given by

$\begin{array}{ll}\begin{array}{c}\text{A}\left(\theta \right)=\left(\begin{array}{cccc}-\left(k_{21}+k_{31}\right)& k_{12}& 0& k_{14}\\ k_{21}& -k_{12}& 0& 0\\ k_{31}& 0& -\left(k_{03}+k_{43}\right)& k_{34}\\ 0& 0& k_{43}& -\left(k_{14}+k_{34}\right)\end{array}\right),\\ \text{B}\left(\theta \right)=\left(\begin{array}{c}D\\ 0\\ 0\\ 0\end{array}\right),\hfill \\ \text{C}\left(\theta \right)=\left(\begin{array}{cccc}\frac{1}{D}& 0& 0& 0\\ \frac{1}{D}& \frac{1}{D}& \frac{1}{D}& \frac{1}{D}\end{array}\right).\hfill \end{array}& (16)\end{array}$

Note that the administration of a bolus dose of size D is now represented as an impulse at time t = 0, rather than a non-zero initial condition, such that ${\text{x}}_{\mathrm{\text{T}}}(0,\theta )={(0,0,0,0)}^{\mathrm{\text{T}}}$.

Taking Laplace transforms of Equations (15), the input-output relation is given by Y(s) = G(s)U(s), where G(s) is the transfer function matrix, given by G(s) = C(θ)(sI − A(θ))⁻¹ B(θ), where I is the 4 × 4 identity matrix. G(s) has two elements, corresponding to the two observed outputs, which are given by

$\begin{array}{ll}\begin{array}{c}{G}_1(s)=\frac{{\varphi }_1+{\varphi }_{2}s+{\varphi }_3{s}^2+s^3}{{\varphi }_4+{\varphi }_{5}s+{\varphi }_6{s}^2+{\varphi }_7{s}^3+s^4}\hfill \\ G_2(s)=\frac{{\varphi }_8+{\varphi }_{9}s+{\varphi }_{10}{s}^2+s^3}{{\varphi }_{11}+{\varphi }_{12}s+{\varphi }_{13}{s}^2+{\varphi }_{14}{s}^3+s^4},\end{array}\hfill & (17)\end{array}$

where the coefficients of s, $\Phi (\theta )={({\varphi }_1(\theta ),{\varphi }_2(\theta ),\dots ,{\varphi }_{14}(\theta ))}^{\mathrm{\text{T}}}$, are nonlinear expressions in the parameters. The coefficients of s, Φ(θ), are given by

$\begin{array}{l}{\varphi }_1(\theta )=k_{12}\left(k_{03}\left(k_{14}+k_{34}\right)+k_{14}{k}_{43}\right)\\ {\varphi }_2(\theta )=k_{03}\left(k_{12}+k_{14}+k_{34}\right)+k_{14}{k}_{43}+k_{12}\left(k_{14}+k_{34}+k_{43}\right)\\ {\varphi }_3(\theta )=k_{03}+k_{12}+k_{14}+k_{34}+k_{43}\\ {\varphi }_4(\theta )={\varphi }_{11}(\theta )=k_{03}{k}_{12}{k}_{31}(k_{14}+k_{34})\\ {\varphi }_5(\theta )={\varphi }_{12}(\theta )=k_{03}((k_{21}+k_{31})(k_{14}+k_{34})+k_{12}(k_{14}+k_{31}\\ +k_{34}))+k_{14}{k}_{21}{k}_{43}+k_{12}(k_{14}(k_{31}+k_{43})+k_{31}(k_{34}+k_{43}))\\ {\varphi }_6(\theta )={\varphi }_{13}(\theta )=k_{14}{k}_{21}+k_{14}{k}_{31}+k_{21}{k}_{34}+k_{31}{k}_{34}+k_{03}(k_{12}\\ +k_{14}+k_{21}+k_{31}+k_{34})+k_{14}{k}_{43}+k_{21}{k}_{43}+k_{31}{k}_{43}+k_{12}\\ (k_{14}+k_{31}+k_{34}+k_{43})\\ {\varphi }_7(\theta )={\varphi }_{10}(\theta )={\varphi }_{14}(\theta )=k_{03}+k_{12}+k_{14}+k_{21}+k_{31}\\ +k_{34}+k_{43}\\ {\varphi }_8(\theta )=k_{03}(k_{12}+k_{21})(k_{14}+k_{34})+k_{14}{k}_{21}{k}_{43}\\ +k_{12}(k_{14}(k_{31}+k_{43})+k_{31}(k_{34}+k_{43}))\\ {\varphi }_9(\theta )=k_{14}{k}_{21}+k_{14}{k}_{31}+k_{21}{k}_{34}+k_{31}{k}_{34}+k_{03}(k_{12}+k_{14}\\ +k_{21}+k_{34})+k_{14}{k}_{43}+k_{21}{k}_{43}+k_{31}{k}_{43}\\ +k_{12}(k_{14}+k_{31}+k_{34}+k_{43}).& (18)\end{array}$

The coefficients Φ(θ) are unique with respect to the input-output relationship represented by the transfer function. Introducing an alternative parameter vector, $\overline{\theta }={({\overline{k}}_{21},{\overline{k}}_{31},{\overline{k}}_{12},{\overline{k}}_{14},{\overline{k}}_{03},{\overline{k}}_{43},{\overline{k}}_{34})}^{\mathrm{\text{T}}}$, and equating $\Phi (\theta )=\Phi (\overline{\theta })$, the resulting set of simultaneous equations is solved for θ using the Solve function in Mathematica (27). The only solution is $\theta =\overline{\theta }$; therefore all of the parameters in θ are structurally uniquely identifiable.

3.1.2. Parameter Estimation

The parameter vector $\theta ={(k_{21},k_{31},k_{12},k_{14},k_{03},k_{43},k_{34})}^{\mathrm{\text{T}}}$ was estimated for each subject using unweighted least squares, by fitting the timecourse data described in section 2.2. The “true” parameter vector for an individual is denoted by θ₀. For an individual subject it is assumed that y_i(t, θ₀), i = 1, 2, is observed with error at measurement times $t_1^{(i)},\dots ,t_{N_i}^{(i)},i=1,2$, where $t_1^{(1)}=t_1^{(2)}=0$. The observed (with error) values of y_i(t, θ₀), i = 1, 2, are now denoted by ${\tilde{y}}_i(t_j^{(i)},{\theta }_0)$ for i = 1, 2 and j = 1, …, N_i. Both outputs y₁ and y₂ were fitted simultaneously, therefore the cost functional for θ is given by

$\begin{array}{ll}J({\theta }_0,\theta )={\displaystyle \sum _{i=1}^2}{J}_i({\theta }_0,\theta ),& (19)\end{array}$

where

$\begin{array}{ll}{J}_i({\theta }_0,\theta )={\displaystyle \sum _{j=1}^{N_i}{\left({\tilde{y}}_i(t_j^{(i)},{\theta }_0)-y_i(t_j^{(i)},\theta )\right)}^2}.& (20)\end{array}$

Differential evolution was implemented using the NonlinearModelFit function in Mathematica (27). The differential evolution algorithm was chosen because there is little information available about the parameters, in particular the parameters k₁₄, k₀₃, k₄₃, and k₃₄. Differential evolution is a stochastic, global minimization algorithm that does not require the user to specify initial guesses for the parameter values (29). All parameters were constrained to be positive. The maximum number of iterations was set to 5,000, which was sufficient for the algorithm to converge in all cases. In differential evolution an initial population of parameter vectors is generated randomly. The algorithm was run for each subject's data with integer seeds for the pseudorandom number generator between 1 and 10; thus 10 estimates for θ were obtained for each subject.

Differential evolution maintains a population of parameter vectors which evolves iteratively. For each new generation of the algorithm, a mutant and trial vector are produced from the current generation and the trial vector is compared with a target vector from the current generation. Either the target or trial vector is selected to move forward to the new generation based on which has the smallest value of the cost function to be minimized. The scaling factor (SF) is used to produce the mutant vector and generally a larger value of SF means a broader search of the parameter space. The crossover probability (CR) is the probability that each element of the mutant vector is used to produce the trial vector, rather than the corresponding element of the target vector. SF and CR were tuned by trial and error for each subject. The settings F = 0.5 and CR = 0.9 were tried initially, as recommended by Storn and Price (29) for faster convergence. For subjects C and F the settings were adjusted to F = 0.7, for a broader search of the parameter space, and CR = 0.95, to speed convergence. The settings for the differential evolution algorithm are given in Table 2.

TABLE 2

Table 2. Settings for differential evolution.

Each run of the algorithm, with a unique seed for the pseudorandom number generator, can produce unique parameter estimates; it is therefore recommended to perform multiple runs with unique, randomly chosen starting populations of parameter vectors (29). The parameter estimates and root mean square error (RMSE) for each run and each subject are tabulated in Table 3. The parameter estimates from multiple runs should be close to one another so that they can be averaged (29, 30); however, in some cases, the different runs give very different parameter estimates, implying that the algorithm has difficulty finding the global minimum and that there may be many local minima. It is therefore not certain that the global minimum has been found for each subject. It is also possible that certain parameters are highly correlated, such that different parameter vectors produce very similar model outputs. This is reflected in the diversity of parameter vectors obtained within subjects using differential evolution.

TABLE 3

Table 3. Parameter values estimated from timecourse data.

In some cases the model parameters are estimated to be zero, or very close to zero, for example k₃₄ for subject B, k₁₂ and k₁₄ for subject C, k₃₄ for subject D, and k₂₁ and k₁₄ for subject E. For each of these subjects the data can be well represented by a reduced model in which either IgG-FcRn binding is irreversible (k₃₄ = 0), there is no transfer from the peripheral compartment to plasma (k₁₂ = 0) or vice versa (k₂₁ = 0), or bound IgG molecules are not recycled into plasma (k₁₄ = 0). This result suggests that the model complexity is not supported by the available data.

The data and the model outputs using the parameter estimates in Table 3 are plotted in Figure 5. In each panel of Figure 5, the model outputs y₁(t) and y₂(t) are plotted for each of the estimated parameter vectors from 10 runs. The model outputs are very similar for all of the estimated parameter vectors for an individual. For some subjects there are small but noticeable differences between the fits, for example: in the first and last 5 days of y₂(t) for subject A; in the first 2 days of y₁(t) for subject B; for all of y₁(t) and the latter part of y₂(t) for subject C; between days 2 and 6 for y₁(t) and the initial 2 days of y₂(t) for subject E; and the first 10 days and final 5 days of y₂(t) for subject F. The similarity between the outputs for the parameter estimates obtained across different runs is shown by the similar values of RMSE within each subject. The model appears to fit the data reasonably well and in some subjects extremely well.

FIGURE 5

Figure 5. Timecourse data [y₁(t) blue circles; y₂(t) red circles] and model fits [y₁(t) blue line; y₂(t) red line] for (A–F) subjects A–F.

The results of the multiple runs of differential evolution show that in many cases, highly different parameter vectors produce very similar model outputs. The spread of the parameter estimates from multiple runs is conveyed using the coefficient of variation (CV), that is, the standard deviation of the estimates of a parameter from 10 runs, divided by the mean of those estimates. The CV is tabulated in Table 4. For some parameters and subjects, the estimates for the parameters have a small CV, for example the first four parameters for subject A and parameter k₁₂ for subject D. In other instances however the CV is much larger, reflecting the highly different estimates obtained for these parameters. The similarly high quality fits produced by diverse parameter vectors implies that, whilst the parameters are structurally identifiable, they are not all practically identifiable for the quality of data that are available.

TABLE 4

Table 4. Coefficients of variation of parameter estimates obtained from 10 runs of differential evolution, for each of subjects A–F.

3.2. Parameter Identification Using Fractional Catabolic Rate Data

Authors who have studied a two-compartment model of IgG metabolism have previously estimated parameters from FCR vs. plasma IgG concentration data (16, 21). In this section we investigate whether it is possible to estimate parameters of the four-compartment model from these data, which are described in section 2.2. In section 2.6 two expressions for the FCR were introduced: the FCR of the tracer (Equation 11) and the FCR of the endogenous IgG in steady state (Equation 14). In practice FCR_T is measured; however it is difficult to obtain a closed form expression for FCR_T. In contrast, we can easily obtain an expression for FCR_E in terms of the model parameters and the quantity of endogenous IgG in plasma, x_1,E, as given by Equation (14). In this section model parameters are estimated by fitting the expression for FCR_E vs. x_1,E Equation (14) to the FCR_T vs. x_1,E data. It is assumed that FCR_E is a good approximation to FCR_T and the parameter estimates are validated in section 3.2.3 using synthetic data.

3.2.1. Structural Identifiability Analysis

The relationship between FCR_E and x_1,E is given by Equation (14). Given that the parameters k_on and v₃ only appear in the model (Equations 3) as the ratio k_on/v₃, we re-write Equation (14), defining ϕ₁ = k_on/v₃, giving

$\begin{array}{l}{\text{FCR}}_{\text{E}}=\frac{1}{2{\varphi }_1{x}_{1,\text{E}}}(k_{31}{\varphi }_1{x}_{1,\text{E}}-k_{14}{\varphi }_1{R}_{\text{tot}}-k_{03}{k}_{14}-k_{03}{k}_{\text{off}}\\ +\sqrt{4k_{03}{k}_{31}\left(k_{14}+k_{\text{off}}\right){\varphi }_1{x}_{1,\text{E}}+{\left(-k_{31}{\varphi }_1{x}_{1,\text{E}}+k_{14}{\varphi }_1{R}_{\text{tot}}+k_{03}\left(k_{14}+k_{\text{off}}\right)\right)}^2}).& (21)\end{array}$

We wish to know whether the parameter vector $\varphi ={({\varphi }_1,k_{31},k_{14},R_{\mathrm{\text{tot}}},k_{03},k_{\mathrm{\text{off}}})}^{\mathrm{\text{T}}}$ is structurally identifiable with respect to the relationship in Equation (21). The structural identifiability problem amounts to determining whether there exists an alternative parameter vector $\overline{\varphi }={({\overline{\varphi }}_1,{\overline{k}}_{31},{\overline{k}}_{14},{\overline{R}}_{\mathrm{\text{tot}}},{\overline{k}}_{03},{\overline{k}}_{\mathrm{\text{off}}})}^{\mathrm{\text{T}}}$ such that ${\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}(x_{1,\mathrm{\text{E}}},\varphi )={\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}(x_{1,\mathrm{\text{E}}},\overline{\varphi })$.

From Equations (13) and (2),

$\begin{array}{ll}{\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}=\frac{I_0}{{\widehat{x}}_1}.& (22)\end{array}$

I₀ is given in terms of ${\widehat{x}}_1$ by the solution of the following quadratic equation, obtained by rearranging the first equation of Equations (2) and setting ϕ₁ = k_on/v₃:

$\begin{array}{l}-{\varphi }_1{I}_0^2+(-k_{03}(k_{14}+k_{\mathrm{\text{off}}})+{\varphi }_1(k_{31}{\widehat{x}}_1-k_{14}{R}_{\mathrm{\text{tot}}}))I_0\\ +k_{03}{k}_{31}(k_{14}+k_{\mathrm{\text{off}}}){\widehat{x}}_1=0.& (23)\end{array}$

Substituting ${\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}{\widehat{x}}_1$ in place of I₀ and setting ${\widehat{x}}_1=x_{1,\mathrm{\text{E}}}$ gives the following quadratic equation in FCR_E:

$\begin{array}{l}-{\varphi }_1{x}_{1,\mathrm{\text{E}}}^2{\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}^2+(-k_{03}(k_{14}+k_{\mathrm{\text{off}}})+{\varphi }_1(k_{31}{x}_{1,\mathrm{\text{E}}}-k_{14}{R}_{\mathrm{\text{tot}}}))\\ x_{1,\mathrm{\text{E}}}{\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}+k_{03}{k}_{31}(k_{14}+k_{\mathrm{\text{off}}})x_{1,\mathrm{\text{E}}}=0.& (24)\end{array}$

Dividing Equation (24) throughout by the coefficient of ${\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}^2$ gives

$\begin{array}{l}{\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}^2+\left(\frac{k_{03}\left(k_{14}+k_{\mathrm{\text{off}}}\right)-k_{31}{\varphi }_1{x}_{1,\mathrm{\text{E}}}+k_{14}{\varphi }_1{R}_{\mathrm{\text{tot}}}}{{\varphi }_1{x}_{1,\mathrm{\text{E}}}}\right){\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}\\ -\frac{k_{03}{k}_{31}\left(k_{14}+k_{\mathrm{\text{off}}}\right)}{{\varphi }_1{x}_{1,\mathrm{\text{E}}}}=0.& (25)\end{array}$

The expression for FCR_E given by Equation (21) is one of the two solutions of Equation (25). We therefore wish to know whether there exists an alternative parameter vector $\overline{\varphi }$ such that,

$\begin{array}{ll}\begin{array}{c}{\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}^2+\left(\frac{k_{03}\left(k_{14}+k_{\mathrm{\text{off}}}\right)-k_{31}{\varphi }_1{x}_{1,\mathrm{\text{E}}}+k_{14}{\varphi }_1{R}_{\mathrm{\text{tot}}}}{{\varphi }_1{x}_{1,\mathrm{\text{E}}}}\right){\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}\hfill \\ -\frac{k_{03}{k}_{31}\left(k_{14}+k_{\mathrm{\text{off}}}\right)}{{\varphi }_1{x}_{1,\mathrm{\text{E}}}}\hfill \\ ={\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}^2+\left(\frac{{\overline{k}}_{03}\left({\overline{k}}_{14}+{\overline{k}}_{\mathrm{\text{off}}}\right)-{\overline{k}}_{31}{\overline{\varphi }}_1{x}_{1,\mathrm{\text{E}}}+{\overline{k}}_{14}{\overline{\varphi }}_1{\overline{R}}_{\mathrm{\text{tot}}}}{{\overline{\varphi }}_1{x}_{1,\mathrm{\text{E}}}}\right)\hfill \\ {\mathrm{\text{FCR}}}_{\mathrm{\text{E}}}-\frac{{\overline{k}}_{03}{\overline{k}}_{31}\left({\overline{k}}_{14}+{\overline{k}}_{\mathrm{\text{off}}}\right)}{{\overline{\varphi }}_1{x}_{1,\mathrm{\text{E}}}}.\hfill \end{array}\hfill & (26)\end{array}$

From the uniqueness of interpolating polynomials (31, p. 98), the coefficients of the quadratic in Equation (25) are unique, therefore the problem amounts to solving the simultaneous equations:

$\begin{array}{ll}\begin{array}{c}\frac{k_{03}(k_{14}+k_{\mathrm{\text{off}}})-k_{31}{\varphi }_1{x}_{1,\mathrm{\text{E}}}+k_{14}{\varphi }_1{R}_{\mathrm{\text{tot}}}}{{\varphi }_1{x}_{1,\mathrm{\text{E}}}}\\ =\frac{{\overline{k}}_{03}({\overline{k}}_{14}+{\overline{k}}_{\mathrm{\text{off}}})-{\overline{k}}_{31}{\overline{\varphi }}_1{x}_{1,\mathrm{\text{E}}}+{\overline{k}}_{14}{\overline{\varphi }}_1{\overline{R}}_{\mathrm{\text{tot}}}}{{\overline{\varphi }}_1{x}_{1,\mathrm{\text{E}}}}\\ -\frac{k_{03}{k}_{31}(k_{14}+k_{\mathrm{\text{off}}})}{{\varphi }_1{x}_{1,\mathrm{\text{E}}}}=-\frac{{\overline{k}}_{03}{\overline{k}}_{31}({\overline{k}}_{14}+{\overline{k}}_{\mathrm{\text{off}}})}{{\overline{\varphi }}_1{x}_{1,\mathrm{\text{E}}}}.\end{array}& (27)\end{array}$

The solution was found using the SolveAlways function in Mathematica. The only solution to Equations (27), for all values of x_1,E, is given by

$\begin{array}{ll}\begin{array}{c}{\overline{k}}_{31}=k_{31}\\ {\overline{k}}_{14}{\overline{R}}_{\mathrm{\text{tot}}}=k_{14}{R}_{\mathrm{\text{tot}}}\\ \frac{{\overline{k}}_{03}({\overline{k}}_{14}+{\overline{k}}_{\mathrm{\text{off}}})}{{\overline{\varphi }}_1}=\frac{k_{03}(k_{14}+k_{\mathrm{\text{off}}})}{{\varphi }_1}.\end{array}& (28)\end{array}$

Therefore, only k₃₁ and the expressions k₁₄R_tot and k₀₃(k₁₄+k_off)/ϕ₁, containing original parameter combinations, are structurally identifiable with respect to the relationship between FCR_E and x_1,E.

3.2.2. Parameter Estimation

Having analyzed the structural identifiability of the expression for FCR_E vs. x_1,E, it becomes clear that we can rewrite the expression in Equation (14) by combining parameters into new structurally identifiable parameters, as follows:

$\begin{array}{l}{\text{FCR}}_{\text{E}}(x_{1,\text{E}},\psi )=\frac{1}{2x_{1,\text{E}}}(k_{31}{x}_{1,\text{E}}-{\psi }_1-{\psi }_2\\ +\sqrt{k_{31}^2{x}_{1,\text{E}}^2+2k_{31}{x}_{1,\text{E}}\left({\psi }_1-{\psi }_2\right)+{\left({\psi }_1+{\psi }_2\right)}^2}),& (29)\end{array}$

where

$\begin{array}{ll}\begin{array}{c}{\psi }_1=\frac{k_{03}{v}_3(k_{14}+k_{\mathrm{\text{off}}})}{k_{\mathrm{\text{on}}}}\hfill \\ {\psi }_2=k_{14}{R}_{\mathrm{\text{tot}}}\hfill \end{array}\hfill & (30)\end{array}$

are uniquely identifiable parameters. ψ₁ and ψ₂ have units of μmol day⁻¹. The parameter vector to be estimated is now ψ = (k₃₁, ψ₁, ψ₂).

It is assumed that Equation (29) is a close approximation to the relationship between the measured FCR_T and x_1,E. Waldmann and Strober (21) provide FCR_T vs. plasma IgG concentration data. The plasma concentrations of endogenous IgG were multiplied by the average plasma volume v₁, from Table 1, in order to obtain the quantity of endogenous IgG in plasma, x_1,E. The data for FCR_T vs. x_1,E were then fitted using the interior point algorithm implemented within the NonlinearModelFit function in Mathematica. The starting value for the minimization was set to 1 for each parameter. The parameter estimates were constrained to be positive.

Since the data were obtained from 41 individuals, the estimated parameter values are assumed to represent the average parameter values within the population. The parameter estimates and their standard errors are provided in Table 5. The fitted expression given by Equation (29) is plotted alongside the data in Figure 6A. The residuals vs. the fitted values are plotted in Figure 6B. On inspection, the model appears to fit the data well. The residuals appear reasonably homoscedastic and there is no obvious autocorrelation.

TABLE 5

Table 5. Parameter estimates from fitting FCR_E expression to FCR_T vs. x_1,E data.

FIGURE 6

Figure 6. (A) Expression for FCR_E vs. x_1,E, given by Equation (29), fitted to FCR_T vs. x_1,E data from Waldmann and Strober (21). (B) Residuals vs. fitted values.

3.2.3. Validation of Parameter Estimates

There are several issues that may cause the estimates of k₃₁, ψ₁, and ψ₂ to be inaccurate. Firstly, the data were obtained from a sample of 41 individuals, each with their own unique parameter vector; this variability is not accounted for by the estimation procedure. Secondly, the parameters were estimated by fitting the expression for FCR_E vs. x_1,E; however the data are for the FCR_T, which is not equivalent to the FCR_E. In addition, the FCR_T is in practice calculated from measurements of radioactivity in plasma and urine; the form of the measurement errors is therefore not clear.

Due to the aforementioned issues, the validity of the parameter estimates obtained in section 3.2.2 was investigated by estimating the parameters from synthetic data. It is assumed that the parameter values in Table 5 are true population parameter values. Data for FCR_T vs. x_1,E were simulated according to the experimental methodology, described by Waldmann and Strober (21). The data were simulated for 100 sets of 41 subjects. The parameter values were then estimated from the synthetic data, generating 100 estimates for ψ = (k₃₁, ψ₁, ψ₂).

In order to simulate the FCR_T data, parameter values are required for all model parameters (see Equations 1), not just k₃₁, ψ₁, and ψ₂. Population parameter values are therefore required for all model parameters in order to randomly generate unique parameter vectors for individual subjects. The population parameter values for k₂₁, k₁₂, k₁₄, k₀₃, and k_off were fixed to the values from the literature in Table 1. The population value of k₃₁ was fixed to the estimated value in Table 5. The population values of R_tot and k_on/v₃ were calculated by substituting the previously fixed parameter values into Equations (30) and solving. In this way, a population parameter vector was found, for which k₃₁, ψ₁, and ψ₂ are equal to their estimated values. Unique parameter values for 41 individuals were randomly generated from a lognormal distribution, with the median given by the population parameter values. The variance was tuned by trial and error in order to replicate the size of the errors seen in the real data. This process was repeated to produce 100 sets of 41 individual parameter vectors and thus 100 sets of FCR_T vs. x_1,E data. Full details of how the synthetic data were generated are provided in the Mathematica code in the Supplementary Material.

The parameter estimates as a proportion of the true parameter values are plotted in Figure 7, showing the spread of the parameter estimates. It is clear from this plot that the parameter k₃₁ is estimated with higher precision than ψ₁ and ψ₂. The sample mean (μ), sample standard deviation (s.d.), bias (b), and variability (v) of the parameter estimates are given in Table 6. The bias is given by

$\begin{array}{ll}b=\mu -p,& (31)\end{array}$

where p is the true value of the parameter. The variability is given by

$\begin{array}{ll}v=\frac{\sqrt{{\mathrm{\text{s.d.}}}^2+b^2}}{p}.& (32)\end{array}$

Variability as defined by Equation (32) has been used by Chen et al. (32) to evaluate the performance of estimation methods when the assumptions relied upon by the methods, in particular relating to noise, are violated. A larger value of v represents a worse performance of an estimation method. The results suggest that k₃₁ has been estimated with a good level of accuracy (v = 0.0735), but that the parameters ψ₁ and ψ₂ were estimated with a higher level of variability. Based on this result, a future study may look at improving experimental design, for example by increasing the number of subjects, in order to improve upon the variability of the estimates of ψ₁ and ψ₂.

TABLE 6

Table 6. Mean, standard deviation, bias, and variability of the estimates of k₃₁, ψ₁, and ψ₂.

FIGURE 7

Figure 7. Parameter estimates for k₃₁, ψ₁, and ψ₂ divided by the true parameter value.

3.3. Simulation of IgG Responses in Multiple Myeloma

It has been shown that parameter estimates obtained using timecourse data are not robust; however, the parameters k₃₁, ψ₁, and ψ₂ may be obtained with reasonably low variability using FCR data. The results from fitting the timecourse data suggest that the model (Equations 1) may be overparameterized with respect to the available data; we therefore ask whether the plasma IgG response can be sufficiently determined using only the parameters k₂₁, k₁₂, k₃₁, ψ₁, and ψ₂.

Firstly we investigate the plasma IgG response given by the full system model (Equations 1), when the parameters k₃₁, ψ₁, and ψ₂ are equal to the values estimated in section 3.2.2. Random values were generated for certain model parameters and the remaining parameter values calculated so that k₃₁, ψ₁, and ψ₂ are equal to their estimated values. Three parameters (not including both R_tot and k₁₄) out of k₀₃, R_tot, k_off, k₁₄, and k_on/v₃ were fixed to randomly generated values and substituted into Equations (30), yielding a linear system of two equations in two unknowns. Equations (30) were then solved for the remaining two parameters. There are seven sets of three parameters from k₀₃, R_tot, k_off, k₁₄, and k_on/v₃, which can be fixed to give the remaining two parameters. Parameters were generated 10 times, as described, for each of these seven sets, giving 70 parameter vectors in total. The randomly generated parameter values were obtained by assuming a lognormal distribution, in order to ensure positivity, with median set to the parameter value from the literature, given in Table 1, and variance 1. The values generated in this way for the parameters k₀₃, R_tot, k_off, k₁₄, and k_on/v₃ are depicted in Figure 8, showing the extremely wide range of parameter values used. The parameter k₃₁ was set to the estimated value given in Table 5. The values for k₂₁ and k₁₂ were set to the values given in Table 1.

FIGURE 8

Figure 8. Parameter values used to simulate IgG responses in multiple myeloma, plotted on a logarithmic scale. The method used to generate the parameter values is described in section 3.3. The model predictions generated using these parameter values are shown in Figure 9A. Despite the large amount of variation in the parameter values, the model predictions for plasma IgG are extremely similar.

In order to simulate the model under realistic clinical conditions, a model for the IgG synthesis rate in multiple myeloma was used, which has been found to predict responses consistent with real patient data (33). The IgG synthesis rate is described by

$\begin{array}{ll}I(t)=(I_0-I_{\infty })exp(-k_{\mathrm{\text{kill}}}t)+I_{\infty }.& (33)\end{array}$

The following parameter values were used to produce the simulation: I₀ = 76 μmol day⁻¹, I_∞ = 26.5 μmol day⁻¹, and k_kill = 0.055 day⁻¹ (33).

A simulation of the responses in all four model compartments is shown in Figure 9A. Each variable is simulated for 70 unique parameter vectors. The predicted trajectories for plasma IgG and peripheral IgG, respectively, are extremely similar for all 70 parameter vectors; however there is some variation in the responses of IgG in intracellular endosomes, particularly the IgG that is not bound to FcRn receptors. The simulation suggests that, under the investigated conditions, the response in the plasma compartment is relatively insensitive to changes in the individual parameters k₀₃, R_tot, k_off, k₁₄, and k_on/v₃, provided that the parameters k₃₁, ψ₁, and ψ₂ are fixed. The maximal difference between any two trajectories for x₁(t) at any simulated time point is 0.2%.

FIGURE 9

Figure 9. (A) Simulation of responses of plasma IgG [(x₁(t)], peripheral IgG [x₂(t)], IgG in endosomes [x₃(t)], and IgG bound to FcRn in endosomes [x₄(t)]. The scenario shown represents a decreasing tumor burden during therapy. Each variable is simulated for 70 unique parameter vectors. (B) Simulation of responses of plasma IgG [x₁(t)], peripheral IgG [x₂(t)], compared for the four-compartment model and the proposed two compartment model. The responses are indistinguishable by inspection for the two models.

The lack of variation within the predicted responses for plasma and peripheral IgG, when parameters k₂₁, k₁₂, k₃₁, ψ₁, and ψ₂ are fixed, suggests that it may be possible to simulate these two variables using a reduced order model based upon the newly derived expression for the FCR (Equation 29). The equations for this model are given by

$\begin{array}{ll}\begin{array}{c}{ẋ}_1(t)=-(k_{21}+f(x_1(t)))x_1(t)+k_{12}{x}_2(t)+k_{14}{x}_4(t)+I(t)\hfill \\ {ẋ}_2(t)=k_{21}{x}_1(t)-k_{12}{x}_2(t),\hfill \end{array}\hfill & (34)\hfill \end{array}$

where

$\begin{array}{l}f\left(x_1\left(t\right)\right)=\frac{1}{2x_1(t)}(k_{31}{x}_1(t)-{\psi }_1-{\psi }_2\\ +\sqrt{k_{31}^2{x}_1{(t)}^2+2k_{31}{x}_1(t)\left({\psi }_1-{\psi }_2\right)+{\left({\psi }_1+{\psi }_2\right)}^2}).& (35)\hfill \end{array}$

The assumption behind this model is that the fractional rate of IgG catabolism is equal to its fractional rate of catabolism at steady state. A simulation of this model, alongside the original four-compartment model, is shown in Figure 9B. The model is simulated with values of k₂₁ and k₁₂ from Table 1 and all other parameter values from Table 5. The responses for x₁(t) and x₂(t) are very similar for the two models and appear overlayed in Figure 9B. The maximal difference between x₁(t) predicted by the two-compartment model and x₁(t) predicted by the four-compartment model, for any of the 70 parameter vectors used and at any simulated time point, is 0.2%. The responses are indistinguishable by inspection for the two models.

The proposed two-compartment model is based upon the assumption that the fractional rate of IgG catabolism is equal to its fractional rate of catabolism in steady state. When the system is in steady state, this assumption is of course true. However, faster dynamics, caused by a rapid change in the IgG synthesis rate, will cause this assumption to progressively weaken. Further study of the proposed model is required to analyse its relationship with the original four-compartment model and to determine under what conditions the proposed model predictions are within an acceptable region of the four-compartment model predictions.

4. Discussion

The motivation behind the research presented in this paper was to investigate a model suitable for predicting IgG responses in patients with IgG multiple myeloma. When producing predictive simulations of a biomedical system, it is important to know the level of confidence in the model parameter values.

There are numerous published models of FcRn-mediated recycling of IgG in the literature, some of which are cited in the Introduction. Most of these models were developed for IgG-based therapeutic monoclonal antibodies and may not be suitable for characterizing endogenous IgG. Those models characterizing endogenous IgG, for example the models of Li et al. (34) and Chen and Balthasar (10), rely upon a mixture of animal and human data for sourcing parameter values.

For example, parameter values provided by Li et al. (34) for endogenous IgG were taken from the literature, apart from the catabolic clearance (corresponding to k₀₃ of the present study), the vascular reflection coefficient (not included in our model), and the recycling rate constant (corresponding to k₁₄ of the present study). These parameter values were obtained by manually varying the parameters within the model until the results showed a mean half-life of 21 days, a mean IgG synthesis rate of 34 mg kg⁻¹ day⁻¹ and a realistic fold reduction in IgG concentration when FcRn is not present. The values used for the half-life and synthesis rate are those obtained from normal human data by Waldmann and Strober (21) and Waldmann and Terry (23).

One of the problems with the approach taken in previous papers is that, whilst parameter values have been found that provide a half-life of 21 days for an IgG synthesis rate of 34 mg kg⁻¹ day⁻¹, it is not clear what would happen to the half-life when the IgG synthesis rate increases or decreases, under the obtained parameter values. This approach is therefore akin to fitting a model to a curve having only one data point. The nonlinear relationship between synthesis or concentration of IgG and its half-life, which is fundamental to the FcRn-IgG recycling system, may therefore not be captured accurately using this approach.

Another issue with this earlier approach is that it requires the parameter values obtained from the literature to be fixed while the remaining values are varied, therefore implicitly assuming complete confidence in the fixed parameter values that were sourced in the literature. One would question what would happen if one or more of these parameter values were inaccurate by, say, 10% or more, what would be the effect on the corresponding values obtained for k₁₄ and k₀₃?

Having considered the models available in the literature and their issues in respect of parameter identifiability, we identified the need for a semi-mechanistic model with parameter values obtained using only in vivo human data. This approach necessitated a simpler model than those available in the literature and previously discussed. The model studied in this paper is therefore missing some of the mechanisms of the more complex models. However, simplified compartmental models can often be derived from complex physiologically-based models by lumping compartments and processes. Lumped models may be adequate for describing processes of interest, for example responses in a central/plasma compartment. Fronton et al. (14) demonstrate the correspondence between a physiologically-based model and several compartmental model structures for IgG. A similar study could be performed using the models presented in this paper in future work.

In this paper, two observed model outputs were considered: the timecourse of the proportion of a dose of IgG remaining in plasma and in the body of an individual subject; and the FCR vs. the quantity of endogenous IgG in plasma, measured in a cohort of subjects with a range of plasma IgG concentrations. We derived mathematical descriptions of these experimental observations based on the underlying model. Structural identifiability analysis was performed with respect to these observations in order to determine which parameters are structurally uniquely identifiable from the available outputs.

In section 3.1 we estimated parameter values using data for the timecourse of an administered dose of radiolabeled IgG in plasma and in the body. We found that all parameters of the linearized model are structurally globally identifiable. Whilst the model is capable of fitting the data well, the results of 10 runs of differential evolution suggest that the parameter estimates are not robust. Highly different parameter vectors, as illustrated by the relative standard deviations of parameter estimates from 10 runs, produce similarly excellent fits to the data. These results suggest that the available data do not support the complexity of the model. A future study may apply a systematic analysis of model sensitivity and parameter correlations, for example using the profile-likelihood method of Raue et al. (35) or generalized sensitivity functions of Thomaseth and Cobelli (36) [extended to multiple output models by Kappel and Munir (37)]. Another potential study for future work could involve estimating model parameters from synthetic timecourse data, to see whether more frequent sampling or a longer observation period provides more stable parameter estimates. However, as highly different parameter values produce similarly excellent fits to the data, the type of data needed for robust parameter estimation is likely to be of a very high quality. As the data are obtained by taking blood samples, there is a practical limitation on the sampling frequency for an individual subject.

The data used were obtained from tracer experiments that were performed between 1963 and 1990. More recent IgG timecourse data are available; however, these data pertain to therapeutic monoclonal antibodies, which can have different kinetics (38). Timecourse data are also available for patients with IgG multiple myeloma, whose serum IgG concentration is monitored during therapy. However, the production rate of IgG in these patients is determined by the status of the disease. Using these data to estimate model parameters would therefore require simultaneous estimation of IgG production parameters. This would require a more complex structural identifiability analysis and may be considered in future work. For these reasons, more recent data were not used in this study.

The structural identifiability of the relationship between FCR_E and the quantity of endogenous IgG in plasma, x_1,E, was analyzed. We found that the parameter k₃₁ and newly defined parameters ψ₁ = (k₀₃v₃(k₁₄ + k_off))/k_on and ψ₂ = k₁₄R_tot are structurally globally identifiable. These new parameters were estimated using least squares estimation. Estimation with synthetic data shows that these parameters can be estimated with a reasonable level of variability. The parameters k₃₁ and ψ₂ are physiologically meaningful: k₃₁ is the rate at which plasma IgG is internalized into intracellular endosomes and ψ₂ is the maximal rate of recycling of IgG from endosomes into plasma. The 95% confidence interval for k₃₁ (0.135–0.174 day⁻¹) is similar to other values reported in the literature [0.13 day⁻¹ (17); 0.18 day⁻¹ (21); 0.16 day⁻¹ (33)]. The 95% confidence interval for ψ₂ (12.3–39.2 μmol day⁻¹) is smaller than previously reported values [68.6 μmol day⁻¹ (16); 103 μmol day⁻¹ (17)]; however it overlaps with the 95% confidence interval (19.1–60.9 μmol day⁻¹) reported by Kendrick et al. (33).

In applications in which the behavior of the variables x₃(t) and x₄(t), representing unbound and bound IgG in intracellular endosomes, respectively, are of great importance, clearly parameter values are required which determine their behavior, including receptor-ligand binding (k_on/v₃, k_off, and R_tot), recycling of bound IgG into plasma (k₁₄) and degradation of unbound IgG (k₀₃). The results presented in this paper suggest that it is not possible to estimate these parameters from the available data that are only based upon measurements in plasma. In section 3.3, it is shown that these parameters can be varied by several orders of magnitude (see Figure 8) whilst having a minimal effect on the plasma IgG response (see Figure 9A). It is possible that the actions of the parameters determining recycling, degradation, association and dissociation can approximately balance each other out with respect to the dynamics in the plasma compartment, even though the responses of IgG in the endosome are affected by changes in these parameter values. For investigations limited to the behavior of IgG in plasma, model reduction using the parameters k₃₁, ψ₁, and ψ₂ could be investigated in future work. A two-compartment model based upon the newly derived expression for the FCR has been proposed in section 3.3. Further analysis of this model is required to determine whether it is suitable for investigating IgG responses under a range of clinical conditions.

In future work the models studied in this paper could be used to simulate plasma IgG responses in clinical applications, such as the bone marrow cancer multiple myeloma, in which malignant plasma cells secrete large quantities of monoclonal Ig (M-protein). It has been suggested that the FcRn-IgG interaction may play a significant role in the detection of M-protein using a recently-developed mass spectrometry-based method (4). It was found that in patients with IgG-producing disease, the test result was more likely to be positive for M-protein after three months than in patients with IgA-producing disease, whereas after 12 months the patients were equally likely to have a positive test result. Mills et al. (4) have suggested that this effect is due to FcRn-mediated recycling extending the half-life of IgG, emphasizing the importance of assessment times of response. FcRn-mediated recycling also plays a role in the pharmacokinetics of the novel monoclonal IgG agent for multiple myeloma, daratumumab. Yan et al. (5) found that the isotype of the patient's M-protein has an effect on drug exposure, with IgG patients having significantly lower daratumumab concentrations than patients with other M-protein types. Yan et al. (5) proposed that competition between the IgG M-protein and IgG-based daratumumab for FcRn receptors is the reason for this phenomenon. These recent studies show the importance of FcRn-mediated recycling of IgG in multiple myeloma and the need for mathematical modeling and simulation of this system. The model studied in this paper could be used in future work to investigate such problems.

There is a trade-off in modeling between model accuracy, which is more often represented in complex physiologically-based pharmacokinetic models, and accuracy of parameter values, which is more easily achieved with simplified compartmental models. At present, there are very few studies available on parameter estimation for models of IgG-FcRn kinetics using human data due to issues of parameter identifiability. This paper not only provides useful parameter estimates and suggests a novel model structure, but also exposes some of the difficulties in achieving this aim. Researchers pursuing physiologically-based models of IgG in the future may find it useful to compare the rate of IgG internalization into endosomes and the maximal rate of IgG recycling in their model with the values that we have estimated from human data [considering the approach of Li et al. (34) discussed above]. Furthermore, our paper shows the level of analysis (including structural identifiability analysis, estimation from synthetic data, for example) required in order to have confidence in parameter estimates obtained and an understanding of their meaning to the model.

5. Conclusion

It is not possible to estimate all of the model parameters robustly; however certain structurally identifiable parameter combinations have been estimated with a good level of variability. Plasma IgG responses, under typical clinical conditions, are insensitive to large changes in many of the model parameters, provided that certain parameters and parameter combinations are fixed. A reduced-order model, based upon the newly derived expression for the FCR, shows potential for simulating plasma IgG responses under clinical conditions.

Author Contributions

FK performed model analyses. FK, MC, and NE wrote the manuscript. SH initiated the work. MC, NE, OB, and SH supervised the work. SH and OB provided discussion on the clinical application of the work. All authors reviewed and approved the final manuscript.

Funding

This research was supported by a Biotechnology and Biological Sciences Research Council (BBSRC) studentship, through the Midlands Integrative Biosciences Training Partnership (MIBTP), and an Engineering and Physical Sciences Research Council (EPSRC) Impact Acceleration Account (IAA) award.

Conflict of Interest Statement

OB and SH were employed by company The Binding Site Limited.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2019.00674/full#supplementary-material

References

1. Hall A, Yates C (eds.). Immunology. 1st ed. Oxford: Oxford University Press (2010).

Google Scholar

2. Junghans RP, Anderson CL. The protection receptor for IgG catabolism is the β2-microglobulin-containing neonatal intestinal transport receptor. Proc Natl Acad Sci USA. (1996) 93:5512–6.

PubMed Abstract | Google Scholar

3. Kumar S, Paiva B, Anderson KC, Durie B, Landgren O, Moreau P, et al. International Myeloma Working Group consensus criteria for response and minimal residual disease assessment in multiple myeloma. Lancet Oncol. (2016) 17:328–46. doi: 10.1016/S1470-2045(16)30206-6

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Mills JR, Barnidge DR, Dispenzieri A, Murray DL. High sensitivity blood-based M-protein detection in sCR patients with multiple myeloma. Blood Cancer J. (2017) 7:1–5. doi: 10.1038/bcj.2017.75

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Yan X, Clemens PL, Puchalski T, Lonial S, Lokhorst H, Voorhees PM, et al. Influence of disease and patient characteristics on daratumumab exposure and clinical outcomes in relapsed or refractory multiple myeloma. Clin Pharmacokinet. (2017) 57:529–38. doi: 10.1007/s40262-017-0598-1

CrossRef Full Text | Google Scholar

6. Ferl GZ, Wu AM, DiStefano JJ. A predictive model of therapeutic monoclonal antibody dynamics and regulation by the neonatal Fc receptor (FcRn). Ann Biomed Eng. (2005) 33:1640–52. doi: 10.1007/s10439-005-7410-3

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Garg A, Balthasar JP. Physiologically-based pharmacokinetic (PBPK) model to predict IgG tissue kinetics in wild-type and FcRn-knockout mice. J Pharmacokinet Pharmacodyn. (2007) 34:687–709. doi: 10.1007/s10928-007-9065-1

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Fang L, Sun D. Predictive physiologically based pharmacokinetic model for antibody-directed enzyme prodrug therapy. Drug Metab Dispos. (2008) 36:1153–65. doi: 10.1124/dmd.107.019182

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Urva S, Yang V, Balthasar J. Physiologically based pharmacokinetic model for T84. 66: a monoclonal anti-CEA antibody. J Pharm Sci. (2010) 99:1582–600. doi: 10.1002/jps.21918

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Chen Y, Balthasar JP. Evaluation of a catenary PBPK model for predicting the in vivo disposition of mAbs engineered for high-affinity binding to FcRn. AAPS J. (2012) 14:850–9. doi: 10.1208/s12248-012-9395-9

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Deng R, Meng YG, Hoyte K, Lutman J, Lu Y, Iyer S, et al. Subcutaneous bioavailability of therapeutic antibodies as a function of FcRn binding affinity in mice. mAbs. (2012) 4:101–9. doi: 10.4161/mabs.4.1.18543

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Yan X, Chen Y, Krzyzanski W. Methods of solving rapid binding target-mediated drug disposition model for two drugs competing for the same receptor. J Pharmacokinet Pharmacodyn. (2012) 39:543–60. doi: 10.1007/s10928-012-9267-z

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ng CM, Loyet KM, Iyer S, Fielder PJ, Deng R. Modeling approach to investigate the effect of neonatal Fc receptor binding affinity and anti-therapeutic antibody on the pharmacokinetic of humanized monoclonal anti-tumor necrosis factor-α IgG antibody in cynomolgus monkey. Eur J Pharm Sci. (2014) 51:51–8. doi: 10.1016/j.ejps.2013.08.033

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Fronton L, Pilari S, Huisinga W. Monoclonal antibody disposition: a simplified PBPK model and its implications for the derivation and interpretation of classical compartment models. J Pharmacokinet Pharmacodyn. (2014) 41:87–107. doi: 10.1007/s10928-014-9349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xiao JJ. Pharmacokinetic models for FcRn-mediated IgG disposition. J Biomed Biotechnol. (2012) 2012:282989. doi: 10.1155/2012/282989

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kim J, Hayton WL, Robinson JM, Anderson CL. Kinetics of FcRn-mediated recycling of IgG and albumin in human: pathophysiology and therapeutic implications using a simplified mechanism-based model. Clin Immunol. (2007) 122:146–55. doi: 10.1016/j.clim.2006.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hattersley JG. Mathematical Modelling of Immune Condition Dynamics: A Clinical Perspective. Ph.D. thesis, University of Warwick (2009).

18. Solomon A, Waldmann T, Fahey J. Metabolism of normal 6.6 s γ-globulin in normal subjects and in patients with macroglobulinemia and multiple myeloma. J Lab Clin Med. (1963) 62:1–17.

PubMed Abstract | Google Scholar

19. Shah DK, Betts AM. Towards a platform PBPK model to characterize the plasma and tissue disposition of monoclonal antibodies in preclinical species and human. J Pharmacokinet Pharmacodyn. (2012) 39:67–86. doi: 10.1007/s10928-011-9232-2

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Hattersley JG, Chappell MJ, Zehnder D, Higgins RM, Evans ND. Describing the effectiveness of immunosuppression drugs and apheresis in the treatment of transplant patients. Comput Methods Programs Biomed. (2013) 109:126–33. doi: 10.1016/j.cmpb.2011.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Waldmann TA, Strober W. Metabolism of immunoglobulins. Progr Allergy (1969) 13:1–110.

PubMed Abstract | Google Scholar

22. Hansen RJ, Balthasar JP. Pharmacokinetic/pharmacodynamic modeling of the effects of intravenous immunoglobulin on the disposition of antiplatelet antibodies in a rat model of immune thrombocytopenia. J Pharm Sci. (2003) 92:1206–15. doi: 10.1002/jps.10364

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Waldmann TA, Terry WD. Familial hypercatabolic hypoproteinemia. A disorder of endogenous catabolism of albumin and immunoglobulin. J Clin Invest. (1990) 86:2093–8. doi: 10.1172/JCI114947

PubMed Abstract | CrossRef Full Text | Google Scholar

24. OriginLab Corporation. OriginPro 2016. Northampton (2016).

Google Scholar

25. Cobelli C, Foster D, Toffolo G. Tracer Kinetics in Biomedical Research. 1st ed. New York, NY: Springer (2002).

Google Scholar

26. Anderson DH. Compartmental Modeling and Tracer Kinetics. 1st ed. Lecture Notes in Biomathematics. Berlin; Heidelberg: Springer-Verlag (1983).

Google Scholar

27. Wolfram Research Inc. Mathematica Version 11.1. Champaign, IL (2017).

28. Bellman R, Åström KJ. On structural identifiability. Math Biosci. (1970) 7:329–39.

Google Scholar

29. Storn R, Price K. Differential evolution– a simple and efficient heuristic for global optimization over continuous spaces. J Glob Optim.. (1997) 11:341–59.

Google Scholar

30. Ghosh S. A differential evolution based approach for estimating minimal model parameters from IVGTT data. Comput Biol Med. (2014) 46:51–60. doi: 10.1016/j.compbiomed.2013.12.014

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Biswal PC. Numerical Analysis. 1st ed. New Delhi: PHI Learning (2008).

Google Scholar

32. Chen KW, Huang SC, Yu DC. The effects of measurement errors in the plasma radioactivity curve on parameter estimation in positron emission tomography. Phys Med Biol. (1991) 36:1183–200.

PubMed Abstract | Google Scholar

33. Kendrick F, Evans ND, Arnulf B, Avet-Loiseau H, Decaux O, Dejoie T, et al. Analysis of a compartmental model of endogenous immunoglobulin G metabolism with application to multiple myeloma. Front Physiol. (2017) 8:149. doi: 10.3389/fphys.2017.00149

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Li L, Gardner I, Dostalek M, Jamei M. Simulation of monoclonal antibody pharmacokinetics in humans using a minimal physiologically based model. AAPS J. (2014) 16:1097–109. doi: 10.1208/s12248-014-9640-5

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, et al. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics. (2009) 25:1923–9. doi: 10.1093/bioinformatics/btp358

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Thomaseth K, Cobelli C. Generalized sensitivity functions in physiological system identification. Ann Biomed Eng. (1999) 27:607–16.

PubMed Abstract | Google Scholar

37. Kappel F, Munir M. Generalized sensitivity functions for multiple output systems. J Inverse Ill-Posed Problems. (2017) 25:499–519. doi: 10.1515/jiip-2016-0024

CrossRef Full Text | Google Scholar

38. Wang W, Wang EQ, Balthasar JP. Monoclonal antibody pharmacokinetics and pharmacodynamics. Clin Pharmacol Ther. (2008) 84:548–58. doi: 10.1038/clpt.2008.170

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Routh EJ. A Treatise on the Stability of a Given State of Motion: Particularly Steady Motion. 1st ed. London: Macmillan (1877).

Google Scholar