Mathematical Modeling of Intravascular Blood Coagulation under Wall Shear Stress.

Increased shear stress such as observed at local stenosis may cause drastic changes in the permeability of the vessel wall to procoagulants and thus initiate intravascular blood coagulation. In this paper we suggest a mathematical model to investigate how shear stress-induced permeability influences the thrombogenic potential of atherosclerotic plaques. Numerical analysis of the model reveals the existence of two hydrodynamic thresholds for activation of blood coagulation in the system and unveils typical scenarios of thrombus formation. The dependence of blood coagulation development on the intensity of blood flow, as well as on geometrical parameters of atherosclerotic plaque is described. Relevant parametric diagrams are drawn. The results suggest a previously unrecognized role of relatively small plaques (resulting in less than 50% of the lumen area reduction) in atherothrombosis and have important implications for the existing stenting guidelines.

1 Introduction

The permeability of the vessel wall with respect to primary procoagulants plays an important role in the initiation of intravascular blood coagulation [1-8]. Within the context of the Virchow’s triad [9] variations in wall permeability contribute to vascular change and, alongside blood flow alterations [10, 11] and abnormalities of blood constituents, have been associated with numerous cases of clinical thrombosis [12].

Shear stress plays a pivotal role in controlling the barrier function of the vessel wall [8, 13–16]. In this paper we aim to investigate how increased shear stresses associated with stenosis—a local narrowing of the vessel wall—contribute to thrombogenicity of atherosclerotic plaques.

Up to now the analysis of blood coagulation under increased shear stress conditions has been primarily concerned with platelet activation—a distinct type of platelet response occurring in shear flows in which the shear rates exceed ∼ 5400s −1 [17-19]. A number of reports however describe clotting induced by shear rates as low as 1000 s −1 [5]—much lower than required for platelet activation. We speculate that some of these incidences can be explained by the shear stress-induced permeability of the vessel wall to procoagulants building up inside the atherosclerotic plaque.

To investigate this scenario we propose a mathematical model of blood coagulation in a vessel with a local stenosis and a source of procoagulants in the surrounding tissue. Importantly, the permeability of the vessel wall to these procoagulants is a function of wall shear stress.

Numerical analysis of the model reveals the existence of two thresholds in the hydrodynamic activation of blood coagulation in the system: one associated with the increase of wall permeability and another one linked to the washing out of procoagulants in the intensified blood flow. Additionally typical scenarios of intravascular clot formation are described and the relevant parametric diagrams of the blood liquid state stability are drawn.

The model also sheds some light onto the relationship between the thrombogenicity of the atherosclerotic plaque and it’s shape and severity. The results obtained are important for critical revision of the indications for stenting procedures currently used in the clinic.

2 Model description

Mathematical modeling of blood coagulation has been a subject of intensive research in the last three decades. Extended mathematical models have been proposed, however criticism for them has been wide-ranging [20, 21]. For this reason the model introduced in this paper is built within the framework of a phenomenological approach in which blood coagulation is treated as a phase transition from a liquid to polymerized and gelled state [22-25]. To supplement the approach we introduce a set of closing relations that allow us to formally treat the polymerized and non-polymerized blood states in a unified manner.

We consider the problem in a two-dimensional approximation assuming a steady blood inflow for simplicity. The thrombogenic properties of the tissue surrounding the vessel are represented as a non-zero concentration of “primary activator” in the vicinity of the vessel wall. The activator can leak into the lumen through a permeable wall. Importantly, the permeability of the vessel is a function of shear stress on the vessel wall.

2.1 Hydrodynamics

The geometry of the problem is shown in Fig 1. The vessel walls are assumed to be rigid. For specificity, the vessel profile is approximated with the following formula: where L denotes vessel length, H denotes plaque height and d corresponds to one half of stenosis width (see Fig 1).

Fig 1

Geometry of the vessel fragment.

L , L and H correspond to vessel length, cross section diameter and plaque height. Γ and Γ denote inlet and outlet boundaries respectively. Γ+ and Γ− refer to upper and lower vessel walls respectively.

The blood flow is described using the modified Navier-Stokes equations [26]: where is blood velocity, t denotes time, is the Hamilton’s operator, p is pressure, and ρ and ν depict blood density and kinematic viscosity, respectively.

