Automated Scale Reduction of Nonlinear QSP Models With an Illustrative Application to a Bone Biology System.

Integrating quantitative systems pharmacology (QSP) into pharmacokinetics/pharmacodynamics (PKPD) has resulted in models that are highly complex and often not amenable to further exploration via estimation or design. Because QSP models are usually depicted using nonlinear differential equations it is not straightforward to apply some model reduction techniques, such as proper lumping. In this study, we explore the combined use of linearization and proper lumping as a general method to simplification of a nonlinear QSP model. We illustrate this with a bone biology model and the reduced model was then applied to describe bone mineral density (BMD) changes due to denosumab dosing. The methodologies used in this study can be applied to other multiscale models for developing a mechanism-based structural model for future analyses.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Proper lumping has been used for order reduction of complicated models. This technique is, however, not straightforward to apply for nonlinear differential equations that are not uncommon in QSP models.

WHAT QUESTION DID THIS STUDY ADDRESS?

Can nonlinear QSP models be simplified by a two‐stage process (1) inductive linearization of the system and (2) application of proper lumping? Application of this process is illustrated with a case example.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

Nonlinear QSP models can be simplified using a two‐stage process (inductive linearization and proper lumping). Both components can be automated, although, in this work, only proper lumping is fully automated. The reduced model was then able to be applied to extrapolate long‐term responses that were not able to be captured by an empirical approach.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

The scale reduction of a QSP model provides a simple mechanistic model, which can then be used to predict long‐term responses from short‐term or middle‐term data.

Quantitative systems pharmacology (QSP) models are increasingly used in drug development to provide a deeper understanding of the mechanisms of action of drugs and their likely effects on the system, or to identify appropriate disease targets in preclinical settings.1, 2 Irrespective of the purpose of development, such models are generally not suitable for estimation purposes due to the large number of states and parameters yielding long runtimes and numerical instability.3, 4

Based on identifying a specific input‐output relationship, however, the system may be reduced to fewer states and parameters that may then be suitable for estimation purposes. Proper lumping is one such technique that has been used for model order reduction with some applications to systems models.5, 6, 7 With this method, the original states of the model are lumped into a reduced number of pseudo‐states, resulting in a model of lower dimensionality but with similar input‐output behavior as well as physiological interpretation. For linear systems, parameter values as well as initial conditions of the lumped states can be obtained directly using the lumping formulae. This is, however, not the case for nonlinear ordinary differential equations (ODEs) that are not uncommon in QSP models.8, 9 In such cases, only initial conditions are provided from lumping formulae,7, 10 therefore, lumping nonlinear systems is not straightforward. We note that solving the (nonlinear) system itself is not usually difficult by using numerical time‐stepping solvers (e.g., Runge–Kutta).

The primary aim of this study was to explore simplification of a nonlinear QSP model. We chose an existing 28‐state bone biology model just as an illustrative example of a nonlinear QSP model. This model links calcium homeostasis and bone biology,11 and was extended for the prediction of lumbar spine bone mineral density (BMD).12 It should be noted that this work is not intended to provide a complete treatment of bone biology or models that have been developed to describe drug effects on bone biology (readers are referred to ref. 13, 14 for more extensive discussions on this therapeutic area). The second aim of this study was, therefore, to investigate if the reduced model provides adequate prediction of responses. Here, we consider the extrapolation with respect to “time” as one of the ways to assess the utility of the model: (i) calibrate the reduced model (estimate model parameters) using short‐term to middle‐term data, then (ii) assess the utility of the model using long‐term data (i.e., extrapolation), in which denosumab (a monoclonal antibody osteoporosis treatment targeting receptor activator of nuclear factor kappa‐B ligand (RANKL)) is used as a test drug. The predictive performance will be compared to two empirical models developed from data over the same period.

Methods

Model order reduction consisted of three steps:

Linearization of the original nonlinear model.

Scale reduction of the linearized original model.

Finalization of the reduced model for parameter estimation.

The model was then applied to data describing the effects of denosumab.

Linearization of the original nonlinear model

An inductive approximation was used to obtain the linearized system from a nonlinear model.15 The method consists of two steps: (i) linearization to create a linear version of the nonlinear ODE; and (ii) then solving the linear ODE. Here, we briefly cover the main concepts of the linearization. QSP as well as pharmacokinetic/pharmacodynamic (PK/PD) models are often defined as a set of nonlinear ODEs of the general form:

where t is time, (t) = (y 1(t), …, y (t))T is an M × 1 vector of response variables (i.e., M states in the system), (t,) is an M × 1 vector, and (t,) is an M x M nonsingular matrix of which each element includes the corresponding parameters (e.g., rate constants). The initial conditions for the ODE are given by 0 (M × 1). The ODE is linearized by replacing the coefficients that are dependent on the current state with a previous iteration of that state [n−1] at time t. The linearization is initialized (when n = 1) with [0](t) (our initial guess for the first iteration in linearization at time t). The initial guess could be set to our initial conditions for the ODE, 0, for all t (an uninformative starting point). The dimensions of [0] depend on the number of ODEs in the nonlinear system (see ref. 15 for further details). The linearization is performed for the solution [n](t) inductively by:

Here, [n−1](t) is a response vector obtained from the previous iteration and, hence, considered to be a known quantity. Because and now depend on a known quantity, [n−1](t), rather than on the current value of [n](t), Eq. (2) now represents a linear (time‐varying) system that can be expressed in the following form:

With a sufficient n, the inductive approximation (Eq. (2)) is considered to form the true solution. We note that the inductive approximation was applied in order to enable the direct use of proper lumping for nonlinear QSP models. Other possible benefits of using the inductive approximation can be found elsewhere.15

Because the PK of denosumab and the QSP bone biology model were expressed with nonlinear ODEs, they were linearized sequentially. A two‐compartment PK model with a quasi‐steady‐state approximation of the target‐mediated drug disposition model was used to describe denosumab PK.16 Model equations are shown in Appendix . The output of interest in this study is the change in BMD from baseline. Two components in the bone biology model, osteoblasts (OBs) and osteoclasts (OCs), directly characterize the change in BMD in the form12:

Therefore, we focused on those two states when inductively approximating the original nonlinear system. The iteration was stopped when the following maximal relative error fell below 10‐3:

where is a vector of response variables of OB and OC, such that = (OB [1], …, OB [; OC [1], …, OC [)T of dimensions n × 2. A matrix exponential solution was used to solve linearized ODEs using appropriately small time steps. All calculations and simulations were performed using MATLAB R2015b or higher (MathWorks, Natick, MA).

Scale reduction of the linearized original model

The technique of proper lumping was used for scale reduction of the linearized original model. The detail of the proper lumping has been described elsewhere.5 Here, we briefly cover the main concepts. Using Eq. (3), the system for the lumped states can be written in the form:

where is the vector of response from lumped states, and are the vector and matrix of parameters in the reduced system, respectively. Then it is known that and can be derived as:

in which is the lumping matrix transforming the original system to reduced models in the form:

+ in Eq. (8) is the Moore‐Penrose inverse, sometimes referred to as pseudo inverse of . Please note that and can take different values at each time as and as in the original system. It is obvious that and cannot be obtained for nonlinear systems. As in this case, and depend on the original response , which cannot be obtained (solved) from the ODEs of lumped states (Eq. (6)). This is why the linearization is required.

A final lumped model was selected using the following composite criterion17:

T 1 and T 2 represent model performance and complexity, respectively, m is the number of states in the reduced model, and the two indices (T 1 and T 2) were weighted with a user‐defined mixing constant . For any given value of weighting the smallest criterion value will provide the best tradeoff between complexity and performance. The search for the lumped model was started with , and this value was finally determined by qualification of the simplified model based on a visual predictive check. The detail of this process is shown in Appendix .

Considering the important role of responding OB,14, 18 it was forced to remain unlumped in addition to the components for the output of interest (i.e., OC and OB). In this model, OB has been separated into two states where the removal rate is either fast (differentiation into osteocytes) or slow (apoptosis), yielding the minimum number of states to be five (corresponding to responding OB (ROB), fast OB (FOB), slow OB (SOB), OC, and a lumped state from the rest), excluding the denosumab PK components. We note that OB is the sum of FOB and SOB. All calculations and simulations were performed using MATLAB R2015b or higher (MathWorks).

Finalization of the reduced model for parameter estimation

The reduced model includes time‐varying composite parameters (derived from Eqs. (7) and (8)) which may vary nonmonotonically over time and, therefore, be difficult to express with typical functions (e.g., exponential, power, etc.). This is not an issue if time‐varying values of the composite parameters are provided as input data (e.g., $INPUT in NONMEM), although one may want to avoid this, especially for parameter estimation. To circumvent this issue, the time‐varying parameters were re‐converted (unlinearized) back to the original nonlinear functions for states which remained unlumped (replace y [n−1] with y [n] then remove [n] from the equation). The conversion to linear time‐invariant parameters using baseline values (at t = t 0) was also applied for states that were lumped together with other states. Through this process all the parameters in the final reduced model are provided as time‐invariant.

Identifiability analysis

Identifiability analysis was performed using an information approach that evaluated structural and deterministic identifiability.19, 20 In addition, an informal heuristic approach was used to evaluate whether further parameters could be estimated by using a sensitivity analysis. This was performed through parameter estimation using published BMD data until 1 year from the initiation of denosumab dosing,21 in which each parameter was assessed for its influence on the objective function value of NONMEM univariately.

Application of the reduced model

The reduced model and two empirical models were trained on the first 12 months of BMD data extracted from a publication describing the treatment effects of denosumab over 48 months in postmenopausal women.21 The extracted data is available in Table . The BMD data, 12 months after the initiation of drug dosing, are often evaluated as primary analyses in typical phase II studies for osteoporosis treatment.22, 23, 24 The three models were then used to predict the BMD response over the next 36 months. The BMD data were digitized using the software application WebPlotDigitizer version 3.12.

Two empirical models were considered, a direct response (with time‐varying maximum effect (Emax)) and a turnover model, and also used for fitting to data until 12 months:

Empirical model 1

Empirical model 2

where C is the serum concentration of denosumab. The Emax is the maximum response, C50 is the serum concentration of denosumab causing 50% of maximum response, k is the rate constant for the onset of denosumab drug effect, R in is the zero‐order rate constant, k out is the first‐order rate constant for BMD, and Imax is the maximum inhibitory effect of denosumab. The number of parameters to be estimated is three in both models because in Eq. (12) R in is derived as k out times 100 (normalized baseline BMD). After parameter estimation, the three models (the reduced and the two empirical models) were then used to extrapolate the long‐term BMD responses over the following 36 months. Parameter estimation was performed using NONMEM version 7.3.0 (ICON Development Solutions).25 Figures were prepared using SAS version 9.3 (SAS Institute Inc, Cary, NC).

As a positive reference, an existing semimechanistic model of BMD26 was also used for fitting to data until 12 months and for extrapolating the long‐term BMD responses over the following 36 months. Model equations are shown in Table .

Results

Linearization was performed for both the PK model for denosumab and the multiscale bone systems model. The results of the linearization of denosumab when dosing 60 mg (approved dose) every 6 months are shown in Figure 1 a. It is seen that the inductive approximations approached the solution provided by the original nonlinear model, and the iteration stopped at n = 47 with a suitably accurate prediction (E S < 0.001 in Figure 1 b). By fixing this result, the profiles of both OC and OB were then inductively linearized. Successive inductive approximations are shown in Figure 1 c for OC and Figure 1 e for OB. Early iterations, however, showed a nonmonotonic decline in E S, which occurs due to the recursive feedback mechanisms in the system (shown in Figure 1 d for OC and Figure 1 f for OB). An average between successive (odd + even) iterations declined monotonically following an exponential decay for both OC and OB (broken lines in Figure 1 d and Figure 1 f, respectively). The iteration finally stopped at n = 20 (E S < 0.001). The runtime for the combined linearized drug and bone biology model was less than 2 minutes.

Figure 1

Inductive linearization (a for denosumab pharmacokinetics (PK), c for osteoclasts (OCs), and e for osteoblasts (OBs)) and its convergence process (b for denosumab PK, d for OCs, and f for OBs). n, number of iterations in inductive approximations; , the error of successive approximations of the inductive solution. The broken lines in d and f show the average between successive (odd + even) iterations.

Two components in Eq. (2) after the inductive linearization (i.e., (t, [) and (t, [)) where n = 20, were directly used for the calculation of parameter values in each reduced model using Eqs. (7) and (8). Figure 2 shows the lumping criterion value (CC(m, α) in Eq. (10)) for each number of states tested (m). When α = 0.5 in CC(m, α), a 7‐state model was selected (shown in Figure 2) and the model predictions are shown in Figure 3 a. This model roughly agreed with the original profile, but clearly underpredicted until 1 year and overpredicted from years 1–2 but was consistent from years 3–4. With α = 0.7 (more contributions of T 1 in Eq. (10)), the 8‐state model was selected, as shown in Figure 2, and the plots of the model predictions are shown in Figure 3 b. This model provided similar predictions of BMD to the original profile, and was considered superior to the seven‐state model over the whole response‐time profile. This eight‐state model was accordingly chosen as a final reduced model. The schematic representation of the reduced model is shown in Figure 4. In this model, four states (ROB, FOB, SOB, and OC) were forced to remain unlumped (see the Methods section), and additional two states (i.e., active transforming growth factor‐beta (TGF‐β) and receptor activator of nuclear factor kappa‐B (RANK)/RANKL complex), also remained unlumped as a result of the automatic search. Other states in the original model were grouped together to the states of either RANK or RANKL. Equations for the eight‐state model were obtained using lumping formulas (Eqs. (7) and (8)) in the form:

Figure 2

Lumping criterion values for each number of states when a weighting factor α is either 0.5 or 0.7.

Figure 3

Visual predictive check for a reduced model of size m = 7 (a), and size m = 8 (b). Each plot shows bone mineral density (BMD) predictions from the reduced (solid line) and original (dashed thin line) models, respectively, and a 95% credible interval (CI; dashed thick lines) versus time after the initiation of 60 mg of denosumab every 6 months. m, number of states in reduced models.

Figure 4

Schematic representation of reduced eight‐state model. Full arrows indicate flows, and broken arrows indicate control mechanisms (+, stimulating; ‐, inhibition). Three original states included in the receptor activator of nuclear factor kappa‐B ligand (“RANKL)‐related state” are RANKL, parathyroid max capacity, and osteoprotegerin‐RANKL complex. A TGF‐β transforming growth factor‐beta; OB, osteoblast; OC, osteoclast; ROB, responding osteoblast.

where L1 and L2 are the lumped states representing RANK and RANKL, respectively, CMX is the RANK‐RANKL complex, and TGF is the active TGF‐β. The parameters R are the zero‐order production rate of X, k is the first‐order rate constant from X to Y, and d is the first‐order degradation rate constant of X. The parameters k int and K SS are the elimination rate constant of the denosumab‐RANKL complex and steady‐state constant for denosumab‐RANKL binding affinity, respectively,16 and C is the serum concentration of denosumab (same variables shown in Eqs. (11) and (12)).

In Eqs. (13) to 20, the parameters with “(t)” are time‐varying and contain some information that is included within nonlinear components of the original model. These time‐varying parameters were reconverted to the original nonlinear functions so the model can be handled easily in parameter estimation. The conversion to linear time‐invariant parameters using baseline values (at t = t 0) was also applied when the original nonlinear function was no longer available (for states that were lumped together with other states). Through this process, the time‐varying feature “(t)” has been removed from all the equations, and Eqs. (16) and (18) were modified using the original nonlinear functions. The final equations were then obtained in the form:

where γ is a sigmoidicity term influencing the steepness of response, a is the maximum anticipated response, ρ is the minimum anticipated response, and δ is the value of the driver that produces the half‐maximal response.11 The parameter values are shown in Table 1 as “before estimation.” The reduced model has in total 37 parameters including Eq. (4).

Table 1

Derived parameter values using lumping formulae (before estimation), and parameter estimates when using bone mineral density data until 1 year from the initiation of denosumab dosing (after estimation)

Parameter	Description	Value
		Before estimation	After estimation (%RSE)
R L1	Production rate of L1	75.0
k OB→L1	Rate constant expressing the effect of OB to the production of L1	55.3
k L2→L1	Rate constant from L2 to L1	160
k CMX→L	Rate constant from CMX to lumped state (L1 or L2)	0.112
d L1	Degradation rate constant of L1	0.970
R L2	Production rate of L2	0.000160	0.00337 (9.1%)
k OB→RANKL	Rate constant expressing the effect of OB to the production of RANKL	0.234
d L2	Degradation rate constant of L2	0.00110
d RANKL	Degradation rate constant of RANKL	0.00290
k int	Elimination rate constant of the denosumab‐RANKL complex	0.00795a
K SS	Steady‐state constant for denosumab‐RANKL binding affinity (ng/ml)	138a	63.4 (63.7%)
k L2→CMX	Rate constant from L2 to CMX	0.0000190
R OC	Production rate of OC	0.00000298
d OC	Degradation rate constant of OC	0.0292b	0.0898 (5.4%)
a 1	Maximum anticipated response of TGF to the degradation of OC	2.18b
ρ1	Minimum anticipated response of TGF to the degradation of OC	0.200b
δ1	Amount of TGF that produces the half‐maximal response to the degradation of OC	16.2b
γ1	Sigmoidicity term for the effect of TGF to the degradation of OC	1b
a 2	Maximum anticipated response of CMX to the degradation of OC	3.80b
ρ2	Minimum anticipated response of CMX to the degradation of OC	0.470b
γ2	Amount of CMX that produces the half‐maximal response to the degradation of OC	0.000013b
γ2	Sigmoidicity term for the effect of CMX to the degradation of OC	3.09b
k OC→TGF	Rate constant expressing the effect of OC to the production of TGF	5.66
d TGF	Degradation rate constant of TGF	0.0298
R ROB	Production rate of ROB	0.000003
a 3	Maximum anticipated response of TGF to the production of ROB	4.18b
ρ3	Minimum anticipated response of TGF to the production of ROB	0.202b
δ3	Amount of TGF that produces the half‐maximal response to the production of ROB	34.0b
γ3	Sigmoidicity term for the effect of TGF to the production of ROB	1.81b
k ROB→OB	Rate constant from ROB to OB	0.003
f	Fraction converting from ROB to SOB	0.06
d FOB	Degradation rate constant of FOB	0.01
d SOB	Degradation rate constant of SOB	0.000001
σ2	Variance of additive residual error for %BMD	–	0.382 (24.0%)

Parameters without RSE were fixed and not estimated.

L1 and L2 are the lumped states representing RANK and RANKL, respectively; CMX, RANK‐RANKL complex; TGF, active TGF‐β; FOB, fast osteoblast; SOB, slow osteoblast; ROB, responding osteoblast; OB, osteoblast; OC, osteoclast; BMD, bone mineral density; RANKL, receptor activator of nuclear factor kappa‐B ligand; RSE, relative standard error. Units for zero‐order production rate and other rate constants are amount/h and/h, respectively.

a Eextracted.14

bExtracted/derived.11

From the identifiability analysis d OC, R , K SS, and k out in Eq. (4) were considered estimable. Other parameters were fixed to the prior values provided by lumping formulae (shown in Table 1 as “before estimation”). An informal heuristic approach derived a consistent result for which estimating more parameters does not improve the model fit.

The extracted mean BMD data used for parameter estimation were obtained from 54, 47, 42, and 47 patients for denosumab 6‐month regimen (dosing every 6 months) groups of 14 mg, 60 mg, 100 mg, and 210 mg, respectively.22 Parameter estimates on the reduced model and the two empirical models, when using data until 1 year, are shown in Table 1 (as “after estimation”) and Table 2, respectively. The estimate of k out in the reduced model was 0.000533/h (24.6%RSE). Figure 5 a shows the fitting results for the three models. All models captured the BMD profile for all doses (administered every 6 months). Two empirical models, however, were not able to capture the increase in BMD after 1 year when dosing 60 mg (approved dose) of denosumab (Figure 5 b). In contrast, the reduced model performed well when used to extrapolate beyond the 1 year period. Similarly, an existing semimechanistic model also performed well (Figures and ). Parameter estimates on this model are also shown in the Table .

Table 2

Parameter estimates from two empirical models when using bone mineral density data until 1 year from the initiation of denosumab dosing

Parameter	Description	Value (%RSE)
Direct response model with time‐varying E max
E max	Maximum response	5.92 (11.8%)
T 1/2	Half‐life for the onset of drug effect (month)a	3.12 (24.0%)
C 50	Serum concentration of denosumab causing 50% of maximum response (ng/mL)	2.05 (30.9%)
σ2	Variance of additive residual error for %BMD	0.328 (29.4%)
Turnover model
I max	Maximum inhibitory effect of denosumab	0.0577 (15.2%)
k out	First‐order rate constant for BMD (/h)	0.000336 (26.5%)
C 50	Serum concentration of denosumab causing 50% of maximum response (ng/mL)	52.5 (67.6%)
σ2	Variance of additive residual error for %BMD	0.387 (26.4%)

BMD, bone mineral density; RSE, relative standard error.

aRate constant k is derived as k = ln(2)/T 1/2.

Figure 5

Fitting (a) and extrapolation (b) results of bone mineral density (BMD) response. A vertical dotted line represents 12 months. CFB, change from baseline.

Discussion

We explored simplification of a nonlinear QSP model by inductively linearizing the system followed by an automatic model order reduction technique. A detailed exploration of PK/PD and/or QSP models to bone biology was not considered as the main aim of this study was to explore a method for scale reduction rather than a model for bone biology. Proper lumping was chosen as a technique for model simplification in which the QSP model was linearized prior to lumping. With the lumping, parameter values as well as initial conditions of the lumped states can be obtained directly using algebraic formulae in Eqs. (7) and (8), as long as the system is linear. Most QSP systems, however, consist of a set of nonlinear ODEs. Standard linearization methods exist for such systems (e.g., a Taylor expansion), but their application is known to typically incur a relatively high degree of error.27 As an alternative, an inductive linearization has been shown to be applicable to nonlinear systems,15 and worked with arbitrary accuracy also for the relatively big system in this study (Figure 1).

After linearization, we applied an automated proper lumping method for model simplification. To find a final lumped model the composite criterion was applied. In the original paper for the criterion, the granularity in m was raised as a limitation17 (i.e., for modest full models (e.g., <10 states)) then the criterion may tend to choose simpler systems. However, for higher dimensional models (e.g., the BMD used in this study), the granularity of m is high and reasonable seven‐state or eight‐state models were appropriately chosen based on the criterion values. The eight‐state model provided an almost indistinguishable BMD profile to the original profile (Figure 3 b). An automatic search resulted in active TGF‐β and the RANK‐RANKL complex to remain unlumped in the final reduced model. This seems reasonable because the TGF‐β is known to play an important role for expressing nonlinear changes in BMD as it mediates the differentiation and apoptosis rates of bone remodeling OC and OB cells.12, 13, 14, 18 Due to the mechanism of action of denosumab (i.e., input in this lumping example) it is known that the RANK‐RANKL complex contributes to the process for BMD change and while the original two states before binding (RANK and RANKL) remain independent they also represent the final composite lumped states in the final reduced eight‐state model. We believe the final reduced eight‐state model (Figure 4) retains the mechanistic nature of the original model.

For extrapolation purposes of long‐term BMD responses, the final reduced model was fit to BMD data until 1 year from the initiation of denosumab dosing. Through structural identifiability analyses parameters directly relating to the input (denosumab as RANKL inhibitor) or output (BMD as a function of OB and OC) in the system were considered estimable. Using BMD data until 1 year, those parameters were precisely estimated with reasonable precision (relative standard error (RSE) <100%). Two typical top‐down empirical PK/PD models (a time‐varying direct response and turnover models 28, 29) were also fit to the 1‐year BMD data in order to extrapolate to the full 4‐year data. We expect the poor performance of these models when extrapolating beyond 1 year (Figure 5 b) because neither model contains the long‐term feedback processes (e.g., TGF‐β) that are present in the original and reduced models that would further modify the trajectory of the response. A potential method to improve predictive performance beyond 1 year would be to consider an additional biomarker of drug response (e.g., markers for bone formation/resorption),26, 30 or a longer training dataset.

An existing semimechanistic model was also used as a positive reference model.26 The model also captured the increase in BMD (see Figures ). This would support the simulation results obtained from the reduced model because both models incorporate the mechanistic nature of bone biology. Predicting long‐term responses from short‐term or middle‐term data is one of the most challenging but important tasks in drug development. If one fails to extrapolate accurately, then bias would affect projected effective doses for future studies (e.g., from phase II to phase III studies). The role of mechanistic models that are amenable to estimation is, therefore, critical for drug development.

Other semimechanistic models have also been reported.13, 14 The former model comprises an integrated systems biology‐continuum micromechanics approach (e.g., bone cell population model) and can be used to simulate bone volume (rather than bone mineral components). The latter model is a simplified version of a full mechanism‐based OB‐OC model31 (a part of the original bone biology model used in this study11). The model expressed the contribution of RANK and TGF‐β to the behavior of OC and OB by their receptor occupancies. This was not present in the reduced model developed in this study. The model, however, has a similar feedback mechanism in the OB‐OC relationship to the reduced model (i.e., RANK occupancy (function of OB) increases the amount of OC and TGF‐β occupancy (function of OC) stimulates the production of OB). Considering this and similar performance of the developed reduced model to the existing one,26 the key benefits of this work are the utility of the scale reduction methods. The method proposed here involves linearization of the original system and then simplification. The latter component can be automated and, therefore, does not require the user to have innate knowledge of the mathematical properties of the system. A tutorial for conducting linearization is under review—we also anticipate that automatic linearization is potentially feasible. We note that linearization is reversible, in that a nonlinear model can be linearized to attend to a particular operation that requires a linear form, and then unlinearized when the original form is desired. Going back to the model itself, the reduced model as well as the existing models13, 14 lack the mechanism of bone mineralization. Secondary mineralization is generally considered a primary mechanism of action of antiresorptive treatments. This feature has been included in another semimechanistic model,26 although the structure may be still empirical. These would suggest the development of mechanistic models, including such mineralization components.

The reduced model was developed by focusing on the increase in BMD (output) after denosumab dosing (input). It is anticipated that the structure of a reduced model will depend on the specific input‐output relationship. For example, when selecting bisphosphonates as the drug modifier, which directly stimulates the degradation of OC,32 the model may be further reduced by letting the states (directly) relating to RANK‐RANKL binding (i.e., L1, L2, and CMX be grouped together or grouped to other states). The reduced model from this work (Figure 4) would still be applicable for that input, although the model includes states that are likely to be of little importance. Finally, the model can be applied to different dose levels as long as the model predictions are based on the nonlinear reduced model rather than the linearized reduced model (to accommodate the impact of potential feedback mechanisms).

A potential additional method to scale reduction would be just doing an identifiability analysis and fixing unidentifiable parameters and then using the original systems model for parameter estimation. It may perform as well (or not) as scale reduction, but is not part of this work. Wendling et al.33 have considered this idea on a physiologically based PK model, but finally proceeded with reducing the system for parameter estimation.10 The main reason for such choices is the consideration that some parameter values will need to be fixed based on in vitro results, which may not reflect in vivo phenomena. By reducing the dimension of the model, this issue can be minimized as the number of parameters that will be required to be fixed will be dramatically reduced. They have also indicated that runtime and numerical instability may be an issue when using the original systems model for parameter estimation.

In conclusion, a nonlinear bone biology model was successfully reduced to an eight‐state model by inductively linearizing the system followed by automatic proper lumping. The reduced model described an increase in BMD after denosumab dosing while maintaining physiological meaning. The method used in this study is automatic, and can be applied directly to other multiscale models for developing a mechanism‐based structural model for future analyses. It is important to note, however, that the reduced model developed in this work was for exploratory purposes to evaluate the model simplification methodology and its potential utility, and would need further evaluation if intended to be used in drug development or clinically.

In addition to the supplements, the code for the linearized full model (automatic lumping code) and the code for the final form of the nonlinear reduced model will be made available with this article and as a download with the original bone model on GitHub. The linearized full model is amenable to future development and the code may be helpful for those who would like to lump the system for any other input‐output‐drug relations of interest.

Funding

No funding was received for this work.

Conflict of Interest

As an Associate Editor for CPT: Pharmacometrics & Systems Pharmacology, Stephen Duffull was not involved in the review or decision process for this article. The authors declared no competing interests for this work.

Author Contributions

C.H. wrote the manuscript. C.H. and S.D. designed the research. C.H. performed the research. C.H. analyzed the data. C.H. and S.D. contributed new reagents/analytical tools.

Figure S1. Fitting results of BMD response when using an existing semimechanistic model.24 A vertical dotted line represents 12 months. CFB, change from baseline.

Click here for additional data file.

Figure S2. Extrapolation results of BMD response when using an existing semimechanistic model.24 A vertical dotted line represents 12 months. CFB, change from baseline.

Click here for additional data file.

Table S1. Parameter estimates when using BMD data until 1 year from the initiation of denosumab dosing with an existing semimechanistic model.

Click here for additional data file.

Table S2. Digitized BMD data used in the manuscript.19

Click here for additional data file.

Appendix S1. Supplementary information.

Click here for additional data file.

Appendix S2. Model code for lumping.

Click here for additional data file.

Appendix S3. Model code for the final nonlinear lumped model.

Click here for additional data file.