Reaction-diffusion waves of blood coagulation

The main function of the coagulation system is terminating bleeding, caused by the injury of the vessel wall, with a fibrin clot. The reaction of fibrin polymerization appears at the final stage of the proteolytic enzymatic cascade where the activated clotting factors act as catalyst for activation of the others (Butenas and Mann, 2001; Orfeo et al, 2005). Mature form of fibrin molecules can aggregate into the long branching fibers and form a complex network which serves as thrombus scaffold. The key enzyme of the coagulation cascade is thrombin as it catalyzes fibrinogen conversion to fibrin and its concentration has a crucial influence on the kinetics of the clot formation (Hemker, 1993; Hemker and Béguin, 1995; Butenas and Mann, 2001). To prevent the spontaneous formation of thrombi the activation reactions are normally regulated by the action of plasma inhibitors (Pieters et al, 1988; Koppelman et al, 1995; Monkovic and Tracy, 1990; Panteleev et al, 2002). The balance between coagulation and anti-coagulation systems is important for the normal organism functioning and any alternations can lead to severe pathological states: thrombosis or, on the contrary, disseminative bleeding (Colman, 2006; Askari and Lincoff, 2010).

The influence of different factors on the coagulation process was studied both experimentally and using multiple theoretical approaches. As compared with the experiment, parameters in the theoretical studies can be varied much easier allowing to detect not only experimentally observed regimes of blood coagulation (Stortelder and Hemker, 1997; Leiderman and Fogelson, 2011; Krasotkina et al, 2000; Ataullakhanov et al, 1998; Bouchnita et al, 2016) but also to suppose their possible variations for the conditions that are hard to reproduce in the experiment (Zarnitsina et al, 2001). Model results also provide data about the possible spatiotemporal distribution of all the blood factors participating in the coagulation cascade, while the main parameter used to measure the dynamics of the clot growth experimentally is fibrin clot density(Ataullakhanov et al, 1998; Krasotkina et al, 2000; Ataullakhanov et al, 2002; Panteleev et al, 2006; Ovasenov et al, 2005). Recently, the hypothesized in multiple mathematical models shape of the thrombin concentration profile in blood plasma was obtained in the direct experiment (Dashkevich et al, 2012).

The biochemical reactions in plasma are normally described with a set of ordinary or partial differential equations (Stortelder and Hemker, 1997; Leiderman and Fogelson, 2011; Krasotkina et al, 2000; Jones and Mann, 1994; Butenas et al, 2004; Lacroix, 2012) while the blood flow is modeled using continuous, discrete or mixed methods. Despite the differences of the modeling approaches, the systems of equations used for the description of the biochemical reactions are rather close, and they can usually be reduced to some simplified models (Zarnitsina et al, 2001; Ataullakhanov et al, 1998; Bouchnita et al, 2016). In Section 2 of the current paper we consider a reaction-diffusion model of the intrinsic pathway of blood coagulation cascade and perform its reduction under certain assumptions. For the reduced system we prove a number of properties of its solutions (Section 3,4) that allow us to replace the analysis of the whole system by the analysis of one equation on the dynamics of thrombin concentration.

Thrombin is the key enzyme of the coagulation cascade as it catalyzes fibrinogen cleavage to fibrin which in turn forms haemostatic clot. Thrombin formation is a result of the prothrombin activation that can be initiated by different blood factors. Normally the process is launched by the tissue factor that is expressed on the surfaces of all human cells except the endothelium and is bared in case of the vessel wall damage (extrinsic pathway), or through the activation of the factor XI after the contact with the foreign surface (contact activation) (Orfeo et al, 2005, 2008; Gailani and Broze, 1991). Both pathways lead to the activation of factor X (Orfeo et al, 2005) that contributes to the thrombin activation. Once the thrombin concentration reaches the threshold value, further prothrombin activation takes place due to the positive feedback loops of the coagulation cascade (intrinsic pathway) (Orfeo et al, 2008, 2005; Panteleev et al, 2006). Thrombin controls activation of factor XI (Gailani and Broze, 1991) and also of factors V (Monkovic and Tracy, 1990) and VIII whose activated forms (Va, VIIIa) increase activity of factors Xa and IXa by formation of the prothrombinase and intrinsic kinase complexes respectively (Butenas et al, 2004; Orfeo et al, 2005; Baugh and Krishnaswamy, 1996; Scandura and Walsh, 1996) (Fig. 1).

