PLoS Computational Biology
Public Library of Science
A three-dimensional phase-field model for multiscale modeling of thrombus biomechanics in blood vessels
Volume: 16, Issue: 4
DOI 10.1371/journal.pcbi.1007709
  • PDF   
  • XML   

Thromboembolism is associated with detachment of small thrombus pieces from the bulk in the blood vessel. These detached pieces, also known as emboli, travel through the blood flow and may block other vessels downstream, e.g., they may plug the deep veins of the leg, groin or arm, leading to venous thromboembolism (VTE). VTE is a significant cause of morbidity and mortality and it affects more than 900,000 people in the United States and result in approximately 100,000 deaths every year. Mechanical interaction between flowing blood and a thrombus is crucial in determining the deformation of the thrombus and the possibility of releasing emboli. In this study, we develop a phase-field model that is capable of describing the structural properties of a thrombus and its biomechanical properties under different blood flow conditions. Moreover, we combine this thrombus model with a particle-based model which simulates the initiation of the thrombus. This combined framework is the first computational study to simulate the development of a thrombus from platelet aggregation to its subsequent viscoelastic responses to various shear flows. Informed by clinical data, this framework can be used to predict the risk of diverse thromboembolic events under physiological and pathological conditions.

Zheng, Yazdani, Li, Humphrey, Karniadakis, and Marsden: A three-dimensional phase-field model for multiscale modeling of thrombus biomechanics in blood vessels


Under physiological conditions, blood clots form at sites of vascular injury to prevent blood loss, but they eventually resolve as the vascular wall heals [1, 2]. The situation can be very different in pathological cases, however, in aortic dissection, which represents a severe injury to the vessel wall that manifests as a delamination that propagates through the media [3], a key clinical question is whether the thrombosis in the false lumen will remain biologically active or be resolved or remodeled as part of the healing process. Other pathological conditions, such as deep vein thrombosis [46], pulmonary embolism [7, 8], and atherothrombosis [9, 10], exhibit excessive and undesirable thrombi within the lumen that can result in partial or complete thrombotic vessel occlusion. Thrombus formation is a multiscale process [1113]: upon a vascular injury or adverse hemodynamics, freely flowing platelets in the blood become activated and form aggregates that cover the thrombogenic area within tens of seconds. Subsequently, a series of biochemical reactions [14, 15], can contribute to the formation of fibrin from the blood-borne fibrinogen. Fibrin monomers form a network that strengthens the platelet aggregates and facilitates thrombus maturation within minutes to hours.

Hemodynamics is essential in the process of thrombus growth [1618]: for example, prolonged residence time of blood cells in regions of weak flow promotes the accumulation of these cells, while high shear rates on the surface of a growing thrombus may alter the morphology of the thrombus and inhibit its expansion into the flow field. If not physiologically degraded, thrombus can begin to undergo a progressive remodeling with replacement of fibrin with collagen fibers [1922], causing a dramatic increase in stiffness, strength, and stability [22, 23]. The mechanical integrity of a thrombus plays an important role in many pathologies. In deep vein thrombosis [46], growth of highly stable thrombus can cause complete vessel occlusion, whereas formation of an unstable thrombus increases the possibility of thromboembolism—a process where pieces of thrombus detach from their original site, transit with blood, and block distal blood vessels—causing life-threatening complications [2426]. Therefore, a comprehensive understanding of the structural constituents of a thrombus and the corresponding mechanical properties is vital to predict thrombus shape and deformation under various hemodynamic conditions and to evaluate the risk of thromboembolism and other pathological consequences.

Extensive experimental studies have increased our understanding of the growth, deformation, and embolization or resolution of a thrombus under physiological and pathological flows [2730]. These studies demonstrated that the deformation of a thrombus and its propensity to generate emboli depend strongly on thrombus mechanical properties and the hemodynamics. Whereas, there is a lack of quantitative information for thrombus deformation and embolization from in vivo studies [31], computational modeling can complement experimental studies by simulating thrombus formation while including sub-processes such as platelet aggregation and coagulation kinetics [3235], as well as thrombus biomechanics in flow [31, 3638]; see also reviews [3942]. There are, however, few studies dedicated to quantifying the integrated process of thrombus initiation and development as well as its deformation and possible embolization under different hemodynamic conditions due to two primary computational challenges: first, thrombus formation is a slow biological process, with the time scale of platelet aggregation on the order of seconds while that of clot remodeling is on the order of days to weeks; second, a constitutive equation describing poro-viscoelastic properties of blood clots with different fibrin concentrations is lacking. Experimental investigations [29, 43, 44] show that a mature thrombus consists of a central core filled with densely packed fibrin and activated platelets. This core is cove with a shell consisting of loosely packed and partially activated platelets, allowing interstitial blood flow through. Hence, a thrombus can be considered as porous medium exhibiting viscoelastic behavior with a fibrous network contributing to both the viscous and the elastic properties [31, 45, 46].

In contrast to previous computational studies, which primarily focused on modeling either platelet aggregation or deformation of an existing thrombus, we present a framework that can capture the integrative process of thrombus formation and its subsequent deformation with possible generation of emboli under dynamic shear flow conditions. We simulate the aggregation of platelets by coupling a spectral/hp element method (SEM) [47] with a force coupling method (FCM) [48], following the strategy in [3, 33]. Once the platelets aggregate becomes stable, we convert the coarse-grained platelet distribution into a three-dimensional (3D) continuum field to estimate the clot volume fraction and continue the simulation of thrombus deformation within a blood flow using a phase-field model (see Fig 1), assuming that further formation of thrombi during the phase-field simulation. The parameters in the phase-field thrombus model are calibrated using existing in vitro experimental data. The connection between FCM particle and phase-field continuum simulations is facilitated by the volume fraction of the thrombus, which is extracted from the FCM simulation as a continuum field; it serves as the initial condition for the subsequent phase-field modeling. As a demonstration of the capability of the proposed framework (FCM + phase-field), we simulate thrombus formation and subsequent interactions with flowing blood in a dilated circular vessel, mimicking a blood vessel with a fusiform aortic aneurysm.

Illustration of the simulation methodology.
Fig 1
(a) A multiscale computational framework to simulate the initiation, growth, and remodeling of a thrombus. (b) Volume fraction Φfcm field for 6 planar cross-sections in the false lumen of a dissecting aortic aneurysm taken along the flow direction from A to F. The platelet Lagrangian distribution in the aggregate obtained from FCM simulation is converted into continuum fields of thrombus volume fraction, which serves as an input for the phase-field method to simulate the interaction between flowing blood and thrombus. Figs. (a) and (b) are adopted from [3].Illustration of the simulation methodology.

This paper is organized as follows: in Results section, we study thrombus biomechanics while varying its permeability and blood flow-induced shear stresses. Subsequently, we calibrate the model parameters that control the permeability and calibrate viscoelastic properties of thrombus using in vitro experimental data. After that, we employ the phase-field model to simulate thrombus deformation in an idealized aneurysm following the initial formation and growth of the platelet aggregate captured by the FCM model. In Discussion section, we discuss the results, limitations, and strengths of the model. In Methods section, we detail the methodology for the proposed 3D framework including the FCM for simulating platelet aggregation and the phase-field method for simulating subsequent interactions between flowing blood and a thrombus.


Effects of thrombus permeability and hemodynamics on thrombus deformation

The relative densities of cross-linked fibrin (fibrous) and platelet aggregates (granular) dictate thrombus permeability. Prior studies show that thrombus permeability and local blood flow-induced shear stresses play important roles in both the deformation of [31] and interstitial flows within thrombus [49]. To examine particular effects of permeability and hemodynamics on thrombus deformation, we first perform 2D and 3D simulations of blood flow passing a semi-spherical porous obstacle that consists of a shell and a core region, (see Fig 2), following the setup used by [31].

Setup of the simulation for modeling thrombus deformation in a 2D rectangular channel (domain Ω = {(x, y)| 0 ≤ x ≤ 6, 0 ≤ y ≤ 2}).
Fig 2
The initial volume fraction distribution of the blood/thrombus is denoted by the phase-field profile, with the core less permeable than the shell. hs and hc are dimensions of the shell and core regions of the thrombus, respectively.Setup of the simulation for modeling thrombus deformation in a 2D rectangular channel (domain Ω = {(x, y)| 0 ≤ x ≤ 6, 0 ≤ y ≤ 2}).

Thrombus deformation in 2D rectangular and 3D cuboid channels

