A Novel Simplified Model for Blood Coagulation: A piecewise dynamical model for thrombin with robust predictive capabilities

TL;DR

This paper introduces a simplified piecewise dynamical model for thrombin activity in blood coagulation, capable of predicting prolonged enzyme activity in specific plasma conditions, aiding patient-specific clotting analysis.

Contribution

It presents a novel, simplified switching model for thrombin dynamics that captures prolonged activity, advancing the understanding of clotting behavior in patient-specific scenarios.

Findings

Model accurately predicts prolonged thrombin activity.

Simplified variables effectively describe complex clotting dynamics.

Potential for studying patient-specific clotting in realistic geometries.

Abstract

Realistic description of patient-specific mechanical properties of clotting dynamics presents a major challenge. Available patient-specific data falls short of robustly characterizing myriads of complex dynamic interactions that happen during clotting. We propose a simplified switching model for a key part of the coagulation cascade that describes dynamics of just four variables. The model correctly predicts prolonged activity of thrombin, an important enzyme in the clotting process, in certain plasma factor compositions. The activity sustains beyond the time which is conventionally considered to be the end of clotting. This observation along with the simplified model is hypothesized as a necessary step towards effectively studying patient-specific properties of clotting dynamics in realistic geometries.

Figures

Tables

Table 1. Table 1: Switching Condition

	[IIa] 2 nM	[IIa] $\ge$ 2 nM
$K_S$	$k_s>0$	0
$K_I$	$k_{i2}>0$	$k_{i1}>0$
$K_P$	0	$k_p>0$

Equations

$$\frac{\partial [Y_i]}{\partial t}+div([Y_i]\mathbf{v})=div(D_{i}grad[Y_i])+G_i,i=1,...,N,$$

$$\frac{d[Y_i]}{dt}=G_i$$

$$s.t.[Y_i](0)$$

$$\mathbf{S}[\mathbf{Y}]$$

$$[Y_i]$$

$$r_i=\frac{k_i{x}_i^{{\eta }_i}}{1+k_i{x}_i^{{\eta }_i}}\text{min}\left(\frac{k_j{x}_j^{{\eta }_j}}{1+k_j{x}_j^{{\eta }_j}},\frac{k_m{x}_m^{{\eta }_m}}{1+k_m{x}_m^{{\eta }_m}}\right)$$

$$\frac{d}{dt}[II]=-K_S-K_P[II][IIa]$$

$$\frac{d}{dt}[IIa]=K_S+K_P[II][IIa]-K_I[IIa][AT]$$

$$\frac{d}{dt}[AT]=-K_I[IIa][AT]$$

$$\frac{d}{dt}[IIa-AT]=K_I[IIa][AT].$$

$$[II]+[IIa]+[IIa-AT]$$

$$[AT]+[IIa-AT]$$

$$u=\frac{1}{M}\sum _{m=1}^M\sum _{i=1}^4{\left(\frac{C_i^{reduced}(m;k)-C_i^{full}(m)}{C_i^{constant}}\right)}^2$$

Full text

1 Introduction

Patient-specific geometry modeling, simulations with more realistic boundary conditions, multiscale models that combine molecular mechanisms with clinical manifestation are some of the open problems discussed in vascular biomechanics [1]. Blood constituents vary drastically in patients and their composition plays a significant role in determining the mechanical properties of the resulting clot [2, 3, 4].

Current methods model generation and depletion of various constituents in blood flow using convection reaction diffusion equations [5],

[TABLE]

where $[Y_{i}]$ is a constituent, $\mathbf{v}$ is velocity, $D_{i}$ are diffusion coefficients, and $G_{i}$ are source terms due to clotting reaction dynamics. Number of constituents $N$ is in the order of thousands; typical models consider between 30 to 40 constituents [6, 5]; and clotting reaction models considering hundreds of constituents have been proposed [7, 8].

1.1 Need for Simple Models

Models accounting for myriads of chemical reactions are considered as a first step towards understanding the complex phenomena of clotting [5]. Further modeling many such constituents are useful in a diagnostic setting where often diseases are reflected as abnormalities in the coagulation cascade [9]. In particular reaction dynamics alone, for example, dynamics of thrombin [10, 11], is known to be associated with many diseases such as acute coronary syndrome [12], acute cerebrovascular disease [13], rheumatoid arthritis [14], etc. Recent studies in this direction include use of machine learning tools to characterize and use such abnormalities for disease classification tasks such as identifying patients prone to heart attacks [15, 16].

However, simplified models for the coagulation cascade are required so that:

the model predictions are more readily verifiable using experiments 2. 2.

augmenting chemical reaction dynamics in other phenomena such as flow and mechanical characterization of clots in patient-specific terms become feasible and readily comprehensible.