As the result of these biochemical reactions, the thrombin wave propagates in the direction from the injury site to the vascular lumen. According to the experimental data, after thrombin reaches substantial concentration the speed of the clot growth does not anymore depend on the way of the initial activation of the coagulation system (Orfeo et al, 2008; Ovasenov et al, 2005) and thrombin wave profile stays constant in time (Butenas et al, 2004; Ataullakhanov et al, 1998; Tokarev et al, 2006). The speed of propagation of the thrombin concentration serves then as one of the important criteria for the validation of the mathematical models of the coagulation system. Moreover, this speed depends strongly on the activity of the blood factors and can be significantly lowered or increased in case of illnesses disturbing the normal functioning of the coagulation system (Colman, 2006; Askari and Lincoff, 2010). This phenomena can be simulated using mathematical models of the whole coagulation cascade (Tokarev et al, 2006; Zarnitsina et al, 1996b; Lobanov and Starozhilova, 2005) or the combination of analytical and numerical approaches as it was done by Pogorelova and Lobanov (2014).

In Section 5 of the current work we propose a new theoretical approach for the estimation of the speed of the thrombin wave propagation. After reducing the system of PDE describing the coagulation cascade to one equation, we build the analytical estimates on the speed of the propagation of the traveling wave solution using two different methods. We then compare the estimates given by analytical formulas with computational values of the speed as well as with the experimental data.

The main biochemical reactions of the coagulation cascade can be described by the following system of partial differential equations:

Here, $U_{i}$ and $V_{i}$ denote the concentrations of the activated and inactivated forms of the $i$-th factor respectively, $C_{1}$ denotes the concentration of the prothrombinase complex and $C_{2}$ corresponds to the concentration of the intrinsic tenase complex. The first term in each equation describes diffusion of the factor in blood plasma. The process of factor activation depends both on concentrations of its inactivated form and on the concentration of the activating factor, thus we consider the factor activation as the bimolecular reactions with the rate constants $r_{i}$. The inhibition of the activated factors is assumed to be the first order reaction with the rate constants $q_{i}$. Reaction of the thrombin activation obeys to the Michaelis-Menten mechanism with Xa and prothrombinase complex playing the role of enzyme and $K_{2m},\;\overline{K_{2m}}$ being the Michaelis constants (Zarnitsina et al, 2001, 1996a; Anand et al, 2006). The reactions of complex formation are considered as the bimolecular reactions. The rate constants of the complex formation are $k_{i}$; as the rates of the complex degradation are much lower, we neglect the corresponding terms. In this work we do not consider the reaction of the fibrinogen conversion to fibrin as it does not influence the other reactions of the cascade and thus can be decoupled from the system. The resulting model consists of fourteen PDEs and it takes into account the main reactions of the intrinsic pathway (1).

To perform further analysis of system (1) we will reduce the number of equations. First, let us apply the assumption of the detailed equilibrium to the complex formation equations and set:

This approximation is justified if the reaction constants $k_{510},h_{510},k_{89},h_{89}$ are sufficiently large. Next, we assume that the concentrations of the inactivated factors remain constant. This assumption implies that these factors are in excess and they do not limit the rate of the corresponding reactions. It allows us to consider the activation reactions as first-order reactions with $k_{i}=r_{i}V_{i}^{0}$. Then, we replace the Michaelis-Menten kinetics by the corresponding first-order reaction assuming that the constants $K_{2m},\;\overline{K_{2m}}$ are sufficiently large. Finally, we replace $P$ in the equation for $T$ using the approximation $P=P_{0}-T$. This relation is exact if $h_{2}=0$ and it gives a good approximation for small $h_{2}$ (Figure 2). The reduced system has the form:

where $\overline{k_{10}}=r_{10}U_{10}^{0},\;\overline{k_{2}}=r_{2}T_{0}$.

A similar model was previously considered by Zarnitsina et al (1996b); Anand et al (2006). The parameter values of this model were verified using the experimental data. We take these values for the current study (Table 1, Appendix C). The thrombin wave propagates with a constant velocity of about $0.05$ mm/min (Fig. 2) that is in a good agreement with the experimental data (Tokarev et al, 2006).

Consider the system of ordinary differential equations:

Its equilibrium points satisfy the following relations:

where $T$ is a solution of the equation $P(T)=0$. Here $P(T)=aT^{4}+bT^{3}+cT^{2}+dT$,

Hence stationary points of system (5) can be found through the stationary points $T^{*}$ of the equation

and equalities (6), (7).

Let us determine the number of positive roots of the polynomial $P(T)$. We set $P(T)=TQ(T)$, where $Q(T)=aT^{3}+bT^{2}+cT+d$. The number of positive roots of $Q(T)$ can be found as follows. First, we consider a function $Q^{\prime}(T)=3aT^{2}+2bT+c$. If it has no zeros, then $Q(T)$ is increasing and has one positive root if and only if $Q(0)<0$. Otherwise, we denote by $T_{1},T_{2}$ the nonzero solutions of the equation $Q^{\prime}(T)=0$: $T_{1,2}=(-b\pm\sqrt{b^{2}-3ac})/(3a)$. Then the polynomial $Q(u)$ has one positive root in one of the cases:

• •

  $T_{1}\leq 0,\;Q(0)<0$,

• •

  $0\leq T_{1}<T_{2},\;Q(0)<0\;{\rm and}\;Q(T_{1})>0,\;Q(T_{2})>0\;{\rm or}\;Q(T_{1})<0$

and it has two positive roots if $0<T_{2},\;Q(0)>0,\;Q(T_{2})<0$.

Stability of stationary points of the system (5) can be determined from the stability of stationary points of equation (8). The following theorem holds (see Appendix A for the proof).

Thus, we can make the following conclusions about the existence and stability of stationary points of the kinetic system of equation (5). It always has the trivial solution $u^{*}=0$. It has one (two) positive solution if and only if the polynomial $P(T)$ has one (two) positive roots. A positive solution $u^{*}$ is stable if and only if $P^{\prime}(T^{*})>0$.

We can now formulate the theorem on the existence of waves for the system (3).

The proof of the theorem follows from the general results on the existence of waves for monotone systems of equation (Volpert, 2008; Volpert et al, 2000). Let us note that conditions on the stability of stationary points follow from the assumption of the Theorem 2 and Theorem 1. We have $P^{\prime}(T^{*})>0$ in both cases since it is the largest root of the polynomial increasing at infinity. The sign of $P^{\prime}(0)$ is negative if there is no other root in between and positive if there is one more root.

Monotone travelling wave solutions of monotone systems are asymptotically stable (Volpert, 2008; Volpert et al, 2000) that gives global stability in the bistable case. In the monostable case the wave is globally stable for the minimal speed $c_{0}$ and stable with respect to small perturbations in a weighted norm for $c>c_{0}$ (Volpert et al, 2000).

The unique wave speed in the bistable case and the minimal wave speed in the monostable case admit minimax representations. Below we will use them for the bistable system as this case is more appropriate for the applications considered in the current work. Indeed, travelling wave solution of system (3) describes propagation of the thrombin concentration in the blood flow due to the reactions of the coagulation cascade. In this system the convergence to the travelling wave solution takes place only if the initial concentrations of blood factors exceed some critical level, otherwise the clot formation does not start because of the action of plasma inhibitors. Described dependency on the initial conditions and stability of zero solution correspond to the bistable case. In the monostable case, on the contrary, a small perturbation would result in the solution converging to the propagating wave. In terms of the coagulation system functioning, monostable case corresponds to the spontaneous disseminated coagulation blocking blood circulation.

