Scale reduction of a systems coagulation model with an application to modeling pharmacokinetic-pharmacodynamic data.

Bridging systems biology and pharmacokinetics-pharmacodynamics has resulted in models that are highly complex and complicated. They usually contain large numbers of states and parameters and describe multiple input-output relationships. Based on any given data set relating to a specific input-output process, it is possible that some states of the system are either less important or have no influence at all. In this study, we explore a simplification of a systems pharmacology model of the coagulation network for use in describing the time course of fibrinogen recovery after a brown snake bite. The technique of proper lumping is used to simplify the 62-state systems model to a 5-state model that describes the brown snake venom-fibrinogen relationship while maintaining an appropriate mechanistic relationship. The simplified 5-state model explains the observed decline and recovery in fibrinogen concentrations well. The techniques used in this study can be applied to other multiscale models.

Venom-induced consumption coagulopathy (VICC) is one of the most important snake envenoming syndromes. It is characterized by activation of the coagulation pathway by procoagulant toxins in the venom, resulting in the consumption of clotting factors such as fibrinogen, factor V, VIII, and X.[1,2,3,4] A comprehensive systems pharmacology model of the coagulation network has been shown to describe the time course of changes in coagulation factors in response to various anticoagulants[5] and Australian elapid envenoming.[6] Brown snakes are the most common cause of serious snakebite in Australia,[7] and their venom contains a potent prothrombin activator. Brown snake envenomation results in VICC with depletion of fibrinogen, factor V and factor VIII.[1,6] The time course of changes in fibrinogen, from any cause, is poorly understood,[8,9,10,11,12] making this an ideal setting to investigate the hypofibrinogenaemic states.

An established mathematical model of the coagulation network[6] is described by 62 ordinary differential equations (ODEs) and 178 parameters. See Supplementary Table S1 for the numbered states of the model with their corresponding species and initial conditions. ODEs for the model are also provided in the Supplementary Data (original_model.m). The values of the parameters in the model can be obtained from our previous publications.[5,6,13] This model consists of multiple inputs, outputs, and modifiers of these input–output relationships. Based on any given set of available data on an input–output relationship, it is possible that while certain states and parameters in the model are relevant to the relationship, others are either less important or may have no influence at all. It is important to differentiate between parameters and mechanisms that would be informed by a new data set and those that would not. Eliminating the parameters and reactions that are not informed by the data would result in a simpler version of the model that is specific to that particular input–output relationship. In doing this, it is important to retain the mechanistic structure even for simpler models.

Model order reduction methods for reducing the number of states and parameters of dynamical systems that are defined by ODEs currently exist.[14] Lumping is a model order reduction method where the original states of the model are lumped or merged to a reduced number of pseudo-states resulting in a reduced number of equations and parameters but with effectively the same or similar input–output behavior. Proper lumping is a special case of lumping where each of the original states contributes to only one of the pseudo-states of the reduced system, thereby forming groups that retain a clear physical interpretation. Note there may be more than one state lumped into a pseudo-state, but no state will be lumped into more than one pseudo-states. The technique of proper lumping has been used to lump physiologically based pharmacokinetic models.[15,16,17,18] An application of proper lumping based on searching through the possible combinations of lumpable states has been described for NF-κB signaling pathways to demonstrate its use and performance.[19]

Even though the number of parameters in the blood coagulation model can be reduced using a model reduction method like proper lumping, it is important that the model is identifiable to get precise estimates of the parameters. Various methods for the assessment of identifiability of linear and nonlinear models have been proposed.[20,21] An informal, unified approach using an information theoretic framework has recently been developed for simultaneous assessment of structural and deterministic identifiability for fixed- and mixed-effects pharmacokinetic–pharmacodynamic (PKPD) models. The approach has been applied to Bateman and Dost[22] models and to a combined parent-metabolite model describing the pharmacokinetics of ivabradine and its N-desmethylated metabolite, S-18982.[23]

Aims

This work aimed to explore a simplification of a blood coagulation systems pharmacology model for use in modeling PKPD data. Specific aims included the following:

To create a simplified model for exploring fibrinogen recovery after envenomation that mechanistically aligns with the coagulation systems pharmacology model

To extract the simplified model for estimation purposes

To assess structural identifiability of the simplified model based on the inputs and outputs available in the data set

To develop a population PKPD model for fibrinogen concentration–time data based on the mechanisms apparent in the simplified model

