key: cord-0841333-k39dlw6u
authors: Ratto, N.; Bouchnita, A.; Chelle, P.; Marion, M.; Panteleev, M.; Nechipurenko, D.; Tardy-Poncet, B.; Volpert, V.
title: Patient-Specific Modelling of Blood Coagulation
date: 2021-03-27
journal: Bull Math Biol
DOI: 10.1007/s11538-021-00890-8
sha: 98844587602d9b165bb2707116cbc55502e27d37
doc_id: 841333
cord_uid: k39dlw6u

Blood coagulation represents one of the most studied processes in biomedical modelling. However, clinical applications of this modelling remain limited because of the complexity of this process and because of large inter-patient variation of the concentrations of blood factors, kinetic constants and physiological conditions. Determination of some of these patients-specific parameters is experimentally possible, but it would be related to excessive time and material costs impossible in clinical practice. We propose in this work a methodological approach to patient-specific modelling of blood coagulation. It begins with conventional thrombin generation tests allowing the determination of parameters of a reduced kinetic model. Next, this model is used to study spatial distributions of blood factors and blood coagulation in flow, and to evaluate the results of medical treatment of blood coagulation disorders.

Schematic representation of the model. It begins with an ODE model for a simplified scheme of coagulation reactions (blue) in order to describe thrombin generation curves and to determine patientspecific parameters (upper left rectangle). The same kinetic model with diffusion is used to describe the propagation of coagulation wave in quiescent plasma (lower left rectangle). A similar model completed by the Navier-Stokes equations and fibrin polymer production describes clot growth in flow (lower right rectangle). Parameters of patient-specific treatment are determined from thrombin generation curves and then used to model clot growth in flow (upper right rectangle) (Color figure online) from very simple containing several equations to very complex with more than one hundred of equations (see Belyaev et al. 2018; Chelle et al. 2018; Hemker 1993; Hemker and Beguin 1995; Hockin 2002; Hoffman et al. 2005; Mann et al. 2003 and the references therein). The main difficulty of complex models is that they contain many reaction rate constants, mainly unknown or with a large discrepancy of their values in the existing literature (Andreeva et al. 2018; Chelle et al. 2018 ). Since the experimental and clinical data are often reduced to thrombin generation curves (thrombin concentration in time), the inverse problem for the determination of the parameters of the model becomes highly under-determined. On the other hand, simplified models with few equations and a relatively small number of parameters imply certain assumptions and approximations (quasi-stationary approximation, detailed equilibrium, etc.). The validity of these assumptions can be difficult to verify. Therefore, clinical applications of kinetic models encounter serious difficulties which are not yet completely overcome. This is even more true for the spatial models, especially, for blood coagulation in flow, aimed to model in vivo clot growth (see, e.g., Anand et al. 2003; Biasetti et al. 2012; Kuharsky and Fogelson 2001; Sequeira et al. 2011; Xu et al. 2010) , because many other factors can influence this process.

Keeping in mind clinical applications as a long-term objective, we develop in this work an approach to patient-specific modelling of blood coagulation (Fig. 1) . We begin with thrombin generation tests for two group of subjects. The first one contains healthy and hemophilic subjects. The second group contains healthy and thrombotic subjects ("Appendix 1"). Fitting thrombin generation curves (TGCs) with a kinetic model, we identify patient-specific parameters. In view of the discussion above, we use a simplified model of three equations with 9 kinetic constants. Data analysis shows that three of them essentially characterize different subgroups of patients, while other constants are less important. The derivation of the simplified model and fitting the patient-specific thrombin generation curves were done in Ratto et al. (2020a) .

The main objective of this work is to use the simplified kinetic model with the patient-specific parameters in order to study spatial models of blood coagulation, clot growth in blood flow, and to develop patient-specific treatment of coagulation disorders. The kinetic ODE model is presented in the next section. The kinetic model augmented by the diffusion term is used to study the propagation of coagulation wave in a quiescent plasma (Sect. 3). We study the existence and stability of such waves, and the existence of pulses, which determine the range of initial conditions providing clot growth. These mathematical results provide a more precise understanding of these processes. We use patient-specific values of parameters to determine the speed of clot growth. In agreement with the previous study (Tokarev et al. 2006) , it is smaller in the hemophilic plasma in comparison with the normal plasma. The introduction of blood flow (Sect. 4) allows us to study clot growth for the three groups of subjects modelling in vivo conditions. The rate of clot growth and its final size are larger for thrombotic subjects, intermediate for healthy ones, and the smallest for hemophilic patients.

Finally, in Sect. 5 we model the treatment of thrombosis with anti-thrombotic drugs (warfarin, heparin). We begin with thrombin generation curves and choose the dosage of treatment which returns TGCs to their normal characterization. The same dosage is then used in the modelling of blood coagulation in flow.

Thrombin generation curves in platelet-poor plasma can be described by various ODE models representing equations of chemical kinetics for the reactions of plasma proteins called blood factors (see literature reviews in Belyaev et al. 2018; Chelle et al. 2018; Tokarev et al. 2019) . Coagulation system is a cascade of proteolytic enzymatic reactions with each level consisting of two processes: zymogen (coagulation factor) activation to the active enzyme (activated coagulation factor), followed by its rapid irreversible trapping by inhibitors always circulating in blood. On each level, the shortliving coagulation factor catalyses the reaction of activation on the next cascade level. The final product of coagulation cascade is fibrin which rapidly polymerizes into a three-dimensional mesh (gel) slowing blood flow and aggregating platelets and other blood cells into the clot.