The last term in Eq (2) aims to account for the effect of thrombus formation on the blood flow. The effect is described in terms of the Darcy’s filtration law, with the multiplier α (M 1, M 2) standing for the filtration resistance of a blood clot (in case one has formed). The values of M 1 and M 2 reflect first and second statistical moments of fibrin polymer mixture (see section 2.3). An explicit dependence of α on M 1 and M 2 is given in section 2.4.

No-slip boundary conditions are assumed for the vessel walls Γ+ and Γ−. The inlet and outlet boundary conditions are given as follows: where L denotes vessel width and V 0 is inlet blood flow velocity in the center of the vessel.

2.2 Blood vessel and activation source

The pro-coagulant factors that enter the lumen from atherosclerotic plaques are mostly products of inflammatory processes that take place inside the plaques [27, 28]. In this work we do not focus on their biochemical nature and describe their effect in terms of a phenomenological “primary activator” of blood coagulation [22-25] denoted as u.

The concentration of primary activator under the vessel wall u 0 is assumed to be non-zero in the vicinity of the stenosed wall Γ−. The “healthy” non-atherosclerotic tissue is assumed to be devoid of primary activator:

The activator is assumed to be “leaking” into the vessel from the underlying tissue. In the mathematical model this is represented by the boundary condition imposed on u on Γ±: where μ is the vessel wall permeability, D denotes diffusion coefficient, operator denotes the space derivative normal to Γ±, u∣Γ is the concentration of the primary activator in the blood flow near the corresponding vessel wall (Γ±).

The permeability μ is assumed to depend on wall shear stress . The dependence is approximated with a piecewise-linear function: where τ 1 and τ 2 are threshold shear stress values that define the degree of vessel wall resistance to stress; μ 1 denotes the permeability of an intact vessel wall, while μ 2 corresponds to the maximal permeability of the vessel.

The concentration of the primary activator u is taken to be zero on the inlet boundary Γ. Zero-gradient conditions are assumed on Γ [29, 30].

2.3 Blood coagulation cascade

The kinetics of blood coagulation reactions is described in the framework of the phenomenological model introduced in [22-25]: Here θ and φ denote concentrations of the activator and the inhibitor of biochemical network of blood coagulation reactions (see [22, 31–33]) and F corresponds to fibrinogen (fibrin precursor) concentration. The components of the coagulation cascade are subject to convection and diffusion as described by the last term in Eqs (11)–(14).

The formation of fibrin polymers from fibrinogen is described in terms of the first two statistical moments of fibrin chain length distribution:

The first and second fibrin moments M 1 and M 2 are defined through the concentration of k-mers of fibrin F as [22, 34]:

M 1 reflects the total amount of fibrin-monomer molecules in polymerized and non-polymerized forms. The ratio M 2/M 1 determines weight-averaged number of fibrin-monomers in polymer molecules of fibrin [35-37]:

Just as Eqs (11)–(14) the polymerization equations take into account diffusion and convection, accounting for the fact that, generally speaking, both depend on the polymer molecule size. This dependency is incorporated into the coefficient of polymer chains transport b and diffusion coefficient D . The explicit expressions for b and D as functions of fibrin-polymer moments are given in the next section.

The values of u, θ, φ, M 1 and M 2 on the boundary Γ are supposed to be equal to zero in all numerical experiments, while F concentration on Γ is assumed to be equal to initial fibrinogen concentration . The zero-gradient conditions were adopted for all chemicals on Γ [29, 30]. Vessel walls are taken to be impermeable for θ, φ, F , M 1 and M 2.

2.4 Special features of polymers movement

To close the mathematical description explicit definition of D , b and α is needed. For the purpose we develop a phenomenological approach that allows us to formally treat the polymerized and non-polymerized states of blood in a unified manner.

To describe dependence of polymer chains diffusion coefficient D on polymer chain size the following asymptotic expression is used: where N is a phenomenological parameter of reptation theory [38, 39]. An asymptotic Eq (19) has very clear limit forms:

when fibrin molecules experience only hydrodynamic friction [39, 40], N ≪ N , for the diffusion coefficient one has: D ∼ 1/N ;

in accordance with reptation theory [38, 39], when N ≫ N and .

The expression for b was obtained in the following way. Single polymer chain in fluid flow experiences the following drag force (we neglect hydrodynamic interaction effects in this work):

