Deriving mechanism‐based pharmacodynamic models by reducing quantitative systems pharmacology models: An application to warfarin
A good understanding of the determinants of drug‐effect size is essential for optimal drug dosing on an individual level. This is of particular interest for the widely used anticoagulant warfarin, as it has a narrow therapeutic window and large interindividual variability (IIV) in drug concentration and effect. A standard measure for the effect of warfarin therapy is the international normalized ratio (INR), a normalized coagulation time. A higher‐than‐desired INR is associated with an increased risk of major bleeding events, whereas with a lower‐than‐desired INR, thromboembolic events cannot effectively be prevented. The large IIV complicates optimal individual dosing, causing more than 10‐fold differences in the dose requirement. Current approaches to dose individualization include regression‐based algorithms (predicting the maintenance dosing) , and pharmacokinetic (PK)/pharmacodynamic (PD) model‐based approaches (predicting the warfarin effect to optimize the dose). , , , , However, a large proportion of the variability observed in warfarin dose requirements is not yet explained by the identified covariates in current approaches. Therefore, identification of further covariates or better dose adaptation after early INR measurements is required, for which PK/PD model‐based approaches are better suited than regression‐based approaches. In PK/PD model‐based approaches, dose adaptation typically relies on updating the model parameters based on INR measurements. How to use biomarkers, such as concentrations of coagulation factor, to further improve the model predictions is not always apparent. In contrast, quantitative systems pharmacology (QSP) models are well suited to identify possible drug targets or useful biomarkers. For warfarin, two QSP models have previously been used to study the treatment effect on the INR. , QSP models can often be used to simulate different scenarios, for example, warfarin treatment as well as envenomation after a snake bite. In the context of analyzing clinical data, however, the complexity of QSP models prevents straightforward parameter estimation. To leverage the knowledge in large‐scale QSP models, it would be desirable to systematically derive small‐scale, mechanism‐based PD models suitable for the analysis of clinical trials in a nonlinear mixed effect or Bayesian statistical context. In this article, we extended the model‐reduction method in Knöchel et al. from pure state elimination to state and parameter elimination, simplification of reactions, and analytic solution of model parts. In addition, model reduction is performed for a diverse virtual population, accounting for the expected variability in real‐world data. To this end, we leverage concepts from parameter identifiability, reaction simplification, and robust model reduction. , The proposed model‐reduction approach maintains a user‐specified threshold on the approximation error of the response for at least 95% of the individuals in a diverse virtual population. In application to a blood coagulation QSP model, we obtained a small‐scale warfarin‐INR model that predicts the INR in terms of the product of three coagulation factors, which are indirectly inhibited by warfarin. Under random variability and genotype heterogeneity, the small‐scale warfarin/INR model maintains a prespecified approximation quality to the original QSP model.   First, we describe the biological background and the blood coagulation QSP model, how it can be used to simulate the INR under warfarin treatment, and how we augmented the model to include variability. Then, we introduce the workflow to reduce the specific scenarios in the warfarin application and finally the general model‐reduction process and how it builds on different reduction methods. The model reduction was implemented in MATLAB 2021a and is accessible from https://doi.org/10.5281/zenodo.7417886 . Biological background Warfarin acts by inhibiting vitamin K epoxide reductase complex 1 ( VKORC1 ), thereby decreasing the rate at which vitamin K hydroquinone (VKH 2 ) is synthesized in the vitamin K cycle. The reduction in VKH 2 decreases the synthesis of important coagulation factors (e.g., II, VII, IX, and X) and thereby the coagulability of the blood. The warfarin effect can be measured by taking a blood sample and performing a prothrombin time (PT) test. The PT test is a typical way to measure an anticoagulant effect, in which coagulation is induced artificially in a blood sample by adding a defined, high amount of tissue factor (TF). The duration until the blood coagulation starts is denoted as PT. The INR is then defined as (1) INR = PT PT ref , where PT ref denotes the PT of a control sample. Simulating the INR with a QSP model The starting point for our analysis was the blood coagulation QSP model from Wajima et al. (see Figure for an illustration). In the QSP model, a one‐compartment PK model with oral absorption is used to simulate the warfarin concentration after multiple dosing. The QSP model does not consider the enantiomers in the racemic mixture separately. The warfarin effect on VKH 2 is modeled via a maximal effect ( E max ) model of the warfarin concentration. The PT, from which the INR is calculated by normalization, is defined by a threshold on the area under the curve (AUC) of fibrin ( F ) (2) PT = min τ ≥ 0 ∫ 0 τ F t d t ≥ δ . The threshold δ = 1500 s·nmol/L had been determined empirically to correspond to a 30% reduction in fibrinogen in Wajima et al., which is physiologically plausible for clotting. It results in a reasonable response of PT msub ≈ 11 s for the reference parameterization in the absence of warfarin. To model the INR under warfarin therapy, the QSP model is used in two scenarios: (i) the in vivo scenario to predict the action of warfarin on the coagulation factors and (ii) the in vitro scenario to predict the PT. The state vector from the in vivo scenario at a given time, divided by three to account for dilution, serves as the initial value for the in vitro scenario; this corresponds to taking a blood sample. The different scenarios can be simulated with different sets of parameter values in the QSP model. The simulation of the INR under treatment of 4 mg daily is illustrated in Figure . The need to repeatedly simulate the QSP model in two disconnected scenarios makes the computation costly, which is especially relevant for parameter estimation. Inclusion of IIV in the QSP model We aim for a PD model that can describe a population's variability and approximates the QSP model also for individuals deviating from the reference (reference denoting the parameterization reported in the original QSP model). To consider the model reduction under variability, we augmented the blood coagulation QSP model to include IIV. Because the reduced model is designed to guarantee error thresholds only for the considered population, we need to generate a diverse enough virtual population to cover a realistic variability. We first considered variability introduced by the genotypes of VKORC1 and cytochrome P450 isoenzyme 2C9 ( CYP2C9 ), by which warfarin is partly metabolized. Relative differences of the warfarin clearance parameter CL dependent on CYP2C9 genotype and of the warfarin sensitivity parameter IC 50 (the half maximal inhibitory concentration) dependent on VKORC1 genotype were adopted from a published PK/PD warfarin model ; see Supplementary Material Section for details. In addition, we considered random IIV on all parameters and initial values q i , independently distributed according to a log‐normal distribution around the reference values q ref (3) q i ∼ log N q ref i , 0.4 2 . A virtual population was generated by deterministically choosing genotypes such that the allele frequencies matched those reported in Hamberg et al. for the Warfarin Genetics study and randomly sampling the parameter values according to Equation . Mathematical model notation The blood coagulation QSP model is defined by a system of ordinary differential equations (ODEs): (4) d x t d t = f x t p , x 0 = x 0 + u . Here, x t ∈ ℝ n denotes the vector of state variables at time t ∈ [ 0,t end ], and p ∈ ℝ d denotes the vector of parameters. The initial value consists of the prestimulus state vector x 0 and the input/stimulus u . The input is the warfarin dose history in the in vivo scenario and the addition of TF in the in vitro scenario. We later also consider the extended parameter vector q = x 0 p ∈ ℝ n + d . The model comprises a system of n = 62 ODEs defined by the function f ; different sets of parameter values allow to simulate in vivo or in vitro settings. The response of interest is defined by either (5) y t = h x t or (6) y = h x ⋅ = min t ≥ 0 x t fulfills some specific condition . Equation  is useful for the in vivo setting, where the function h would be the determination of the INR dependent on the current state vector x ( t ) (via the in vitro scenario). Equation  is useful for the in vitro setting, where the response is the PT as defined in Equation . Workflow of reducing the scenarios considering IIV Ultimately, we are interested in a model reduction for the combined scenario, as seen in Figure on the right. To this end, we first reduced the in vitro scenario separately (but with an ensemble of virtual blood samples obtained from the in vivo scenario) and then reduced the in vivo scenario with the reduced in vitro model as response h . The reduction workflow specific to the combined scenario for warfarin treatment is summarized in Figure . A virtual population including covariate‐explained and random IIV was generated for the in vivo scenario as described in . We simulated with a fixed warfarin regimen of 4 mg daily for 30 days. To constrain the INR values to a clinically relevant range, we used a reduced dose of 1 mg daily for individuals for which the fixed dosing would lead to a steady‐state INR above 4. From the virtual population, we generated a diverse ensemble of blood samples for the in vitro scenario by virtually sampling blood at days t = 1, 4, 11, 30. The timepoints were chosen to give a good representation of the different stages of warfarin treatment in addition to reflecting the IIV. General automatic model‐reduction procedure The goal of the model reduction is to yield a model as simple as possible while ensuring that its response y red approximates the response y QSP of the original QSP model with a user‐defined maximal approximation error. We determined the approximation quality dependent on a randomly sampled virtual population. To be more robust regarding the realization of the random parameters and to account for possible unphysiological individuals, we propose to require the error threshold to hold for only 95% of the population. Therefore, we accept a reduced model if for at least 95% of the population, the maximal relative error is below a predefined threshold (here chosen to be 10%), that is, (7) where Q 0.95 denotes the 95% sample quantile, and the superscript ( i ) refers to the i th individual. We have chosen the maximal relative error, but other error measures can also be used. As a postprocessing step, we suggest to analyze the individuals for which the threshold is not attained. If they appear physiologic but have an error only slightly larger than the threshold, the reduced model might still be deemed acceptable for the population. If the model is not deemed acceptable, a higher quantile can be used to ensure that the previously excluded but critical individuals are accounted for. In our model‐reduction procedure, we differentiate between (i) model order reduction , in which the number of states/ODEs is reduced; and (ii) simplification of the functional form of reaction rates/ODEs; see Figure for an illustration. Both are fully automatized using the MATLAB Symbolic Toolbox in the model simplification. Next, we describe both steps in detail. Model order reduction For model order reduction , we employed the method proposed in Knöchel et al. In the model order reduction, each state variable is either classified as Environmental (env), that is, its dynamic is deemed unimportant, and it is approximated by a constant equal to its initial value; Negligible (neg), that is, considered completely unimportant and set constant to zero; or Dynamic (dyn), that is, its dynamics are modeled by an ODE, as in the original QSP model. An iterative approach is used to determine the classification of the states, with classifications of already assessed states constituting an intermediate reduced model that is updated after each further state is assessed. Each state variable is classified depending on the impact this would have on the approximation quality of the reduced model. If setting the state negligible or environmental meets the error criterion, the classification with the smaller error is accepted, otherwise, the state is classified dynamic. The states were ordered from lowest to highest importance using the sensitivity‐based input–response (ir) indices from Knöchel and considered for reduction in that order. If the r th state is the response and the i th state the input, the ir‐index (8) ir k t * = 1 t end ∫ t * t end S r , k t , t * 2 d t 1 2 ⋅ ∣ S k , i t * , t 0 ∣ is defined as the product of two terms: the left term represents the impact of the k th state on the output and is averaged over the remaining time interval, and the right term represents the impact of the input on the k th state. The larger the ir‐index, the more important the dynamics of the k th state variable for the ir relationship. If the response is not a state of the system but given by a response function h , the first sensitivity needs to be replaced by a term dependent on h . For the general definition and the sensitivity S n , m t 2 t 1 of the n th state at time t 2 to the m th state at time t 1 , see equations  and  in Supplementary Material Section . We performed the model order reduction with two important extensions compared with Knöchel et al. : (i) for the error criterion, we considered a virtual population (see Equation ) instead of a single reference parametrization; and (ii) we extended the definition of the ir‐indices to event‐type response functions. For the application to the in vitro scenario, the ir‐indices of the k th state with event‐type response are defined as (9) ir k t * = f x ref PT ref r − 1 ∣ S r , k PT ref t * ∣ ⋅ ∣ S k , i t * t 0 ∣ , where the subscript r refers to the AUC of fibrin and the subscript i to TF. For details of the derivation, see Supplementary Material Section . Model simplification In the model simplification step, we further simplified the ODEs of the dynamic states. This included parameter reduction and simplification of the functional form of the reaction rates. For parameter reduction, we measured the importance of a parameter [ p ] j for the response by the parameter sensitivity (10) P j = 1 t end ∫ 0 t end ∂ h x ⋅ t ∂ p j 2 d t 1 2 . We first ordered the parameters from the lowest to the highest parameter sensitivity. Then, a parameter was neglected if after the neglection the threshold in Equation  was still attained. This procedure is similar to the model order reduction procedure, in which the states were iteratively considered for neglection or elimination, ordered by their ir‐indices. Neglection is tested by setting the parameter to zero, or if zero‐neglection does not meet the error bound, to infinity. The parameter reduction can simplify the reactions of the remaining state variables significantly, as reaction rates drop out if containing a multiplicative factor that is neglected by setting to zero or a divisor that is neglected by setting to infinity. Note that the same holds for reaction rates in which a neglected state is a multiplicative factor. This applies to all ODEs in which reaction rates with the respective parameter or state are part. A typical source of nonlinearity in QSP models are Michaelis–Menten reaction rates (11) r = V max ⋅ S K M + S , that include linear dependence on [S] or constant behavior. If S ≫ K M over the time span of interest, then K M can be set to zero in the parameter reduction as this introduces only a slight error, and the reaction rate thus becomes constant, simplifying the ODEs in which the rate is part. To simplify ODEs with Michaelis–Menten kinetics also in the case that S ≪ K M , we considered the remaining reaction rates for Taylor approximation in x 0 , which linearizes the Michaelis‐Menten kinetics if S 0 = 0 . Note that if one reaction rate is part of multiple ODEs, they are considered separately, ODE‐wise; however, this did not occur in our application. As with the other reduction steps, the simplification by Taylor approximation was only realized if this did not violate the error criterion (Equation ). After the model‐reduction procedure, we evaluated if a reduced model could be solved analytically, and in the in vitro scenario we conducted a postprocessing step to obtain an analytic solution for the INR.   Warfarin acts by inhibiting vitamin K epoxide reductase complex 1 ( VKORC1 ), thereby decreasing the rate at which vitamin K hydroquinone (VKH 2 ) is synthesized in the vitamin K cycle. The reduction in VKH 2 decreases the synthesis of important coagulation factors (e.g., II, VII, IX, and X) and thereby the coagulability of the blood. The warfarin effect can be measured by taking a blood sample and performing a prothrombin time (PT) test. The PT test is a typical way to measure an anticoagulant effect, in which coagulation is induced artificially in a blood sample by adding a defined, high amount of tissue factor (TF). The duration until the blood coagulation starts is denoted as PT. The INR is then defined as (1) INR = PT PT ref , where PT ref denotes the PT of a control sample.   INR with a QSP model The starting point for our analysis was the blood coagulation QSP model from Wajima et al. (see Figure for an illustration). In the QSP model, a one‐compartment PK model with oral absorption is used to simulate the warfarin concentration after multiple dosing. The QSP model does not consider the enantiomers in the racemic mixture separately. The warfarin effect on VKH 2 is modeled via a maximal effect ( E max ) model of the warfarin concentration. The PT, from which the INR is calculated by normalization, is defined by a threshold on the area under the curve (AUC) of fibrin ( F ) (2) PT = min τ ≥ 0 ∫ 0 τ F t d t ≥ δ . The threshold δ = 1500 s·nmol/L had been determined empirically to correspond to a 30% reduction in fibrinogen in Wajima et al., which is physiologically plausible for clotting. It results in a reasonable response of PT msub ≈ 11 s for the reference parameterization in the absence of warfarin. To model the INR under warfarin therapy, the QSP model is used in two scenarios: (i) the in vivo scenario to predict the action of warfarin on the coagulation factors and (ii) the in vitro scenario to predict the PT. The state vector from the in vivo scenario at a given time, divided by three to account for dilution, serves as the initial value for the in vitro scenario; this corresponds to taking a blood sample. The different scenarios can be simulated with different sets of parameter values in the QSP model. The simulation of the INR under treatment of 4 mg daily is illustrated in Figure . The need to repeatedly simulate the QSP model in two disconnected scenarios makes the computation costly, which is especially relevant for parameter estimation.   IIV in the QSP model We aim for a PD model that can describe a population's variability and approximates the QSP model also for individuals deviating from the reference (reference denoting the parameterization reported in the original QSP model). To consider the model reduction under variability, we augmented the blood coagulation QSP model to include IIV. Because the reduced model is designed to guarantee error thresholds only for the considered population, we need to generate a diverse enough virtual population to cover a realistic variability. We first considered variability introduced by the genotypes of VKORC1 and cytochrome P450 isoenzyme 2C9 ( CYP2C9 ), by which warfarin is partly metabolized. Relative differences of the warfarin clearance parameter CL dependent on CYP2C9 genotype and of the warfarin sensitivity parameter IC 50 (the half maximal inhibitory concentration) dependent on VKORC1 genotype were adopted from a published PK/PD warfarin model ; see Supplementary Material Section for details. In addition, we considered random IIV on all parameters and initial values q i , independently distributed according to a log‐normal distribution around the reference values q ref (3) q i ∼ log N q ref i , 0.4 2 . A virtual population was generated by deterministically choosing genotypes such that the allele frequencies matched those reported in Hamberg et al. for the Warfarin Genetics study and randomly sampling the parameter values according to Equation .   The blood coagulation QSP model is defined by a system of ordinary differential equations (ODEs): (4) d x t d t = f x t p , x 0 = x 0 + u . Here, x t ∈ ℝ n denotes the vector of state variables at time t ∈ [ 0,t end ], and p ∈ ℝ d denotes the vector of parameters. The initial value consists of the prestimulus state vector x 0 and the input/stimulus u . The input is the warfarin dose history in the in vivo scenario and the addition of TF in the in vitro scenario. We later also consider the extended parameter vector q = x 0 p ∈ ℝ n + d . The model comprises a system of n = 62 ODEs defined by the function f ; different sets of parameter values allow to simulate in vivo or in vitro settings. The response of interest is defined by either (5) y t = h x t or (6) y = h x ⋅ = min t ≥ 0 x t fulfills some specific condition . Equation  is useful for the in vivo setting, where the function h would be the determination of the INR dependent on the current state vector x ( t ) (via the in vitro scenario). Equation  is useful for the in vitro setting, where the response is the PT as defined in Equation .   IIV Ultimately, we are interested in a model reduction for the combined scenario, as seen in Figure on the right. To this end, we first reduced the in vitro scenario separately (but with an ensemble of virtual blood samples obtained from the in vivo scenario) and then reduced the in vivo scenario with the reduced in vitro model as response h . The reduction workflow specific to the combined scenario for warfarin treatment is summarized in Figure . A virtual population including covariate‐explained and random IIV was generated for the in vivo scenario as described in . We simulated with a fixed warfarin regimen of 4 mg daily for 30 days. To constrain the INR values to a clinically relevant range, we used a reduced dose of 1 mg daily for individuals for which the fixed dosing would lead to a steady‐state INR above 4. From the virtual population, we generated a diverse ensemble of blood samples for the in vitro scenario by virtually sampling blood at days t = 1, 4, 11, 30. The timepoints were chosen to give a good representation of the different stages of warfarin treatment in addition to reflecting the IIV.   The goal of the model reduction is to yield a model as simple as possible while ensuring that its response y red approximates the response y QSP of the original QSP model with a user‐defined maximal approximation error. We determined the approximation quality dependent on a randomly sampled virtual population. To be more robust regarding the realization of the random parameters and to account for possible unphysiological individuals, we propose to require the error threshold to hold for only 95% of the population. Therefore, we accept a reduced model if for at least 95% of the population, the maximal relative error is below a predefined threshold (here chosen to be 10%), that is, (7) where Q 0.95 denotes the 95% sample quantile, and the superscript ( i ) refers to the i th individual. We have chosen the maximal relative error, but other error measures can also be used. As a postprocessing step, we suggest to analyze the individuals for which the threshold is not attained. If they appear physiologic but have an error only slightly larger than the threshold, the reduced model might still be deemed acceptable for the population. If the model is not deemed acceptable, a higher quantile can be used to ensure that the previously excluded but critical individuals are accounted for. In our model‐reduction procedure, we differentiate between (i) model order reduction , in which the number of states/ODEs is reduced; and (ii) simplification of the functional form of reaction rates/ODEs; see Figure for an illustration. Both are fully automatized using the MATLAB Symbolic Toolbox in the model simplification. Next, we describe both steps in detail. Model order reduction For model order reduction , we employed the method proposed in Knöchel et al. In the model order reduction, each state variable is either classified as Environmental (env), that is, its dynamic is deemed unimportant, and it is approximated by a constant equal to its initial value; Negligible (neg), that is, considered completely unimportant and set constant to zero; or Dynamic (dyn), that is, its dynamics are modeled by an ODE, as in the original QSP model. An iterative approach is used to determine the classification of the states, with classifications of already assessed states constituting an intermediate reduced model that is updated after each further state is assessed. Each state variable is classified depending on the impact this would have on the approximation quality of the reduced model. If setting the state negligible or environmental meets the error criterion, the classification with the smaller error is accepted, otherwise, the state is classified dynamic. The states were ordered from lowest to highest importance using the sensitivity‐based input–response (ir) indices from Knöchel and considered for reduction in that order. If the r th state is the response and the i th state the input, the ir‐index (8) ir k t * = 1 t end ∫ t * t end S r , k t , t * 2 d t 1 2 ⋅ ∣ S k , i t * , t 0 ∣ is defined as the product of two terms: the left term represents the impact of the k th state on the output and is averaged over the remaining time interval, and the right term represents the impact of the input on the k th state. The larger the ir‐index, the more important the dynamics of the k th state variable for the ir relationship. If the response is not a state of the system but given by a response function h , the first sensitivity needs to be replaced by a term dependent on h . For the general definition and the sensitivity S n , m t 2 t 1 of the n th state at time t 2 to the m th state at time t 1 , see equations  and  in Supplementary Material Section . We performed the model order reduction with two important extensions compared with Knöchel et al. : (i) for the error criterion, we considered a virtual population (see Equation ) instead of a single reference parametrization; and (ii) we extended the definition of the ir‐indices to event‐type response functions. For the application to the in vitro scenario, the ir‐indices of the k th state with event‐type response are defined as (9) ir k t * = f x ref PT ref r − 1 ∣ S r , k PT ref t * ∣ ⋅ ∣ S k , i t * t 0 ∣ , where the subscript r refers to the AUC of fibrin and the subscript i to TF. For details of the derivation, see Supplementary Material Section . Model simplification In the model simplification step, we further simplified the ODEs of the dynamic states. This included parameter reduction and simplification of the functional form of the reaction rates. For parameter reduction, we measured the importance of a parameter [ p ] j for the response by the parameter sensitivity (10) P j = 1 t end ∫ 0 t end ∂ h x ⋅ t ∂ p j 2 d t 1 2 . We first ordered the parameters from the lowest to the highest parameter sensitivity. Then, a parameter was neglected if after the neglection the threshold in Equation  was still attained. This procedure is similar to the model order reduction procedure, in which the states were iteratively considered for neglection or elimination, ordered by their ir‐indices. Neglection is tested by setting the parameter to zero, or if zero‐neglection does not meet the error bound, to infinity. The parameter reduction can simplify the reactions of the remaining state variables significantly, as reaction rates drop out if containing a multiplicative factor that is neglected by setting to zero or a divisor that is neglected by setting to infinity. Note that the same holds for reaction rates in which a neglected state is a multiplicative factor. This applies to all ODEs in which reaction rates with the respective parameter or state are part. A typical source of nonlinearity in QSP models are Michaelis–Menten reaction rates (11) r = V max ⋅ S K M + S , that include linear dependence on [S] or constant behavior. If S ≫ K M over the time span of interest, then K M can be set to zero in the parameter reduction as this introduces only a slight error, and the reaction rate thus becomes constant, simplifying the ODEs in which the rate is part. To simplify ODEs with Michaelis–Menten kinetics also in the case that S ≪ K M , we considered the remaining reaction rates for Taylor approximation in x 0 , which linearizes the Michaelis‐Menten kinetics if S 0 = 0 . Note that if one reaction rate is part of multiple ODEs, they are considered separately, ODE‐wise; however, this did not occur in our application. As with the other reduction steps, the simplification by Taylor approximation was only realized if this did not violate the error criterion (Equation ). After the model‐reduction procedure, we evaluated if a reduced model could be solved analytically, and in the in vitro scenario we conducted a postprocessing step to obtain an analytic solution for the INR.   For model order reduction , we employed the method proposed in Knöchel et al. In the model order reduction, each state variable is either classified as Environmental (env), that is, its dynamic is deemed unimportant, and it is approximated by a constant equal to its initial value; Negligible (neg), that is, considered completely unimportant and set constant to zero; or Dynamic (dyn), that is, its dynamics are modeled by an ODE, as in the original QSP model. An iterative approach is used to determine the classification of the states, with classifications of already assessed states constituting an intermediate reduced model that is updated after each further state is assessed. Each state variable is classified depending on the impact this would have on the approximation quality of the reduced model. If setting the state negligible or environmental meets the error criterion, the classification with the smaller error is accepted, otherwise, the state is classified dynamic. The states were ordered from lowest to highest importance using the sensitivity‐based input–response (ir) indices from Knöchel and considered for reduction in that order. If the r th state is the response and the i th state the input, the ir‐index (8) ir k t * = 1 t end ∫ t * t end S r , k t , t * 2 d t 1 2 ⋅ ∣ S k , i t * , t 0 ∣ is defined as the product of two terms: the left term represents the impact of the k th state on the output and is averaged over the remaining time interval, and the right term represents the impact of the input on the k th state. The larger the ir‐index, the more important the dynamics of the k th state variable for the ir relationship. If the response is not a state of the system but given by a response function h , the first sensitivity needs to be replaced by a term dependent on h . For the general definition and the sensitivity S n , m t 2 t 1 of the n th state at time t 2 to the m th state at time t 1 , see equations  and  in Supplementary Material Section . We performed the model order reduction with two important extensions compared with Knöchel et al. : (i) for the error criterion, we considered a virtual population (see Equation ) instead of a single reference parametrization; and (ii) we extended the definition of the ir‐indices to event‐type response functions. For the application to the in vitro scenario, the ir‐indices of the k th state with event‐type response are defined as (9) ir k t * = f x ref PT ref r − 1 ∣ S r , k PT ref t * ∣ ⋅ ∣ S k , i t * t 0 ∣ , where the subscript r refers to the AUC of fibrin and the subscript i to TF. For details of the derivation, see Supplementary Material Section .   In the model simplification step, we further simplified the ODEs of the dynamic states. This included parameter reduction and simplification of the functional form of the reaction rates. For parameter reduction, we measured the importance of a parameter [ p ] j for the response by the parameter sensitivity (10) P j = 1 t end ∫ 0 t end ∂ h x ⋅ t ∂ p j 2 d t 1 2 . We first ordered the parameters from the lowest to the highest parameter sensitivity. Then, a parameter was neglected if after the neglection the threshold in Equation  was still attained. This procedure is similar to the model order reduction procedure, in which the states were iteratively considered for neglection or elimination, ordered by their ir‐indices. Neglection is tested by setting the parameter to zero, or if zero‐neglection does not meet the error bound, to infinity. The parameter reduction can simplify the reactions of the remaining state variables significantly, as reaction rates drop out if containing a multiplicative factor that is neglected by setting to zero or a divisor that is neglected by setting to infinity. Note that the same holds for reaction rates in which a neglected state is a multiplicative factor. This applies to all ODEs in which reaction rates with the respective parameter or state are part. A typical source of nonlinearity in QSP models are Michaelis–Menten reaction rates (11) r = V max ⋅ S K M + S , that include linear dependence on [S] or constant behavior. If S ≫ K M over the time span of interest, then K M can be set to zero in the parameter reduction as this introduces only a slight error, and the reaction rate thus becomes constant, simplifying the ODEs in which the rate is part. To simplify ODEs with Michaelis–Menten kinetics also in the case that S ≪ K M , we considered the remaining reaction rates for Taylor approximation in x 0 , which linearizes the Michaelis‐Menten kinetics if S 0 = 0 . Note that if one reaction rate is part of multiple ODEs, they are considered separately, ODE‐wise; however, this did not occur in our application. As with the other reduction steps, the simplification by Taylor approximation was only realized if this did not violate the error criterion (Equation ). After the model‐reduction procedure, we evaluated if a reduced model could be solved analytically, and in the in vitro scenario we conducted a postprocessing step to obtain an analytic solution for the INR.   We applied our model‐reduction approach to a QSP model of the effect of warfarin on blood coagulation (see also Figures and for illustrations). Model‐reduction results for in vitro and in vivo scenarios The first step to obtaining a reduced model for the original QSP model with 62 ODEs and 174 parameters is to apply the model order reduction. The resulting reduced order models for the in vitro and in vivo scenarios, each involving six ODEs, are shown in Figure . The states of the reduced in vitro model are colored blue, and the states of the reduced in vivo model are colored orange, with the dynamic states (those modeled by an ODE) in darker shades and the environmental states (set constant to their initial value) in lighter shades. As a result of the parameter reduction, 8 and 13 parameters remained in the reduced in vitro model and reduced in vivo model, respectively, compared with 174 parameters each in the original QSP model. Of the 13 parameters in the in vivo model, three are synthesis rates and are determined via parameter interdependencies by the prestimulus concentrations as stated in Supplementary Material Section , leaving 10 actual parameters. After checking for linearization, only linear reactions remained in the reduced in vitro model, whereas the reduced in vivo model was not further simplified. In the virtual ensemble of blood samples, the reduced in vitro model approximates the original QSP model with the INR as response in accordance with the error criterion (Equation ) in ≥99.8% of cases. In the virtual population, the reduced in vivo model approximates the QSP model with the INR equation as response in accordance with the error criterion (Equation ) in 100% of the cases. As the reduced in vitro model consisted only of linear ODEs, it allowed for an analytic solution; see Supplementary Material Section for solutions for all states. Specifically, the analytic solution for the concentration of fibrin, on which the PT definition (Equation ) directly depends, is given by (12) F t = II 0 ⋅ VII 0 ⋅ X 0 ⋅ Fg 0 ⋅ c 1 ⋅ q t + p t ⋅ exp ‐ c 2 t , where c 1 , c 2 are positive constants, Fg 0 is the initial fibrinogen concentration, and p ( t ) and q ( t ) are cubic and linear polynomials as a function of the in vitro time t . Recall that the initial concentrations II 0 , VII 0 , and X 0 in a blood sample were obtained from an individual's concentrations, II i t * , VII i t * , X i t * at some in vivo time t *. Notably, the fibrin concentration depends on the three warfarin‐dependent coagulation factor concentrations only via their product II 0 ⋅ VII 0 ⋅ X 0 . The results presented in this section were obtained using the automatic model‐reduction procedure. The insights gave rise to the manual approximation presented in the next section. An algebraic INR equation As the INR depends on fibrin via Equation , it can only depend on the coagulation factors via their product II 0 ⋅ VII 0 ⋅ X 0 . We can thus represent the individually normalized INR γ as a function of the individually normalized factor concentration product to characterize the effect of warfarin. Plotting the normalized INR against the normalized product in a log–log plot in Figure , we find that the individual INR is well approximated in the most relevant INR range of 2 to 3 by (13) INR i t INR i 0 = II i t II i 0 ⋅ VII i t VII i 0 ⋅ X i t X i 0 γ , with the exponent chosen to be γ = −0.1975. The reduced in vivo model and the INR equation can be combined to yield a small‐scale warfarin/INR model as the reduced model for the combined scenario. Small‐scale warfarin/ INR model has good approximation quality under IIV The small‐scale warfarin/INR model (see Figure ) simulates the combined scenario by accounting for the effect of warfarin on the coagulation Factors II, VII, and X via inhibition of VKH 2 and translating the relative reduction in the product of the coagulation factor concentrations into an increase in the INR. The small‐scale warfarin/INR model consists of 11 parameters (10 from the in vivo model and the exponent in the INR equation). We evaluated the INR approximation quality of the small‐scale warfarin/INR model for a virtual population including genotype‐induced and unexplained random IIV. While during the model reduction only either the in vitro or the in vivo scenario was considered at a time, we now evaluate the approximation quality of the small‐scale warfarin/INR model to the QSP model for the combined scenario. Figure (top) shows the INR simulated with the QSP model versus the small‐scale warfarin/INR model for the virtual population. The very good approximation quality shows that the small‐scale warfarin/INR model robustly predicts the INR in a heterogeneous population. In addition, we assessed the approximation in INR‐time profiles for different genotypes of CYP2C9 and VKORC1 . The comparisons in Figure (bottom) show simulations of fixed warfarin dosing (4 mg daily) for different genotypes. The small‐scale warfarin/INR model approximates the QSP model well for all genotypes. The poorest approximation quality is observed for larger INRs ( CYP2C9*3/*3 simulation in Figure [bottom left]); however, the relative error is still below 10%. Moreover, in clinical practice, the dose would be reduced when such a high INR is observed, and simulating the same individual with a smaller dose resulted in a much improved approximation. The small‐scale warfarin/INR model allows to determine the warfarin effect on the INR computationally much more efficiently, as all additional simulations of the QSP model in the in vitro scenario are avoided. Biomarker proposal to predict steady‐state INR early In the small‐scale warfarin/INR model, the INR is calculated from the relative concentrations of the coagulation Factors II, VII, and X; therefore, it is natural to consider them in the search of useful biomarkers. Assume that the relative reductions from their pretreatment values II 0 , VII 0 , X 0 to their steady‐state values II ss , VII ss , X ss are related via (14) This allows to predict the steady‐state INR ss from the INR Equation  by (15) INR ss = INR 0 ⋅ a ⋅ b ⋅ r ss ⋅ r ss ⋅ r ss γ . Of note, in the QSP model, it is assumed that a = b = 1. Any of the three factors can be used to determine r ss once it is in steady state. Because Factor VII has the shortest half‐life (~ 6 h), it will adapt fastest, much faster than Factors II (~ 69 h) and X (~ 39 h), suggesting to use (16) r ss = VII ss VII 0 to predict the steady‐state INR value from measurements of VII . Importantly, Factor VII measurements already account for the interindividual differences in the vitamin K cycle and in the PK and thus allow to assess their impact on steady‐state INR. At early timepoints, measuring the Factor VII concentrations and calculating the expected steady‐state INR using Equations  and  should be more informative for adapting the warfarin dose than only INR measurements.   The first step to obtaining a reduced model for the original QSP model with 62 ODEs and 174 parameters is to apply the model order reduction. The resulting reduced order models for the in vitro and in vivo scenarios, each involving six ODEs, are shown in Figure . The states of the reduced in vitro model are colored blue, and the states of the reduced in vivo model are colored orange, with the dynamic states (those modeled by an ODE) in darker shades and the environmental states (set constant to their initial value) in lighter shades. As a result of the parameter reduction, 8 and 13 parameters remained in the reduced in vitro model and reduced in vivo model, respectively, compared with 174 parameters each in the original QSP model. Of the 13 parameters in the in vivo model, three are synthesis rates and are determined via parameter interdependencies by the prestimulus concentrations as stated in Supplementary Material Section , leaving 10 actual parameters. After checking for linearization, only linear reactions remained in the reduced in vitro model, whereas the reduced in vivo model was not further simplified. In the virtual ensemble of blood samples, the reduced in vitro model approximates the original QSP model with the INR as response in accordance with the error criterion (Equation ) in ≥99.8% of cases. In the virtual population, the reduced in vivo model approximates the QSP model with the INR equation as response in accordance with the error criterion (Equation ) in 100% of the cases. As the reduced in vitro model consisted only of linear ODEs, it allowed for an analytic solution; see Supplementary Material Section for solutions for all states. Specifically, the analytic solution for the concentration of fibrin, on which the PT definition (Equation ) directly depends, is given by (12) F t = II 0 ⋅ VII 0 ⋅ X 0 ⋅ Fg 0 ⋅ c 1 ⋅ q t + p t ⋅ exp ‐ c 2 t , where c 1 , c 2 are positive constants, Fg 0 is the initial fibrinogen concentration, and p ( t ) and q ( t ) are cubic and linear polynomials as a function of the in vitro time t . Recall that the initial concentrations II 0 , VII 0 , and X 0 in a blood sample were obtained from an individual's concentrations, II i t * , VII i t * , X i t * at some in vivo time t *. Notably, the fibrin concentration depends on the three warfarin‐dependent coagulation factor concentrations only via their product II 0 ⋅ VII 0 ⋅ X 0 . The results presented in this section were obtained using the automatic model‐reduction procedure. The insights gave rise to the manual approximation presented in the next section.   INR equation As the INR depends on fibrin via Equation , it can only depend on the coagulation factors via their product II 0 ⋅ VII 0 ⋅ X 0 . We can thus represent the individually normalized INR γ as a function of the individually normalized factor concentration product to characterize the effect of warfarin. Plotting the normalized INR against the normalized product in a log–log plot in Figure , we find that the individual INR is well approximated in the most relevant INR range of 2 to 3 by (13) INR i t INR i 0 = II i t II i 0 ⋅ VII i t VII i 0 ⋅ X i t X i 0 γ , with the exponent chosen to be γ = −0.1975. The reduced in vivo model and the INR equation can be combined to yield a small‐scale warfarin/INR model as the reduced model for the combined scenario.   INR model has good approximation quality under IIV The small‐scale warfarin/INR model (see Figure ) simulates the combined scenario by accounting for the effect of warfarin on the coagulation Factors II, VII, and X via inhibition of VKH 2 and translating the relative reduction in the product of the coagulation factor concentrations into an increase in the INR. The small‐scale warfarin/INR model consists of 11 parameters (10 from the in vivo model and the exponent in the INR equation). We evaluated the INR approximation quality of the small‐scale warfarin/INR model for a virtual population including genotype‐induced and unexplained random IIV. While during the model reduction only either the in vitro or the in vivo scenario was considered at a time, we now evaluate the approximation quality of the small‐scale warfarin/INR model to the QSP model for the combined scenario. Figure (top) shows the INR simulated with the QSP model versus the small‐scale warfarin/INR model for the virtual population. The very good approximation quality shows that the small‐scale warfarin/INR model robustly predicts the INR in a heterogeneous population. In addition, we assessed the approximation in INR‐time profiles for different genotypes of CYP2C9 and VKORC1 . The comparisons in Figure (bottom) show simulations of fixed warfarin dosing (4 mg daily) for different genotypes. The small‐scale warfarin/INR model approximates the QSP model well for all genotypes. The poorest approximation quality is observed for larger INRs ( CYP2C9*3/*3 simulation in Figure [bottom left]); however, the relative error is still below 10%. Moreover, in clinical practice, the dose would be reduced when such a high INR is observed, and simulating the same individual with a smaller dose resulted in a much improved approximation. The small‐scale warfarin/INR model allows to determine the warfarin effect on the INR computationally much more efficiently, as all additional simulations of the QSP model in the in vitro scenario are avoided.   INR early In the small‐scale warfarin/INR model, the INR is calculated from the relative concentrations of the coagulation Factors II, VII, and X; therefore, it is natural to consider them in the search of useful biomarkers. Assume that the relative reductions from their pretreatment values II 0 , VII 0 , X 0 to their steady‐state values II ss , VII ss , X ss are related via (14) This allows to predict the steady‐state INR ss from the INR Equation  by (15) INR ss = INR 0 ⋅ a ⋅ b ⋅ r ss ⋅ r ss ⋅ r ss γ . Of note, in the QSP model, it is assumed that a = b = 1. Any of the three factors can be used to determine r ss once it is in steady state. Because Factor VII has the shortest half‐life (~ 6 h), it will adapt fastest, much faster than Factors II (~ 69 h) and X (~ 39 h), suggesting to use (16) r ss = VII ss VII 0 to predict the steady‐state INR value from measurements of VII . Importantly, Factor VII measurements already account for the interindividual differences in the vitamin K cycle and in the PK and thus allow to assess their impact on steady‐state INR. At early timepoints, measuring the Factor VII concentrations and calculating the expected steady‐state INR using Equations  and  should be more informative for adapting the warfarin dose than only INR measurements.   We introduce a model‐reduction method that allows deriving a small‐scale warfarin/INR model from a blood coagulation QSP model. The small‐scale warfarin/INR model calculates the INR via a linear relationship in the log–log space from concentrations of Factors II, VII, and X. The INR prediction of the small‐scale warfarin/INR model approximates the INR prediction of the QSP model to within 10% for more than 99% of a diverse virtual population. Without clinical data, but solely based on the small‐scale warfarin/INR model, we identified Factor VII as a possible biomarker. The model‐reduction method, including model order reduction and model simplification, is fully automatic and can be automatically applied to other QSP models as well. The INR equation was derived manually using the automatically obtained analytic equation for fibrin. Application of the method requires defining an input (e.g., the dosing history) and an output function (e.g., the effect) for the model as well as generating a virtual population; together they define the scenario in which the model can be applied. Different scenarios typically require different reduced models , ; in this article, we concentrated on a scenario modeling warfarin treatment and subsequent INR measurements. Importantly, it is based on the common PT test and a range of INR values below 4; the virtual population includes polymorphisms in CYP2C9 ( *1 , *2 , *3 ) as well as polymorphisms in VKORC1 (A, G). Many structurally different warfarin PK/PD models are available in the literature, , , , , , which makes it difficult to judge which of them to use. A QSP model as a starting point, however, makes the underlying processes and assumptions explicit, so they can be discussed and tested. Similar to our model, existing empirical warfarin PK/PD models account for factor concentration‐time courses (although not always explicitly) mostly via an E max model for inhibition. They subsequently translate decreased factor concentrations into increased INRs via an INR equation. Our mechanism‐based INR representation reinforces the interpretation that the model components in Hamberg et al. (the two parallel transit chains) represent relative concentrations of Factors II and VII. Moreover, our INR equation derivation shows that the INR equations represent the in vitro processes in the PT test. The small‐scale warfarin/INR model is also comparable with empirical PK/PD models in terms of numbers of structural parameters (11 parameters compared with 12 parameters in Hamberg et al., 11 parameters in Xue et al., and nine parameters in Ohara et al. ). For a better comparison, we excluded parameters related to covariate effects, as the models include a different number of covariates. A linear relationship in the log‐space between INR and coagulation factor concentrations, as we have identified in our model, has previously been observed under stable therapy of acenocoumarol, another vitamin K antagonist. Factor VII has been proposed as a biomarker in warfarin treatment before. , Also for the vitamin K antagonist acenocoumarol, the variability in steady‐state INR was well explained by only Factor VII concentrations. Based on our low‐scale warfarin‐INR PD model, we can offer a simple equation to improve the early prediction of steady‐state INR based on steady‐state Factor VII concentrations. Notably, the simple equation helps to make the assumptions (e.g., about the factor's relative reduction from pretreatment to steady‐state values) explicit, under which we expect this approach to give good results. The impact of violation of assumptions can be easily examined, and assumptions might be weakened, for example, we do not require the factor's relative reduction to be the same, but only that the ratio is known. To assess the impact of dose individualization based on the biomarker relationship identified in this analysis, the feasibility of performing these measurements in clinical practice and cost‐effectiveness remains to be evaluated. To reduce the QSP model to a low‐scale PD model enabling parameter estimation and accounting for parameter variability, we combined multiple reduction approaches: state and parameter reduction, , reaction reduction, and robust model reduction. , The combination of different reduction approaches is essential to reduce QSP models with different properties. Lumping , , and time scale separation can easily be included in our approach but were not used in this analysis as they did not substantially improve the model reduction. Another method to reduce complex rational rate expressions is term‐based identifiability analysis. Because the warfarin example included only linear or Michaelis–Menten–type reaction rates, the simpler approach of reaction reduction by first‐order Taylor approximation was used and resulted in a very good approximation. In contrast to reducing the blood coagulation QSP model using a neural network approach, the model‐reduction approach presented in this article is fully parametric, enabling biomarker identification. Variability was introduced into the QSP model by considering different CYP2C9 and VKORC1 genotypes and randomly distributed parameters according to a log‐normal distribution with a 40% coefficient of variation. We assumed parameter variability to be uncorrelated due to the lack of more detailed knowledge on the correlation structure. This, however, is not a required feature of the virtual population, and correlated parameter variability can be included dependent on the state of knowledge. The population and input should be chosen such that the parameter sets represent therapeutically relevant scenarios because the reduced model guarantees an error threshold only for the considered population. If the population includes unrealistic parameter sets or inputs (e.g., dosing history), this might unnecessarily impair the reducibility. In the model reduction, the approximation quality was assessed in a virtual population of 1000 individuals. The resulting virtual population‐based reduction approach to account for variability is similar as in Dokoumetzidis and Aarons. However, instead of restricting the approximation error only for the population mean, we focused on individual prediction and chose to ensure an acceptable approximation for at least 95% of the virtual population. Using the 95% threshold, we addressed the existence of possible unphysiological parameter combinations and thus unphysiological responses in the virtual population, which is a known problem in the automatic generation of virtual populations (see, e.g., Duffull and Gulati ). Consequently, we suggest to a posteriori examine the characteristics of the excluded virtual individuals to avoid excluding critical but uncommon individuals. In our case, the obtained small‐scale warfarin/INR model actually attains the error threshold for more than 99% of the population. We judged the excluded individuals (5 of 1000) as uncritical because they are still relatively well approximated, with errors between 10% and 13%. QSP and physiologically based PK models have previously been used to predict individual outcomes, for example, the INR or drug exposure. , A small‐scale model that predicts the response well for a diverse population also enables dose adaptation using full Bayesian updating. The reduced model, together with the included random unexplained IIV, could either directly act as a prior to estimate individual posterior parameters or first have parameters reestimated for given data. By systematically deriving mechanism‐based PD models from QSP models, we bridge the gap between mechanistic and empirical modeling and make a step toward exploiting QSP models to guide dose adaptation within model‐informed precision dosing.   U.F., J.K., C.K., and W.H. designed the research. U.F. performed the research. U.F., J.K., and W.H. analyzed the data. U.F., J.K., C.K., and W.H. wrote the manuscript.   Funding provided by the graduate research training program PharMetrX: Pharmacometrics & Computational Disease Modeling, Berlin/Potsdam, Germany. Funded by the Deutsche Forschungsgemeinschaft (German Research Foundation)–Projektnummer 491466077.   C.K. and W.H. report research grants from an industry consortium (AbbVie Deutschland GmbH & Co. K.G., AstraZeneca, Boehringer Ingelheim Pharma GmbH & Co. K.G., Grúnenthal GmbH, F. Hoffmann‐La Roche Ltd., Merck KGaA, and SANOFI) for the PharMetrX program. In addition, C.K. reports research grants from the Innovative Medicines Initiative‐Joint Undertaking (“DDMoRe”) and Diurnal Ltd. C.K. report grants from the Federal Ministry of Education and Research within the Joint Programming Initiative on Antimicrobial Resistance Initiative. J.K. is an employee of AstraZeneca and owns stock in AstraZeneca. U.F. declared no competing interests for this work.   Appendix S1. Click here for additional data file.   