ORIGINAL RESEARCH article
Sec. Molecular Innate Immunity
Volume 10 - 2019 | https://doi.org/10.3389/fimmu.2019.00674
Parameter Identification for a Model of Neonatal Fc Receptor-Mediated Recycling of Endogenous Immunoglobulin G in Humans
- 1School of Engineering, University of Warwick, Coventry, United Kingdom
- 2Department of Research and Development, The Binding Site Group Limited, Birmingham, United Kingdom
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.
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 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
where x1(t), x2(t), x3(t), and x4(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, kij, represent the rate of material flow from compartment j to compartment i, with the convention that 0 represents the environment outside the system. kon and koff are the receptor-ligand binding constants of IgG and FcRn. We denote the volumes of plasma, the peripheral compartment and the endosomes by v1, v2, and v3, respectively. We assume a constant total (bound and unbound) quantity of FcRn, Rtot (6). This means that the quantity of unbound FcRn is represented by [Rtot − x4(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. States and parameters of four-compartment model of IgG metabolism, with parameter values sourced in the literature.
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 (Rtot − x4) represents unbound FcRn receptors and has been eliminated from the model equations. The arrow, labeled k14, from the IgG-FcRn complex compartment (x4) to the unbound FcRn receptor compartment (Rtot − x4), 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) = I0, the system has a stable equilibrium point given by
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. (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 xi(t) = xi,T(t) + xi,E(t), where xi,T(t) and xi,E(t) denote the quantities in μmol in compartment i of radiolabeled and endogenous IgG, respectively, gives
IE (μ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 x1,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 I0 = IE. In summary, the initial conditions are given by
where 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 y1(t)] and in the body [denoted by y2(t)] during the experiment. The observation functions are thus given by
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
where x1,T(t), x2,T(t), x3,T(t), and x4,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 k34 and k43 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).
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−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 IE = 15 µmol day−1 was used; however the linearized model was still valid for D = 1 µmol when comparatively very small values of IE 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. 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 x1,T(t) and x2,T(t), and x3,T(t) and x4,T(t), respectively.
2.6. 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 x3,T(t) and x1,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 FCRT(t) is not constant, as shown by the dependence on time in Equation (8). A simulation of FCRT(t) during the experiment is shown in Figure 4. After around day 5, for the particular parameter values used, FCRT(t) approaches a steady state value, which is denoted here by FCRT,∞:
Solving Equations (6) gives
where Aij and λj (j = 1, …, 4) are expressions in terms of the model parameters and |λ1| > |λ2| > |λ3| > |λ4|. After sufficient time, xi, T(t) can be approximated by Ai4 exp(λ4t); thus, FCRT,∞ is given by
Figure 4. Simulation of FCRT(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, FCRT(t) is equal to the rate of change of radiolabeled IgG in all compartments, divided by the quantity of radiolabeled IgG in plasma:
From Equation (12), it can be seen that FCRT(t) is equal to the slope of the observation y2(t) divided by y1(t), showing how FCRT(t) can be obtained from the observations y1(t) and y2(t). In practice, the experimenter obtains a value for FCRT(tN), where tN is a time toward the end of the experiment, such that FCRT(tN) can be assumed a close approximation of FCRT,∞. Henceforth, the quantity obtained from experiments, FCRT(tN), is referred to simply as FCRT.
It is also possible to derive an expression for the FCR with respect to the endogenous IgG, FCRE. If the endogenous IgG is assumed to remain in steady state, then from the definition of the FCR,
where and are the quantities of IgG in compartments 1 and 3 in steady state, given by Equations (2). Substituting the expression for from Equations (2) into Equation (13), eliminating I0 in favor of using the first equation of Equations (2), and setting , gives the following expression for the FCRE in terms of the quantity of IgG in plasma, x1,E:
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 y1(t) and y2(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 as
where and 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
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 .
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, , 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, , 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.
3.1.2. Parameter Estimation
The parameter vector 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 yi(t, θ0), i = 1, 2, is observed with error at measurement times , where . The observed (with error) values of yi(t, θ0), i = 1, 2, are now denoted by for i = 1, 2 and j = 1, …, Ni. Both outputs y1 and y2 were fitted simultaneously, therefore the cost functional for θ is given by
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 k14, k03, k43, and k34. 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 k34 for subject B, k12 and k14 for subject C, k34 for subject D, and k21 and k14 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 (k34 = 0), there is no transfer from the peripheral compartment to plasma (k12 = 0) or vice versa (k21 = 0), or bound IgG molecules are not recycled into plasma (k14 = 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 y1(t) and y2(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 y2(t) for subject A; in the first 2 days of y1(t) for subject B; for all of y1(t) and the latter part of y2(t) for subject C; between days 2 and 6 for y1(t) and the initial 2 days of y2(t) for subject E; and the first 10 days and final 5 days of y2(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. Timecourse data [y1(t) blue circles; y2(t) red circles] and model fits [y1(t) blue line; y2(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 k12 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. 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 FCRT is measured; however it is difficult to obtain a closed form expression for FCRT. In contrast, we can easily obtain an expression for FCRE in terms of the model parameters and the quantity of endogenous IgG in plasma, x1,E, as given by Equation (14). In this section model parameters are estimated by fitting the expression for FCRE vs. x1,E Equation (14) to the FCRT vs. x1,E data. It is assumed that FCRE is a good approximation to FCRT and the parameter estimates are validated in section 3.2.3 using synthetic data.
3.2.1. Structural Identifiability Analysis
The relationship between FCRE and x1,E is given by Equation (14). Given that the parameters kon and v3 only appear in the model (Equations 3) as the ratio kon/v3, we re-write Equation (14), defining ϕ1 = kon/v3, giving
We wish to know whether the parameter vector 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 such that .
From Equations (13) and (2),
I0 is given in terms of by the solution of the following quadratic equation, obtained by rearranging the first equation of Equations (2) and setting ϕ1 = kon/v3:
Substituting in place of I0 and setting gives the following quadratic equation in FCRE:
Dividing Equation (24) throughout by the coefficient of gives
The expression for FCRE 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,
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 x1,E, is given by
Therefore, only k31 and the expressions k14Rtot and k03(k14+koff)/ϕ1, containing original parameter combinations, are structurally identifiable with respect to the relationship between FCRE and x1,E.
3.2.2. Parameter Estimation
Having analyzed the structural identifiability of the expression for FCRE vs. x1,E, it becomes clear that we can rewrite the expression in Equation (14) by combining parameters into new structurally identifiable parameters, as follows:
are uniquely identifiable parameters. ψ1 and ψ2 have units of μmol day−1. The parameter vector to be estimated is now ψ = (k31, ψ1, ψ2).
It is assumed that Equation (29) is a close approximation to the relationship between the measured FCRT and x1,E. Waldmann and Strober (21) provide FCRT vs. plasma IgG concentration data. The plasma concentrations of endogenous IgG were multiplied by the average plasma volume v1, from Table 1, in order to obtain the quantity of endogenous IgG in plasma, x1,E. The data for FCRT vs. x1,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.
Figure 6. (A) Expression for FCRE vs. x1,E, given by Equation (29), fitted to FCRT vs. x1,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 k31, ψ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 FCRE vs. x1,E; however the data are for the FCRT, which is not equivalent to the FCRE. In addition, the FCRT 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 FCRT vs. x1,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 ψ = (k31, ψ1, ψ2).
In order to simulate the FCRT data, parameter values are required for all model parameters (see Equations 1), not just k31, ψ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 k21, k12, k14, k03, and koff were fixed to the values from the literature in Table 1. The population value of k31 was fixed to the estimated value in Table 5. The population values of Rtot and kon/v3 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 k31, ψ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 FCRT vs. x1,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 k31 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 k31 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.
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 k31, ψ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 k21, k12, k31, ψ1, and ψ2.
Firstly we investigate the plasma IgG response given by the full system model (Equations 1), when the parameters k31, ψ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 k31, ψ1, and ψ2 are equal to their estimated values. Three parameters (not including both Rtot and k14) out of k03, Rtot, koff, k14, and kon/v3 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 k03, Rtot, koff, k14, and kon/v3, 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 k03, Rtot, koff, k14, and kon/v3 are depicted in Figure 8, showing the extremely wide range of parameter values used. The parameter k31 was set to the estimated value given in Table 5. The values for k21 and k12 were set to the values given in Table 1.
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
The following parameter values were used to produce the simulation: I0 = 76 μmol day−1, I∞ = 26.5 μmol day−1, and kkill = 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 k03, Rtot, koff, k14, and kon/v3, provided that the parameters k31, ψ1, and ψ2 are fixed. The maximal difference between any two trajectories for x1(t) at any simulated time point is 0.2%.
Figure 9. (A) Simulation of responses of plasma IgG [(x1(t)], peripheral IgG [x2(t)], IgG in endosomes [x3(t)], and IgG bound to FcRn in endosomes [x4(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 [x1(t)], peripheral IgG [x2(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 k21, k12, k31, ψ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 by
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 k21 and k12 from Table 1 and all other parameter values from Table 5. The responses for x1(t) and x2(t) are very similar for the two models and appear overlayed in Figure 9B. The maximal difference between x1(t) predicted by the two-compartment model and x1(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.
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 k03 of the present study), the vascular reflection coefficient (not included in our model), and the recycling rate constant (corresponding to k14 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 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 k14 and k03?
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 FCRE and the quantity of endogenous IgG in plasma, x1,E, was analyzed. We found that the parameter k31 and newly defined parameters ψ1 = (k03v3(k14 + koff))/kon and ψ2 = k14Rtot 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 k31 and ψ2 are physiologically meaningful: k31 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 k31 (0.135–0.174 day−1) is similar to other values reported in the literature [0.13 day−1 (17); 0.18 day−1 (21); 0.16 day−1 (33)]. The 95% confidence interval for ψ2 (12.3–39.2 μmol day−1) is smaller than previously reported values [68.6 μmol day−1 (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 x3(t) and x4(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 (kon/v3, koff, and Rtot), recycling of bound IgG into plasma (k14) and degradation of unbound IgG (k03). 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 k31, ψ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 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.
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.
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.
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2019.00674/full#supplementary-material
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Keywords: biological systems, lumped-parameter systems, immunoglobulin G, neonatal Fc receptor, parameter estimation, structural identifiability
Citation: Kendrick F, Evans ND, Berlanga O, Harding SJ and Chappell MJ (2019) Parameter Identification for a Model of Neonatal Fc Receptor-Mediated Recycling of Endogenous Immunoglobulin G in Humans. Front. Immunol. 10:674. doi: 10.3389/fimmu.2019.00674
Received: 10 September 2018; Accepted: 12 March 2019;
Published: 08 April 2019.
Edited by:Latha P. Ganesan, The Ohio State University, United States
Reviewed by:Jayajit Das, The Ohio State University, United States
William C. L. Stewart, The Research Institute at Nationwide Children's Hospital, United States
Copyright © 2019 Kendrick, Evans, Berlanga, Harding and Chappell. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Michael J. Chappell, email@example.com