There are two pathways of coagulation activation (Hoffman et al. 2005; Mann et al. 2003) . Activation by the extrinsic pathway begins when blood comes in contact with tissue factor (TF): this transmembrane protein is expressed by the majority of cells except those normally being in contact with blood. TF binds with FVIIa, which circulates in tiny amounts (1% of total FVII), making FVIIa able to cleave FIX and FX to their activated forms. Intrinsic activation pathway (Kondratovich et al. 2002) is initiated by the contact of blood with any "foreign" surface. Upon adsorption on this surface, FXII becomes activated due to conformational changes and then stimulates its own formation both autocatalytically and by activating prekallikrein to kallikrein: Fig. 2 The main activation reactions of the intrinsic pathway of the coagulation cascade. Thrombin (IIa) catalyses activation of factors V, VIII, XI; factors XIa and IXa catalyse activation of factors IX and X, respectively; factors VIIIa and Va form active complexes with factors IXa and Xa, respectively, and further increase thrombin production. Thrombin accelerates fibrin (F) production from fibrinogen (Fg). Fibrin polymer Fp forms the clot kallikrein activates its cofactor high molecular weight kininogen and FXII. Generated FXIIa activates FXI, FXIa activates FIX, and FIXa activates FX.

Both pathways unite at the activation of FX to FXa. FXa cleaves prothrombin (FII) to thrombin (FIIa), the central coagulation enzyme. In addition to cleavage of fibrinogen to fibrin, thrombin controls at least three positive feedback loops activating FV, FVIII, and FXI which are located above in the cascade. Two of these loops lead to the activation of cofactors FVa and FVIIIa which bind with FXa and FIXa, respectively, forming prothrombinase and intrinsic tenase complexes having activities 10 4 − 10 5 times larger than free enzymes have. Therefore, upon the initial activation by any pathway, local thrombin concentration increases in a dramatically nonlinear manner leading to full fibrinogen conversion to fibrin.

We consider the model of coagulation cascade shown schematically in Fig. 2 . The corresponding ODE system can be reduced to a simplified system of three equations for prothrombin P, thrombin T , and activated factor X , u:

Complete kinetic system of equations, and the method of reduction using physiologically based approximations are presented in (Ratto et al. 2020a ). The constant k 1 in Eq.

(3) characterizes the initiation stage, while all other positive terms of these equations correspond to the propagation stage, and two negative terms in Eqs.

(2), (3) to the termination stage. Reprinted with permission from Ratto et al. (2020a) Let us recall that system (2.1)-(2.3) is derived from a more complete model describing the coagulation cascade with individual reactions between the blood factors (Ratto et al. 2020a ). Due to this reduction, the coefficients of the simplified model represent combinations of reaction rates and initial concentrations of the more complete model. Furthermore, reaction rate constants can be sensitive to some other factors, which are not considered in the complete model but which can depend on individual patients. For example, thrombin activity depends on the concentration of Na+ (Huntington 2008) , etc. Thus, the coefficients of the simplified model of three equations can be patient-dependent.

The concentration of thrombin as a function of time (thrombin generation curve, TGC) during thrombin generation test is considered as an important characterization of blood coagulation process (Hemker 1993; Hemker and Beguin 1995; Mann et al. 2003) . Thrombin generation curves are mainly characterized by lag time t lagbeginning of explosive growth of thrombin concentration, time to maximum of thrombin concentration t max (also called time to peak, TtP), the maximal value of thrombin concentration T max , and endogenous thrombin potential (ETP)-the area under the curve (Fig. 3, left) . Typical TGCs for normal and hemophilia subjects are shown in Fig. 3 .

Various kinetic models are used to model thrombin generation curves (see, e.g., Belyaev et al. 2018; Chelle et al. 2018; Hemker 1993; Hemker and Beguin 1995; Tokarev et al. 2019 and the references therein). The interest of model (1)-(3) is that it contains a relatively small number of parameters, it includes the main blood factors, and it has a clear biophysical structure including the initiation, prolongation and inhibition stages of the coagulation cascade. Furthermore, this model is derived from a more complete model (Tokarev et al. 2006) , and it is validated by comparison with the Hockin model (Hockin 2002) considered as one of the benchmark models.

Analysis of data from Chelle et al. (2018) shows that there is a linear correlation between T max and ETP. The remaining three independent parameters t lag , t max , and ETP can be used to identify normal and hemophilic subjects (Ratto et al. 2020a ). If we Normal and hemophilic subjects are separated in the parameter space (k 2 , k 6 , k 9 ) characterizing thrombin generation curves in the case of platelet poor plasma. Hemophilic patients are shown by circles, healthy subjects by crosses. Reprinted with permission from Ratto et al. (2020a) use only one or two parameters, this separation is not so precise, and the two groups partially overlap. It should be noted that this method does not allow the separation of hemophilia A and B subjects. However, this result can depend on the organization of thrombin generation tests, physiological parameters and clinical conditions of the group of patients. In the second database considered in this work ("Appendix 1"), thrombotic and normal subjects can be segregated by T max .

Fitting of TGCs by system (1)-(3) allows the identification of healthy and hemophilia subjects by the parameters of this system. They are clearly separated in the space of three parameters k 2 , k 6 , k 9 (Fig. 4) , while the dependence on other parameters is less important. The details of the fitting procedure and parameter dependence can be found in Ratto et al. (2020a) ; Ratto (2020) . For the convenience of reading, it is shortly presented in the appendix. Let us note that parameters k 2 and k 6 represent combinations of original reaction rate constants and factor concentrations (Ratto et al. 2020a ). However, they mainly characterize the activation of factor X by thrombin and thrombin self-activation, respectively, while parameter k 9 describes the action of antithrombin.

