Simulation-based simplification of target-mediated drug disposition model of denosumab
Yu Fu1, Ye Yao1, Peiming Ma2, Xuan Zhou2, Wei Lu1, Tianyan Zhou1*           
1. Beijing Key Laboratory of Molecular Pharmaceutics and New Drug Delivery System, Department of Pharmaceutics, School of Pharmaceutical Sciences, Peking University Health Science Center, Beijing 100191, China
2. GSK Research and Development, Shanghai201203, China
Abstract: Target-mediated drug disposition (TMDD) model is one of the main modeling theories for studying nonlinear pharmacokinetics (PK) of monoclonal antibodies. However, there are too many parameters in full TMDD model to be estimated based on limited clinical data, leading to instability of the final model. In the present study, we analyzed the predictive ability and applicability of a simplified quasi-steady state (QSS) model with the assumption that the total target concentration was a constant parameter during treatment with monoclonal antibody in clinical data modeling. Based on the parameters of a published TMDD model of denosumab, simulations were performed at population and individual levels. Then, a simplified TMDD model, QSS model, was used to examine the effects of hypotheses, in which the total receptor concentration was constant or variable on model fit and stability of parameter estimation. Both simulations at the population level and model fit results of simulated individual data showed that at the therapeutic doses, the total receptor concentration had little influence on changes in drug concentration, and the model with constant total receptor concentration had the same predictive power. The validated hypothesis could be applied to clinical trial design and selection of the optimal PK model in the development of monoclonal antibodies.   
Keywords: Target-mediated drug disposition model; Monoclonal antibody; Nonlinear pharmacokinetics; Denosumab; Simulation
CLC number: R943.5                Document code: A                 Article ID: 10031057(2018)1176710
1. Introduction
Pharmacokinetic (PK) characteristics of some bio-macromolecular drugs (such as monoclonal antibodies) may be affected when they bind with high affinity to their pharmacological targets (such as receptors), mainly manifested as nonlinear elimination of the drugs[1]. Levy et al. have first proposed the concept of target-mediated drug disposition (TMDD) to describe this nonlinear elimination in 1994. TMDD refers to the phenomenon that the binding of drugs with their pharmacological targets significantly affects their disposition[2]. TMDD is one of the main causes of  nonlinear drug elimination. The TMDD theory was first applied by Bauer et al. to the population PK/PD modeling of human antibody hu1124[3]. Subsequently, Mager and Jusko have systematically constructed the framework of TMDD model, including differential equations for the amounts of drug, target and drug-targetcomplex. It quantitatively describes the potential biologicalprocesses of the target, and it is commonly referred to as full TMDD model[4].
TMDD model provides a theoretical basis for quantification of the in vivo physiological processes of monoclonal antibodies, and it can be used duringPK and pharmacodynamic (PD) modeling of various macromolecular drugs with nonlinear PK behavior[24].The application of TMDD model usually requires extensive data information, including drug concentration, target concentration, and/or drug-target complex concentration data. In addition, the frequency and time range of sampling are also crucial for characterizing the dynamics of drug-target binding process[5]. Data information in each phase reflects the characteristics of one or more specific parameters in TMDD model. The data can support the establishment of full TMDD model only when collected with sufficient sampling frequency and appropriate time scale[6].
However, in most cases, it is quite difficult to collect concentration data of target or drug-target complex in clinical practice, which is limited by the methodology of target detection or other ethical reasons. A stable full TMDD model cannot be established when drug concentration is the only abundant data available[4,7,8]. Under this condition, reasonable assumptions can be made during drug-target binding process based on drug concentration data, such as rapid equilibrium or rapid binding hypothesis, and then full TMDD model is simplified to its approximate model, such as QE (quasi equilibrium) or QSS (quasi steady state) and so on. These approximate models accord with the mechanism of TMDD, and meanwhile exhibit robust predictive ability, so that they can be applied to characterize nonlinear PK behavior of macromolecules[911].
In QSS model, it is assumed that the free drug, target and complex are in a quasi-steady state, when the bindingrate of drug and target is comparable to the rate of dissociation and clearance of complex[12]. In this study, denosumab, a human monoclonal immunoglobulin G2 antibody against bone metastasis in solid malignancy, was selected as a model drug to test the stability of parameter estimation of QSS approximate model in clinical practice[13]. Simulation analysis was performed mainly based on population PK model parameter estimates of denosumab in literatures. Then, simulated individual data were fit by simplified QSS models in order to investigate the applicability and feasibility of the simplification methodology concerning whether total receptor concentration could be a constant.
2. Methods
2.1. Model structure
A population PK model has been established by Gibiansky et al. in 2012 based on concentration data of denosumab collected from 14 clinical trials[13]. The PK data was described using a two-compartment model with first-order absorption, linear clearance and QSS approximation of TMDD, in which the total concentration of receptor, Rtot, was expressed as a function of time. The model structure is shown in Figure 1, and the equations of the model are as follows:
where D represents the amount of drug in the depot compartment; Ka is the absorption rate constant; C and Vc represent the free drug (i.e. unbound drug) concentrationand apparent volume of distribution in the central compartment, respectively; Ctot is a sum of unbound and bound drug; T is the free drug concentration in the peripheral compartment; R and P indicate the unbound receptor and drug-receptor complex concentrations, respectively; Rtot is the sum of R and P, and ke is linear clearance rate constant of free drug (mainly through catabolism or renal filtration). Binding and dissociation rate constants of drug-receptor complex are expressed as kon and koff, respectively. kct and ktc are distribution rate constants between the two compartments. ksyn and kdeg represent the intrinsic production and degradation rate constants of receptor, respectively. kmet is the rate of complex’s metabolism. The parameter KSS refers to the steady state of complex. 
Figure 1. Two-compartment TMDD model structure with extravascular administration (s.c.).D represents the free drug amount in the depot compartment. C and T are the free drug concentrations in the central and peripheral compartments. R and P indicate the unbound receptor and drug-receptor complex concentrations. The definitions of the rate parameters are stated in the text. 
2.2. Theoretical analysis of Rtot
When the metabolisms of target and drug-target complex share the same metabolic enzyme or pathway, kmet might sometimes be equal to kdeg[12]. In this case, there is merely linear elimination of Rtot, i.e. in equation (6), so that . Since the function is continuous and derivable everywhere, the rate of change in Rtot gets to zero, which means that Rtot remains the same as its initial value, which isksyn/kdeg, and becomes a constant parameter.
However, the elimination rates of complex andreceptor are not the same in reality, i.e.kdegkmet. Underthis condition, the equation (6) is obviously not equal to 0, and Rtot changes with time. Whether this variation is significant enough to affect drug concentration in the central compartment, C, and whether turnover of Rtot can be predicted only based on drug concentration data are still worth investigating. Here the range of Rtot and its influence on C will be discussed.
Initially, all the receptors in the body are unbound before drug administration under the condition thatthe drug is not exogenous. According to the initial value of R, we got Rtot(0) = R (0) = ksyn/kdeg.
Once the drug is given, it rapidly binds to its target along with increasing C. According to the parameters estimated by Gibiansky, the elimination of unbound receptor (kdeg = 0.00116) is much slower than that of complex (kmet = 0.0112), suggesting that Rtot is mainly removed in the form of complex (P) rapidly afteradministration. Therefore, when C is much higher than KSS in equation (6), the rate of change in Rtot is approximate to (ksynkmetRtot).
Then, Rtot will tend to a new steady-state level, where Rtot,SS = ksyn/kmet, which also appears to be its limit. Since kdeg was much smaller than kmet in this study, receptor concentration was decreased after drug administration. Finally, the maximum and minimum of Rtot were as follows:
Rtot,max = ksyn/kdeg              (8)
Rtot,min = ksyn/kmet              (9) 
2.3. Simulation
Based on the parameter values estimated in literature, we simulated the population prediction of drug concentration in the central compartment, as well as total concentration of receptor. Briefly, 1000-time simulations of individual PK data of denosumab were executed and fitted by Rtot constant model, Rtot variable model and classic two-compartment model, respectively.
In addition, the three models were quantitatively evaluated using 2 log-likelihood (2LL), Akaike information criterion (AIC) and Bayesian information criterion (BIC). The formula for each evaluation method was as follows:
AIC = 2LL + 2k                                                                  (11)
BIC = 2LL + ln(n)*k                                                            (12)
where n is the number of data, and k is the number of parameters in the model. A lower 2LL, AIC or BIC value means a better model.
All the simulations were performed using R software with mlxR package. Population analysis was carried out in NONMEM (Version 7.3.0) with PsN (Version 4.2.0) and a gFORTRAN compiler (Version 4.8.2) using Stochastic Approximation Expectation-Maximization (SAEM) method.
3. Results
3.1. Simulations at population level
To investigate the changing pattern of Rtot and its effect on TMDD model, dynamics of denosumab and its receptor were simulated based on the population typical values of PK parameters reported in literatures. Drug concentration in the central compartment and total concentration of receptor during 200 d aftermultiple subcutaneous doses (120 mg once every 28 d) or a single subcutaneous dose (120 mg) were illustrated in Figure 2. In panel A, C reached a steady state after four or five dosing cycles (112–140 d), which was characterized by similar peak and trough concentrations. Moreover, no significant nonlinear behavior was observed. Rtot in panel B was rapidly and dramatically decreased after drug administration, reached its steady state in about 10 d and remained below 0.1 mg/L ever after. However, following a single dose, clearance of denosumab was mainly linear during the first 100 d, and showed apparent nonlinearity after 150 d (panel C). As shown in panel D, Rtot was decreased rapidly in about 10 d and remained stable at a low level until it began to rebound quickly after 150 d. Based on the above simulations, limits of Rtot were mathematically derived from model equations and then incorporated into model simplification process.
The maximum and minimum Rtot following both single and multiple doses of 120 mg denosumab accorded with Rtot,max and Rtot,min in equation (8–9), as shown in Figures 2E and 2F. Following an initial decline, Rtot mostly stayed near Rtot,min for 150 d after single administration (Fig. 2E) and during the whole time course with a dosing interval of 4 weeks (Fig. 2F). Simulations of the QSS model for denosumab above indicated that the total concentration of receptor was not dramatically fluctuated after a single dose, and it was always close to its lower bound derived from equation (9) when denosumab was given Q4W. Instead of estimating turnover parameters of receptor, the Rtot constant assumption appeared to be sufficient to fit clinical data with limited information.
Simulations based on a differential equation for Rtot suggested that Rtot did not change much. Therefore, we assumed that Rtot was a constant parameter which did not change over time after drug administration. PK profiles after a single dose were predicted by setting the value of Rtot in the original TMDD model according to the maximum and minimum total concentrations of receptor derived in 2.2. As shown in Figures 3A and 3B, when Rtot was a variable (indicated as red solid line) or equal to Rtot,min (green solid line), C versus time curves overlaid for the first 150 d, which was longer than the actual period with PK observations (18 weeks). The deviation of the blue line from the other two suggested that PK behavior of denosumab was significantly different as Rtot reached its maximum. In this case, receptor-mediated nonlinear elimination of the antibodywas so fast that denosumab was soon cleared. In addition,most of the concentrations within 100 d were higherthan 10*KSS, and with C/(KSS + C) over 0.9, the premise that C is much higher than KSS was met.
Figure 2. Simulations in the central compartment and total concentration of receptor on a semi-logarithmic scale. (A) Drug concentration following multiple doses of 120 mg Q4W. (B) Total concentration of receptor following multiple doses. (C) Drug concentration following a single dose of 120 mg. (D) Total concentration of receptor following a single dose. (E) Simulated total receptor concentration versus time following a single dose (120 mg) on a semi-logarithmic scale. (F) Simulated total receptor concentration versus time following multiple doses (120 mg Q4W) on a semi-logarithmic scale. The red dashed lines indicate Rtot,min = ksyn/kmet, and the green indicate Rtot,max = ksyn/kdeg.
Similar trends were observed after multiple doses. Identical PK profiles were predicted for Rtot variable and Rtot minimum hypothesis, while there was slight difference when Rtot was equal to Rtot,max. Replacing the differential equation of Rtot with its minimum value had no impact on PK simulation of denosumab for the current dosage and model parameters. The precondition that C is much higher than KSS met for both regimens.
On the other hand, when kmet is smaller than kdeg, unbound receptor will decrease after administration, while Rtot will gradually increase along with the accumulationof complex. In this case, the maximum and minimum values of Rtot are:
Rtot,max = ksyn/kmet                       (13)
Rtot,min = ksyn/kdeg                  (14)
PK profiles of denosumab were simulated with different Rtot after exchanging the values of kmet and kdeg, i.e., kmet = 0.00116 day1 and kdeg=0.0112 day1. As shown in Figures 3C and 3D, the three lines were similar following a single dose and even identical after multiple doses, indicating the feasibility of Rtot constant hypothesis. 
Figure 3. Simulated denosumab concentration in the central compartment versus time on a semi- logarithmic scale when kmet = 0.0112 day1 and kdeg = 0.00116 day1, with different Rtot assumptions after a single dose (A) and multiple doses of 120 mg Q4W (B). Simulated denosumab concentration in the central compartment versus time on a semi- logarithmic scale when kmet = 0.00116 day1 and kdeg = 0.0112 day1, with different Rtot assumptions after a single dose (C) and multiple doses of 120 mg Q4W (D). The red solid lines show PK profiles with Rtot changing with time. For the green and blue solid lines, Rtot is a constant with a value of ksyn/kmet and ksyn/kdeg respectively. The red dashed lines indicate 10*KSS 
3.2. Simulation at individual level and model fitting
In order to further investigate the influence of inter-individual variability (IIV) and error on model selection, PK data of denosumab in 1000 individuals were simulated based on the Rtot variable equation, IIV of the parameters and the error model in literatures using mlxR software package. Figure 4 illustrates PK profiles following multiple doses of 120 mg Q4W during 24 weeks. 
Figure 4. Simulated PK data of 1000 individuals following multiple doses of 120 mg Q4W on a semi-logarithmic scale, based on the TMDD model with variable Rtot in literature. The solid line represents the median denosumab concentrations, and the dashed lines are the 2.5th and 97.5th percentiles. The area between the dashed lines means the 95% confidence interval of the simulations. 
Subsequently, the simulated individual data were fitted by Rtot constant model, Rtot variable model and classic two-compartment model separately using NONMEM software. The goodness-of fit plots in Figures S1–S3 suggested that individual predictions of the three models matched the observations well. Table 1 lists the model parameter estimates. Rtot variable model had two more parameters kdeg and IIV in kdeg, compared with Rtotconstant model. Besides, the relative standard error (RSE) of IIV in kdeg in Rtot variable model and that of IIV in Q in two-compartment model were extremely large. 
Table 1. Estimated parameters of the three models.
Table 2 shows the evaluation results of the three models according to –2LL, AIC, and BIC criteria. Rtot constant model with the lowest scores was considered the best to fit the simulated individual data.
Table 2. –2LL, AIC and BIC results of the three models.
4. Discussion
QSS model is one of the most commonly used TMDD approximate models, with predictive power almost comparable to that of full TMDD model. In the traditional QSS model, the total concentration of receptor, Rtot, changes with time after dosing, and it is described by a differential equation[14]. However, in this study, we observed that Rtot remained stable during most of the time, based on the simulations with PK parameters of denosumab in literatures[13]. Therefore, we hypothesized that Rtot could be set as a constant parameter, which did not change in the course of administration. The maximum and minimum Rtot were derived from the model and set separately as the value of Rtot to simulate and compare drug concentration-time curves. The results showed that the model with a fixed Rtot within this range could possibly fit the data, in which Rtot actually changed. The drug concentration, C, was higher than 10*KSS for most of the initial period after administration, which was consistent with the premise assumption that C was much larger than KSS.
When fitting the individual data simulated from the original TMDD model with a time-varying Rtot,goodness of fit in DV versus IPRED plots demonstratedthat the three models all had good descriptive and predictive abilities. The RSEs of parameters estimated in Rtot constant model were all within 30%, indicating the precision of estimation. However, the RSE of IIV in kdeg in Rtot variable model and that of IIV in Q intwo-compartment model were extremely large, indicating that the two models failed to accurately estimate theseparameters. Drug concentration data alone are insufficient to estimate parameters related to target turnover in QSS model, such as variability of kdeg. If prediction of the kinetics of target through modeling is necessary, preclinical studies should be turned to. Models based on preclinical data can help to confirm physiological mechanisms, set clinical PK parameters, and determine elimination pathways of drug and target[12].
Model evaluations according to –2LL, AIC and BIC standards all showed that Rtot constant model was the best to fit this set of data, followed by two-compartment model and Rtot variable model. The three standards examine not only the precision of model prediction, but also the complexity of the model. The more complicated the model structure is, and the more parameters it contains, the higher the score goes[15]. Therefore, though all the three models succeeded in fitting the nonlinear PK data, Rtot constant model with the fewest parameters turned out to be the most stable one during estimation process.
A priori identifiability is a necessary feature of a good mathematical model system and a prerequisite for parameter estimation[16]. The exact arithmetic rank (EAR) method has been used to test the a priori identifiability of the parameters of full TMDD model and its approximate models (QSS, QE and Michaels-Menten equation). The results show that these models all have a prioriidentifiability in structure[17]. On the other hand, a posterioriidentifiability, also called actual identifiability, means that the model parameters can be estimated through specific experimental design and data analysis[18]. In this simulation-based study, parameters of receptorturnover couldn’t be identified by the data simulated in our preset conditions. Meanwhile, PK profiles of denosumab proved to be not quite sensitive to the parameters of receptor turnover. According to the sensitivity analysis of TMDD model predictions towards model parameter fluctuations conducted by Abraham et al., drug concentration is highly sensitive to perturbation of the linear clearance rate parameter and the value of the dissociation constant[19].
However, there are differences in biological parameters, as well as mode and degree of binding to pharmacological targets between different drugs. The differences in parameters of different drugs may be up to hundreds of times[20]. Our study was based only on one model drug and we did not yet explore PK profiles of other monoclonal antibodies or other drugs demonstrating TMDD. Additional model drugs should be included in further study.
5. Conclusions
The applicability and feasibility of a simplification method with constant total target concentration, Rtot, in establishment of TMDD model were explored through simulation. Rtot was set to its maximum and minimum values, respectively, and the predicted population PK profiles of denosumab under both conditions were similar to the time course of drug concentration, when Rtot changed over time. Rtot constant model proved to be the best when fitting simulated individual PKdata compared with Rtot variable model and classic two-compartment model. Rtot constant model had the best predictive power and practicability for the given TMDD model of denosumab. Collectively, it had the potential to describe the PK profiles of other monoclonal antibodies.
[1] Dirks, N.L.; Meibohm, B. Population pharmacokinetics of therapeutic monoclonal antibodies. Clin. Pharmacokinet. 2010, 49, 633–659.
[2] Levy, G. Pharmacologic target-mediated drug disposition. Clin. Pharmacol. Ther.1994, 56, 248–252.
[3] Bauer, R.J; Dedrick, R.L; White, M.L; Murray, M.J. Garovoy MR. Population pharmacokinetics and pharmacodynamics of the anti-CD11a antibody hu1124 in humansubjects with psoriasis. J. Pharmacokinet. Biopharm.1999, 27, 397–420.
[4] Mager, D.E; Jusko, W.J. General pharmacokinetic model for drugs exhibiting target-mediated drug disposition. J. Pharmacokinet. Pharmacodyn. 2001, 28, 507–532.
[5] Aston, P.J; Derks, G.; Raji, A.; Agoram, B.M.; van der Graaf, P.H. Mathematical analysis of the pharmacokinetic-pharmacodynamic (PKPD) behaviour of monoclonal antibodies: predicting in vivo potency. J. Theor. Biol. 2011, 281, 113–121.
[6] Peletier, L.A. Gabrielsson J. Dynamics of target-mediated drug disposition: characteristic profiles and parameter identification. J. Pharmacokinet. Pharmacodyn. 2012, 39, 429–451.
[7] Mager, D.E.; Jusko, W.J. Receptor-mediated pharmacokinetic/pharmacodynamic model of interferon-beta 1a in humans. Pharm. Res. 2002, 19, 1537–1543.
[8] Meijer, R.T.; Koopmans, R.P.; ten Berge, I.J.; Schellekens,P.T. Pharmacokinetics of murine anti-human CD3 antibodies in man are determined by the disappearance of target antigen. J. Pharm. Exp. Ther. 2002, 300, 346–353.
[9] Mager, D.E.; Krzyzanski, W. Quasi-equilibrium pharmacokineticmodel for drugs exhibiting target-mediated drug disposition. Pharm. Res. 2005, 22, 1589–1596.
[10] Woo, S.; Krzyzanski, W.; Jusko, W.J. Target-mediated pharmacokinetic and pharmacodynamic model of recombinant human erythropoietin (rHuEPO). J. Pharmacokinet. Pharmacodyn. 2007, 34, 849–868.
[11] Ma, P. Theoretical considerations of target-mediated drug disposition models: simplifications and approximations.Pharm. Res. 2012, 29, 866–882.
[12] Gibiansky, L.; Gibiansky, E. Target-mediated drug disposition model: approximations, identifiability of model parameters and applications to the population pharmacokinetic-pharmacodynamic modeling of biologics. Expert. Opin. Drug Metabolism Toxicol.2009, 5, 803–812. 
[13] Gibiansky, L.; Sutjandra, L.; Doshi, S.; Zheng, J; Sohn W; Peterson MC; et al. Population pharmacokinetic analysisof denosumab in patients with bone metastases from solid tumours. Clin. Pharmacokinet. 2012, 51, 247–260.
[14] Gibiansky, L.; Gibiansky, E.; Kakkar, T.; Ma, P. Approximations of the target-mediated drug dispositionmodel and identifiability of model parameters. J. Pharmacokinet. Pharmacodyn.2008, 35, 573–591.
[15] Pan, W. Akaike’s information criterion in generalized estimating equations. Biometrics.2001, 57, 120–125.
[16] Jacquez, J.A; Greif, P. Numerical parameter identifiability and estimability: integrating identifiability, estimability, and optimal sampling design. Math. Biosci.1985, 77, 201–227.
[17] Eudy, R.J; Riggs, M.M.; Gastonguay, M.R. A Priori Identifiability of Target-Mediated Drug Disposition Models and Approximations. AAPS J.2015, 17, 1280–1284.
[18] Lamberton, T.O; Condon, N.D.; Stow, J.L.; Hamilton, N.A. On linear models and parameter identifiability in experimental biological systems. J. Theor. Biol.2014, 358, 102–121.
[19] Abraham, A.K.; Krzyzanski, W.; Mager, D.E. Partial derivative-based sensitivity analysis of models describingtarget-mediated drug disposition. AAPS J. 2007, 9, E181–189.
[20] Gabrielsson, J.; Peletier, L.A,; Hjorth, S. In vivo potency revisited-Keep the target in sight. Pharmacol. Ther.2018, 184, 177–188.
付毓1, 姚烨1, 马培敏2, 周绚2, 卢炜1, 周田彦1*
1. 北京大学医学部 药学院 分子药剂学与新释药系统北京市重点实验室, 北京 100191
2. 葛兰素史克研发中心, 上海 201203
摘要: 靶点介导的药物处置模型(Target-mediated drug disposition, TMDD)是研究单克隆抗体非线性PK的主要模型理论之一。然而全量TMDD模型参数较多,有限的临床数据并不能支持模型估算全部参数。因此,本研究旨在通过公式推导和模型仿真探讨TMDD模型的简化方法,提升TMDD模型在有限临床数据建模中的适用性。基于一项关于地诺单(denosumab)TMDD模型研究的结果及其模型参数,在群体水平上和个体水平上进行仿真分析。然后,利用准稳态近似模型,在模型拟合和参数估计的稳定性上,检验受体总浓度的两种假设的影响。结果表明,在治疗剂量下, 总受体浓度对药物浓度变化的影响很小,而假设恒定总受体浓度的模型具有相同的预测能力。该简化方法可应用于单抗类药物研发中合理实验设计的制定和最优PK模型的选择。 
关键词: 靶点介导的药物处置模型; 单克隆抗体; 非线性药物动力学; 地诺单抗;仿真 
Received: 2018-05-15,Revised: 2018-05-29, Accepted: 2018-09-30.
*Corresponding author. Tel.: +86-010-82801717, E-mail:

本作品采用知识共享署名-非商业性使用 4.0 国际许可协议进行许可。