This force makes the chain to move with the speed:

Thus b is:

It was assumed that in blood the prevailing mechanism of fibrin molecules nucleation is a heterogeneous one and that microparticles and some of the formal blood elements serve as nucleation centers [41]. If the concentration of heterogeneous centers of nucleation is equal to n 0 ([n 0] = [cm −3]), the distance between the centers may be estimated as .

Assuming fibrin polymer chains to be Gaussian, characteristic size of the coil may be estimated as [40]: where l 0 is the characteristic size of single fibrin-monomer molecule, and K is number of fibrin monomers in Kuhn’s segment. l 0 can be estimated from the fibrin molecule volume: .

The transition between diluted and semidiluted solution takes place when and when the distance between centers of nucleation is equal to R. This gives the value [38]:

It is easy to see that the coefficient of polymer chains transport by flow b is essentially decreased if N > N . At the same time it is known that entanglement effects are essential when . Hence, it seems to be reasonable to suppose that the magnitudes of N and are of the same value:

Coefficient α determines the degree of slowing down the velocity of blood when large polymers or gel are formed. In this case α is a filtration resistance inversely proportional to filtration permeability in the Darcy’s law. In the case of dilute polymer solution there is almost no slowing effect, α is given by: where ξ represents characteristic size of polymer gel/solution. It is known that for diluted solutions b goes to unity. In accordance with expression (26) the value of filtration resistance drops to zero. The value of ξ reflects average pore size in gel and it is equal to the distance between polymer molecules in solution.

In dilute solution the distance between polymer molecules is more than molecules size ξ ≫ R and the average concentration of monomers in the solution, M 1, is much smaller then the local concentration of monomers in one polymer M : M ≫ M 1. In mature gel opposite conditions are satisfied: ξ < R, M < M 1. General form of expression for ξ may be given as [36, 38]: where power m must be chosen to make ξ independent on N .

Local concentration of fibrin-monomer in the Gaussian coil is estimated as: where N /N is a number of moles of monomers in one polymer (N is Avogadro number) and R 3 estimates the volume that polymer occupies.

Substituting Eq (28) to Eq (27) and choosing m = 1 [38] gives: This gives for α : where

Equations describing polymerization in terms of statistical moments may have singular solutions (particularly, when M 2 blows up) [42]. In present work we focus on the early stages of fibrin gel formation in the blood flow. It was assumed that gel is worth to be considered as sufficiently “mature” when the mean weight-averaged number of fibrin-monomers in polymer molecules N exceeds the value corresponding to the semidiluted solution condition by at least two orders of magnitude (). Research of more condensed gel evolution remains outside the scope of this work’s objectives.

In the mathematical model presented in this work, if N exceeded “mature gel” value () in some region in the vessel, then M 2 at that region is assumed to be equal to its maximal value .

2.5 Initial conditions and numerical parameter values

At the initial moment t = 0 all variables but F were assumed to be equal to zero in the interior part of calculation domain while F was assumed to be equal to . The and p fields are taken to correspond to stationary flow under given boundary conditions.

Numerical simulations were performed assuming L = 7.5 cm and L = 1 cm. The values of other parameters used in numerical calculations are listed in Table 1 (see also [24, 25]).

Table 1

Parameter values.

Parameter	Value	Refs.	Parameter	Value	Refs.
α	3.33 ⋅ 10−2 s −1	[22, 23]	k b	1.67 ⋅ 10−3 s −1	[22, 23]
θ 0	5 nM	—“—	n 0	1010 cm −3	[43]
χ 1	8.33 ⋅ 10−4 s −1	—“—	Fg0	9 ⋅ 103 nM	[22, 23]
γ	8.33 ⋅ 10−2 (nM ⋅ s)−1	—“—	D u	3 ⋅ 10−7 cm 2/s	—“—
β	2.5 ⋅ 10−5 s −1	—“—	D φ	3 ⋅ 10−7 cm 2/s	—“—
c	5 nM	—“—	D θ	3 ⋅ 10−7 cm 2/s	—“—
ϵ g	1.66 ⋅ 10−6 s −1	—“—	D g	3 ⋅ 10−7 cm 2/s	—“—
φ 0	0.05 nM	—“—	D	3 ⋅ 10−7 cm 2/s	—“—
χ 2	0.35 nM	—“—	k d	1.66 ⋅ 10−6 s −1	—“—
k g	5 ⋅ 10−6 (nM ⋅ s)−1	—“—	k u	1.66 ⋅ 101 s −1	—“—
k p	2.5 ⋅ 10−4 (nM ⋅ s)−1	—“—	k r	1.67 ⋅ 10−2 s −1	—“—
τ 1	10 dyn/cm 2	[15, 44, 45]	μ 1	2 ⋅ 10−12 cm/s	[24, 25]
τ 2	20 dyn/cm 2	[15, 44, 45]	u 0	100 nM	[24, 25]
ν	5 ⋅ 10−2 cm 2/s	[43]	ρ	1 g/cm 3	[43]
K	10	[46]	v 0	3 ⋅ 10−18 cm 3	[46, 47]