Results

Simplification of the blood coagulation systems pharmacology model

The original 62-state blood coagulation systems pharmacology model was simplified to a 5-state model () that described the brown snake venom–fibrinogen input–output relationship. An in silico brown snakebite followed by an in silico antivenom administration at 4 h resulted in a similar consumption-recovery profile for fibrinogen using the 5-state and 62-state models. During this process, some parameter(s) were adjusted heuristically by modifying their values manually to refine the predictions of the different structure of the simplified model. The values of the parameter(s) were adjusted (manually) based on how well the predictions from the simplified model compared with the original model.

compares the model predictions from the original 62-state model with the observed fibrinogen concentrations post–brown snakebite; compares the model predictions from the lumped 5-state model with the observed fibrinogen concentrations post–brown snakebite; compares the model predictions from the 62-state and 5-state models; and shows the percent relative differences in the predictions from the original 62-state and lumped 5-state models. It was observed that out of the five lumped states formed, fibrinogen and the brown snake venom absorption and plasma states were the ones that were initially left unlumped; the fourth and fifth states being factor IIa and a single lumped state that was a group of the remaining 58 states. summarizes the five lumped states with their initial conditions and their corresponding states of the original model that were lumped. Lumping the states further worsened the fibrinogen profile significantly.

Extraction of the simplified model

ODEs for the unlumped states in the 5-state lumped model were explicitly written by eliminating the reactions from the original ODEs that were lumped based on the proper lumping technique. Because state 5 was a lumped state that was formed as a result of merging of various states from the original model, its ODE had to be explicitly written as if it had been an unlumped state. The clotting factor that was most relevant to the brown snake venom–fibrinogen relationship represented the fifth state. With this in mind, the fifth state was therefore assumed to be dominated by factor II. The extraction of the ODEs of the simplified model from the ODEs of the original model resulted in reduction of the total number of parameters to 11 compared to 178 in the original model. Because the ODEs were extracted from the ODEs of the original model, this version of the model was referred as the “extracted model.” The following are the ODEs of the extracted model with their respective initial conditions: Eq. 1, ODE for state 1—venom absorption state; Eq. 2, ODE for state 2—venom plasma state; Eq. 3, ODE for state 3—fibrinogen; Eq. 4, ODE for state 4—factor IIa; Eq. 5, ODE for state 5—lumped state.

Note that in these equations, the form of the model reflects a time-varying system rather than a nonlinear Michaelis–Menten process. Hence, the values of vmax and Km do not necessarily reflect their implicit meaning in standard enzyme kinetics (see Wajima et al.[5]).

Identifiability of the simplified model

Assessment of identifiability of the extracted model using Population OPTimal (POPT) design software found that fixing the two parameters, pLumped_state and dLumped_state, resulted in both the predefined conditions for identifiability being met. Thus, 9 parameters out of the total 11 parameters in the extracted model were identifiable. The remaining two unidentifiable parameters, pLumped_state and dLumped_state, arose from the lumped state and were required to be fixed at 147 nmol/l/h and 0.01/h , respectively. These values were assumed to be equal to the values of pII and dII in the original model.

Modeling the fibrinogen concentration–time data using the simplified model

The decline and eventual recovery of fibrinogen after brown snake envenomation was described by the 5-state model. Testing additive, exponential, and combined error models for residual unexplained variability showed that the combined error model described the data the best. Comparison of the models with between-subject variabilities on one structural parameter or more showed that the model with between-subject variabilities on dIIa, , , dFg, and Fvenom provided the best description of the data. The parameter estimates for the final model are shown in . Individual fits obtained from the final model are shown in Supplementary Figures S1–S3. It can be seen that the model was able to capture the trend in the data for all individuals. shows the corresponding visual predictive check using the parameter estimates from the final model. The visual predictive check showed that the model explained the observed data well. In comparison to the 5-state model, the original 62-state model was only used to predict into the full data set, and this was not a model-fitting process but rather an overlay of the model predictions over the data. There are, therefore, no individual fits from the original model that can be shown.

Discussion

This study explores the simplification of a complex and complicated systems pharmacology model and its application to modeling PKPD data. A comprehensive model of the blood coagulation network was simplified and used to model data from a large study of snake envenomations. The simplified version of the model was shown to describe the time course of changes in fibrinogen concentrations following brown snake envenoming, which was comparable to that with the original version of the model.