First, we compare simulation results between a 2D thrombus model based on a vector potential ψ (see Eq 1 in Supporting information) and the 3D model based on the deformation gradient tensor F (see Eq 5) implemented in the present work. The initial setup of the 3D simulation is similar to that for the 2D simulation shown in Fig 2 except that we extend the 2D geometry into the z-direction with a thickness of 1 and impose a periodic boundary condition in the z-direction. The geometric and mechanical parameters for both simulations are: the thickness of the shell hs = 0.8, the thickness of the core hc = 0.6, the elastic shear modulus of the thrombus λe = 0.5, the density ratio of thrombus to blood ρ2/ρ1 = 2, the viscosity ratio of thrombus to blood η2/η1 = 1, and the maximum inlet velocity 0.75. As shown in Fig 3, under a steady flow, for both the 2D and 3D simulations, the shell region of the thrombus deforms more toward the downstream direction of the blood flow than the core region since the shell is more permeable to the blood flow than the core. Also, Fig 3(a)–3(d) show that the phase-field and velocity contours obtained from the 3D simulation at cross-section z = 0.5 are similar to those from the 2D simulation, although they are computed using different mathematical formulations. We also plot the time history of the streamwise u-velocity at three different locations as well as the L2 norm of the u -velocity field for 2D and 3D simulations in Fig 4, which suggests that the simulation results for the 2D and 3D models are consistent. These results confirm that the general 3D formulation can reduce to the 2D formulation when the third direction is homogeneous with periodic boundary conditions.

Comparison of simulation results between 2D and 3D models.
Fig 3
The 3D channel has periodic boundaries in the z-direction. The phase-field and streamwise u-velocity are viewed in an xy-plane at time T = 0.96, specifically at z = 0.5 in the 3D simulation. First row: (a) 2D (b) 3D phase-field ϕ, with ϕ = 1 for the flowing blood and ϕ = 0 for the thrombus; second row: (c) 2D (d) 3D u-velocity fields around the deformed thrombus. The shell permeability κs = 10000, the core permeability κc = 1e-4, and the elastic shear modulus for the thrombus λe = 0.5.Comparison of simulation results between 2D and 3D models.
Comparison of velocity magnitudes computed from 2D and 3D rectangular models; see Figs 2 and 3.
Fig 4
(a) Three points A, B, C for 2D and 3D velocity comparison. (b)-(d) Time-history of velocity u at three different axial locations in the channel, and (e) L2 norm of velocity u computed from 2D and 3D models.Comparison of velocity magnitudes computed from 2D and 3D rectangular models; see Figs 2 and 3.

Effect of thrombus permeability on the flow field in a circular vessel

Here we examine thrombi with different permeabilities. The computational domain is a circular vessel with a unit diameter and length of 10. A semi-spherical thrombus composed of a core and shell is initially placed in the lumen (see Fig 5). The flow profile at the inlet of the vessel is set to be parabolic with a maximum value of 0.75 and a zero-Neumann boundary condition is prescribed at the outlet. All relevant model parameters are given in the first row of Table 1. To isolate the effect of the permeability κs, we do not consider elastic energy in the permeability tests (i.e., the thrombus is rigid). The core permeability κc = 1e − 4 is fixed for these simulations. Fig 5 shows the axial velocity field for different values of shell permeability κs of thrombus in xy and xz planes; it is observed that, blood velocity outside the thrombus increases significantly as shell permeability decreases, suggesting that the microstructure of the thrombus plays an important role on its neighboring flow field.

Velocity field in a circular vessel with a thrombus inside.
Fig 5
(a) 3D view of simulation setup. Second row: Axial velocity field w shown in the xz−plane cut at y = 0. Third row: a cross-sectional xy−plane cut through the center of thrombus. Here, hs = 0.5, hc = 0.3 and κc = 1e − 4 for all cases. The permeability of the thrombus varied as: (b) κs = 1, (c) κs = 0.01, (d) κs = 0.005 at time T = 0.9.Velocity field in a circular vessel with a thrombus inside.
Table 1
Parameters used for thrombus deformation simulations in a circular vessel in non-dimensional units.

Effect of hemodynamics on the deformation of a thrombus in a circular vessel

Next, we investigate the effect of hemodynamics on the deformation of a thrombus in a circular vessel. The same simulation setup as the permeability tests is used except that the elastic energy of the thrombus is included in the total energy relation. The parameters used in these simulations are given in the second row of Table 1. Fig 6 shows that as the blood velocity increases, the deformation of the shell region becomes more dramatic whereas the deformation of the core region is not significant due to the increased flow velocity. These observations are consistent with prior experimental observations [29, 50] and numerical studies [31].

Thrombus shapes in a vessel under different flow conditions.
Fig 6
The geometric parameters for thrombus shape are hs = 0.5 and hc = 0.3. Phase-field profile at different blood velocity: (a) w¯=0.25, (b) w¯=1.0, (c) w¯=1.6 shown at time T = 0.91. Phase-field plots shown on xz−plane at y = 0 in the first row and on xy−plane in the second row. w¯: Average velocity in z-direction. ϕ = 1 denotes flowing blood and ϕ = 0 denotes solid thrombus.Thrombus shapes in a vessel under different flow conditions.

Calibration of thrombus permeability and viscoelastic properties

In the previous section, we did not consider the connection between thrombus local volume fraction ϕ and its permeability κ or the connection between ϕ and its elastic properties λe when simulating the deformation of an idealized thrombus. In this section, we define these connections and calibrate the thrombus permeability and elastic properties by comparing our simulation results with existing experimental data.

Calibration of the permeability of thrombus model using fibrin gel data

In this section, we use the fibrin gel data reported in [49] to calibrate the permeability of the thrombus model κ(ϕ ) because of the limited experimental data on thrombi. As shown in Fig 7(a) fluid from a reservoir flows through a fibrin gel in a permeation chamber. The pressure is measured before and after the chamber, and the flow rate is calculated by measuring the displacement of the air-buffer interface in the tubing at regular time intervals. The permeability κ is then calculated using Darcy’s law v = −κΔP/η, where v is the interstitial fluid velocity, η is the viscosity of the percolating fluid, and ΔP is the pressure drop.

Calibration of thrombus permeability κ.
Fig 7
(a) Experimental setup for calibrating thrombus permeability, a typical experimental setup contains a reservoir connected to the permeation chamber via a connecting tube. (b) 3D view of the simulation setup. (c-f) Four initial structures of the thrombus models with VF = 0.33, 0.45, 0.57, 0.66, respectively. (g) Axial-velocity u in the vessel with thrombus VF = 0.66. (h) Phase-field ϕ in the vessel with thrombus VF = 0.66. (i) κ(ϕ ) as a function of the bulk thrombus VF in logarithmic scale. Red line represents the simulation results of the phase-field model. Black dots represent the experimental data as reported in [49]. Blue line represents the relation expressed by Davie’s Equation in [49].Calibration of thrombus permeability κ.

To mimic this experiment, we place a thrombus in the center of a straight cuboid channel filling the channel, as shown in Fig 7(b). The computational domain is Ω = {(x, y, z)| 0 ≤ x ≤ 6, 0 ≤ y ≤ 2, 0 ≤ z ≤ 1}, and no-slip boundary conditions are applied to the upper and lower walls. Blood flow with a parabolic velocity profile is imposed at the inlet of the channel and, a zero-Neumann boundary condition is imposed at the outlet. We use the correlation κ(ϕ)/af2=[16Ψf1.5(1+56Ψf3)]1 [49], also known as Davies’ equation, to calculate the permeability. We vary the initial thrombus volume fraction, and simulate nine cases to calculate the pressure drop of the flow passing through thrombus. Four initial simulation setups are shown in Fig 7(c)–7(f). Parameters are density ratio ρ2/ρ1 = 1, viscosity ratio η2/η1 = 2, length L = 0.76mm, and width H = 0.25mm. The maximum inlet velocity vmax is 1mm/s, which gives a Reynolds number Re = ρ1vmaxH/η1 = 0.02. Furthermore, the average volume fraction (VF) is calculated by ∫Ω[(1 − ϕ)/2]dv/Vtotal . The velocity field and the phase-field contour for the case of VF = 0.66 are plotted in Fig 7(g) and 7(h), respectively, which shows a notable deformation of the thrombus. Our simulation results in Fig 7(i) show that the permeability decreases significantly with the increased VF. We also observe a good agreement between our simulation results and both the experimental data and the empirical fit (Davies’ equation) [49].

Calibration of thrombus elastic shear modulus λe