3 Numerical methods

Eqs (2)–(14), (15) and (16) were solved numerically by means of finite volume method [48, 49] with the aid of splitting technique [30, 50]. To discretize convective terms in Eqs (11)–(14), (15) and (16) upwind differencing was used [30, 48]. The Laplacian term was discretized by the means of “over-relaxed correction” technique [49]. To solve equations of fluid motion Eqs (2) and (3) PISO method was used [49, 51] with linear approximation of convective terms. Chemical kinetics ODE’s were solved by backward differencing technique [30, 52].

The program code was written with the help of OpenFOAM [53] and SUNDIALS [52] open-source libraries. Numerical simulations were carried out on unstructured meshes consisting of both triangles and quadrangles. The boundary layers in the vicinity of the vessel walls were meshed by quadrangles whereas the other parts of mesh were meshed by triangles. The meshes were generated using the open-source program SALOME [54] and adaptively refined in the near-wall area and along the separatrix line (see section 4.3). The results were visualized using ParaView [55] and gnuplot [56] open-source programs.

4 Results

Eqs (2), (3) and (11)–(16) with boundary conditions (1)–(6) and (9) and initial conditions described in section 2.5 were solved numerically. Generally speaking, the system considered is governed by two dozens of dimensionless “true” parameters easily determined by scaling procedures [57]. We have focused our analysis on the dependence of the solution on following four non-dimensional parameters related to blood flow intensity, wall permeability and plaque shape:

The rest of the parameters were held fixed as in Table 1.

4.1 Threshold-like activation of intravascular thrombus formation

Of primary importance is mapping the blood coagulation regimes in the system. We define coagulation in terms of fibrin gelation occurring anywhere in the vessel. Mathematically, this corresponds to N = N (x, y, t) reaching in the flow domain at a finite time. In contrast, the no thrombi formation regime is defined such that holds true everywhere in the vessel.

Numerical analysis of Eqs (2), (3) and (11)–(16) shows that both regimes — with and without blood coagulation — are observed in the system within the explored range of parameters. Fig 2 shows the parametric diagram of the regimes in the parameter space. Label “I” is used to denote the stationary regimes with no coagulation in the vessel while “II” marks the regimes with thrombi formation.

Fig 2

Parametric diagram of blood coagulation system regimes in the parameter space.

Re is the Reynolds number and is the non-dimensional vessel wall permeability. Label “I” is used to denote stationary regimes with no coagulation in the vessel while “II” marks the regimes with thrombi formation. Also shown are two typical values Re 1, Re 2 marking the boundary of the coagulation regime for a given and —the permeability of the vessel wall below which clotting does not occur under any hydrodynamic conditions. Finally, Re = Re indicates the lowest Reynolds number at which the shear stress reaches τ 2 and the local permeability of the vessel reaches . h = 0.6, d = 0.4.

Note the characteristic “tongue” shape of the boundary between regions I and II. This shape implies that the range of Reynolds number values at which blood coagulation occurs is limited both from above and below: for any two thresholds of hydrodynamic activation of blood coagulation exist. ( denotes the value of vessel wall permeability such that if clotting does not occur under any hydrodynamic conditions (see Fig 2).)

The two thresholds can be naively interpreted as follows. Let us denote the two threshold Re values, for a given , as Re 1, Re 2 (see Fig 2). If Re < Re 1, the value of wall shear stress is small and primary activator u does not leak into the blood flow in sufficient quantities to induce coagulation. If Re > Re 2, convective flow washes the coagulation substances away. If Re 1 < Re < Re 2, thrombus formation starts.

