Dynamic Modeling of the Human Coagulation Cascade Using Reduced Order Effective Kinetic Models

In this study, we present a novel modeling approach which combines ordinary differential equation (ODE) modeling with logical rules to simulate an archetype biochemical network, the human coagulation cascade. The model consisted of five differential equations augmented with several logical rules describing regulatory connections between model components, and unmodeled interactions in the network. This formulation was more than an order of magnitude smaller than current coagulation models, because many of the mechanistic details of coagulation were encoded as logical rules. We estimated an ensemble of likely model parameters (N = 20) from in vitro extrinsic coagulation data sets, with and without inhibitors, by minimizing the residual between model simulations and experimental measurements using particle swarm optimization (PSO). Each parameter set in our ensemble corresponded to a unique particle in the PSO. We then validated the model ensemble using thrombin data sets that were not used during training. The ensemble predicted thrombin trajectories for conditions not used for model training, including thrombin generation for normal and hemophilic coagulation in the presence of platelets (a significant unmodeled component). We then used flux analysis to understand how the network operated in a variety of conditions, and global sensitivity analysis to identify which parameters controlled the performance of the network. Taken together, the hybrid approach produced a surprisingly predictive model given its small size, suggesting the proposed framework could also be used to dynamically model other biochemical networks, including intracellular metabolic networks, gene expression programs or potentially even cell free metabolic systems. Processes 2015, 3 179


Introduction
Developing mathematical models of biochemical networks is a significant facet of systems biology.Modeling approaches differ in their degree of detail, where the choice of approach is often determined by prior knowledge, or model requirements [1].Ordinary differential equation (ODE) models are common tools for modeling biochemical systems because of their ability to capture dynamics and encode mechanism.However, ODE models typically come with difficult (or sometimes impossible) parameter identification problems.For example, Gadkar et al. showed that even with near-perfect information, it was often impossible to identify all the parameters in typical signal transduction models [2].However, it is not clear whether we actually need precise estimates for all model parameters.Bailey suggested more than a decade ago, that achieving qualitative or even quantitative understanding of biological systems should not require complete structural and parametric knowledge [3].Since Bailey's complex biology with no parameters hypothesis, Sethna showed that model performance is typically sensitive to only a few parameters, a characteristic seemingly universal to multi-parameter models referred to as sloppiness [4].Thus, reasonable predictions may be possible, despite parametric uncertainty, if a few critical parameters are well-defined.For example, Tasseff et al., showed in a model of Retinoic acid (RA) induced differentiation of HL-60 cells, that correct predictions were possible even when 75% of the parameters were known only to an order of magnitude [5].Perhaps more importantly, ODE models require significant mechanistic information, thereby limiting their utility in poorly understood systems, or conversely explode in size when considering multiple pathways or subsystems.Toward this challenge, logical modeling is an emerging paradigm that encodes causal relationships between model components using quasi-mechanistic non-linear transfer functions [6].Logical models are highly flexible, and despite their simplicity, they have captured rich behaviors in a variety of systems important to human health [7][8][9].However, modeling complex dynamics with logical models is challenging.Thus, there is an unmet need for a third approach which combines ODEs and logical models, where ODEs could encode mechanistic information, while missing or incomplete mechanistic knowledge can be approximated using a logical approach.
In this study, we developed a hybrid approach which combined ODE modeling with logical rules to model a well studied biochemical network, the human coagulation system.Coagulation is an archetype proteolytic cascade involving both positive and negative feedback [10][11][12].Coagulation is mediated by a family proteases in the circulation, called factors and a key group of blood cells, called platelets.The central process in coagulation is the conversion of prothrombin (fII), an inactive coagulation factor, to the master protease thrombin (FIIa).Thrombin generation involves three phases, initiation, amplification and termination [13,14].Initiation requires a trigger event, for example vessel injury, which leads to the activation of factor VII (FVIIa).Two converging pathways, the extrinsic and intrinsic cascades, then process and amplify this initial coagulation signal.The extrinsic cascade is generally believed to be the main mechanism of thrombinogenesis in the blood [15][16][17].Initially, thrombin is produced upon cleavage of prothrombin by fluid phase activated factor X (FXa), which itself has been activated by Tissue Factor/activated factor VII (TF/FVIIa) [10].Picomolar amounts of thrombin then activate the cofactors factors V and VIII (fV and fVIII) and platelets, leading to the formation of the tenase and prothrombinase complexes on activated platelets.These complexes amplify the early coagulation signal by further activating FXa, and directly converting prothrombin to thrombin.There are several control points in the cascade that inhibit thrombin formation, and eventually terminate thrombin generation.Tissue Factor Pathway Inhibitor (TFPI) inhibits FXa formation catalyzed by TF/FVIIa, while antithrombin III (ATIII) neutralizes several of the proteases generated during coagulation, including thrombin.Thrombin itself also inadvertently plays a role in its own inhibition; thrombin, through interaction with thrombomodulin, protein C and endothelial cell protein C receptor (EPCR), converts protein C to activated protein C (APC) which attenuates the coagulation response by proteolytic cleavage of fV/FVa and fVIII/FVIIIa.Termination occurs after either prothrombin is consumed, or thrombin formation is neutralized by inhibitors such as APC or ATIII.
Previous coagulation models have typically been formulated as systems of nonlinear ordinary differential equations, using mass action or more complex kinetics, to describe the rates of biochemical conversions [18][19][20][21][22]. Mechanistic ODE coagulation models from our laboratory [23,24] were built upon the earlier studies of Jones and Mann [25], Hockin et al. [26], and later Butenas et al. [27] who developed and then subsequently refined highly mechanistic coagulation models.Other laboratories have also expanded upon Hockin et al. for example by exploring the intrinsic pathway, the role of stochastic fluctuations in coagulation [28], and the dynamics of thrombin mediated clot formation [29] and fibrinolysis [30].Other aspects of coagulation have also been modeled, such as platelet biochemistry [31], multi-scale models of clot formation [32,33], and transport inside clots [34].However, these previous studies were largely based upon extensive mechanistic knowledge.This is possible because blood, while enormously complex, can be systematically interrogated.Other systems, such as intracellular signaling networks, are much more difficult to experimentally interrogate.Towards this unmet need, we formulated a hybrid modeling approach which combines ODEs and logical rules to model biochemical processes for which a complete mechanistic understanding is missing.We tested this approach by modeling the human coagulation cascade.Others have also constructed reduced order human coagulation models.Recently, Papadopoulos and co-workers constructed a phenomenological mathematical model for thrombin generation [35].Using four ordinary differential equations and six parameters, they derived an expression describing the temporal evolution of thrombin generation in a variety of cases.The reduced order Papadopoulos model showed good agreement with experimental data, and underscored that model reduction is possible even for complex positive feedback systems like coagulation.However, the Papadopoulos model was focused on thrombin generation, and had a lesser emphasis on the influence of physiological inhibitors such as ATIII or the protein C pathway.In this study, we focused on building a reduced order coagulation model that included the physiological inhibitors using a hybrid strategy.The hybrid model consisted of only five differential equations augmented with several logical rules.Thus, the model was more than an order of magnitude smaller than comparable purely ODE models in the literature.We estimated the model parameters from in vitro extrinsic coagulation data sets, in the presence of ATIII, with and without the protein C pathway.We then compared the model predictions with thrombin data sets, for both normal and hemophilic coagulation, that were not used for model training.Once validated, we performed flux and sensitivity analysis on the model to estimate which parameters were critical to model performance in several conditions.The reduced order hybrid approach produced a surprisingly predictive coagulation model, suggesting this framework could potentially be used to model other biochemical networks important to human health.Once activated, thrombin catalyzes its own activation (amplification step), as well as its own inhibition via the conversion of protein C to activated protein C (APC).APC and tissue factor pathway inhibitor (TFPI) inhibit initiation and amplification, while antithrombin III (ATIII) directly inhibits thrombin.All inhibition steps and trigger-induced initiation were modeled using a rule-based approach.Likewise, the dependence of amplification on other coagulation factors was also modeled using a rule-based approach.The abundance of the highlighted species (in the dashed boxes) was governed by an ordinary differential equation.All other species were assumed to be constant.