In our model, the elastic shear modulus λe depends on the local volume fraction of a thrombus, which is equivalent to the phase-field variable ϕ in this study. There are few rheometry experiments that used oscillatory shear deformations to measure the viscoelastic properties of thrombus as a function of fibrin concentration, although they did not cover a wide range of fibrin concentration. In [51], a small oscillatory shear deformation with amplitude of γ0 = 0.01 and frequency f = 0.5 Hz was used to quantify the viscoelasticity of a fibrin gel at low concentrations, which translates to a small thrombus volume fraction. The oscillatory deformation is sufficiently small to probe the linear viscoelastic regime using a Kelvin-Voigt model that decomposes shear stress into elastic and viscous contributions. Furthermore, the storage G′ and loss G″ moduli were quantified for the fibrin gel by the experiment. Both moduli may be combined into a phase angle δ = tan−1(G″/G′) that describes the phase shift between the imposed strain and the measured stress. The two moduli are related by

where f is the frequency of the oscillations, η1 is the fluid viscosity, and λs is the relaxation time defined as λs = η2e with η2 the viscosity of the fibrin gel. Eq 2 describes how the storage and loss moduli are calculated, i.e.,
where Ar is the area of the torque-displacement loop, δ is the phase angle, Ts is the torque, and d is the lateral displacement, following the definition in [52].

Fig 8(a) shows the experimental setup for the oscillatory shear test for the viscoelastic properties of a thrombus sample [53], where the diameter of the rheometer disk is 50 mm and the gap width is 300 μm . Fig 8(b) and 8(c) show the computational domain, which is a 3D channel with Ω = {(x, y, z)| 0 ≤ x ≤ 6, 0 ≤ y ≤ 1, 0 ≤ z ≤ 1}. The upper plate moves at a sinusoidal velocity v = 0.2 sin(1.2πt), with a periodic boundary condition for velocity in the z-direction and a no-slip boundary condition on the lower plate. In our simulations, the characteristic velocity and length scales (required for non-dimensionalization of the equations) are 3.33 × 10−3m/s and 3 × 10−4m, respectively. Zero-Neumann boundaries are imposed for the phase-field variable and for the deformation related variable ψ in 2D, specifically, we impose n ⋅ ∇ϕ = 0, n ⋅ ∇(Δϕ) = 0, and n ⋅ ∇ψ = 0. The thrombus is considered to be in the middle of the channel and three different initial volume fractions VF = 0.3419, 0.5129, and 0.6335 are examined. At the frequency of 0.5 Hz, the velocity of the upper plate is estimated to be 6.67 × 10−4m/s with the fluid kinematic viscosity η1 = 1e-6 Pas and density ρ1 = 1000 kg/m3 . We compute the displacement at the thrombus interface, and the shear stress (including the viscous and the elastic contributions) acting on the surface where the thrombus contacts the upper plate. All parameters for the simulations are listed in Table 2.

Calibration of thrombus elastic shear modulus λe.
Fig 8
(a) Experimental setup for a typical oscillatory shear test. (b) 3D view of the simulation setup. (c) 2D view of the simulation setup. Third and Fourth rows: Stress-displacement loop for different λe. Quantities are all in non-dimensional units. Third row (d-g): VF (volume fraction) = 0.3419; Fourth row (h-k): VF = 0.5129. Points are data and lines are the fits. Displacement and stress are nondimensionalized with the characteristic velocity and length scales 3.33 × 10−3m/s and 3 × 10−4m, respectively.Calibration of thrombus elastic shear modulus λe.
Table 2
Parameters used for oscillatory shear simulations in dimensional units.
Frequency (Hz)Velocity (m/s)ρ2ρ1η2η1λe (Pa)

We construct the stress-displacement loops as shown in Fig 8(d)–8(k). For the same volume fraction shown in each row, the shear stress increases as the elastic shear modulus increases. Moreover, for the same elastic shear modulus on each column, the shear stress increases with increasing volume fraction. Our simulation results in Fig 9 show both G′ and G″ increase as the concentration of fibrinogen increases.