Note that the results imply that blood coagulation can be induced in the system via both increasing and decreasing the blood flow intensity — depending on the initial state of the system.

In addition to mapping the boundary between the two coagulation regimes we have explored the scaling behavior of the system in the vicinity of the boundary. Numerical simulations show that the following holds true for the nucleation time T* (given constant Re): Here the nucleation time T* is the non-dimensional time it takes for N to reach in least one point of the vessel in the simulation, is the specific value of , located on the border between zones “I” and “II”, and C 1 is a value independent of . Also, in the vicinity of the left portion of the boundary the following scaling law is valid (given constant ): where Re is the specific value of Re, located on the border between zones “I” and “II”, and C 2 is the value independent on Reynolds number (Re).

In other words, it appears that in the vicinity of the liquid state stability boundary the clot nucleation time grows up to infinity (see Eqs (36) and (37)). This result is reminiscent of the results of the theory of first-order transitions where similar scaling laws exist connecting the extent of supersaturation with the nucleation time [41, 58, 59].

4.2 Effects of stenosis geometry

We further investigated the coagulation regimes by exploring the influence of stenosis geometry on thrombus formation. Fig 3 shows the parametric diagrams in the (h, Re) space for two values of plaque width . As in Fig 2 zone “I” corresponds to stable liquid states (no coagulation) and zone “II” corresponds to thrombus formation.

Fig 3

Parametric diagrams of blood coagulation regimes in the (h, Re) parameter space for two different values of , the non-dimensional plaque diameter.

Region “I” denotes regimes in which there is no coagulation, region “II” represents regimes of macroscopic thrombus formation. (a): , ; (b): , .

The results suggest that for relatively wide stenoses (, Fig 3b) the “thrombogenic” zone “II” in parametric plane is a single-connected domain, while for a narrow stenosis (, Fig 3a) zone “II” is split into two separate areas. There exists a critical value of for which zone “II” consist of two areas connected in 1 point only (not shown in Fig 3). This means that in parametric space the surface dividing zones “I” and “II” is saddle-like.

Interestingly, the results of the mapping also suggest that the range of Reynolds numbers in which macroscopic thrombus formation occurs in general is wider for relatively small plaques (0.3 < h < 0.4) than it is for large plaques (0.8 < h < 0.9), irrespective of .

4.3 Typical scenarios

We also explored the spatial characteristics of blood coagulation in the system under consideration. Numerical simulations suggest that at least two different types of scenarios can be distinguished.

Successive stages of the two different scenarios are shown in Figs 4 and 5. The figures show the distribution of the weight-average number of fibrin monomers in polymer chains N in the vessel. Red color corresponds to clots (), while green shades represent microthrombi with different lengths of polymer chains ().

Fig 4

Successive stages of a solid thrombus formation.

a—thrombus nucleation, b—formation of fibre-like fibrin structure, c—fibre-like structure thickening, d—solid thrombus. Shown are the color maps of N —the weight-average number of monomers in fibrin-polymer in the vessel, with red areas representing regions of fibrin gel formation (). Streamlines are plotted to visualize the flow, and the separatrix, which divides the core of the flow from the recirculation zone, is shown with a dashed line. Parameters used in the simulations are: Re = 130, h = 0.5, , . Note that only a fragment of the vessel closest to the plaque is depicted.

Fig 5

Successive stages of a friable floating fibrin structure formation.

a—thrombus nucleation, b—formation of fibre-like fibrin structure, c—thick and friable floating fibrin structure. Shown are the color maps of N —the weight-average number of monomers in fibrin-polymer in the vessel, with red areas representing regions of fibrin gel formation (). Streamlines are plotted to visualize the flow, and the separatrix is shown with a dashed line. Parameters used in these simulations are: Re = 180, h = 0.6, , . Note that only a fragment of the vessel closest to the plaque is depicted.

Notably, the early stages of coagulation are the same across the explored range of parameters: the nucleation of a macroscopic thrombus occurs near the reattachment point of the recirculation zone that forms behind the stenosis. Then the stage of macroscopic fibre-like structure formation comes (see Figs 4b and 5b). The direction of the growth of the fibre-like structure is determined by the separatrix line, which divides the core of the flow from the recirculation zone.