Thrombin generation curves for platelet-rich plasma are qualitatively similar to those for platelet-poor plasma (PPP). As before, there is a clear separation of healthy and hemophilia subjects in the space of three parameters T t P, E T P, T lag (not shown). In the case of healthy subjects, the presence of platelets basically acts to increase ETP. In the case of hemophilia patients, the presence of platelets essentially widens the distribution of points (subjects) in the parameter space.

In the case of PRP, the corresponding simplified model contains one more equation for the concentration of activated platelets (Ratto 2020) . In order to simplify the presentation, we consider in the remaining part of the paper only the model (1)-(3) for the PPP.

Blood coagulation can be described as a reaction-diffusion wave (see, for example, Guria and Guria 2015; Panfilov 2019; Pogorelova and Lobanov 2014; Zarnitsina et al. 2001) . Wave propagation is based on self-amplifying feedback loop of thrombin production. Thrombin (and other blood factors) diffuses ahead of the reaction front, initiates thrombin production there, and so on, providing self-sustained propagation of the coagulation reactions in space (Sect. 3.1). In general, the speed of reaction-diffusion waves can be positive, zero or negative depending on the values of parameters. Positive speed of the coagulation wave corresponds to clot growth. We will see below that it can be zero or negative for hemophilia patients (Sect. 3.3). In this case, the clot does not form. The propagation of the coagulation wave also depends on the initial thrombin concentration which should be large enough in order to overcome some threshold level. This critical condition is determined by the pulse solution discussed below in Sect. 3.2. In the physiological conditions, the coagulation wave is initiated at the blood vessel wall, e.g. in the case of injury, and it propagates inside the vessel forming the blood clot. The wave propagation is influenced by the blood flow (Sect. 4).

In order to study the patient-dependent properties of such waves, we will consider the reaction-diffusion system corresponding to the kinetic system (1)-(3):

where we suppose that the diffusion coefficients of the three species are equal to each other. We set k 1 = 0 in Eq. (6) in order to describe propagation of reaction-diffusion waves of blood coagulation sufficiently far from the vessel wall. The case where this parameter is positive corresponds to the initiation of coagulation cascade in the bulk during thrombin generation tests.

Consider system (4)-(6) for all real x. If k 9 = 0, then the sum of Eqs. (4) and (5) gives ∂ P/∂t + ∂ T /∂t = 0. Therefore, if the initial conditions P(x, 0) and T (x, 0) for the corresponding concentrations satisfy the equation

where T 0 is some constant, then a similar equality holds of all positive t, P(x, t) + T (x, t) ≡ T 0 . Hence, system (4)-(6) can be reduced to the system

If k 9 is different from 0 but small enough, then system (7), (8) approximates system (4)-(6). Set

Therefore, (7), (8) is a monotone system characterized by the applicability of the maximum principle. Hence, we can use the results on the existence and stability of waves developed for such systems (Volpert and Volpert 1990 ). Let us begin with the analysis of stationary points. From the equation

and substitute into equation F 1 (T , u) = 0:

This equation has solution T = 0, and, according to numerical estimates, up to two positive solutions T * and T * . The existence of two positive solutions corresponds to normal physiological conditions with two stable points of the kinetic system. The first one T = 0, u = 0 corresponds to the case without blood coagulation, while the second one T = T * , u = u * corresponds to the production of thrombin and activated factor X in the coagulation cascade. The intermediate positive stationary point is unstable. Let us recall that reaction-diffusion wave is a solution of system (7),

where c is the wave speed and the functions θ(ξ ), ω(ξ ) satisfy the equations

and the conditions at infinity

Theorem 1 Suppose that Eq. (9) has two positive solutions, T * and T * , T * < T * . Then, problem (10)-(12) has a monotonically (component-wise) decreasing solution for a unique value of c. This solution is globally asymptotically stable, that is, for any monotonically decreasing initial condition T (x, 0), u(x, 0) with limits (12) at infinity, solution T (x, t), u(x, t) of system (7), (8) satisfies the following convergence

where h is some real number. The wave speed c is given by the following minimax representation:

where M is a set of all smooth monotonically decreasing functions w = (w 1 , w 2 ) with limits (12) at infinity.

The proof of this theorem follows from the general results obtained for monotone reaction-diffusion systems (Volpert and Volpert 1990 ).

A positive stationary solution v = (v 1 , v 2 ) of system (7), (8) with 0 limits at infinity is called a pulse solution. It is a solution of the following problem:

We will use the notation F = (F 1 , F 2 ). The existence and instability of a pulse solution are determined by the following theorem. (14) The proof of the existence of solutions is based on the Leray-Schauder method. It is similar to the proof for more complete models of blood coagulation (Galochkina et al. 2018; Ratto et al. 2020b ). The spectral properties of the operator L follow from Volpert and Volpert (2020) . Indeed, since the equation Lu = 0 has a nonzero solution u(x) = v (x), then 0 is an eigenvalue of the operator L and the corresponding eigenfunction w (x) is not positive. Therefore, 0 is not the principal eigenvalue, and there exists a positive eigenvalue of the operator L.

The pulse solution separates initial conditions of problem (7), (8) for which the solution of this problem grows or decays in time.

of the Cauchy problem for system (7), (8) satisfies the inequalities (7), (8) satisfies the inequalities