Finally, let us note that in the Theorem 2 we consider only the case of a single positive root of the polynomial and the case of two positive roots. If $P(T)$ has three positive roots the system would be monostable with a stable intermediate stationary point. While this case is interesting from the point of view of wave existence and stability, it is less relevant for the modelling of blood coagulation, and we will not discuss it here.

In order to simplify the presentation, we describe the method of reduction for the system of two equations:

where $\varepsilon$ is a small parameter, $\dfrac{\partial f}{\partial v}>0$ and the system (9),(10) is bistable. If we multiply the second equation by $\varepsilon$ and take a formal limit as $\varepsilon\to 0$, then we have $v=\dfrac{a}{b}u$, and the first equation can be rewritten as follows:

Let us recall that the value of the speed $c=c_{\varepsilon}$ in system (9), (10) and $c=c_{0}$ for the scalar equation (11) are unknown, and in general they are different from each other. We will demonstrate that $c_{\varepsilon}\to c_{0}$ as $\varepsilon\to 0$:

Singular perturbations of travelling waves are extensively studied by Volpert (2008). Here we present another method of proof based on the estimates of the wave speed. This method is simpler and gives not only the limiting value of the speed for $\varepsilon=0$ but also the estimates of the speed value for any positive $\varepsilon$. In the following section we describe the approach in details and construct the wave speed estimates for the system (9), (10).

We get the following estimates from the minimax representation of the wave speed in the bistable case (Volpert et al, 2000) :

where

$\rho=(\rho_{1},\rho_{2})$ is an arbitrary test function continuous together with its second derivatives, monotonically decreasing (component-wise) and having the same limits at infinity as the wave, $\rho(+\infty)=0,\;\rho(-\infty)=u^{*}$.

Let us choose the following test functions:

where $u_{0}$ is the solution of the equation (11). Neglecting the second-order terms with respect to $\varepsilon$, we get:

where

and $c_{0}$ is the value of the speed in (11). Next,

where

Hence from (14), (15) we obtain the estimate

where $c_{0}$ is the wave propagation speed for (11), the functions $\varphi(x),\;\psi(x)$ are bounded. The proof of Theorem 3 follows from this estimate.

If the reaction rate constants in the equations of system (3) for the variables $U_{9}$, $U_{10}$, $U_{5}$ and $U_{8}$ are sufficiently large, then we can replace these equations by the following algebraic relations (cf. Section 4.1):

Then instead of system (3) we obtain the following system of two equations:

Similarly, we can reduce this system to the single equation:

We realize this reduction in two steps in order to compare the one-equation model with the system (3) as well as the intermediate model of two equations (17). Numerical simulations show that for the values of parameters in the physiological range, all three models give the wave speed of the same order of magnitude (Fig. 3). The two equation model (17) gives a better approximation of the model (3) than the single equation (18). However, the latter demonstrates the same parameter dependence of the wave speed as other models. Taking into account the complexity of the initial model (3), the approximation provided by one equation is acceptable. In what follows we obtain analytical formulas for the wave speed for the one equation model.

In dimensionless variables

we write equation (18) in the following form:

where:

Analysis of the rate constant values allows us to further simplify the equation. As $M_{3}\gg 1$ we can approximate equation (20) by the following equation:

Let us note that in the term $u_{1}^{2}(1+M_{2}u_{1})$ the first component corresponds to the prothrombin activation by the factor Xa and the second one corresponds to the prothrombin activation by the [Va, Xa] complex. Since during the propagation phase the rate of activation by prothrombinase complex is several orders of magnitude higher than the activation by Xa itself, we can neglect the first component. Thus, applying the assumption of the detailed equilibrium for the second equation, we finally obtain the following equation for the thrombin concentration:

where:

The equation (22) can be rewritten in the more general form:

Travelling wave solution of (24) satisfies the equation:

Here we will present two approximate analytical methods to determine the wave speed.