Formulation of Reduced Order Coagulation Models
We developed a reduced order extrinsic coagulation model to test our hybrid modeling approach (Figure 1).The core of our model was based upon the earlier work of Ismagilov and coworkers [36][37][38][39], where we added initiation, factor dependence, and specific inhibition terms to the earlier simplified model.A trigger event initiates thrombin formation (FIIa) from prothrombin (fII) through a lumped initiation step.This step loosely represents the initial activation of thrombin by activated FXa.Once activated, thrombin catalyzes its own formation (amplification step), and inhibition via the conversion of protein C to activated protein C (APC).Antithrombin III (ATIII) inhibits amplification, while APC and tissue factor pathway inhibitor (TFPI) potentially inhibit both initiation and amplification.All initiation and inhibition processes, as well as the dependence of amplification upon other coagulation factors, was approximated using our rule-based approach (Figure 2).Individual regulatory contributions to the activity of pathway enzymes were integrated into control coefficients (v's) using an integration rule (min/max).These control coefficients then modified the rates of model processes at each time step.Hill-like transfer functions 0 ≤ f (Z) ≤ 1 quantified the contribution of components upon a target process.Components were either individual inhibitor or activator levels or some function of levels, e.g., the product of factor levels.In this study, Z corresponded to the abundance of individual inhibitors or activators, with the exception of the dependence of amplification upon specific coagulation factors (modeled as the product of factors).When a process was potentially sensitive to multiple inputs, logical integration rules were used to select which transfer functions influenced the process at any given time.In our proof of concept model, we used a winner takes all strategy; the maximum or minimum transfer function was selected at any given time step.However, other integration rules are certainly possible.Taken together, while the reduced order coagulation model encodes significant biological complexity, it is highly compact (consisting of only five differential equations).Thus, it will serve as an excellent proof of principle example to study the reduction of a highly complex human subsystem.

Identification of Model Parameters Using Particle Swarm Optimization
A critical challenge for any dynamic model is the estimation of kinetic parameters.We estimated kinetic and control parameters simultaneously from eight in vitro time-series coagulation data sets with and without the protein C pathway.The residual between model simulations and experimental measurements was minimized using particle swarm optimization (PSO).A population of particles (N = 20) was initialized with randomized kinetic and control parameters and allowed to search for parameter vectors that minimized the residual.However, not all parameters were varied simultaneously.We partitioned the parameter estimation problem into two subproblems based upon the biological organization of the training data; (i) estimation of parameters associated with thrombin formation in the absence of the protein C pathway and (ii) estimation of parameters associated with the protein C pathway.Only those parameters associated with each subproblem were varied during the optimization procedure for that subproblem, e.g., thrombin parameters were not varied during the protein C subproblem.The PSO procedure was run for 20 generations for each subproblem, where each generation was 1200 iterations.The best particle from each generation was used to generate the particle population for the next generation.We rotated the subproblems, starting with subproblem 1 in the first generation.
The experimental training data for parameter estimation was reproduced from the experiments of Butenas and co-workers [40].In these experiments thrombin generation was initiated by FVIIa-TF using mean plasma concentrations of coagulation proteins and inhibitors.To prepare FVIIa-TF, TF (0.5 nmol/L) was relipidated into 400 µmol/L of phospholipid vesicles (PCPS) by incubation in 20 mmol/L HEPES, 150 mmol/L NaCl, and 2 mmol/L CaCl 2 pH 7.4 (HBS/Ca 2+ ) for 30 min at 37 • C. The relipidated TF was incubated with 10 pmol/L factor VIIa for 20 min to allow the formation of FVIIa-TF.Factors V, VIII and thrombomodulin (Tm) (when protein C activation is required) were added to FVIIa-TF complex.Thrombin generation was then initiated by adding equal volumes of this mixture with a mixture containing prothrombin, factor IX and factor X, TFPI, AT-III and protein C (added when required), protein S (added when required) and factor XI (added when required).In the training data, 5 pmol/L FVIIa-TF was used along with 200 µmol/L of phospholipid vesicles (PCPS) to initiate thrombin generation.All other the coagulation factors and inhibitors i.e., factors X, IX, V, and VIII, prothrombin, TFPI and AT-III, protein C and protein S (when applicable) were at their mean plasma concentrations.Reduced order coagulation model parameters were estimated using particle swarm optimization (PSO) with the protein C pathway as a function of prothrombin.Only APC pathway parameters were allowed to vary in these simulations keeping the parameters estimated without protein C pathways constant.Solid lines denote the simulated mean value of the thrombin profile for N = 20 independent particles, points denote experimental data.The shaded region denotes the 99% confidence estimate of the mean simulated thrombin value (uncertainty in the model simulation).(A-C) depict training data and results for 150%, 100% and 50% of physiological prothrombin levels in the presence of the protein C pathway.The experimental training data was reproduced from the study of Butenas et al. [40,41].Thrombin generation in these experiments was initiated using 5 pmol/L FVIIa-TF in the presence of 200 µmol/L of phospholipid vesicles (PCPS).As depicted in (A-C) the prothrombin levels were at 150%, 100% and 50% of their physiological concentration in the presence of protein C pathway.All factors and control proteins in these experiments were at their physiological concentration unless otherwise denoted.
The reduced order coagulation model captured the role of initial prothrombin abundance, and the decay of the thrombin signal following from ATIII activity (Figure 3).However, we systematically under-predicted the thrombin peak and the strength of ATIII inhibition in this training data set.On the other hand, with fixed thrombin parameters, we captured peak thrombin values and the decay of the thrombin signal (at least for the 150% fII case) in the presence of both ATIII and the protein C pathway (Figure 4).Lastly, we were unable to capture global differences in initiation time across separate data sets with a single ensemble of model parameters.These differences likely resulted from normal experimental variability.For example, different thrombin generation experiments within the training data (at the same physiological factor levels) had significantly different initiation times (data not shown).However, the inability to globally capture initiation time also highlighted a potential shortcoming of the initiation module within the model.To capture the variability in initiation time across training data sets, we included a constant time-delay parameter (T D ) for each data group.The delay parameter was constant within a data set, but allowed to vary across training data sets.Introduction of the delay parameter allowed the model to simulate multiple training data sets using a single ensemble of model parameters.Taken together, the model identification results suggested that our hybrid approach could reproduce a panel of thrombin generation data sets in the neighborhood of physiological factor and inhibitor concentrations.However, it was unclear whether the reduced order model could predict new data, without updating the model parameters.

Validation of the Reduced Order Coagulation Model
We tested the predictive power of the reduced order coagulation model with validation data sets not used during model training.Two validation data sets were used, thrombin generation for various prothrombin and ATIII concentrations with the protein C pathway, and thrombin generation in normal versus hemophilic plasma in the presence of the protein C pathway.Lastly, we compared the qualitative output of the model to rFVIIa addition in the presence of hemophilia.The hemophilia case was an especially difficult test as it was taken from a different study which used a plasma-based in vitro assay involving platelets instead of phospholipid vesicles (PCPS).All kinetic and control parameters were fixed for the validation simulations.The only globally adjustable parameter T D , was fixed within each validation data set but allowed to vary between data sets.The reduced order model predicted the thrombin generation profile for ratios of prothrombin and ATIII in the absence of the protein C pathway (Figure 5).Simulations near the physiological range (fII,ATIII) = (100%, 100%) or (125%, 75%) tracked the measured thrombin values (Figure 5B,C).On the other hand, predictions for factor levels outside of the physiological range (fII,ATIII) = (50%, 150%) or (150%, 50%), while qualitatively consistent with measured thrombin values, did show significant deviation from the measurements (Figure 5A,D).Likewise, simulations of thrombin generation in normal versus hemophilia (missing both fVIII and fIX) were consistent with measured thrombin values (Figure 6).We modeled the dependence of thrombin amplification on factor levels using a product rule (Z = f V × f X × f V III × f IX), which was then was integrated using a min integration rule into the control variable governing amplification.Thus, in the absence of fVIII or fIX, the amplification control variable evaluated to zero, and the only thrombin produced was from initiation (Figure 6B).However, the decay of the thrombin signal was underpredicted in the normal case (Figure 6A), while the activated thrombin level was overpredicted in hemophilia simulations, although thrombin generation was far less than normal (Figure 6B).Taken together, the reduced order model performed well in the physiological range of factors, even with unmodeled components such as platelet activation in the hemophilia data set.(A-D) prediction results for (FII,ATIII): (50%, 150%), (100%, 100%), (125%, 75%) and (150%,50%) of physiological prothrombin and ATIII levels in the absence of the protein C pathway.The experimental validation data was reproduced from the study of Butenas et al. [40].Thrombin generation in these experiments was initiated using 5 pmol/L FVIIa-TF in the presence of 200 µmol/L of phospholipid vesicles (PCPS).As depicted in (A-D) the prothrombin and ATIII levels were at (50%, 150%), (100%, 100%), (125%, 75%) and (150%, 50%) of their physiological concentrations in the absence of protein C pathway.All factors and control proteins were at their physiological concentration unless otherswise denoted.The model ensemble predicted a direct correlation between thrombin generation and rFVIIa addition in hemophilia (Figure 7).In the current model, we cannot distinguish between different initiation sources, e.g., TF/FVIIa versus rFVIIa, as we have only a single lumped initiation source (trigger).Thus, we simulated the addition of rFVIIa in hemophilia by removing fVIII and fIX from the model, and modulating the initial level of trigger.
Simulations with a baseline level of trigger were consistent with the previous hemophilia simulations, where the only thrombin produced was from initiation (Figure 7A, 1× trigger).However, as we increased the trigger strength, the thrombin profile began to approximate normal coagulation, showing a pronounced peak albeit with a slower peak time (Figure 7B, t * * > t * ).Further increases in trigger strength resulted in decreased thrombin peak time and increased maximum thrombin values (Figure 7A, 50× trigger).Thus, for large trigger values (200× trigger), the hemophilic thrombin profile approximated normal coagulation, where peak thrombin was achieved shortly after administration and 95% of the thrombin was gone by 20 min after initiation.We performed flux analysis to understand how the reduced order coagulation model balanced initiation, amplification and termination of thrombin generation for normal and hemophilic coagulation.Analysis of the reaction flux through the reduced order network for thrombin generation in normal, hemophilia and rFVIIa-treated hemophilia identified three distinct operational modes (Figure 8).We calculated the flux through four lumped reactions, initiation, amplification, thrombin-induced APC generation and total thrombin inhibition (including both APC and ATIII action).Directly after the addition of a trigger (e.g., TF/FVIIa or rFVIIa), the lumped initiation flux was the largest for all three cases.However, within a few minutes enough thrombin was generated by the initiation mechanism to induce the amplification stage.During amplification, thrombin catalyzes its own formation and inhibition by generating activated protein C (APC), a potent inhibitor of the coagulation cascade.For normal coagulation, amplification and thrombin inhibition were the dominate reactions by 6 min after initiation (Figure 8, left).After 10 min, the dominate reaction had shifted to thrombin inhibition (both ATIII and APC action).In hemophilia (missing both fVIII and fIX), the amplification reaction did not occur, and thrombin was produced only by initiation (Figure 8, center).Initiation was quickly inhibited by APC, and the thrombin level stabilized (eventually decaying at longer times because of ATIII activity).Lastly, when 50× trigger was used to induce thrombin formation in hemophilia (absence of fVIII/fIX), initiation mechanisms dominated for up to 6 min following initiation (Figure 8, right).Similar to hemophilia alone, no amplification occurred in the 50× trigger+hemophilia case, and the rate of thrombin generation was extinguished by the combined action of ATIII and APC.Taken together, the hybrid modeling approach captured the transition between the modes of thrombin generation, as well as the role that inhibitors play in attenuating the thrombin generation rate.Thus, the transfer function approach encoded the inhibitory logic of this cascade in the absence of specific mechanism.

Global Sensitivity Analysis of the Reduced Order Coagulation Model
We conducted a global sensitivity analysis to estimate which parameters controlled the performance of the reduced order model.We calculated the sensitivity of the time to maximum thrombin (peak time) and the thrombin exposure (area under the thrombin curve) for different levels of prothrombin, and protein C (Figure 9).Globally, 41% of the parameters shifted in importance between the (fII,PC) = (50%, 0%) and (150%, 100%) cases for the peak thrombin time (Figure 9A).The majority of these shifts involved the interaction between increased prothrombin and the protein C pathway, while only 5% were directly associated with increased prothrombin alone.The rate constant for thrombin amplification was the most important parameter controlling the peak thrombin time.While this parameter was differentially important for different prothrombin levels, and in the presence or absence of the activated protein C pathway, it was consistently the most sensitive parameter in the model.The saturation constant governing thrombin amplification was the second most important parameter, followed by the initiation control gain parameter.Other important parameters influencing the thrombin peak time included the control gain for activated protein C formation, and the rate constant controlling ATIII inhibition of thrombin activity.On the other hand, only 27% of the model parameters were differentially sensitive between the (fII,PC) = (50%, 0%) and (150%, 100%) cases for thrombin exposure (Figure 9B).Of these parameters, all of the shifts were associated with the interplay between thrombin formation and the protein C pathway.The rate constant controlling ATIII inhibition was the most important parameter controlling the thrombin exposure.While this parameter was less important in the presence of protein C for 150% prothrombin levels, it was significantly above all other parameters.Similar to the peak time, for 150% prothrombin, the control gain for activated protein C formation was differentially important along with the rate constant controlling amplification.However, the amplification parameter was much less important for thrombin exposure vs. peak time.

Discussion
In this study, we developed a reduced order model of the human coagulation cascade.We modeled coagulation because it is well studied, has a complex architecture, and has an abundance of experimental data available for model identification and validation.However, coagulation was just a proof of concept test of our approach.The proposed hybrid framework could also be used to dynamically model other biochemical networks, including intracellular metabolic networks, gene expression programs or potentially even cell free metabolic systems.The model consisted of five differential equations augmented with several logical rules describing regulatory connections between model components and unmodeled interactions in the network.We estimated model parameters from in vitro extrinsic coagulation data sets, in the presence of ATIII, with and without the protein C pathway.To estimate parameters, the residual between model simulations and experimental measurements was minimized using particle swarm optimization (PSO).However, not all of the model parameters were uniquely identifiable, given the training data.Instead, we estimated an ensemble of likely parameter sets (N = 20) from eight in vitro time-series coagulation data sets with and without the protein C pathway.Ensemble approaches have been used previously for other signal transduction models [43][44][45][46][47], and for metabolic models [48] to estimate the impact of poorly constrained parameter values or poorly understood network structure on simulation performance.Thus, ensemble approaches are common in the dynamic modeling community.However, a unique feature of the current study is the direct connection between our particle swarm approach, and the parameter ensemble; each particle in our swarm uniquely corresponded to a parameter set in our ensemble.Thus, by constraining particles to operate in different parameter regions, giving each particle a different parameter combination to explore, or perhaps even suppling a different model formulation to each particle we can effectively traverse through complex parameter and model spaces.We validated the ensemble using thrombin data sets taken from multiple laboratories for a variety of experimental conditions not used during training.The ensemble predicted thrombin trajectories for conditions not used for model training, including thrombin generation for normal and hemophilic coagulation in the presence of platelets (a significant unmodeled component).We then used flux analysis to understand how the network operated in a variety of conditions, and global sensitivity analysis to identify which parameters controlled the performance of the network.Flux analysis showed the logical rules formulation encoded the transitions between initiation, amplification and termination of thrombin generation.Sensitivity analysis suggested that the amplification rate constant was more important to the time to peak thrombin, while the ATIII inhibition constant controlled thrombin exposure.Taken together, the proposed hybrid framework produced a surprisingly predictive model, suggesting this approach could be used to effectively model other biochemical networks important to human health.
Malfunctions in coagulation can have potentially fatal consequences.Aggressive clotting involved with Coronary Artery Diseases (CADs), collectively accounts for 38% of all deaths in North America [49].Coagulation management during surgery can also be challenging, particularly with the increase in clinical use of antithrombotic drugs [50].Insufficient coagulation due to genetic disorders such as hemophilia can also result in recurrent bleeding.The coagulation factors VIII (fVIII) and IX (fIX) are deficient in Hemophilia A and B, respectively [51][52][53].People with mild hemophilia have 5%-40% of the normal clotting factor levels while severe hemophiliacs have <1% [53].Hemophilia can be controlled with regular infusions of the deficient clotting factors.However, clotting factor replacement sometimes leads to the formation of fVIII and fIX inhibitors in vivo [54].Alternatively, recombinant factor VIIa (rFVIIa) has been used to treat bleeding disorders [55,56] including hemophilia with and without factor VIII/IX inhibitors [57].However, rFVIIa requires frequent administration (every 2-3 h), and many questions remain about its mechanism of action, its effective dosage [54], and its overall utility for the treatment trauma-associated hemorrhage [58].In this study, we did not model rFVIIa-induced coagulation directly.Rather, we modeled a general trigger which initiated the extrinsic coagulation cascade.Since we identified the model using TF/FVIIa, inherent to our rFVIIa simulations (and the rate constant governing initiation) was the presence of TF.However, even with this complication, the model generated potentially useful insight into the rFVIIa mechanism of action, and its possible shortcomings especially for the treatment of hemophilia.The addition of rFVIIa directly activated thrombin through the initiation pathway.However, no amplification of the thrombin signal occurred without fVIII or fIX.Thus, the peak thrombin signal was lower than normal coagulation, the peak thrombin time was longer, and thrombin generation was eventually inhibited by the combined action of ATIII and the protein C pathway.However, as the dose of rFVIIa increased, the peak thrombin time decreased (eventually saturating around 200×nominal trigger), and the peak thrombin value increased such that the thrombin profile resembled normal coagulation.Butenas et al. performed an extensive in vitro study of rFVIIa-induced thrombin generation under normal and hemophilic conditions [59].They found qualitatively similar trends, namely rFVIIa restored normal coagulation (even in the absence of TF) for large enough rFVIIa doses, although rFVIIa-induced coagulation in hemophilia (even for large rFVIIa doses) lagged the normal profile.These results suggest that rFVIIa administration alone might not be able to initiate normal coagulation in recurrent bleeding, unless the dosage is well above a critical threshold.However, defining this threshold, which is likely patient specific, is difficult as there is tremendous patient to patient variability even with a normal coagulation phenotype [60].
The hybrid model simulations of thrombin formation were qualitatively similar to other published coagulation models.Many mathematical models of coagulation dynamics have been built upon the Hockin-Mann model [26].For example, Brummel-Ziedins and co-workers incorporated the Protein C (PC) pathway into the original Hockin-Mann network to investigate thrombin generation in cases of familial PC deficiency [61].Simulations using this model showed that PC mutations caused elevated thrombin levels without changing the initiation time or the slope of thrombin generation.This trend was qualitatively captured by our reduced order model.For example, we showed decreased peak thrombin concentration in the presence of PC pathway and similar thrombin generation slopes, although the initiation time was different as these were different experimental trials.Danforth et al. simulated normal thrombin generation using 5 pM tissue factor with all other factors at their mean physiological level [60].The initiation time in this simulation was approximately 4.4 min.When predicting the normal thrombin generation curve using 0.2 nM of TF/FVIIa, we showed an initiation time between 4-5 min in a platelet based system, and an initiation time of 3-4 min in a PCPS system with 5 pM TF/FVIIa, although these times were largely dictated by the time delay parameter TD.Thus, while were not able to explicitly predict initiation time, we did bracket the initiation time predicted by the Danforth study.Mitrophanov et al. showed, in a study exploring the mode of action of rfVIIa [62], that increasing amounts of rFVIIa accelerated the maximal slope of thrombin generation, and both the peak thrombin and initiation times.While we failed to capture the effect of rFVIIa on initiation time, we correctly predicted that changes in the rFVIIa concentration affected the maximal thrombin slope and the propagation phase.Lastly, detailed mechanistic coagulation models, for example the model by Luan and co-workers [24] or the recent model by Mitrophanov et al. [30], often contain hundreds of proteins and interactions.Thus, it is unlikely that the reduced order hybrid model will replicate all the rich behavior of these other models.However, for qualitatively different cases such as normal versus hemophilia, the hybrid approach gave similar results.For example, akin to the hybrid model, Luan et al. modeled the normal and hemophilia data from Allen [42].However, unlike the previous Luan et al. study, we used this data for validation rather than training.Hybrid model simulations of the Allen et al. data set were surprisingly consistent with the Luan et al. model.For example, the amplification of a normal thrombin signal were comparable, and in hemophilia both correctly predicted decreased thrombin amplification.Taken together, the hybrid model performance was similar to other full scale mechanistic models in the literature, although we consistently failed to predict initiation times across data sets.
The performance of the proof of principle coagulation model was impressive given its limited size.However, there are several issues that could be further explored.First, the prediction of initiation time should be investigated.We were able to estimate initiation time within a data set, but unable to predict initiation time across independent data sets.This suggested that we should update the initiation module to distinguish between different triggers, e.g., TF/FVIIa versus rFVIIa alone, and to include key biological milestones such as FXa activation (a prerequisite to thrombin formation).Next, there are several additional biological modules that could be added to the core model presented here.First, we could include thrombin-induced platelet activation and the role of activated platelets in amplification.We captured thrombin generation data in the presence of platelets, however, the initial shape of the activation curve and the time-scale of activation was not always consistent with the data.Platelets are activated by thrombin through the cleavage of the extracellular domain of protease-activated receptors (PARs) on the platelet surface.Once activated, platelets play an important role in amplification, and are key mediators of the positive feedback driving amplification.Thus, this biology is a potentially important component of an expanded model.We should also add the intrinsic pathway to the model.The intrinsic pathway is triggered by contact activation of the plasma protease factor XI (fXI) by negatively charged surfaces and by thrombin and upstream factors such as activated plasma protease factor XII (FXIIa) [63,64].Activated platelets may also release polyphosphate which directly activates fXII [65].Arguably a minor player in acute bleeding, contact activation could also be important in other wound healing contexts.Finally, to make the model more clinically relevant, we should include the biochemical processes responsible for clot formation and clot dissolution (fibrinolysis).Clot formation is driven by thrombin activity, while fibrinolysis is driven by plasmin activity, a key enzyme that cleaves fibrin (one of the main materials in a clot).Similar to coagulation, fibrinolysis is managed by several activating and inhibitory factors which control the balance between clot formation and dissolution.Tissue plasminogen activator (t-PA) and urokinase activate plasmin, along with contact pathway factors such as fXIa.On the other hand, thrombin activatable fibrinolysis inhibitor (TAFI) inhibits the degradation of fibrin by plasmin.Also, similar to coagulation, there is considerable fibrinolysis and contact pathway data sets that can be used to train the model.Lastly, the choice of max/min integration rules or the particular form of the transfer functions could be generalized to include other rule types and functions.Theoretically, an integration rule is a function whose domain is a set of transfer function inputs, and whose range is v ∈ [0, 1].Thus, integration rules other than max/min could be used, such as the mean or the product, assuming the range of the transfer functions is always f ∈ [0, 1].Alternative integration rules such as the mean might have different properties which could influence model identification or performance.For example, a mean integration rule would be differentiable, which allows derivative-based optimization approaches to be used.The particular form of the transfer function could also be explored.We choose a Hill-like function because of its prominence in the systems and synthetic biology community.However, the only mathematical requirement for a transfer function is that it map a non-negative continuous or categorical variable into the range f ∈ [0, 1].Thus, many types of transfer functions are possible.

Formulation and Solution of the Model Equations
We used ordinary differential equations (ODEs) to model the time evolution of proteins (x i ) in our reduced order coagulation model: where R denotes the number of reactions, M denotes the number of protein species in the model.The quantity r j (x, , k) denotes the rate of reaction j.Typically, reaction j is a non-linear function of biochemical species abundance, as well as unknown kinetic parameters k (K × 1).The quantity σ ij denotes the stoichiometric coefficient for species i in reaction j.If σ ij > 0, species i is produced by reaction j.Conversely, if σ ij < 0, species i is consumed by reaction j, while σ ij = 0 indicates species i is not connected with reaction j.The system material balances were subject to the initial conditions x (t o ) = x o , which were specified by the experimental setup.Each reaction rate was written as the product of two terms, a kinetic term (r j ) and a control term (v j ) that could depend upon many regulatory transfer functions: We used multiple saturation kinetics to model the reaction term rj : where k max j denotes the maximum rate for reaction j, i denotes the enzyme abundance which catalyzes reaction j, and K js denotes the saturation constant for species s in reaction j.The product in Equation ( 3) was carried out over the set of reactants for reaction j (denoted as m − j ).The control term v j depended upon the combination of factors which influenced the activity of enzyme i.For each enzyme, we used a rule-based approach to select from competing control factors (Figure 2).If an enzyme was activated by m metabolites, we modeled this activation as: where 0 ≤ f ij (Z) ≤ 1 was a regulatory transfer function that calculated the influence of metabolite i on the activity of enzyme j.Conversely, if enzyme activity was inhibited by m metabolites, we modeled this inhibition as: Lastly, if an enzyme had both m activating and n inhibitory factors, we modeled the control term as: where: The quantities j + and j − denoted the sets of activating and inhibitory factors for enzyme j.If a process has no modifying factors, we set v j = 1.There are many possible functional forms for 0 ≤ f ij (Z) ≤ 1.However, in this study, each individual regulatory transfer function took the form: where Z j denotes the abundance of the j factor (e.g., metabolite abundance), and k ij and η are control parameters.k ij was the species gain parameter, while η was a cooperativity parameter (similar to a Hill coefficient).Applying the general framework to the reduced coagulation network resulted in five ordinary differential equations: where x = (f II, F IIa, P C, AP C, AT III) T .The terms r * v * in the balance equations denote corrected kinetic expressions for initiation, amplification and inhibition processes.The rate of initiation rinit was modeled as: where k init , K init,f II are the rate and saturation constants governing initiation, respectively.The rate of initiation was modified by v init , the control parameter governing initiation.Initiation was sensitive to the level of trigger (activator) and TFPI (inhibitor): where the transfer functions f took the form of Equation (9).The rate of thrombin amplification was given by: ramp = k amp (x 2 ) x 1 K amp,f II + x 1 (17) where k amp , K amp,f II denote the rate and saturation constants governing amplification, respectively.The amplification control term, which modified amplification rate, was modeled as a combination of multiple inhibition terms and one activation term: where ) is an activating term, we included it in the min integration rule; the factors in Z amp were essential for amplification (if any of these factors was missing the amplification reaction would not occur).Thus, the factors in Z amp were required components, a classification that we implemented by the min selection rule.The rate activated protein C formation was given by: rapc = k AP C,f ormation (T M ) x 3 K f ormation,P C + x 3 (19) where k AP C,f ormation and K f ormation,P C denote the rate and saturation constants governing activated protein C formation, respectively and T M denotes the thrombomodulin abundance.We modeled the control term which governed APC formation as a single thrombin-dependent activation term: Lastly, we included direct irreversible inhibition of FIIa by ATIII: rinh,ATIII = k AT III,inhibition (x 5 x γ 2 ) ( where γ was estimated to be γ = 1.26.For ATIII inhibition of FIIa, the control variables v inh,AT III was taken to be unity.The model equations were encoded using the Python programming language and solved using the ODEINT routine of the SciPy module [66].The model files can be downloaded from http://www.varnerlab.org.

Estimation of Model Parameters From Experimental Data
Model parameters were estimated by minimizing the difference between simulations and experimental thrombin measurements (squared residual): where xj (τ ) denotes the measured value of species j at time τ , x j (τ, k) denotes the simulated value for species j at time τ , and ω j (τ ) denotes the experimental measurement variance for species j at time τ .The outer summation is with respect to time, while the inner summation is with respect to state.We minimized the model residual using Particle swarm optimization (PSO) [67].PSO uses a swarming metaheuristic to explore parameter spaces.A strength of PSO is its ability to find the global minimum, even in the presence of potentially many local minima, by communicating the local error landscape experienced by each particle collectively to the swarm.Thus, PSO acts both as a local and a global search algorithm.For each iteration, particles in the swarm compute their local error by evaluating the model equations using their specific parameter vector realization.From each of these local points, a globally best error is identified.Both the local and global error are then used to update the parameter estimates of each particle using the rules: where (θ 1 , θ 2 , θ 3 ) are adjustable parameters, L i denotes the local best solution found by particle i, and G denotes the best solution found over the entire population of particles.The quantities r 1 and r 2 denote uniform random vectors with the same dimension as the number of unknown model parameters (K ×1).In thus study, we used (θ 1 , θ 2 , θ 3 ) = (1.0,0.05564, 0.02886).The quality of parameter estimates was measured using goodness of fit (model residual).The particle swarm optimization routine was implemented in the Python programming language.All plots were made using the Matplotlib module of Python [68].

Global Sensitivity Analysis of Model Performance
We conducted a global sensitivity analysis, using the variance-based method of Sobol, to estimate which parameters controlled the performance of the reduced order model [69].We computed the total sensitivity index of each parameter relative to two performance objectives, the peak thrombin time and the area under the thrombin curve (thrombin exposure).We established the sampling bounds for each parameter from the minimum and maximum value of that parameter in the parameter set ensemble.We used the sampling method of Saltelli et al. [70] to compute a family of N (2d + 2) parameter sets which obeyed our parameter ranges, where N was the number of trials, and d was the number of parameters in the model.In our case, N = 10,000 and d = 22, so the total sensitivity indices were computed from 460,000 model evaluations.The variance-based sensitivity analysis was conducted using the SALib module encoded in the Python programming language [71].

Figure 1 .
Figure1.Schematic of the connectivity of the reduced order coagulation model.A trigger compound, e.g., TF/FVIIa initiates thrombin production (FIIa) from prothrombin (fII).Once activated, thrombin catalyzes its own activation (amplification step), as well as its own inhibition via the conversion of protein C to activated protein C (APC).APC and tissue factor pathway inhibitor (TFPI) inhibit initiation and amplification, while antithrombin III (ATIII) directly inhibits thrombin.All inhibition steps and trigger-induced initiation were modeled using a rule-based approach.Likewise, the dependence of amplification on other coagulation factors was also modeled using a rule-based approach.The abundance of the highlighted species (in the dashed boxes) was governed by an ordinary differential equation.All other species were assumed to be constant.

Figure 2 .Figure 3 .Figure 4 .
Figure 2. Schematic of rule-based effective control laws.Traditional enzyme kinetic expressions, e.g., Michaelis-Menten or multiple saturation kinetics are multiplied by an enzyme activity control variable 0 ≤ v j ≤ 1.Control variables are functions of many possible regulatory factors encoded by arbitrary transfer functions of the form 0≤ f j (Z) ≤ 1.At each simulation time step, the v j variables are calculated by evaluating integration rules such as the max or min of the set of transfer functions f 1 , . . ., f n influencing the activity of enzyme E j .

Figure 5 .
Figure 5. Reduced order coagulation model predictions versus experimental data for normal coagulation.The reduced order coagulation model parameter estimates were tested against data not used during model training.Simulations of different levels of prothrombin and ATIII were compared with experimental data in the absence of the protein C pathway.Solid lines denote the simulated mean value of the thrombin profile for N = 20 independent particles, points denote experimental data.The shaded region denotes the 99% confidence estimate of the mean simulated thrombin value (uncertainty in the model simulation).(A-D) prediction results for (FII,ATIII): (50%, 150%), (100%, 100%), (125%, 75%) and (150%,50%) of physiological prothrombin and ATIII levels in the absence of the protein C pathway.The experimental validation data was reproduced from the study of Butenas et al.[40].Thrombin generation in these experiments was initiated using 5 pmol/L FVIIa-TF in the presence of 200 µmol/L of phospholipid vesicles (PCPS).As depicted in (A-D) the prothrombin and ATIII levels were at (50%, 150%), (100%, 100%), (125%, 75%) and (150%, 50%) of their physiological concentrations in the absence of protein C pathway.All factors and control proteins were at their physiological concentration unless otherswise denoted.

Figure 6 .Figure 7 .
Figure 6.Reduced order coagulation model predictions versus experimental data with and without coagulation factors VIII (fVIII) and IX (FIX).The reduced order coagulation model parameter estimates were tested against data not used during model training.Simulations of normal thrombin formation with ATIII and the protein C pathway were compared with thrombin formation in the absence of fVIII and fIX.Solid lines denote the simulated mean value of the thrombin profile for N = 20 independent particles, points denote experimental data.The shaded region denotes the 99% confidence estimate of the mean simulated thrombin value (uncertainty in the model simulation).(A,B) prediction results for normal thrombin generation and thrombin generation in hemophilia.All factors and control proteins were at their physiological concentration unless others noted.Coagulation was initiated with 0.2 nmol/L FVIIa.The experimental validation data was reproduced from the study of Allen et al.[42].

Figure 8 .
Figure 8. Reaction flux distribution as a function of time for thrombin generation under normal (left), hemophilia (center) and rFVIIa treated hemophilia (right).Reaction flux was calculated for each particle at T = 0, 4, 6, 8, 10, 12, 14 min after the initiation of coagulation.Reaction fluxes were calculated for each particle in the parameter ensemble (N = 20).Blue colors denote low flux values while red colors denote high flux values.
Points denote the mean total sensitivity value, while the area around each point denotes the uncertainty in the sensitivity value.The gray dashed line denotes the 45 • degree diagonal, if sensitivity values are equal for different conditions they will lie on the diagonal.Sensitivity values significantly above or below the diagonal indicate differentially important model parameters.The radius of the shaded region around each total sensitivity value was the maximum uncertainty in that value estimated by the Sobol method.