The proof of this theorem is presented in "Appendix 4". Thus, Theorems 2 and 3 provide the conditions of clot growth, that is, of growth of solution of system (7), (8): the wave speed c in problem (10)-(12) is positive, and the initial condition is greater than the pulse solution. The second condition implies that the quantities of thrombin and activated factor X produced at the initiation stage should be larger than the threshold determined by the pulse solution.

Let us recall that in the physiological conditions clot growth starts at the damaged wall of the blood vessel. The clot grows perpendicular to the wall inside the vessel. In this case, we consider system (7), (8) on the half-axis x ≥ 0 with the homogeneous Neumann boundary condition. Similar results on the existence and stability of solutions can be obtained for the problem on the half-axis.

Approximation of system (4)-(6) by system (7), (8) allows us to study existence and stability of waves and pulses, but it modifies the distribution of thrombin concentration during wave propagation. Thrombin distribution is not monotone for system (4)-(6) (Fig. 5, left) . However, if the value of the constant k 9 is small enough, the wave speeds for these two systems are close to each other.

Minimax representation (12) allows the estimation of the wave speed c approximating the experimental data (Galochkina et al. 2017a, b) . In this work, we carry out Wave speed for healthy subjects (circles) and hemophilia patients (crosses) (right) in µm/s. There are three subgroups of hemophilia patients: without wave propagation (zero speed), with a slow propagating coagulation wave, with a fast propagating wave numerical simulations of system (4)-(6) in order to determine the wave speed for the values of parameters corresponding to healthy subjects and to hemophilic patients. Numerical implementation of this system is based on an implicit-explicit finite difference method where the diffusion operator is considered implicitly while the reaction terms explicitly with respect to time. Thomas algorithm is used to inverse the tridiagonal matrix corresponding to the diffusion operator. Accuracy of numerical simulations is controlled by decreasing the space and time steps.

The values of parameters in numerical simulations are determined by fitting TGCs (Sect. 2). The wave speed for healthy and hemophilic subjects is shown in Fig. 5 (right). For most of the hemophilic patients, the wave propagation is not observed. According to the analysis of the previous section, we conclude that either the wave speed is not positive in this case or the initial thrombin production is not sufficient. The initial conditions in numerical simulation are taken sufficiently large. Therefore, we conclude that the wave speed is not positive, and we set it 0 in Fig. 5 (right) .

For some hemophilic patients, the wave speed is about 0.2 µm/s. This values is less than the wave speed for the majority of normal subjects. Finally, there are several patients with a high value of speed. The characterization of TGCs for these patients is also close to normal. We can assume that the lack of factors VIII or IX, specific for hemophilia A and B, is compensated in this case by some other factors. The first two groups of hemophilic patients are separated from normal subjects with respect to the wave speed.

Thus, we have identified three sub-groups of hemophilic subjects for which either the clot does not grow at all (zero wave speed), or it grows slowly (small positive wave speed), or it grows normally (high positive wave speed).

In this section, we will consider blood coagulation in platelet-free plasma taking into account the influence of blood flow. Platelet-free plasma is a Newtonian fluid, and fluid motion is described by the Navier-Stokes equations ∂v ∂t

where v = (v x , v y ) is the velocity vector, p is pressure, ν is the kinematic viscosity, ρ is fluid density, and K (F p ) is the hydraulic permeability of the fibrin clot which depends on the concentration F p of fibrin polymer (Wufsus et al. 2013) :

Equations for the concentrations of prothrombin, thrombin and activated factor X are now considered with convective terms

Let us recall that fibrin clot is formed by fibrin polymer due to the conversion of fibrinogen into fibrin catalysed by thrombin (Figs. 1, 2) (see, for example, Blomback 1996; Gailani and Renne 2007; Mosesson 2005 and the references therein). Therefore, the previous equations should be completed by the equations for fibrinogen F g , fibrin F, and fibrin polymer F p :

The reaction term in Eq. (24) describes consumption of fibrinogen during its conversion into fibrin with thrombin as an enzyme. Reaction terms in Eq. (25) characterize production of fibrin and its consumption for the production of fibrin polymer. It does not diffuse and is not convected [see Eq. (26)]. (26) is considered in two-dimensional rectangular domain 0 ≤ x ≤ L, 0 ≤ y ≤ H corresponding to the interior part of blood vessel (lumen, Fig. 6 ). Direction of blood flow is from left to right; the damaged part of the blood vessel wall is at the bottom. The clot grows from the damaged wall inside the vessel.

We consider no-slip boundary conditions for the flow velocity v at the upper and lower vessel walls b and t , and parabolic velocity profiles at the lateral walls l and r :