Estimating viscoelastic properties of thrombus using the phase-field model.
Fig 9
(a) Storage (G′ and (b) loss moduli (G″) as a function of fibrinogen concentration cFbg at different λe computed from simulation results.Estimating viscoelastic properties of thrombus using the phase-field model.

Fig 10 shows the results of calibrating the relaxation time for the rheometry experiment. The dotted green line shows the relaxation time computed based on experiments corresponding to initial fibrinogen concentrations of 1, 3 and 6 mg/mL, where G=102.6*log10(cFbg)+log10(10). In addition, we estimate G=101.7*log10(cFbg)+log10(0.7), and the relaxation time is calculated using λs = [G″/(2πf) − η1]/G ′ [51]. The solid lines show the relaxation time from extrapolation and extension from our simulations, where we use a linear fit for the simulation data (only the two points corresponding to VF = 0.3419 and 0.5219) and extrapolate to the concentration range of 1 − 6 mg/mL . By comparing our simulation results with experimental data (see Fig 10), we select λe = 0.44 Pa. As a result, we set the relaxation time as λs(ϕ)=100.7172*log10(cFbg)1.1140 in the following 3D simulation.

Relaxation time as a function of fibrinogen concentration cFbg corresponding to different λe.
Fig 10
Simulation results are represented by solid lines (data are fitted and extrapolated) and experimental result is represented by the green dotted line [51].Relaxation time as a function of fibrinogen concentration cFbg corresponding to different λe.

Simulation of thrombus formation and deformation in an idealized aneurysm

In the case of an abdominal aortic aneurysm (AAA), the low shear rate region inside the most enlarged region combined with a long shear history experienced by the entering platelets can promote activation and aggregation of the platelets and thus thrombus initiation [54, 55]. Thrombus, in turn, causes further enlargement of the AAA and possible thromboembolic events [56, 57]. To demonstrate the capability of the proposed framework in predicting thrombus deformation and possible embolism, in this section, we model the initiation and development of a thrombus in an idealized AAA as well as the subsequent interaction between the thrombus and blood flow.

Modeling thrombus formation with FCM

To avoid solving directly the spatio-temporal “thrombus deposition potential” [54, 55], we compare three realistic platelet deposition sites for illustrative purposes, namely, on the proximal and distal neck and the bottom of the aneurysm, respectively, as shown by the blue patches in Fig 11. Initially, passive platelets (particles) are distributed uniformly inside the aneurysm (Fig 11(a)). When the passive platelets interact with the deposition sites, they become triggered and change to an activated state after a delay time τact = 0.1 − 0.3s [33, 48]. Upon activation, the platelets begin to adhere at deposition sites and to other activated platelets, then modeled as pseudo-platelets (i.e., platelets and associated fibrin) growing in size to an effective radius reff (as illustrated by the schematic in Fig 11(b)). This coarse-grained approach reduces the computational cost by simulating fewer platelets than the physiological concentration, following our prior work [3]. The Reynolds number of the flow is Re = 181.9, which is calculated based on blood viscosity 3.77 × 10−3Pas and blood density 1060 kg/m3 . A generic pulsatile blood flow is prescribed at the inlet for which the average velocity profile within one cardiac cycle is shown in Fig 11(c).

Geometry and boundary condition for simulating thrombus formation in a fusiform aortic aneurysm.
Fig 11
(a),(b) Passive platelet particles are initially distributed throughout the computational domain, including within the aneurysm, while new platelets are inserted into the proximal domain and outgoing platelets are deleted from distal domain. Upper inset: passive platelets as spherical particles of radius rp = 1.5 μm. Pseudo-platelets change their size once they come into contact with the deposition sites and become activated, catalyzing the conversion of blood borne fibrinogen into semi-solid fibrin. Lower inset: activated platelets with increased effective radius reff = 60rp. (c) Pulsatile flow profile imposed at the inlet; a zero pressure condition is prescribed at the outlet.Geometry and boundary condition for simulating thrombus formation in a fusiform aortic aneurysm.

Once the simulation begins, the platelets close to the deposition sites become activated and adhered to the wall. As the simulation continues, increased numbers of activated platelets adhere to the aggregates. We find that platelet aggregates become stable after simulating 15 cardiac cycles. As shown in Fig 12(a), platelet aggregates are mainly located around the platelet deposition sites at the neck and bottom of the aneurysm, consistent with findings from prior computational modeling [54]. Next, as illustrated in Fig 12(b), we convert the Lagrangian distribution of the platelet aggregates to volume fractions that can be used as an initial condition for the subsequent phase-field modeling of interactions between the thrombus and blood flow with possible embolization.

Platelet aggregation modeled by the FCM and the initial conditions for the phase-field model.
Fig 12
(a) Side and top views of the simulated platelet aggregates in the aneurysm after 15 cardiac cycles. (b) The corresponding input for the phase-field simulation, where the top row is the central cross-section of the vessel and the bottom row is the view from the top.Platelet aggregation modeled by the FCM and the initial conditions for the phase-field model.

Multiphase continuum modeling of a thrombus

After the initial phase of platelet aggregation modeled by FCM, we simulate the interaction of blood and thrombus under the same flow conditions using the phase-field model. We assume that platelets will no longer accumulate on the thrombus or wall. The permeability κ(ϕ) calculated by the Davies’ equation and the thrombus viscoelastic properties λe(ϕ) estimated in Results section are used to represent the biomechanical properties of the thrombus.

Fig 13 shows a sequence of snapshots of the thrombus shape under the pulsatile flow in different cross-sectional planes, where the detachment of small pieces of thrombus due to the interaction with blood flow is observed at different axial locations. There is mild pinch-off at both the upstream x ≈ 5cm and downstream x ≈ 8.5 cm of the thrombus, while moderate thrombus contraction occurs proximal to the aneurysm at x ≈ 5.5 cm, 6.5 cm, and 7.5 cm. Thrombus volume fraction VF is low close to the center of the vessel and distal to the aneurysm higher in the deeper regions of the aneurysm. We also perform simulations with a steady inflow boundary condition to investigate the effect of pulsatility on the results. As shown in Fig 6 in Supporting information, our simulations show that a thrombus under steady flow does not release emboli, which is different from the case of pulsatile flow shown in Fig 13. This difference suggests that pulsatility of the blood flow could be important in thromboembolic events secondary to intraluminal thrombus in AAAs.

Physiologic phase-field modeling of thrombus deformation in an idealized aneurysm.
Fig 13
Contours of thrombus volume fraction in the aneurysm shown using the phase-field variable at (a) T = 0s (black lines indicate planar cross-sections at x = 5, 5.5, 6.5, 7.5, and 8.5 cm), (b) T = 0.92s, (c) T = 1.53s, and (d) T = 2.88s. Rows 3-6: phase-field contours plotted at different x locations (e) x = 5cm (f) x = 5.5 cm (g) x = 6.5 cm (h) x = 7.5 cm (i) x = 8.5 cm varying with time. Third row: T = 0.25s; Fourth row: T = 1.17s; Fifth row: T = 1.87s; Sixth row: T = 2.88s.Physiologic phase-field modeling of thrombus deformation in an idealized aneurysm.


In this work, we present a multiscale framework that combines FCM with a phase-field method to model both the initial formation of a thrombus through platelet aggregation and fibrin deposition and subsequent interactions between the flowing blood and the thrombus. Such modeling is challenging because of the multiphysics and multiscale biomechanics. Neither the popular ALE approach for fluid-structure interactions nor the level-set method have been applied to model these aspects of a thrombus due to the difficulty in capturing the complicated rheology at the interface of the thrombus and blood flow. The phase-field method, on the other hand, can capture these effects naturally because it is fully-Eulerian and able to capture the interfacial rheology and the shear stress on the thrombus surface. A recent numerical study [31] simulated thrombus deformation with a 2D phase-field model, but did not consider how the viscoelastic properties of the thrombus affect its deformation. In this work, we implemented a 3D phase-field model to investigate effects of the permeability and viscoelastic properties of a thrombus on its deformation in a flow field. Our simulations show that the permeability and the elastic shear modulus are essential in determining the responses of the thrombus. Our simulations of the rheology and flow-induced shape changes of the thrombus were consistent with results in a prior computational study [31] and experimental observations [29, 50], suggesting the potential of this approach in predicting the risk of thrombotic and thromboembolic events under physiologic and pathological conditions. Indeed, even though we used the volume fraction of a thrombus to define the phase-field variable ϕ , which leads to a finite thickness interface in our simulations that may compromise the original phase-field assumption of a thin smooth transition layer, our assumption simplified the model and yielded a reasonable congruence with available additional data [49].

A unique feature of the proposed framework is the combination of a phase-field method with FCM by converting the volume fraction of platelet aggregation into a phase-field variable, which enables consistent modeling of thrombus initiation, development, and deformation under blood flow. When simulating thrombus formation in an idealized aortic aneurysm, our FCM simulations suggest that the thrombus does not grow uniformly; moreover, the volume fractions of the blood and thrombus also vary spatially. Portions of thrombus with lower volume fractions VF have a larger propensity to deform toward the downstream direction of the blood flow. Small pieces even can detach from the bulk thrombus and flow downstream, similar to the process of thrombus embolization.

We note that thrombus can form via different mechanisms that lead to different compositions and stability [58, 59]. The composition of thrombus in AAAs is even more complex. Intraluminal thrombus may contain layers with different compositions, representing blood cells that were captured over multiple time scales [60]. Complexities similarly arise in intramural thrombus in the aortic dissection, with some regions demonstrating remodeling of fibrin-rich thrombus into collagen-rich fibrotic tissue [21]. In this work, we considered the formation of a thrombus in an idealized AAA, in which the low shear rate region inside the aneurysm combines with a long shear history experienced by the platelets to promote local platelet deposition [54, 55]. The aim of our simulations was thus not to predict the differential formation of a thrombus in aneurysms, but rather to demonstrate the capability of the proposed framework to simulate the formation of thrombus that initiated via platelet activation and deforms due to both external and interstitial blood flow. Finally, we calibrated both the permeability and viscoelasticity of the model using data from fibrin gels given the lack of data on the viscoelastic properties of platelet-rich thrombus. These parameter values can be refined when new experimental data become available.

Taken together, we developed a 3D phase-field model for simulating blood-thrombus interactions. Calibration of model parameters was achieved by comparing the simulation results with existing in vitro experimental data for the poro-viscoelastic behavior. The proposed 3D phase-field model can be combined with FCM to simulate the process from platelet activation and deposition to subsequent deformation while exposed to steady or pulsatile blood flow. Guided by patient-specific clinical data, such as lesion geometry and the local blood flow rate, this multiscale framework has the potential to predict the risk of thrombotic-embolic events, which are responsible for significant morbidity and mortality. Furthermore, this framework can be further expanded by combining it with models for thrombus remodeling [6163] to improve our understanding of how changes of biomechanical properties of a thrombus during its maturation affect its role in both physiology and pathophysiology.


Platelet transport and aggregation

We simulate the transport of platelets within flowing blood and their aggregation on deposition sites using FCM integrated with SEM [33]. FCM is used to describe the motion of platelets and their (bi-directional) interactions with the flowing blood, which is determined using SEM to solve the flow field on a fixed Eulerian grid. This multiscale approach has been successfully used in modeling platelet aggregation in venules [48], stenotic channels [33], and intramural aortic dissections [3]. In our simulations, platelets exhibit three different states, namely passive, triggered, or activated. In passive or triggered states, platelets have their physiological radius rp = 1.5 μm and are non-adhesive. Once platelets in the passive state interact with activated platelets or a deposition site, they switch to a triggered state and then activate after a delay time τact = 0.1 − 0.3 s [48]. After activation, pseudo-platelets of platelets embedded in a local fibrin network) grow to an effective radius reff ∼ 60rp . This approach allows us to use fewer platelets than the physiological concentration when modeling thrombus formation in large computational domains [3].

The margination of platelets due to collisions with blood cells in flowing blood has been studied extensively, both experimentally and theoretically, for straight channels and idealized vessels [6468]. To avoid the added complexity of simulating blood cells explicitly in the blood flow, we account for platelet margination by using a master profile for platelet distribution within the aorta, quantified in in [64] when platelet particles are inserted into the inlet of a vessel. The interactions amongst the platelets as well as between the platelets and the deposition sites are then described by a Morse potential (attractive interactions between platelets) and an exponential repulsion potential (exclusive effects of the platelet particles) (see S5 Text). The adhered platelet particles are considered to be ‘stationary’ or part of the thrombus when their moving distance within one cardiac cycle is less than 1/100 of their diameter. The interaction forces are shear-rate dependent; the force values were calibrated using data from four independent studies, including two in vivo [3, 27] and two in vitro [69, 70] experiments, which measured platelet aggregation at different shear rates. More detailed information of the FCM model can be found in S5 Text and in [3, 33].