Use of proper lumping as a model order reduction technique in this study showed that a 62-state blood coagulation network model could be significantly reduced to a 5-state model that described the fibrinogen consumption and recovery after envenomation from a brown snake comparable to the original model. Brown snake venom–fibrinogen serves as an example of an input–output relationship to show the application of proper lumping to reduce the number of states in a systems pharmacology model. The technique, however, can be applied to any other input–output relationship.

The ODEs of the reduced system had to be written explicitly, and the resulting extracted model only describes the brown snake venom–fibrinogen relationship. Lumped models are specific to the input–output being studied and other lumped models will be required from the same systems model for different input–output relationships. Use of an Information Theoretic Approach to assess the simplified model resulted in the list of parameters that were identifiable and thus could be estimated precisely. The structural identifiability method used in this work was based on that proposed by Shivva et al.[23] This method is a local identifiability method based on an information theoretical framework. This could be termed as a “pragmatic” identifiability solution.

In this study, the half-life of fibrinogen was estimated to be 1.5 days. This compares to a value of 1 day from the original 62-state model (as per Wajima et al.).[5] Also, this value is similar to another study of taipan bites in which the half-life was found to be 1 day.[24] Others have found longer values of half-life between 4 and 5 days,[25] and it may need to be explored whether there are other components of snake venom that may eliminate fibrinogen via other mechanisms or there are other processes that get stimulated following snake envenoming. The half-life of brown snake venom was estimated as equal to 55 min. There are no previous studies of brown snake venom in humans that have investigated the half-life of the venom. However, in other snake species like red-bellied black snake, the half-life has been reported to be ~8–12 h (G.K. Isbister, personal communication). This was only able to be determined because the antivenom was not given in many cases of red-bellied black snakebites. In addition, the half-life refers to the presence of venom as a whole. The half-life of brown snake venom estimated in this study refers in particular to the activity of the prothrombin activator in the venom that is Xa:Va like. However, it is possible that the half-life of the functional activity of the procoagulant toxin is different to the half-life measured as the presence of whole venom in serum.

The use of lumping principles has been illustrated previously for systems models in the area of PKPD. These works have concentrated on physiologically based pharmacokinetic models, and only limited attention has been placed specifically on systems models describing pharmacological effects. The studies on physiologically based pharmacokinetic models characterized the pharmacokinetics of a homologous series of barbiturates in rat[17] and for compounds like 1,3-butadiene[15] that are largely used in the production of plastics and synthetic rubber. Proper lumping principles have also been applied to a systems biology model describing the signaling pathways of NF-κB.[19] These studies, however, only discuss the illustration of the lumping process. The work carried out in our study appears to be the first time that a systems pharmacology model has been simplified to model PKPD data.

It might be theoretically possible to perform an identifiability analysis on the original model and then fix unidentifiable parameters. However, this model (i) would be cumbersome to use and computationally intensive to model and (ii) would not be able to be described simply in any report or future work. The theoretical benefit of the methods presented here is the observation that a model can be “extracted” from the original model to accommodate potentially any input–output relationship that is of interest. In addition, the use of mechanistic elements in models provides for greater predictive performance to settings outside of which the model was created. Using the simplification techniques discussed in this study, any multiscale mechanistic model that includes the necessary input–output relationship of interest can be used as a starting point for model extraction and then building. The resulting simplified model preserves the mechanism as seen in the original systems model. This allows for rapid input–output model building that effectively eliminates lengthy structural model selection processes. However, it should be noted that model order reduction techniques are locally dependent on the parameter space. Care should be taken to ensure that the performance of the final reduced model meets the intended needs.

To conclude, a simplified version of a multiscale model was created for exploring fibrinogen recovery after snake envenomation. The simplified model mechanistically aligned with the coagulation systems pharmacology model and was extracted for estimation purposes. The simplification process is multistaged and could benefit from some level of automation. A population PKPD model for fibrinogen concentration–time data based on the mechanisms apparent in the simplified model was developed and was able to describe the recovery of fibrinogen following brown snake envenomation without the need for further structural model development.

Methods

The technique of proper lumping based on a previously published method by Dokoumetzidis and Aarons[19] was used to simplify the blood coagulation systems pharmacology model. The work of Dokoumetzidis and Aarons consisted of two components (i) a formal method for lumping and (ii) an automated algorithm for lumping a systems model. The formal method of lumping demonstrated the use of a lumping matrix to simplify a model with X numbers of states to a model with Y numbers of states (where Y < X). In our work, we used their formal method of using the lumping matrix, but instead of using their automated algorithm, we manually undertook the process to control more easily what states were lumped.