1.2 Challenges in Modeling the Coagulation Cascade

The chemical reaction kinetics particularly poses many challenges. Models describing kinetics using $M$ reactions are of the form,

[TABLE]

where $S_{ij}$ are stoichiometric coefficients ($i=1,...,N,j=1,...,M$), $r_{j}$ are reaction rates, $\mathbf{k}_{j}$ are parameters of the reaction rates, and $[Y_{i}^{0}]$ are the initial plasma composition.

Concentration of proteins involved in coagulation and rates of reactions vary by orders of magnitude (see Figure 1). Essentially, picomoles of trigger results in the formation of hundreds of nanomoles of certain enzymes. This in turn results in macroscopic formation of clots. Typically stochastic methods [17] are used to properly account for low concentrations of species [18, 19].

The nonlinear chemical kinetics is modeled using reaction rates that typically has quadratic terms or Hill-like terms [20] . The models are very stiff and solution trajectories are unstable in many directions [21]. The rates involve negative feedback loop or cycles in the reaction cascade [22]. This poses challenges in numerically coupling chemical reactions with flow simulations.

Moreover, there is uncertainty in the parameters of the reaction kinetics models. It is hard to measure most of these protein factors. Many rate parameters are inferred indirectly rather than being directly measured. Models used to simulate clot in small two dimensional regions ($\sim$100 $\mu$m) consider the dynamics of many reactants [23]. These are inappropriate for simulations in realistic 3 dimensional flow conditions in arteries, say, in order to study atheroscelerosis or thrombosis [23]. Reduced dimensional simplified chemical kinetics models will help to advance the quality of patient specific simulations in realistic geometries.

1.3 Relevant Literature on Reduced Order Models

Papadopoulos et al. [23] suggested a phenomenological model for thrombin dynamics. Based on the mechanism of thrombin dynamics, they propose a simplified thrombin dynamics model using four reactions. The reactions include dynamics of thrombin, prothrombin, platelets, and activated platelets. Using the assumption of fast platelet activation, they derive analytical expression for thrombin generation. These are similar to thrombin generation functions prescribed by Hemker et al. [24]. The model essentially fits patient specific thrombin generation profiles and the effect of the plasma factor composition and inhibitors on the dynamics of thrombin were not emphasized.

This motivated Sagar et al. [25] to come up with a dynamical model for thrombin generation using a hybrid strategy. The strategy combines differential equations and several logical rules to model thrombin dynamics. They design their approach to model systems where mechanistic insights are poor and experimental interrogation is difficult. This results in reduced order model that has rates of the product of Hill-like terms and transfer functions that act as the logical rules, i.e.,

[TABLE]

where $r$ is the reaction rate, $k$ and $\eta$ are parameters, and $x$ pertains to protein concentration or activity. Though the model shows good performance for thrombin dynamics, the transparency in the mechanistic models such as [24] and [23] is lost, i.e., the functionality of the model parameters and effect of their changes is not evident. We seek a middle ground between the two simplified models where we find a dynamical model that makes use of the mechanistic knowledge of blood coagulation and is able to account for changes in plasma factor composition.

1.4 Proposed Model Properties

We suggest a simple phenomenological model for thrombin dynamics based on chemical kinetics:

We model the stoichiometry and dynamics of certain important and easily measurable chemical species. The model is based on the classically viewed initiation, propagation, and termination of thrombin dynamics. Hence the chemistry involved in the simulations offer physiological insight. 2. 2.

We model the initiation phase and the propagation/termination phase of thrombin dynamics separately using a switching criteria for the rate parameters. The switch separates the slow initiation phase from the latter which is orders of magnitude faster. 3. 3.

The model we propose shows varied responses. The functionality of the parameters are evident and different aspects of thrombin dynamics are easily alterable. 4. 4.

A good model should be able to capture the necessary rich behavior of the phenomena as well as generalize well in order to predict important qualitative and quantitative responses. The model accounts for the physiological effect of antithrombin and is able to predict certain important changes in thrombin dynamics due to changes in prothrombin and antithrombin concentration. 5. 5.

The model is dynamical in the sense that rates for the model species could be calculated given the current values of the model variables. In other words, we essentially model $G_{i}$ for the important species in equation 1 so that the model could be extended to account for transport.

The simplified model proposed here is expected to make patient-specific mechanical characterization of clots in realistic geometries feasible.

2 Thrombin Generation in the Extrinsic Pathway

2.1 Full Model