To compute the volume fraction of the formed thrombus VF (i.e., the local volume occupied by the thrombus divided by the volume of blood), we first use a continuum representation of the adhered platelets from FCM simulations to evaluate the local volume fraction of pseudo-platelet particles as Ψfcm(x,t)=n=1NVpnΔ(xYn(t)), where Vpn is the volume of each pseudo-platelet, x is the position of the background Eulerian grid, and Yn is the position of each FCM particle. Knowing the local concentration of fibrinogen and platelets and using the computed Ψfcm field, we estimate the thrombus volume fraction as VF = [Ψf(cFbg)+ Ψp(cplat)] Ψfcm, where Ψf is the volume fraction of fibrin and Ψp is the volume fraction contribution of platelets [3, 49]. Ψf can be computed by an empirical relation Ψf = cFbg/[ρFbg(0.015 log(cFbg) + 0.13)], where cFbg (mg/mL) is the concentration of fibrinogen and ρFbg = 1.4 g/mL is the density of a single fibrinogen molecule [49]. Ψp can be evaluated directly based on each platelet’s volume and the local platelet number density [49].

Modeling clot deformation using the phase-field method

Various computational methods have been implemented to model interactions between fluids and viscoelastic solids, including the arbitrary Lagrangian Eulerian formulation (ALE), and the immersed boundary method as well as, level set, volume-of-fluid, and phase-field methods [7173]. An advantage of phase-field approaches in modeling multiphase systems is that they describe the system by characterizing the free energy. In this way, different physical/biological phenomena can be considered through an appropriate modification of a unified set of governing equations for the free energy formulated over the entire computational domain [71]. The phase-field method has been used to address interactions between flowing blood and thrombus as well as between the flowing blood and a biofilm [31, 38]. These studies were limited, however, to 2D simulations for predefined thrombus and biofilm shapes. Other studies elucidated the roles of thrombus permeability on its deformation and embolization [31], but did not include effects of spatial variability of the volume fraction of the thrombus. In the current study, we extend the prior work by performing new 3D simulations of interactions between the thrombus and flowing blood, while considering both the spatially variable permeability and viscoelastic properties of a thrombus. We believe that this approach provides more biologically realistic and physically consistent results.

The key idea of the phase-field approach is to use a phase-field variable ϕ to describe a two-phase system, with the interface between two different phases modeled by a thin smooth transition layer [71]. In our simulation of a vascular lumen Ω containing a thrombus, ϕ = 1 represents the blood, 0 < ϕ < 1 represents a mixture of blood and thrombus, and ϕ = 0 represents thrombus only. The phase-field equation is derived by minimizing the total energy Etot of the blood-thrombus system, which is the sum of the kinetic energy, cohesive energy of the mixture, and the elastic energy of the fibrin network and platelets in a thrombus. The cohesive energy Ecohesive of the mixture is adopted from Cahn & Hilliard [74]. We describe the viscoelastic fibrin network in thrombus using a Kelvin-Voigt type model where a neoHookean relation is used to describe the elastic behavior and a dashpot is used to describe the viscous behavior [51, 75].

The total free energy for the blood-thrombus system can thus be written by adding the cohesive and elastic energies, respectively, as

where λ is the mixing energy density, λe is the elastic shear modulus, h is a characteristic length scale of the interface thickness, tr is the trace of a matrix and I is the identity matrix. Relative to an Eulerian reference frame, the deformation gradient tensor F(x,t)=x^/X satisfies the following evolution equation [75] (it can be shown that ∇ ⋅ F = 0 if ∇ ⋅ F0 = 0 at time t = 0):
where u is the velocity, x^ and X are the current and original locations of a material point. In prior 2D models [31, 38], an auxiliary vector ψ was introduced to calculate the deformation gradient tensor F through F = ∇ × ψ and the corresponding equation for ψ is ψt+u·ψ=0. To improve the computational efficiency in 2D, we solve for ψ = [ψ1, ψ2] instead of the four components of Eq 1 in Supporting information.

The governing equations for the 3D phase-field model are summarized as follows

where u, p(ϕ), ρ(ϕ), η(ϕ), κ(ϕ) are the velocity, pressure, mass density, dynamic viscosity and permeability, respectively, with the phase-field variable ϕ ∈ [0, 1]. Furthermore, γ is a constant interfacial mobility of the phase-field, τ is a relaxation parameter, and g1(ϕ) is the derivative of the double well potential ϕ2(ϕ − 1)2/2h2 . Eq 5(a) and 5(b) account for linear momentum balance and mass balance, respectively. Eq 5(c) calculates the deformation gradient tensor F and Eq 5(d) computes the phase-field variable ϕ. The free energy potential μ1 for the blood-thrombus system is expressed by Eq 5(e), which is derived by taking the derivative of the total free energy Etotal with respect to ϕ.

Following Dong et al . [71], our model and algorithms have several attractive features. The governing equations are solved on the entire domain with a fixed Eulerian grid and all key variables, such as the velocity u, the phase-field ϕ and the deformation gradient tensor F are decoupled. The resulting coefficient mass and stiffness matrices are constant and thus can be computed at the beginning of each simulation to increase the computational efficiency. Furthermore, the 4th order phase-field equation (i.e ., Cahn-Hilliard equation [74]) is reformulated into two second-order Helmholtz equations, which can be solved by the spectral or finite element method. We use an entropy viscosity method to stabilize the hyperbolic Eq 5(c) for F (or vector ψ in 2D) [76].

The surface tension parameter σ is related to the mixing energy density λ in Eq 5 using λ=322σh. We tested different values of the surface tension (λ changes accordingly) with other parameters fixed, and found that the effect of σ on the mechanical properties of the thrombus is negligible. Therefore, we assume a weak surface tension force in the momentum Eq 5(a) compared to the viscous and elastic forces.



RW Colman. Hemostasis and thrombosis: basic principles and clinical practice. Lippincott Williams & Wilkins; 2006.


JC Chapin, KA Hajjar. . Fibrinolysis and the control of blood coagulation. Blood Reviews. 2015;29(1):, pp.17–24. , doi: 10.1016/j.blre.2014.09.003


A Yazdani, H Li, MR Bersi, P Achille, J Insley, JD Humphrey, et al. Data-driven modeling of hemodynamics and its role on thrombus size and shape in aortic dissections. Scientific Reports. 2018;8(1):, pp.2515, doi: 10.1038/s41598-018-20603-x


EE Weinmann, EW Salzman. . Deep-vein thrombosis. New England Journal of Medicine. 1994;331(24):, pp.1630–1641. , doi: 10.1056/NEJM199412153312407


AW Lensing, P Prandoni, MH Prins, H Büller. . Deep-vein thrombosis. The Lancet. 1999;353(9151):, pp.479–485.


PA Kyrle, S Eichinger. . Deep vein thrombosis. The Lancet. 2005;365(9465):, pp.1163–1174.


S Sevitt, N Gallagher. . Venous thrombosis and pulmonary embolism. A clinico-pathological study in inju and burned patients. British Journal of Surgery. 1961;48(211):, pp.475–489. , doi: 10.1002/bjs.18004821103


EF Mammen. . Pathogenesis of venous thrombosis. Chest. 1992;102(6):, pp.640S–644S. , doi: 10.1378/chest.102.6_supplement.640s


M Yasaka, T Yamaguchi, M Shichiri. . Distribution of atherosclerosis and risk factors in atherothrombotic occlusion. Stroke. 1993;24(2):, pp.206–211. , doi: 10.1161/01.str.24.2.206


S Endo, N Kuwayama, Y Hirashima, T Akai, M Nishijima, A Takaku. . Results of urgent thrombolysis in patients with major stroke and atherothrombotic occlusion of the cervical internal carotid artery. American Journal of Neuroradiology. 1998;19(6):, pp.1169–1175.


Z Xu, N Chen, MM Kamocka, ED Rosen, M Alber. . A multiscale model of thrombus development. Journal of the Royal Society Interface. 2008;5(24):, pp.705–722.


A Belyaev, J Dunster, J Gibbins, M Panteleev, V Volpert. . Modelling thrombosis in silico: frontiers, challenges, unresolved problems and milestones. Physics of Life Reviews. 2018.


DA Fedosov. . Hemostasis is a highly multiscale process: Comment on Modeling thrombosis in silico: Frontiers, challenges, unresolved problems and milestones by AV Belyaev et al. Physics of Life Reviews. 2018;. , doi: 10.1016/j.plrev.2018.06.017


N Mackman, RE Tilley, NS Key. . Role of the extrinsic pathway of blood coagulation in hemostasis and thrombosis. Arteriosclerosis, Thrombosis, and Vascular biology. 2007;27(8):, pp.1687–1693. , doi: 10.1161/ATVBAHA.107.141911


MS Chatterjee, WS Denney, H Jing, SL Diamond. . Systems biology of coagulation initiation: kinetics of thrombin generation in resting and activated human blood. PLoS Computational Biology. 2010;6(9):, pp.e1000950, doi: 10.1371/journal.pcbi.1000950