There are two possible formulations of flow conditions, with a given velocity at the entrance of the vessel or with a fixed pressure difference between the entrance and the exit. In the first case, the constant a is given; in the second case, it is unknown and chosen in such a way that p 1 − p 0 = δ p, where p 1 is the average pressure at the boundary l , p 0 at the boundary r , and δ p is a given number. The difference between these two conditions becomes essential if the clot size (height and length) is large enough. Then, its resistance influences velocity distribution. If the flow velocity at the entrance and, consequently, total flow rate are given, then flow velocity above the clot accelerates. Flow washes out blood factors from the clot and decelerates its grows. Complete vessel occlusion for a dense clot is impossible in this case. If the pressure difference is given and not the total flow rate, then flow velocity above the clot can either increase or decrease depending on its width (Bouchnita et al. , 2020 . Decelerating flow velocity promotes further clot growth which can lead to the complete vessel occlusion. We consider here the case of given pressure difference because it is more realistic from the physiological point of view. Since prothrombin and fibrinogen circulate in blood flow and provide necessary compounds for the formation of fibrin clot, we set their concentrations equal to some given constants at the entrance of the vessel and we consider no-flux boundary conditions for these concentrations at all other boundaries:

The no-flux boundary conditions for thrombin and fibrin are considered at all boundaries: Finally, simplifying the coagulation cascade, we assume that there is a flux of activated factor X from the damaged wall:

Here n is the outer normal vector, and is the total boundary.

Problem (19)-(30) is solved numerically. The details of the numerical implementation are presented in "Appendix 2". For numerical simulations, we use patient-specific parameters determined from thrombin generation curves for healthy, hemophilic, and thrombotic subjects ("Appendix 3"). The corresponding flow velocity and fibrin concentrations are shown in Figs. 7 and 8 at the moment of time when the clot stops growing and does not change any more. We see from Fig. 9 that clot growth stops faster for hemophilic and thrombotic subjects than for normal subjects. In the case of thrombotic subjects, the clot reaches the opposite vessel wall resulting in complete occlusion.

In the case of hemophilic plasma, the final clot size is small, and it does not practically influence the flow velocity. The intermediate clot size is reached for the healthy subject, and flow perturbation becomes visible. Finally, complete vessel occlusion occurs for thrombotic plasma, and flow velocity is close to 0. There is still some remaining flow due to clot permeability.

Let us note that clot shape is determined as the level line where the concentration F p (x, y, t) of fibrin polymer reaches 0.5 of the maximal concentration. This method gives a good accuracy since the gradient of this concentration is very sharp. The clot height is considered as the maximal distance from the vessel wall to this level line in the direction perpendicular to the wall.

Dynamics of clot growth is shown in Fig. 9 . For thrombotic subjects, clot reaches the upper vessel wall leading to complete occlusion. In the case of normal subjects, clot grows slower, and it will reach its final size later (not shown). Clot growth rate is even slower for hemophilic subjects. Thus, the difference in thrombin generation curves for normal, thrombotic, and hemophilic subjects manifests itself also in the case of blood coagulation in flow. This result can be expected, but it should be verified because spatial distribution of blood factors and blood flow influences clot growth. Moreover, the initiation of coagulation cascade occurs uniformly in the whole volume in thrombin generation tests, while it occurs at the vessel wall in the case of blood coagulation in flow.

Patient characterization in the space of three parameters (Sect. 2) allows us to develop treatment protocols which transfer thrombotic and hemophilic patients to the healthy zone. Such treatment of hemophilic patients was modelled in Ratto et al. (2020a) . Here, we model the action of anti-thrombotic drugs heparin and warfarin. Since characterization of TGCs in two experimental setups is different, we use here virtual thrombotic patients obtained by a random change of the coefficients k 2 , k 5 (increased), and k 9 (decreased). In order to model the experimental TGCs for thrombotic patients in the second experimental group, it is sufficient to increase the value of ETP keeping two other parameters the same as for healthy subjects.

Heparin increases the concentration of antithrombin. In order to take it into account, we introduce factor α in the corresponding coefficients, αk 4 and αk 9 , with α ≥ 1. Warfarin influences the activity of thrombin and of factor X. We describe its action by decreasing the initial concentrations of prothrombin, P 0 , and of factor X denoted by u 0 . We introduce the multiplicative factor β, 0 ≤ β ≤ 1 such that the initial concentrations of prothrombin and of factor X become, respectively, β P 0 and βu 0 . Then, β = 1 corresponds to the case without treatment, and β = 0 to the hypothetical case of complete elimination of these substances. Hence, the value of β measures the application of warfarin.

Since we model the action of heparin through the coefficients k 4 and k 9 , we analyse how these coefficients influence thrombin generation curves. Figure 10 shows their dependence on the coefficient k 9 . The dependence of k 4 is qualitatively similar (not shown). As expected, ETP decreases, while lag time and time to maximum increase. A similar behaviour of these parameters is observed for decreasing coefficient β (warfarin). However, the relations between these parameters are quantitatively different for the dependence on α in comparison with the dependence on β. Therefore, from the Figure 11 illustrates the action of warfarin on a virtual thrombotic patient characterized by the parameters of thrombin generation curves. This treatment does not bring the corresponding trajectory to the healthy zone. Similarly, treatment by heparin does not allow us to move TGC to the healthy zone (not shown). The joint action of these drugs is shown in Fig. 12 for three different virtual patients. A proper combination of these drugs chosen for each individual patient modifies the corresponding TGCs in such a way that they are located in the healthy zone.

When we know the parameters of treatment, we can now apply them to model blood coagulation in flow (Fig. 13 ).

Theoretical modelling of blood coagulation provides a better understanding of this process in various experimental and physiological conditions. Kinetics curves in thrombin generation tests are described by kinetic ODE models (Chelle et al. 2018; Dunster and King 2017; Hoffman et al. 2005; Mann et al. 2003 ) (see also Belyaev et al. 2018; Panfilov 2019; Tokarev et al. 2019 for a more complete literature review). Clot growth in a quiescent plasma is investigated with reaction-diffusion equations (Guria and Guria 2015; Pogorelova and Lobanov 2014; Zarnitsina et al. 2001) , clot growth in blood flow is studied by various continuous models (Belyaev et al. 2015; Bouchnita et al. 2020; Sequeira et al. 2011; Bodnar and Sequeira 2008 ; Sequeira and Bodnar 2014), particle Fig. 12 and hybrid models (Tosenberger et al. 2013 (Tosenberger et al. , 2016 Xu et al. 2010 ). These studies concern either platelet-poor plasma, or platelet-rich plasma, or the whole blood.

These theoretical studies have some limitations. Models with more complete kinetics contain too many equations and unknown constants, especially for the platelet-rich plasma. Moreover, these constants are patient-specific with possibly large variations (see Andreeva et al. 2018; Chelle et al. 2018 ) and the references therein). Models with simpler kinetics impose some assumptions which can be difficult to verify (Tokarev et al. 2019) . Specific features of clot growth can depend on physiological conditions and can be influenced by many different factors, such as inflammation or hypoxia . Some important questions remain open even at the qualitative level of understanding, such as clot formation during heart fibrillation (Karim et al. 2020 ). Let us also note that disseminated blood coagulation in small lung arteries is one of the main causes of mortality during the coronavirus disease (Ackermann 2020; Zuo 2020) .

These limitations also concern the application of theoretical modelling in clinical practice, mainly because of the large variation of patient-specific parameters and physiological conditions (Andreeva et al. 2018; Bouchnita et al. 2017b; Chelle et al. 2018) . In this work, we develop a modelling approach taking into account patient-specific parameters. We begin with simplified kinetic models and determine their parameters from the experimental thrombin generation curves. As we discussed above, simplified kinetics implies some assumptions on the model, but it would be impossible to reliably determine parameters of more complete models fitting TGCs. At the next stages of modelling, the same kinetics equations and individual-based parameters are used to study clot growth in quiescent plasma and in blood flow, and to evaluate treatment of hemophilia and thrombosis.

Thrombin generation curves can be efficiently characterized by three parameters, time to peak, lag time, and endogenous thrombin potential (ETP). In particular, hemophilic and healthy subjects are well separated in the space of three parameters, both for platelet-poor plasma and platelet-rich plasma (Ratto et al. 2020a; Ratto 2020) . Comparison of thrombotic and healthy subjects shows that a single parameter (ETP or thrombin maximum) can be sufficient to distinguish them. However, the experimental set-up for thrombin generation tests is different here in comparison with the first one. Therefore, data analysis can also depend on the laboratory conditions. Furthermore, this method does not allow the separation of hemophilia A and B patients.

Numerical simulations show that clot growth in quiescent plasma essentially differs for healthy and hemophilic subjects. The speed of wave propagation (clot growth) for the majority of healthy subjects is in the interval 0.25-0.7 µm/s, while for hemophilic patients, the wave does not propagate at all or its speed is about 0.2 µm/s. This result corresponds to the theoretical analysis establishing the conditions of clot growth: the wave speed should be positive and the initial condition should be large enough. Positive wave speed provides the existence of a pulse solution which determines the threshold for the initial conditions.

It is interesting to note that 5 hemophilic patients have the speed of wave propagation in the normal range, and 3 of them even at the upper limit of the normal range, in spite of the fact that these patients belong to the hemophilic zone in the (k 2 , k 6 , k 9 )-space (Fig. 4, right) . This means that patient characterization with TGCs is not identical to clot growth rate. The latter is essential for the determination of the final clot size in blood flow. Hence, some hemophilic patients can have sufficient clot formation, while application of anti-hemophilic drugs can move them to the thrombotic zone. Furthermore, the dependence of TGC separation on the parameter k 9 , characterizing the action of antithrombin, can be an additional indication of some compensatory mechanisms possibly existing for hemophilia patients. This question requires further investigations.

The necessity of patient-specific hemophilia treatment is also confirmed by thrombin generation tests with added factors VIII or IX lacking in hemophilic plasma. Characterization of TGCs shows the efficiency of this treatment (Ratto et al. 2020a ), but the dosage should be calculated for each individual patient. On the other hand, similar analysis shows that treatment with TFPI does not bring hemophilic TGCs to the healthy zone (Ratto et al. 2020a) .

Treatment of thrombotic patients seems to be even more delicate and patientdependent. We analysed the action of heparin and warfarin. Taken separately, they may not bring TGCs to the healthy zone. Their combination calculated for individual patients can be more efficient from the point of view of modifying thrombin generation curves in such a way that their characteristics become similar to those for healthy subjects.

Numerical simulations of clot growth in blood flow show essential difference between healthy, hemophilic, and thrombotic subjects with respect to clot growth speed and its final size. In order to have a more precise characterization of these groups, more complete studies are needed for a larger cohort of patients and taking into account additional factors, such as platelet aggregation and protein C, which influence the dynamics of clot growth.

Summarizing, we expect that the approach developed in this work can help in patient-specific diagnosis and treatment of coagulation disorders. Further investigations are required to fill the gap between theoretical modelling and clinical applications. to the centre's follow-up protocol so that no additive venipunctures were performed. There are 32 healthy persons and 29 subjects. The blood samples were processed by double centrifugation in order to obtain platelet-free plasma. The study was approved by the ethical committee of the centre (Table 2) . 

We denote the time step by dt, and, at time step n, we denote by v n the flow velocity. The space step is the same for the x−axis and the y−axis and is denoted by h.

Step 1: Nonlinear term. We compute an intermediate velocity

where

In order to avoid convection dominated problems, we use an upwind method to discretize (v n x ) + x , (v n x ) + y , (v n y ) + x and (v n y ) + y .

Step 2: Viscosity terms. We treat the viscosity terms with an implicit method. To compute the second intermediate velocity (v * *

x , v * * y ), we solve the equations:

in which subscripts x x and yy denote the discretization of the second derivative of v * * x and v * * y .

Step 3: Pressure correction. Finally, we compute the velocity v n+1 at time n + 1 by correcting the second intermediate velocity

where the pressure gradient is obtained implicitly by solving the Poisson problem:

To solve problem (40), we use the SOR method, which converges faster than the Jacobi and Gauss-Seidel methods. The SOR method is an iterative method; hence, we repeat the iteration:

until it converges. The constant ω is the relaxation factor. We denote by N x × N y the size of the grid mesh used for the computation. Hence, in order to have a fast convergence, we set: 

) is a pulse solution, (ii) If the initial condition T (x, 0), u(x, 0) of the Cauchy problem for system (7), (8) satisfies the inequalities T (x, 0) ≥ṽ k 1 (x), u(x, 0) ≥ṽ k 2 (x), x ∈ R, for any k, then T (x, t) → T * , u(x, t) → u * as t → ∞, uniformly on every bounded interval. (iii) If the initial condition T (x, 0), u(x, 0) of the Cauchy problem for system (7), (8) satisfies the inequalities