Figure 2 shows a schematic of the key elements of the extrinsic pathway involved during clotting. We consider the extrinsic pathway because hemostasis occurs due to tissue factor initiation. Further, we simply simplify thrombin dynamics. Given that flow properties affect and are affected by fibrin formation, such a simple thrombin dynamics model factors out the two phenomena. Reactions in the intrinsic pathway could be readily accounted for due to the phenomenological aspect of initiation in the proposed model. For simulations of the full model, we use the extrinsic pathway developed by Hockin et al. [6]. A schematic of the model used is shown in Figure 3.

2.2 Key Ideas Used for Simplification

We exploit the following ideas for model simplification:

Dynamical Parametrization:

Patient specificity in blood coagulation is usually modeled using variations in the initial conditions for the differential equations. This leads to a high-dimensional state space and specification of dynamics in that space of many constituents as described earlier. However, we simplify the state of the system using just four important species and then model the dynamics of those four constituents. The effect of other constituents are parametrized in the rates describing the dynamics of the former, i.e., the effect of patient specificity of thrombin dynamics is described using variations in the rate parameters of the simplified model. We show that training the model for a specific choice of initial condition in the original state space (physiological mean) is able to predict qualitative responses of changes in prothrombin and antithrombin concentration in the reduced state space.

Switching:

We separate the parameters in the initiation and propagation/termination phase of thrombin dynamics using a switching model. Such switching models are commonly used to model complex nonlinear systems [26]. Switching models are known to result in simplified models for even chaotic dynamics. For example, piecewise affine models are known to model chaotic behavior of many nonlinear systems well [27]. The proposed model switches from initiation to propagation/termination phase based on the concentration of thrombin. Such threshold-based responses have been used in models describing initiation [23] and platelet activation [5] in the coagulation cascade. The hybrid (switching) model by Makin and Narayanan [28] intends to develop a comprehensive model rather than a simplified model.

2.3 Simplified Model

We use the traditional simplification of the thrombin generation cascade and describe kinetics for prothrombin, thrombin, antithrombin, and thrombin-antithrombin using the following set of reactions:

Thrombin Initiation: Tissue factor activates prothrombin to form thrombin.

{chemmath}

II \reactrarrow0pt2.5cm$K_{S}$Initiation IIa 2. 2.

Thrombin propagation: Given that sufficient amount of thrombin (2 nM) is activated, clotting propagates via a different set of reactions.

{chemmath}

II \reactrarrow0pt2.5cm$K_{P}$Propagation IIa 3. 3.

Thrombin inhibition: Finally normal hemostasis requires that thrombin generation is controlled.

{chemmath}

IIa + AT \reactrarrow0pt2.5cm$K_{I}$Termination IIaAT

Corresponding to the reactions, we define the following stoichiometric model:

[TABLE]

where we model initiation using the rate constant $K_{S}$, propagation using the rate constant $K_{P}$, and inhibition using the rate constant $K_{I}$. In order to be stoichiometrically consistent, our [IIa] is the sum of both forms of thrombin in the full 34 variable model (the variables used in our model are drawn in continuos lines in the full model simulation as seen in Figure 5). [IIa-AT] in our model is the sum of the two antithrombin complex formed due to inhibition111there are other antithrombin complexes formed in the full model but they are 3 orders of magnitude smaller than [IIa-AT] and we neglect them..

Among the four species modeled, methods to measure or infer thrombin, ATIII, and thrombin antithrombin [29, 30, 31, 32] exist. The rates in equation 4 are such that they satisfy a stoichiometric constraint, i.e.,

[TABLE]

$c_{1}$ and $c_{2}$ are constants defined using the initial conditions. These dependence relations are crucial while designing experiments to properly observe the system. For example, calibrated automated thrombogram assay is inadequate to make proper measurements with deficient ATIII [33] due to fluorogenic substrate depletion. In the results section, we will discuss the possibility of sustained thrombin activity when ATIII is deficient leading to the depletion.

We propose the switching rules in Table 1 that changes the response of the model during initiation and propagation/termination. Essentially, thrombin propagation occurs if [IIa] crosses a threshold. For normal clotting, rate of propagation is expected to be orders of magnitude higher than that of rate of initiation. We also use two different inhibition rate constant $k_{i1}$ and $k_{i2}$ so that rate of inhibition could be separately tuned during the two phases.

3 Results and Discussion

3.1 Estimation of Model Parameters

In the simulations of the full model, clotting was initiated with 5 pM and the plasma factor composition was set to physiological mean values [6]. We used the data from the full model to fit parameters for the reduced model. We note that thrombin dynamics of the full model has been corroborated with experimental data [6]. Particle swarm optimization [34] was used for parameter estimation and we obtained one set of parameters for the physiological mean composition.

We used the sum of squared differences of the normalized concentration profiles between the common species in the full model and the reduced model as the objective function,

[TABLE]

