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.

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.

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 (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 IgGproducing 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)(7)(8)(9)(10)(11)(12)(13)(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.

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 bẏ where x 1 (t), x 2 (t), x 3 (t), and x 4 (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 −1 . 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 1 , v 2 , and v 3 , 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 4 (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. When the production rate of IgG is constant, I(t) = I 0 , the system has a stable equilibrium point given bŷ A stability analysis for this equilibrium point is provided in the Supplementary Material.

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.
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).

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, giveṡ I E (µmol day −1 ) 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 0 = I E . In summary, the initial conditions are given by wherex 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 1 (t)] and in the body [denoted by y 2 (t)] during the experiment. The observation functions are thus given by

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 bẏ 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 34 and k 43 are given by 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).

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 −3 and 7 · 10 −3 µ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 −1 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.

Fractional Catabolic Rate
We recall that the FCR (µmol day −1 ) 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 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 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.
Solving Equations (6) gives where A ij and λ j (j = 1, . . . , 4) are expressions in terms of the model parameters and |λ 1 | > |λ 2 | > |λ 3 | > |λ 4 |. After sufficient time, x i,T (t) can be approximated by A i4 exp(λ 4 t); thus, FCR T,∞ is given by The expressions for A 34 and A 14 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.
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: (12) From Equation (12), it can be seen that FCR T (t) is equal to the slope of the observation y 2 (t) divided by y 1 (t), showing how FCR T (t) can be obtained from the observations y 1 (t) and y 2 (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, wherex 1 andx 3 are the quantities of IgG in compartments 1 and 3 in steady state, given by Equations (2). Substituting the expression forx 3 from Equations (2) into Equation (13), eliminating I 0 in favor ofx 1 using the first equation of Equations (2), and settingx 1 = x 1,E , gives the following expression for the FCR E in terms of the quantity of IgG in plasma, x 1,E : 3. RESULTS

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.

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 1 (t) and y 2 (t), given by Equations (4-6). The unknown parameter vector is given by 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 aṡ T and y(t, θ ) = y 1 (t), y 2 (t) 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 (16) 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 x T (0, θ) = (0, 0, 0, 0) 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(θ )) −1 B(θ ), where I is the 4 × 4 identity matrix. G(s) has two elements, corresponding to the two observed outputs, which are given by where the coefficients of s, (θ ) = (φ 1 (θ ), φ 2 (θ), ..., φ 14 (θ )) T , are nonlinear expressions in the parameters. The coefficients of s, (θ ), are given by The coefficients (θ ) are unique with respect to the input-output relationship represented by the transfer function. Introducing an alternative parameter vector, θ = (k 21 , k 31 , k 12 , k 14 , k 03 , k 43 , k 34 ) T , and equating (θ) = (θ ), the resulting set of simultaneous equations is solved for θ using the Solve function in Mathematica (27). The only solution is θ = θ ; therefore all of the parameters in θ are structurally uniquely identifiable.