WS Nesbitt, E Westein, FJ Tovar-Lopez, E Tolouei, A Mitchell, J Fu, et al. A shear gradient–dependent platelet aggregation mechanism drives thrombus formation. Nature Medicine. 2009;15(6):, pp.665, doi: 10.1038/nm.1955


DL Bark, AN Para, DN Ku. . Correlation of thrombosis growth rate to pathological wall shear rate during platelet accumulation. Biotechnology and Bioengineering. 2012;109(10):, pp.2642–2650. , doi: 10.1002/bit.24537


DM Wootton, DN Ku. . Fluid mechanics of vascular systems, diseases, and thrombosis. Annual Review of Biomedical Engineering. 1999;1(1):, pp.299–329. , doi: 10.1146/annurev.bioeng.1.1.299


V Fineschi, E Turillazzi, M Neri, C Pomara, I Riezzo. . Histological age determination of venous thrombosis: a neglected forensic task in fatal pulmonary thrombo-embolism. Forensic Science International. 2009;186(1-3):, pp.22–28. , doi: 10.1016/j.forsciint.2009.01.006


M Nosaka, Y Ishida, A Kimura, T Kondo. . Time-dependent organic changes of intravenous thrombi in stasis-induced deep vein thrombosis model and its application to thrombus age determination. Forensic Science International. 2010;195(1-3):, pp.143–147. , doi: 10.1016/j.forsciint.2009.12.008


A Schriefl, M Collins, D Pierce, GA Holzapfel, L Niklason, J Humphrey. . Remodeling of intramural thrombus and collagen in an Ang-II infusion ApoE-/- model of dissecting aortic aneurysms. Thrombosis Research. 2012;130(3):, pp.e139–e146. , doi: 10.1016/j.thromres.2012.04.009


YU Lee, A Lee, J Humphrey, M Rausch. . Histological and biomechanical changes in a mouse model of venous thrombus remodeling. Biorheology. 2015;52(3):, pp.235–245. , doi: 10.3233/BIR-15058


S Emelianov, X Chen, M O’donnell, B Knipp, D Myers, T Wakefield, et al. Triplex ultrasound: elasticity imaging to age deep venous thrombosis. Ultrasound in Medicine & Biology. 2002;28(6):, pp.757–767.


WH Geerts, JA Heit, GP Clagett, GF Pineo, CW Colwell, FA Anderson, et al. Prevention of venous thromboembolism. Chest. 2001;119(1):, pp.132S–175S. , doi: 10.1378/chest.119.1_suppl.132s


RH White. . The epidemiology of venous thromboembolism. Circulation. 2003;107(23_suppl_1):, pp.I–4.


SZ Goldhaber, H Bounameaux. . Pulmonary embolism and deep vein thrombosis. The Lancet. 2012;379(9828):, pp.1835–1846.


N Begent, G Born. . Growth rate in vivo of platelet thrombi, produced by iontophoresis of ADP, as a function of mean blood flow velocity. Nature. 1970;227:, pp.926–930. , doi: 10.1038/227926a0


B Furie, BC Furie. . Thrombus formation in vivo. Journal of Clinical Investigation. 2005;115(12):, pp.3355, doi: 10.1172/JCI26987


JD Welsh, TJ Stalker, R Voronov, RW Muthard, M Tomaiuolo, SL Diamond, et al. A systems approach to hemostasis: 1. The interdependence of thrombus architecture and agonist movements in the gaps between platelets. Blood. 2014;124:, pp.1808–1815. , doi: 10.1182/blood-2014-01-550335


OV Kim, Z Xu, ED Rosen, MS Alber. . Fibrin networks regulate protein transport during thrombus development. PLoS Computational Biology. 2013;9(6):, pp.e1003095, doi: 10.1371/journal.pcbi.1003095


S Xu, Z Xu, OV Kim, RI Litvinov, JW Weisel, M Alber. . Model predictions of deformation, embolization and permeability of partially obstructive blood clots under variable shear flow. Journal of the Royal Society Interface. 2017;14(136):, pp.20170441.


A Tosenberger, F Ataullakhanov, N Bessonov, M Panteleev, A Tokarev, V Volpert. . Modelling of platelet–fibrin clot formation in flow with a DPD–PDE method. Journal of Mathematical Biology. 2016;72(3):, pp.649–681. , doi: 10.1007/s00285-015-0891-2


A Yazdani, H Li, JD Humphrey, GE Karniadakis. . A general shear-dependent model for thrombus formation. PLoS Computational Biology. 2017;13(1):, pp.e1005291, doi: 10.1371/journal.pcbi.1005291


H Lei, DA Fedosov, GE Karniadakis. . Time-dependent and outflow boundary conditions for dissipative particle dynamics. Journal of Computational Physics. 2011;230(10):, pp.3765–3779. , doi: 10.1016/


W Wang, TG Diacovo, J Chen, JB Freund, MR King. . Simulation of platelet, thrombus and erythrocyte hydrodynamic interactions in a 3D arteriole with in vivo comparison. PLoS One. 2013;8(10):, pp.e76949, doi: 10.1371/journal.pone.0076949


W Wang, JP Lindsey, J Chen, TG Diacovo, MR King. . Analysis of early thrombus dynamics in a humanized mouse laser injury model. Biorheology. 2014;51(1):, pp.3–14. , doi: 10.3233/BIR-130648


M Tomaiuolo, TJ Stalker, JD Welsh, SL Diamond, T Sinno, LF Brass. . A systems approach to hemostasis: 2. Computational analysis of molecular transport in the thrombus microenvironment. Blood. 2014;124:, pp.1816–1823. , doi: 10.1182/blood-2014-01-550343


G Tierra, JP Pavissich, R Nerenberg, Z Xu, MS Alber. . Multicomponent model of deformation and detachment of a biofilm under fluid flow. Journal of The Royal Society Interface. 2015;12(106):, pp.20150045.


K Leiderman, A Fogelson. . An overview of mathematical modeling of thrombus formation under flow. Thrombosis Research. 2014;133:, pp.S12–S14. , doi: 10.1016/j.thromres.2014.03.005


SL Diamond. . Systems analysis of thrombus formation. Circulation Research. 2016;118(9):, pp.1348–1362. , doi: 10.1161/CIRCRESAHA.115.306824


MH Flamm, TV Colace, MS Chatterjee, H Jing, S Zhou, D Jaeger, et al. Multiscale prediction of patient-specific platelet function under flow. Blood. 2012;120(1):, pp.190–198. , doi: 10.1182/blood-2011-10-388140


Z Xu, O Kim, M Kamocka, ED Rosen, M Alber. . Multiscale models of thrombogenesis. Wiley Interdisciplinary Reviews: Systems Biology and Medicine. 2012;4(3):, pp.237–246. , doi: 10.1002/wsbm.1160


DB Cines, T Lebedeva, C Nagaswami, V Hayes, W Massefski, RI Litvinov, et al. Clot contraction: compression of erythrocytes into tightly packed polyhedra and distribution of platelets and fibrin. Blood. 2014;123(10):, pp.1596–1603. , doi: 10.1182/blood-2013-08-523860


TJ Stalker, JD Welsh, M Tomaiuolo, J Wu, TV Colace, SL Diamond, et al. A systems approach to hemostasis: 3. Thrombus consolidation regulates intrathrombus solute transport and local thrombin activity. Blood. 2014;124:, pp.1824–1831. , doi: 10.1182/blood-2014-01-550319


GB Thurston. . Viscoelasticity of human blood. Biophysical Journal. 1972;12(9):, pp.1205–1217. , doi: 10.1016/S0006-3495(72)86156-3


TH van Kempen, WP Donders, FN van de Vosse, GW Peters. . A constitutive model for developing blood clots with various compositions and their nonlinear viscoelastic behavior. Biomechanics and Modeling in Mechanobiology. 2016;15(2):, pp.279–291. , doi: 10.1007/s10237-015-0686-9


G Karniadakis, S Sherwin. Spectral/hp element methods for computational fluid dynamics. 3rd edOxford University Press; 2013.


IV Pivkin, PD Richardson, G Karniadakis. . Blood flow velocity effects and role of activation delay time on growth and form of platelet thrombi. Proceedings of the National Academy of Sciences. 2006;103(46):, pp.17164–17169.


A Wufsus, N Macera, K Neeves. . The hydraulic permeability of blood clots as a function of fibrin and platelet density. Biophysical Journal. 2013;104(8):, pp.1812–1823. , doi: 10.1016/j.bpj.2013.02.055


TJ Stalker, EA Traxler, J Wu, KM Wannemacher, SL Cermignano, R Voronov, et al. Hierarchical organization in the hemostatic response and its relationship to the platelet signaling network. Blood. 2013;123:, pp.1875–1885.