where $m$ denotes time points and $i$ denotes model species. Mean concentration of prothrombin and antithrombin were used as normalization constants $C_{i}^{constant}$. Comparison of the reduced and full model simulation for the physiological mean initial composition is shown in Figure 6.

3.2 Parameter Study of the Simplified Model

We show the effect of parameters by changing them one at a time. Clot time depends exponentially on $K_{S}$ (seen in Figure 7). Rate constant $k_{i2}$ also controls the clot time Figure 8. For certain combinations of $K_{S}$ and $k_{i2}$ it takes more than 1200 seconds for clot initiation. Both the parameters together offer more control over dynamics of clot initiation.

Figures 9 and 10 show the effect of changes in the rate constants $K_{P}$ and $k_{i1}$ respectively. The parameters offer a wide range of thrombin generation rates. Similar to the initiation, there are certain values of $K_{P}$ (for a given value of $k_{i1}$) and vice versa where thrombin generation is too low. These two parameters together offer control over simulating a wide range of thrombin propagation.

3.3 Prediction of Variation in Prothrombin and Antithrombin

Finally, we check the qualitative response of the model predictions for variations in initial prothrombin and antithrombin concentrations. As seen in Figure 11, higher values of prothrombin are able to predict more thrombin generation. This has been observed in experiments [35]. Thrombin rates during termination could be improved using better training data and using different reaction rates. Similarly, lower values of antithrombin are able to predict higher thrombin generation.

One of the most important predictions of this model is that thrombin termination appears to halt at non-zero values (Figure 12). Such sustained activity is also in observed in experiments when there is too much prothrombin compared to antithrombin [35]. Inhibitors such as activated protein C may need to modeled in order to account for oscillations observed in such sustained activity. In this model, the reaction essentially runs out of the inhibitor [AT] when initiated with a certain plasma factor composition. In such a scenario, other phenomena like diffusion and convection will control the extent of clotting. For example, when there is less inflow of antithrombin concentration, as in the case of stasis, we would expect more clotting due to the presence of excess active thrombin. The significance of antithrombin deficiency in simulations accounting for diffusion has already been reported [36]. For a given value of initial prothrombin, and given that thrombin generation proceeds normally due to sufficient rates, the simplified model quantifies the amount of antithrombin required for thrombin generation to terminate.

Further, active thrombin could be propagated downstream and could potentially cause clotting elsewhere. This observation and model prediction on sustained activity of thrombin are hypothesized to play a necessary role towards effectively studying clotting in realistic geometries.

3.4 Limitations and Improvements

The simplified model proposed here, when coupled with fluid flow and transport models, is expected to make patient-specific mechanical characterization of clots in realistic geometries feasible. The parameters of the model need to be altered in order to better simulate thrombin dynamics on a wider range of plasma factor compositions. Figure 15 shows simulations from parameters identified for acute coronary syndrome population mean composition [12] and a hemophilia patient [11] (plasma factor composition of patient C in Figure 3. (B)). The hemophilia patient simulations suggest the form of rate functions as well as the threshold value (2 nM [IIa]) could be modeled better.

For the mean plasma composition, Figure 15 shows the full model thrombin rates for the simplified model at different times during clotting plotted against the corresponding time points of the full model. The reaction rates could be improved to better align the simulation data with the $x=y$ line. We also note that the objective used for simplification in equation 6 does not explicitly penalize deviations in rates. Figure 15 shows thrombin rates for the simplified model at different points in the state space along a particular trajectory of the full model plotted against the corresponding rates of the full model. The alignment is better in this case compared to figure 15. This suggests that the simplified model takes a slightly different path compared to the full model. We are currently working towards modifying the dynamics and the objective so that the rates could be better simulated. We are also currently testing the performance of the model while accounting for the effects of diffusion.

Augmenting the simplified model with fibrin dynamics would allow patient- and event-specific characterization of thromboelastography experiments [37, 38]. This would be useful to better characterize clotting during surgeries [39]. This work could be further improved by coupling with simplified platelet aggregation models. To conclude, the drastically simplified model proposed here is a novel and a fertile step and subsequent progress has potential applications such as virtual flow-diverter treatment planning [40], predicting thrombotic risk accounting for flow [41], factors like injuries in tandem [42] and stenosed arteries [43].

4 Conclusion

We proposed a simplified model for thrombin dynamics based on the stoichiometry of certain important chemical species. The model, using switching of parameters based on a threshold, describes dynamics of thrombin akin to jump starting a car. The simplified model fits corresponding aspects of the full model well and different features of thrombin dynamics are easily alterable. The model is able to predict certain important changes in thrombin dynamics due to changes in prothrombin and antithrombin concentration. This prediction using the simplified model with the potential for clinical manifestation in hypercoagulable diseases is hypothesized as a necessary step in patient-specific simulations of clotting in realistic geometries.