for any k, then T (x, t) → 0, u(x, t) → 0 as t → ∞, uniformly on the whole axis (Tables 5) .

Proof Consider a sequence h k such that h k → 0 as k → ∞ and set v k (x) = max(v(x), v(x + h k )),v k (x) = min(v(x), v(x + h k )).

Then, condition (i) is satisfied.

when thrombin concentration approaches to zero. Consider a sequence x n corresponding to the values of some concentration at the moments of time t n . Then, we set y 0 = x 0 , y n+1 = αx n+1 + (1 − α)x n . Therefore, the first value of the sequence remains the same, while all other values are obtain as a linear interpolation between two neighbouring points. The value of α determines the rate of smoothing. We use α = 0.9 (weak smoothing) in the beginning of the thrombin generation curve and α = 0.2 at the tail of the curve where the perturbations are stronger. The smoothing algorithm is used twice, and the fitting algorithm described below is applied to the smoothed curves.

For a set of parameters K = (k 1 , . . . , k 9 ) of the reduced model (1)-(3), we denote by K the corresponding solution. We aim to compute the gradient of K with respect to K . A numerical approach can be implemented, since the gradient of K can be computed considering small variations of K . For the data set D(t i ), i = 1, . . . , n determined at time t = t i , we introduce the error function

where · denotes Euclidian norm. Depending on the available data, D(t i ) can be a scalar or a vector variable, and K corresponds here to the same components of the vector. The gradient with respect to K of the error function is as follows:

where T denotes the transposition. The gradient method converges but the error function E(K ) possesses several local minima, so that the result can depend on the choice of the initial approximation. We use Monte Carlo algorithm to find several starting points, and then, we apply the gradient method for each one of them. We keep the set of parameters that minimize the error, and we stop the process if the relative errors for the ETP, the lag time, the maximum and the time to peak are less or equal to 10 %.

The relative error in the data fitting is given in the following Since we know how the characterization of thrombin generation curves depends on each parameter (see, for example, Fig. 10 for k 9 ), we can estimate the error in the determination of parameter values. For the parameter k 9 , we use the dependence of time to peak and maximum of thrombin (not shown). The relative error for k 9 is estimated by ±5%. The dependencies of ETP and of lag time give larger error estimates for which the value of time to peak will be outside the admissible interval. Therefore, we should use the smallest estimate.

In further estimates, we restrict ourselves to the parameters k 2 and k 6 since other parameters are less important for the characterization of thrombin generation curves. We obtain 30% for k 2 (lag time) and 10% for k 6 (maximum of thrombin). Let us note that the worst estimate 30% remains acceptable compared to the estimates from the other works (cf. Chelle et al. 2018) .

Pulmonary vascular endothelialitis, thrombosis, and angiogenesis in Covid-19

A model incorporating some of the mechanical and biochemical factors underlying clot formation and dissolution in flowing blood

Mathematical modelling of platelet rich plasma clotting. Pointwise unified model

Boundary conditions at a naturally permeable wall

Threshold of microvascular occlusion: injury size defines the thrombosis scenario

Modeling thrombosis in silico: frontiers, challenges, unresolved problems and milestones

An integrated fluid-chemical model toward modeling the formation of intra-luminal thrombus in abdominal aortic aneurysms

Fibrinogen and fibrin -proteins with complex roles in hemostasis ans thrombosis

Numerical simulation of the coagulation dynamics of blood

Conditions of microvessel occlusion for blood coagulation in flow

Modeling of the effects of IL-17 and TNF-on endothelial cells and thrombus growth

A mathematical model to quantify the effects of platelet count, shear rate, and injury size on the initiation of blood coagulation under venous flow conditions

Evaluation and calibration of in silico models of thrombin generation using experimental data from healthy and haemophilic subjects

Numerical solution of the Navier-Stokes equation

Mathematical modelling of thrombin generation: asymptotic analysis and pathway characterization

Intrinsic pathway of coagulation and arterial thrombosis

Reaction-diffusion waves of blood coagulation

Traveling wave solutions in the mathematical model of blood coagulation

Initiation of reaction-diffusion waves of blood coagulation

Spatial aspects of blood coagulation: two decades of research on the self-sustained traveling wave of thrombin

Flow effects on coagulation and thrombosis

Thrombin generation, an essential step in haemostasis and thrombosis

Thrombin generation in plasma: its assessment via the endogenous thrombin potential

Rethinking the coagulation cascade

A model for the stoichiometric regulation of blood coagulation

How Na+ activates thrombin-a review of the functional and structural data

The left atrial appendage in humans: structure, physiology, and pathogenesis

Spatiotemporal dynamics of contact activation factors of blood coagulation

Surface-mediated control of blood coagulation: the role of binding site densities and platelet deposition