TH van Kempen, AC Bogaerds, GW Peters, FN van de Vosse. . A constitutive model for a maturing fibrin network. Biophysical Journal. 2014;107(2):, pp.504–513. , doi: 10.1016/j.bpj.2014.05.035


M Puig-de Morales-Marinkovic, KT Turner, JP Butler, JJ Fberg, S Suresh. . Viscoelasticity of the human blood cell. American Journal of Physiology-Cell Physiology. 2007;293(2):, pp.C597–C605. , doi: 10.1152/ajpcell.00562.2006


P Sousa, J Carneiro, R Vaz, A Cerejo, F Pinho, M Alves, et al. Shear viscosity and nonlinear behavior of whole blood under large amplitude oscillatory shear. Biorheology. 2013;50(5-6):, pp.269–282. , doi: 10.3233/BIR-130643


L Mountrakis, E Lorenz, A Hoekstra. . Where do the platelets go? A simulation study of fully resolved blood flow through aneurysmal vessels. Interface Focus. 2013;3(2):, pp.20120089, doi: 10.1098/rsfs.2012.0089


P Di Achille, G Tellides, J Humphrey. . Hemodynamics-driven deposition of intraluminal thrombus in abdominal aortic aneurysms. International Journal for Numerical Methods in Biomedical Engineering. 2017;33(5):, pp.e2828.


SS Hans, O Jareunpoon, M Balasubramaniam, GB Zelenock. . Size and location of thrombus in intact and ruptured abdominal aortic aneurysms. Journal of Vascular Surgery. 2005;41(4):, pp.584–588. , doi: 10.1016/j.jvs.2005.01.004


JS Wilson, L Virag, P Di Achille, I Karšaj, JD Humphrey. . Biochemomechanics of intraluminal thrombus in abdominal aortic aneurysms. Journal of Biomechanical Engineering. 2013;135(2):, pp.021011, doi: 10.1115/1.4023437


Sen-Chowdhry S, Gordon RJ. One weapon, two blows in the war against the thrombus; 2014.


Dongjune KimZ LChristopher Bresette, DN Ku. . Occlusive thrombosis in arteries. APL Bioeng. 2019;3(3):, pp.041502, doi: 10.1063/1.5115554


J Tong, GA Holzapfel. . Structure, mechanics, and histology of intraluminal thrombi in abdominal aortic aneurysms. Annals of Biomedical Engineering. 2015;43(7):, pp.1488–1501. , doi: 10.1007/s10439-015-1332-5


J Humphrey, K Rajagopal. . A constrained mixture model for growth and remodeling of soft tissues. Mathematical Models and Methods in Applied Sciences. 2002;12(03):, pp.407–430.


I Karšaj, JD Humphrey. . A mathematical model of evolving mechanical properties of intraluminal thrombus. Biorheology. 2009;46(6):, pp.509–527. , doi: 10.3233/BIR-2009-0556


JD Humphrey, GA Holzapfel. . Mechanics, mechanobiology, and modeling of human abdominal aorta and aneurysms. Journal of Biomechanics. 2012;45(5):, pp.805–814. , doi: 10.1016/j.jbiomech.2011.11.021


P Aarts, S Van Den Broek, GW Prins, G Kuiken, JJ Sixma, RM Heethaar. . Blood platelets are concentrated near the wall and blood cells, in the center in flowing blood. Arteriosclerosis: An Official Journal of the American Heart Association, Inc. 1988;8(6):, pp.819–824.


H Zhao, ES Shaqfeh. . Shear-induced platelet margination in a microchannel. Physical Review E. 2011;83(6):, pp.061924.


DA Reasor Jr, M Mehrabadi, DN Ku, CK Aidun. . Determination of critical parameters in platelet margination. Annals of Biomedical Engineering. 2013;41(2):, pp.238–249. , doi: 10.1007/s10439-012-0648-7


A Jordan, T David, S Homer-Vanniasinkam, A Graham, P Walker. . The effects of margination and cell augmented platelet diffusivity on platelet adhesion in complex flow. Biorheology. 2004;41(5):, pp.641–653.


R Zhao, MV Kameneva, JF Antaki. . Investigation of platelet margination phenomena at elevated shear stress. Biorheology. 2007;44(3):, pp.161–177.


E Westein, AD van der Meer, MJ Kuijpers, JP Frimat, A van den Berg, JW Heemskerk. . Atherosclerotic geometries exacerbate pathological thrombus formation poststenosis in a von Willebrand factor-dependent manner. Proceedings of the National Academy of Sciences. 2013;110(4):, pp.1357–1362.


M Li, DN Ku, CR Forest. . Microfluidic system for simultaneous optical measurement of platelet aggregation at multiple shear rates in whole blood. Lab on a Chip. 2012;12(7):, pp.1355–1362. , doi: 10.1039/c2lc21145a


S Dong, J Shen. . A time-stepping scheme involving constant coefficient matrices for phase-field simulations of two-phase incompressible flows with large density ratios. Journal of Computational Physics. 2012;231(17):, pp.5788–5804.


X Zheng, H Babaee, S Dong, C Chryssostomidis, G Karniadakis. . A phase-field method for 3D simulation of two-phase heat transfer. International Journal of Heat and Mass Transfer. 2015;82:, pp.282–298.


X Zheng, G Karniadakis. . A phase-field/ALE method for simulating fluid–structure interactions in two-phase flow. Computer Methods in Applied Mechanics and Engineering. 2016;309:, pp.19–40.


JW Cahn, JE Hilliard. . Free energy of a nonuniform system. I. Interfacial free energy. The Journal of Chemical Physics. 1958;28(2):, pp.258–267.


FH Lin, C Liu, P Zhang. . On hydrodynamics of viscoelastic fluids. Communications on Pure and Applied Mathematics. 2005;58(11):, pp.1437–1471.


JL Guermond, R Pasquetti, B Popov. . Entropy viscosity method for nonlinear conservation laws. Journal of Computational Physics. 2011;230(11):, pp.4248–4267.

18 Oct 2019

Dear Dr Karniadakis,

Thank you very much for submitting your manuscript 'A Three-dimensional Phase-field Model for Multiscale Modeling of Thrombus Biomechanics in Blood Vessels' for review by PLOS Computational Biology. Your manuscript has been fully evaluated by the PLOS Computational Biology editorial team and in this case also by independent peer reviewers. The reviewers appreciated the attention to an important problem, but raised some substantial concerns about the manuscript as it currently stands. While your manuscript cannot be accepted in its present form, we are willing to consider a revised version in which the issues raised by the reviewers have been adequately addressed. We cannot, of course, promise publication at that time.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

Your revisions should address the specific points made by each reviewer. Please return the revised version within the next 60 days. If you anticipate any delay in its return, we ask that you let us know the expected resubmission date by email at Revised manuscripts received beyond 60 days may require evaluation and peer review similar to that applied to newly submitted manuscripts.

In addition, when you are ready to resubmit, please be prepared to provide the following:

(1) A detailed list of your responses to the review comments and the changes you have made in the manuscript. We require a file of this nature before your manuscript is passed back to the editors.

(2) A copy of your manuscript with the changes highlighted (encouraged). We encourage authors, if possible to show clearly where changes have been made to their manuscript e.g. by highlighting text.

(3) A striking still image to accompany your article (optional). If the image is judged to be suitable by the editors, it may be featured on our website and might be chosen as the issue image for that month. These square, high-quality images should be accompanied by a short caption. Please note as well that there should be no copyright restrictions on the use of the image, so that it can be published under the Open-Access license and be subject only to appropriate attribution.

Before you resubmit your manuscript, please consult our Submission Checklist to ensure your manuscript is formatted correctly for PLOS Computational Biology: Some key points to remember are:

- Figures uploaded separately as TIFF or EPS files (if you wish, your figures may remain in your main manuscript file in addition).

- Supporting Information uploaded as separate files, titled Dataset, Figure, Table, Text, Protocol, Audio, or Video.

- Funding information in the 'Financial Disclosure' box in the online system.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see here

We are sorry that we cannot be more positive about your manuscript at this stage, but if you have any concerns or questions, please do not hesitate to contact us.


Alison Marsden

Associate Editor

PLOS Computational Biology

Jason Haugh

Deputy Editor

PLOS Computational Biology

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact immediately:


Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Attached