Parameter Estimation
The parameter vector θ = (k 21 , k 31 , k 12 , k 14 , k 03 , k 43 , k 34 ) 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 θ 0 . For an individual subject it is assumed that y i (t, θ 0 ), i = 1, 2, is observed with error at measurement times j , θ 0 ) for i = 1, 2 and j = 1, . . . , N i . Both outputs y 1 and y 2 were fitted simultaneously, therefore the cost functional for θ is given by where 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 14 , k 03 , k 43 , and k 34 . 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.
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.
In some cases the model parameters are estimated to be zero, or very close to zero, for example k 34 for subject B, k 12 and k 14 for subject C, k 34 for subject D, and k 21 and k 14 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 34 = 0), there is no transfer from the peripheral compartment to plasma (k 12 = 0) or vice versa (k 21 = 0), or bound IgG molecules are not recycled into plasma (k 14 = 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 1 (t) and y 2 (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 2 (t) for subject A; in the first 2 days of y 1 (t) for subject B; for all of y 1 (t) and the latter part of y 2 (t) for subject C; between days 2 and 6 for y 1 (t) and the initial 2 days of y 2 (t) for subject E; and the first 10 days and final 5 days of y 2 (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.
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 12 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.

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.

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 3 only appear in the model (Equations 3) as the ratio k on /v 3 , we re-write Equation (14), defining φ 1 = k on /v 3 , giving We wish to know whether the parameter vector φ = φ 1 , k 31 , k 14 , R tot , k 03 , k off 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 φ = φ 1 , k 31 , k 14 , R tot , k 03 , k off T such that From Equations (13) and (2), I 0 is given in terms ofx 1 by the solution of the following quadratic equation, obtained by rearranging the first equation of Equations (2) and setting φ 1 = k on /v 3 : Substituting FCR Ex1 in place of I 0 and settingx 1 = x 1,E gives the following quadratic equation in FCR E : Dividing Equation (24) throughout by the coefficient of FCR 2 E gives 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 φ such that, Frontiers in Immunology | www.frontiersin.org 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: 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 Therefore, only k 31 and the expressions k 14 R tot and k 03 k 14 + k off /φ 1 , containing original parameter combinations, are structurally identifiable with respect to the relationship between FCR E and x 1,E .

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: where (30) are uniquely identifiable parameters. ψ 1 and ψ 2 have units of µmol day −1 . The parameter vector to be estimated is now ψ = k 31 , ψ 1 , ψ 2 .
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 1 , 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.

Validation of Parameter Estimates
There are several issues that may cause the estimates of k 31 , ψ 1 , and ψ 2 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 31 , ψ 1 , ψ 2 .
In order to simulate the FCR T data, parameter values are required for all model parameters (see Equations 1), not just k 31 , ψ 1 , and ψ 2 . 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 21 , k 12 , k 14 , k 03 , and k off were fixed to the values from the literature in Table 1. The population value of k 31 was fixed to the estimated value in Table 5. The population values of R tot and k on /v 3 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 31 , ψ 1 , and ψ 2 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 31 is estimated with higher precision than ψ 1 and ψ 2 . 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 where p is the true value of the parameter. The variability is given by 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 31 has been estimated with a good level of accuracy (v = 0.0735), but that the parameters ψ 1 and ψ 2 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 ψ 1 and ψ 2 .

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 31 , ψ 1 , and ψ 2 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 21 , k 12 , k 31 , ψ 1 , and ψ 2 . Firstly we investigate the plasma IgG response given by the full system model (Equations 1), when the parameters k 31 , ψ 1 , and ψ 2 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 31 , ψ 1 , and ψ 2 are equal to their estimated values. Three parameters (not including both R tot and k 14 ) out of k 03 , R tot , k off , k 14 , and k on /v 3 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 03 , R tot , k off , k 14 , and k on /v 3 , 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 03 , R tot , k off , k 14 , and k on /v 3 are depicted in Figure 8, showing the extremely wide range of parameter values used. The parameter k 31 was set to the estimated value given in Table 5. The values for k 21 and k 12 were set to the values given in Table 1.
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 The following parameter values were used to produce the simulation: I 0 = 76 µmol day −1 , I ∞ = 26.5 µmol day −1 , and k kill = 0.055 day −1 (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 03 , R tot , k off , k 14 , and k on /v 3 , provided that the parameters k 31 , ψ 1 , and ψ 2 are fixed. The maximal difference between any two trajectories for x 1 (t) at any simulated time point is 0.2%.
The lack of variation within the predicted responses for plasma and peripheral IgG, when parameters k 21 , k 12 , k 31 , ψ 1 , and ψ 2 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 bẏ 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 21 and k 12 from Table 1 and all other parameter values from Table 5. The responses for x 1 (t) and x 2 (t) are very similar for the two models and appear overlayed in Figure 9B. The maximal difference between x 1 (t) predicted by the two-compartment model and x 1 (t) predicted by the fourcompartment 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.

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 IgGbased 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 03 of the present study), the vascular reflection coefficient (not included in our model), and the recycling rate constant (corresponding to k 14 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 −1 day −1 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 −1 day −1 , it is not clear what would happen to the halflife 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 14 and k 03 ?
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 31 and newly defined parameters ψ 1 = (k 03 v 3 k 14 + k off )/k on and ψ 2 = k 14 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 31 and ψ 2 are physiologically meaningful: k 31 is the rate at which plasma IgG is internalized into intracellular endosomes and ψ 2 is the maximal rate of recycling of IgG from endosomes into plasma. The 95% confidence interval for k 31 (16); 103 µmol day −1 (17)]; however it overlaps with the 95% confidence interval (19.1-60.9 µmol day −1 ) reported by Kendrick et al. (33).
In applications in which the behavior of the variables x 3 (t) and x 4 (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 3 , k off , and R tot ), recycling of bound IgG into plasma (k 14 ) and degradation of unbound IgG (k 03 ). 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 31 , ψ 1 , and ψ 2 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 FcRnmediated 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 physiologicallybased 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.

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.