The dynamics of thrombin formation

Fibrinogen and fibrin structure and functions

Antithrombin III: critical review of assay methods. Significance of variations in health and disease

Activated protein C assays: A review

Two subpopulations of thrombin-activated platelets differ in their binding of the components of the intrinsic factor X-activating complex

Influence of enzymatic reactions on blood coagulation autowave

Modélisation de la coagulation sanguine : propagation d'onde et application à des sujets cliniques

Clustering of thrombin generation test data using a reduced mathematical model of blood coagulation

Existence of pulses for a reaction-diffusion system of blood coagulation

Blood coagulation dynamics: mathematical modeling and stability results

Blood coagulation simulations using a viscoelastic model

Spatial dynamics of contact-activated fibrin clot formation in vitro and in silico in haemophilia B: effects of severity and Ahemphil B treatment

Mathematical modeling of thrombin generation and wave propagation: from simple to complex models and backwards

Modelling of thrombus growth in flow with a DPD-PDE method

Modelling of platelet-fibrin clot formation in flow with a DPD-PDE method

Applications of the rotation theory of vector fields to the study of wave solutions of parabolic equations

Spectrum of elliptic operators and stability of travelling waves

Platelets in inflammation and thrombosis

The hydraulic permeability of blood clots as a function of fibrin and platelet density

A multiscale model of venous thrombus formation with surface-mediated control of blood coagulation cascade

Dynamics of spatially nonuniform patterning in the model of blood coagulation

Prothrombotic autoantibodies in serum from patients hospitalized with COVID-19

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

The work has been supported by the RUDN University Strategic Academic Leadership Program (VV). This research was performed within the framework of the Development Program of the Interdisciplinary Scientific and Educational School of Lomonosov Moscow State University "Photonic and Quantum technologies. Digital medicine" (DN).

We use the database from the clinical study number NCT02300519 on clinicaltrials.gov and named "Thrombin Generation Numerical Models Validation in Haemophilic Case". It contains the thrombin generation data for PPP (platelet-poor plasma) patients. For each patient, we have the evolution of thrombin concentration over time. The initial concentrations of factor II (prothrombin), V, VII, VIII, IX, X, XI, XII, Fg (fibrinogen), AT (antithrombin), and TFPI are measured. For the hemophilic patients, the concentrations of TFPI antibodies are also given. The data include a set of 40 healthy subjects and 86 hemophilia patients. We have the TGC for both PPP and PRP samples. We focus here on PPP samples, but a similar study can be done with PRP samples (Table 1) .

We present here the data for subject with thrombosis tendency. The results presented here came out from an observational study conducted in women undergoing caesarean section at the Research Center for Obstetrics, Gynecology and Perinatology (Moscow, Russian Federation). All participants had elective caesarean and a absence of a history of psychiatric diseases (including alcohol-and drug-induced and of trauma or surgical treatment in the 90 days before the caesarean section). Recruitment occurred at the date of hospitalization and not less than 24 h before the caesarean section. Blood samples were obtained during routine venipunctures after the caesarean delivery due

Numerical implementation used in this work is similar to the one presented in Bouchnita et al. (2017a) . The reaction-diffusion system of equations is solved using a finite difference method. The blood flow is driven by the pressure difference p 1 − p 0 = δ p, where δ p is a given number. We use the Chorin projection method for the implementation of the Navier-Stokes equations, and the successive over-relaxation (SOR) method for the pressure Poisson problem.We now be briefly summarize the implementation. We begin by describing the Chorin projection method (Beavers and Joseph 1967; Chorin 1968) used to solve the Navier-Stokes equations for the incompressible fluid. The computation occurs in three steps and uses some intermediate velocities. The Navier-Stokes equations for the incompressible fluid read:

For the numerical implementation, we use the following value of parameters (Figs. 7 and 8):For the fluid, we use the following values: Modelling with treatment (Fig. 13, Tables 3, 4 and 5).

Theorem 3 There exist sequences of functionsṽ k (x) = (ṽ k 1 (x),ṽ k 2 (x)) andv k (x) = (v k 1 (x),v k 2 (x)) such that:Let us now verify (ii). Functionṽ k (x) is a lower function for system (7), (8). This means that solution T (x, t), u(x, t) with initial conditionṽ k (x) is a growing function of t for any x fixed. Since T * , u * is a stationary solution of this system, andṽ k 1 (x) < T * , v k 2 (x) < u * , then the same inequalities hold for the solution:Hence, there exist the limitsBy virtue of (45),Next, we prove thatIndeed, suppose that this is not the case. Then, there are two possible cases. In the first one, the functions T 0 (x) and u 0 (x) are monotonically increasing or decreasing. This means that system (10), (11) has a solution for c = 0. This contradicts Theorem 2 stating that c > 0. If at least one of these functions is not monotone, then we setFor h sufficiently small,is an upper function for system (7), (8), then the solution with this initial condition is monotonically decreasing with respect to t. We obtain a contradiction with convergence (46). Thus, (47) is proved. In order to prove (ii), it remains to note that for any initial condition greater than v k (x) for all x, the solution remains greater than T (x, t), u(x, t) and, consequently, converges to T * , u * uniformly on every bounded set. The proof of (iii) is similar. The theorem is proved.

For the convenience of reading, we present here a shortened description of the fitting algorithm used in Ratto et al. (2020a) ; Ratto (2020) in order to approximate thrombin generation curves by system of equations (1)-(3).

At the first stage of data treatment, experimental thrombin generation curves are smoothed because the relative value of their perturbation become sufficiently high