Reviewer #2: The authors present a framework to integrate the Force Coupling Method (as a means to model platelet transport in fluid) and a Phase-field modeling technique (as a means to model a thrombus as a poro-elastic structure) to investigate thrombus formation and deformation. To calibrate parameters within their phase-field model, they compare to data from in vitro experiments of clot permeability and rheology. The goal is to use this framework in the future to help determine risk of thrombotic events like embolization in patient-specific geometries. The framework presented is an extension from their recent paper and adds phase field models developed by a different group. The ideas in this paper are very interesting and suited for the journal, but there are many details missing related to the platelet aggregation portion of the model and so that section of the paper is difficult to fully and fairly review. Further, the calibrations for permeability and rheology were based on studies with fibrin rich clots instead of platelet rich clots, which leaves doubt about the relevance for cell-rich clots. See detailed comments below.

Major comments:

1. It is unclear to the reviewer how the platelets (modeled by FCM) are adhering and aggregating, and how their “growth” is considered dynamically in time and in the fluid. In section 2.1 (platelet transport and aggregation) there is no description of how platelets aggregate or adhere to the wall. At what point do they begin to “grow”? How is the spatial distribution of the new larger platelets determined from the distribution of the smaller platelets, especially near the wall? Are the larger platelets stationary? Finally, during the growth stage, how do the authors update the fluid equations surrounding the platelets? Due to lack of details provided on the thrombus formation process makes difficult to fully assess and properly review the study.

2. The usefulness of going from FCM to a PF model is the heterogeneity of the clot – then the effects of heterogeneity on the mechanics of embolization can be studied; the structure of the thrombus that initially forms is clearly extremely important. Without describing how the platelets are aggregating and which ones have actually aggregated, it is not clear that the phase field model is a correct representation of the clot itself, but rather just a mass of platelets that are flowing through the fluid, in a non-aggregated state. It is possible that this is distinguished somehow but it is not made clear to the reviewer.

3. It is unclear to the reviewer where fibrinogen fits into this model. In the methods section, the concentration of fibrinogen is mentioned, and then again in the rheological study section it is used, but there is no evolution equation for fibrinogen for it to exist or interact in the thrombus or at the wall for adhesion.

4. The function g1(phi) represents a double well potential when phi is in [-1,1] with minima at +/-1; what is the rationale for choosing phi in [0,1] in this model?

5. The calibration of the permeability of the thrombi were based on experiments with fibrin clots and the empirical formula of Davies is based on fibrous material. The referenced paper (Wufsus et al.) included studies on platelet rich clots where the permeabilities were fit well with a different empirical formula (Ethier). This suggests that the permeabilities of the thrombi modeled in the current study could be significantly underestimated. The authors should instead fit to the platelet-rich clot data or justify their choice for using the fibrin clots instead of platelet-rich ones.

6. Similarly, the calibrations for the shear modulus are based on rheology experiments performed with fibrin gels while the clots in this study are assumed to be primarily platelets. The section describing how the volume fraction is calculated based on fibrinogen concentration is very confusing. It is not clear what the parameters are in equation 5 or how they relate to the ‘VF’. The final sentence in that section relates the relaxation time to a fibrinogen concentration but there is no fibrinogen in the model. This section should be rewritten more carefully and explicitly state what the outcome is.

7. What is the bases for having the AA completely filled with platelets as the initial placement? It seems that one of the benefits of modeling the initial thrombus formation using FCM is that it will lead to more interesting and physiologically relevant thrombus formations. To keep inline with the journal criteria for significant biological insight, it would enhance this study to show a simulation with clots that have grown in the AA that initially was empty.

8. The claim is that the timescale of the phase-field modeling covers hours of time (Figure 1), but this claim is not justified by the simulations provided in this paper. Can the authors comment on this?

Minor comments:

1. The sensitivity of the surface tension parameter was tested for fixed values of h (and other parameters). The sensitivity was said to be negligible, but it was not clear what tests were used to come to that conclusion. Is the parameter sensitive as h changes? I don’t believe h was ever specified either.

2. In the clot deformation in 2D/3D section, how thick was the z-direction and were the results sensitive to that? As for the two densities and viscosities, are those for inside and outside the thrombus and which is where?

3. “calibration” would be a better description than “validation” for the permeability model, as it is truly calibration of a model representation.

4. Units. In most of the results sections, the units were left off of the parameters. It would be helpful to at least see units on the permeability

5. In the Modeling thrombin formation with FCM section, the term deposition sites is used, but in the figure, initiation sites is used. Also the reference to figure 11b on the last line of the first paragraph should reference 11c

6. What is the upstream distribution of platelets during the dynamic simulation?


Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biologydata availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes


PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Submitted filename: PCOMPBIOL-D-19-01313_Review.pdf

4 Jan 2020

Submitted filename: Review_for_phase_field.pdf

3 Feb 2020

Dear Professor Karniadakis,

We are pleased to inform you that your manuscript 'A Three-dimensional Phase-field Model for Multiscale Modeling of Thrombus Biomechanics in Blood Vessels' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch within two working days with a set of requests.

Reviewer #1's constructive comments should also be considered, and the text revised where appropriate.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology. 

Best regards,

Alison Marsden

Associate Editor

PLOS Computational Biology

Jason Haugh

Deputy Editor

PLOS Computational Biology


Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The revision has been improved considerably. However, there are further comments that should be addressed before publication.

Multiple places: The authors claim white clot is primarily made of fibrin and platelets, which is not consistent with early (Cadroy 89) and recent descriptions (e.g. Jackson 07 or Kim 19). White clot is composed of vWF and platelet. Fibrin is in red clots that may form on top of white clots after the occlusion causes stagnation of the blood.

Question 10: It was understood that the phase-field parameter ϕ is defined in terms of volume fractions. However, the authors use ϕ = 1 to denote blood flow and ϕ = 0 for thrombus. This is not consistent with the custom that the authors use for ϕ_f to denote the volume fraction of fibrin (solid), etc. I suggest defining ϕ = 0 as blood flow and thrombi as ϕ = 1 in the phase field model.

Question 12: the authors validate the permeability calculation using existing fibrin gel data. Fibrin gel is not a thrombus. It is OK to validate the model this way given the limited experimental data. However, the author should not oversell it by calling it "validation of the permeability of the thrombus". This is only a validation of the permeability calculation using fibrin gel data.

Question 13: It is understood that eventually an equilibrium configuration of the clot forms. However, this equilibrium is dynamic, meaning there is thrombi formation and rupture (or platelet adhering and detaching). This answer still does not address the fact the platelet aggregation may keep forming even during so-called "reconfiguration" process when emboli occurs. I suggest the authors state that a simplification/assumption is that they neglect the further formation of thrombi during the phase-field simulation.

Reviewer #2: The authors have answered all the reviewer questions and handled the concerns.


Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biologydata availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: None


PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

21 Apr 2020


A Three-dimensional Phase-field Model for Multiscale Modeling of Thrombus Biomechanics in Blood Vessels

Dear Dr Karniadakis,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Laura Mallard

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom | Phone +44 (0) 1223-442824 | | @PLOSCompBiol three-dimensional phase-field model for multiscale modeling of thrombus biomechanics in blood vessels&author=Xiaoning Zheng,Alireza Yazdani,He Li,Jay D. Humphrey,George E. Karniadakis,Alison Marsden,Jason M. Haugh,Alison Marsden,Jason M. Haugh,Alison Marsden,Jason M. Haugh,Alison Marsden,&keyword=&subject=Research Article,Biology and Life Sciences,Anatomy,Body Fluids,Blood,Platelets,Medicine and Health Sciences,Anatomy,Body Fluids,Blood,Platelets,Biology and Life Sciences,Physiology,Body Fluids,Blood,Platelets,Medicine and Health Sciences,Physiology,Body Fluids,Blood,Platelets,Biology and Life Sciences,Cell Biology,Cellular Types,Animal Cells,Blood Cells,Platelets,Physical Sciences,Physics,Classical Mechanics,Deformation,Physical Sciences,Physics,Classical Mechanics,Damage Mechanics,Deformation,Biology and Life Sciences,Anatomy,Body Fluids,Blood,Blood Flow,Medicine and Health Sciences,Anatomy,Body Fluids,Blood,Blood Flow,Biology and Life Sciences,Physiology,Body Fluids,Blood,Blood Flow,Medicine and Health Sciences,Physiology,Body Fluids,Blood,Blood Flow,Physical Sciences,Materials Science,Material Properties,Permeability,Medicine and Health Sciences,Hematology,Blood Coagulation,Platelet Aggregation,Research and Analysis Methods,Simulation and Modeling,Biology and Life Sciences,Biochemistry,Proteins,Fibrin,Medicine and Health Sciences,Vascular Medicine,Vascular Diseases,Aneurysms,