The formal method of using the lumping matrix to reduce the number of states of a model, as explained previously by Dokoumetzidis and Aarons,[19] is described briefly below. In our work, the resulting simplified model was termed the “lumped” model. The 62-state systems pharmacology model was referred to as the “original” model. Supplementary Table S1 summarizes the 62 states in the original model with their corresponding species and initial conditions. ODEs for the original model are also provided in the Supplementary Data (original_model.m). The values of the parameters in the model can be obtained from our previous publications.[5,6,13]

The 62-state model was written in terms of the concentrations of the species in their respective states given by vector

where y is the vector of states of length 62 and is time, an independent variable.

A lumping matrix M consisting of 0s and 1s and of dimension n×n (n < n) was used to transform the vector of states y of dimension n×1 to a vector of lumped states y of dimension n×1:

The number of rows in the lumping matrix M was equivalent to the number of final lumped states intended to be obtained as a result of proper lumping. The number of columns in M was equal to the number of states in the original model which was 62 in this study. As an example, if the original 62-state model was being lumped to a 35-state model, dimensions of y would be 35 × 1 and that of would be 62 × 1. This would give dimensions of M equal to 35 × 62, which would be used for the transformation of y to y.

Inverse transformation of Eq. 6 is given by:

Here, is the Moore–Penrose pseudoinverse of the lumping matrix M. Thus, given the 62-state coagulation model and the lumping matrix M, which transforms the vector of states y to a vector of lumped states y such that y = My, the aim was to determine f (y), where . The following steps were followed to determine f (y):

Thus, the main lumping formula that was used to simplify the 62-state coagulation model was (as per Li and Rabitz[26]):

For linear systems, the lumping matrix M can be used to transform the parameter vector of the original model (K) to the parameter vector of the lumped model (K). The lumping formula used gives the lumped system explicitly with the values of the parameters of the lumped system. The lumped model can then be used for estimation purposes. For time-varying linear or nonlinear ODE systems, the lumping formula does not provide adjusted parameter values for the lumped model. But it gives the initial conditions of the lumped states that are formed as a result of merging of various states. The lumping process uses the following relation to give the initial conditions of the lumped states:

Here, and are, respectively, the vectors of initial conditions for the original and lumped models. Thus, in the case of time-varying or nonlinear systems, the lumped model produces outputs based on the reduced number of states and modified initial conditions but with the parameters from the original model.

For the time-varying nonlinear coagulation model, the above process was used to simplify the coagulation model up until the point when the lumped model started to give divergent results to the original model using the parameter values of the original model. After that point, the lumping process was not able to simplify the coagulation model further using the same set of parameter values as the original model. Further simplification of the model was only possible heuristically by adjusting values of parameter(s) manually. The values of the parameter(s) were adjusted (manually) based on how well the predictions from the simplified model compared with the original model.

Based on the input–output relationship considered in this study, states that were left unlumped were fibrinogen and the brown snake venom absorption and plasma compartments. For the remaining 59 lumpable states, each state was lumped with another arbitrarily chosen state, and the lumping matrix M was constructed accordingly. This search for the combination of states was conducted in a stepwise manner looking through a single series of lumping steps. If we intended to lump our original 62-state model to just 52 states, there will be 1026 possible combinations. This total number of combinations was based on the binomial coefficient:

Ultimately, an automated global search algorithm is desirable that operates within a predefined set of constraints. But this is not yet available. As the number of possible combinations of states was extremely high in this study, all possible combinations were not tested in this work. The lumping method used here can, therefore, be considered as a local solution.

The initial conditions of the lumped states obtained from the lumping formula and the parameters of the original model were used to simulate the time course of fibrinogen post–brown snakebite. The simulated time courses were then compared among the lumped and the original models. To simulate the time course of fibrinogen in a complete VICC scenario post–brown snakebite, a virtual brown snake venom dose of 0.0075 nmol/l, as in a previous study,[6] was used. All complete VICC patients were given antivenom in the observed snakebite data set,[1] and its administration in the model was set at 4 h after the bite as the average found in the brown snake venom data. Simulations were carried out using MATLAB R2011a. The simulated fibrinogen concentrations from the lumped model were compared with those from the original model at each time point to assess for loss of predictive performance. The lumped model was considered as giving a similar output to the original unlumped version if the percentage relative difference was less than 20% for all time points. The value of 20% was chosen arbitrarily. For differences greater than 20%, the state was unlumped from the pseudo-state being considered at that point and other combinations assessed. This process was stopped when no further combination of states was possible that resulted in a similar time course of fibrinogen post-bite among original and lumped models.