Further behavior depends on the parameter values. The fibre either successively thickens into a solid thrombus in the recirculation zone (see Fig 4d) or underlays the formation of a floating friable structure (see Fig 5c). Interestingly, the floating structure does not have a sharp border and is characterized by a downstream cloud of mircothrombi.

5 Discussion

In this paper we develop a mathematical model to describe the clinically important case of thrombus formation in the vicinity of an atherosclerotic plaque. The clotting is induced by pro-coagulants leaking into the vessel from the sclerotic tissue through the permeable vessel wall. In contrast to earlier work in the field (see, e.g. [11, 60–67]) the dependence of the vessel wall permeability on shear stress is taken into account. The model sheds some light on the role shear stress plays in thrombogenicity of the vessel when the shear rates are too small to induce platelet activation, as soon as in the present work the inequation was valid at all times.

Within the scope of the model blood coagulation is treated as a phase transition from a liquid to polymerized and gelled state. To simplify the analysis across the different states we develop a phenomenological approach that allows us to formally treat the polymerized and non-polymerized blood in a unified manner. This is done through the closing expressions (19), (22) and (30) for D , b and α , implying that these coefficients can be described as functions of statistical moments of fibrin polymer distribution M 1 and M 2. The expressions (19), (22) and (30) are based on the scaling approach [38] and the perturbation theory technique of asymptotic expansions [68].

Numerical analysis of the introduced phenomenological model reveals the existence of two types of thresholds in hydrodynamic activation of blood coagulation (see Fig 2). One is associated with the increase of wall permeability to procoagulants above the critical level due to increased shear stresses in the vicinity of stenosis. The second threshold is concerned with the intensity with which the procoagulants are washed out from the atherosclerotic region.

The model also sheds some light onto the link between the risk of atherothrombosis and the associated stenosis’ shape and severity. Interestingly, plaques with severe stenosis (h > 0.6) induced clotting only at relatively low Reynolds number (Re < 200), while smaller plaques are characterized by a wider range of Reynolds numbers at which macroscopic thrombus formation occurs (Fig 3).

The result suggests that the most thrombogenic plaques are not necessarily those associated with the most severe occlusion and that the worldwide-accepted indications for stenting vessels should be critically re-assessed [69-71].

Finally, we describe two possible scenarios of blood coagulation in the vicinity of an atherosclerotic plaque. The first scenario corresponds to localized thrombus formation in the recirculation zone behind the stenosis (see Fig 4). The second scenario corresponds to the formation of a friable fibrin polymer structure flattering in the flow behind the plaque (see Fig 5). Which of these scenarios is realized in the system depends upon the blood flow conditions. Low Reynolds numbers usually correspond to solid thrombus formation. In contrast, intensive flows result in the formation of a flattering clot.

Importantly, the friable floating clots are characterized by loose long “tails” of microthrombi downstream of the stenosis (see Fig 5c). One can argue that this microthrombi “dust” may lead to microcirculatory disorders in organs distal from the plaque.

It is worth mentioning that in the case of loose thrombus formation, fibrin generation proceeds inside the thrombus. Thrombin generation continues inside the clot structure. It is natural to suppose that the more intense is thrombin generation the denser microthrombi clouds are formed downstream. Thereupon, an old conclusion may be repeated: to estimate thrombogenic perspectives there is more important to asses endogenous thrombin generation potential than time of clot formation (traditional measured by prothrombin or APTT tests) [72]. Endogenous thrombin potential controls the capability of blood for intensive generation of fibrin microthrombi (irrelative of whether they will be condensed on macrothrombus or washed out downstream) [73].

The floating fibre-like thrombi attached to stenosis are routinely observed in our in vitro arrangement for studying blood coagulation [74, 75]. Interestingly, these frail clots should flatter in the pulsating blood flow and thus be acoustically detectable. Presumably, in the past medical doctors could diagnose them by stethoscoping the sclerosed vessels.

It should be noted that the mathematical description adopted in this work has a few important limitations. In particular, we confine our analysis to stationary flows and neglect the dependence of blood viscosity upon fibrin polymerization. We hope to address these limitations in our future work. Yet, we believe that the results described in this paper will serve as a valuable guide for subsequent work in studying atherothrombosis.