Because the aim of this study was to reduce the number of states as well as parameters, it was important to extract the ODEs of the lumped model to model the fibrinogen concentration–time data. ODEs of the lumped model were extracted from the ODEs of the original model by eliminating the “unwanted” reactions that did not have any influence on the fibrinogen profile based on the proper lumping technique. The ODEs of the lumped states (or pseudo-states) that were formed as a result of merging of various states from the original model had to be explicitly written as if they had been unlumped states. The input–output model with the reduced number of states and parameters was referred to as the “extracted model.”

The structural identifiability of the extracted input–output model was assessed using an Information Theoretic Approach.[23] POPT design software was used to carry out the structural identifiability analysis. A venom dose of 0.0075 mg and a generic study design () for the sampling times (1, 2, 3, 4, 5, 6, 8, 10, 16, and 24 h) were assumed for the analysis. A predefined condition, as per Shivva et al.,[23] had to be met for a model to be structurally identifiable. According to the condition, for all values of design variables (), where number of time points of observation is greater than the number of parameters and all design variables have a unique value, the log of the determinant of the Fisher information matrix should approach infinity as the associated noise, the log of the residual variance, approaches zero. If the predefined condition was not met, different sets of parameters, individually and in combination with each other, were fixed and the model reassessed.

Modeling the fibrinogen concentration-time data using the simplified model

The data. The data from the Australian Snakebite Project consisted of timed clotting factor concentrations in a large cohort of patients with VICC resulting from Australian elapid envenomation. There were 138 snakebite cases recruited to the Australian Snakebite Project from over 100 hospitals in Australia between January 2004 and May 2008. Each patient contributed only one case. summarizes the patient demographics. The data used for this analysis consisted of timed fibrinogen concentrations in patients with complete VICC resulting from brown snake envenomation. All patients with complete VICC had been administered antivenom. The data for many patients were sparse, and uncertainty in the venom “dose” was a key feature of the data. The limit of quantification for fibrinogen was 0.2 g/l. The data are shown in .

Data analysis using a population approach. A full population approach was carried out to analyze the fibrinogen data using NONMEM version 7.2 (using first-order conditional estimation with interaction).[27] Beal's method M3[28] was used for handling data that were below the limit of quantification. The extracted structurally identifiable simplified coagulation model was used as the structural model to describe the fibrinogen concentration-time data. The compartmental structure of the model was not subject to further modeling. The dose of brown snake venom was set to 0.0075 mg for all patients. Additive, exponential, and combined error models were tested for the residual unexplained variability. Residuals were assumed to be normally distributed with a mean of zero and variance σ2. Between-subject variability on the structural parameters was modeled assuming an exponential structure and was included on the bioavailable fraction, where the fixed-effects value was fixed to 1.

Any parameters determined to be unidentifiable were fixed at the values in the original 62-state model and thus not estimated. Model comparison was based on the likelihood ratio test with a decrease in the objective function value of at least 3.84 for one parameter difference between nested models to determine statistical significance. Additional model selection criteria included reduction of random residual unexplained variability and the stability of the model. A percentile visual predictive check to evaluate the final model was performed by simulating 1,000 replicates from the model and comparing the 10th, 50th, and 90th percentiles of the observed data and the prediction intervals derived from the simulated data graphically. Censoring of the data was used with the simulated data for the visual predictive check in the same way as the real data.

Author Contributions

A.G. wrote the manuscript; G.K.I. and S.B.D. designed the research; A.G. performed the research; and A.G. analyzed the data.

Conflict of Interest

The authors declared no conflict of interest. As an Associate Editor for CPT:PSP, S.B.D. was not involved in the review or decision process for this paper.

Study Highlights

Table 1

Numbered states of the lumped 5-state coagulation model with their initial conditions and the corresponding states of the original 62-state model

Table 2

Parameter estimates for the final model for fibrinogen

Table 3

Patient demographics