ORIGINAL RESEARCH article
An integrated fluid-chemical model toward modeling the formation of intra-luminal thrombus in abdominal aortic aneurysms
- 1 Department of Solid Mechanics, School of Engineering Sciences, KTH Royal Institute of Technology, Stockholm, Sweden
- 2 Mechanics Division, National Institute of Metrological Research, Turin, Italy
- 3 Department of Molecular Medicine and Surgery, Center for Molecular Medicine, Karolinska Institutet, Stockholm, Sweden
Abdominal Aortic Aneurysms (AAAs) are frequently characterized by the presence of an Intra-Luminal Thrombus (ILT) known to influence their evolution biochemically and biomechanically. The ILT progression mechanism is still unclear and little is known regarding the impact of the chemical species transported by blood flow on this mechanism. Chemical agonists and antagonists of platelets activation, aggregation, and adhesion and the proteins involved in the coagulation cascade (CC) may play an important role in ILT development. Starting from this assumption, the evolution of chemical species involved in the CC, their relation to coherent vortical structures (VSs) and their possible effect on ILT evolution have been studied. To this end a fluid-chemical model that simulates the CC through a series of convection-diffusion-reaction (CDR) equations has been developed. The model involves plasma-phase and surface-bound enzymes and zymogens, and includes both plasma-phase and membrane-phase reactions. Blood is modeled as a non-Newtonian incompressible fluid. VSs convect thrombin in the domain and lead to the high concentration observed in the distal portion of the AAA. This finding is in line with the clinical observations showing that the thickest ILT is usually seen in the distal AAA region. The proposed model, due to its ability to couple the fluid and chemical domains, provides an integrated mechanochemical picture that potentially could help unveil mechanisms of ILT formation and development.
Focal dilatations of the abdominal aorta, known as Abdominal Aortic Aneurysms (AAAs), are frequently observed in the aging male population (Fleming et al., 2005) and their rupture is fatal in up to 90% of cases (Upchurch and Schaub, 2006). Rupture is prevented through elective surgical intervention. Clearly, the ability to predict lesions that are at risk is required to optimize medical and economic outcomes.
AAAs are thought to be the end results of an irreversible pathological remodeling of the extracellular matrix due to an excessive proteolytic activity (Choke et al., 2005; Michel et al., 2011); in particular, loss of structural integrity of the major ground substances elastin and collagen. A thin or thick Intra-Luminal Thrombus (ILT) is a tissue found in nearly all AAAs large enough to indicate risk of rupture (Hans et al., 2005). While the thin ILT is not explored very well, the thick ILT (Folkesson et al., 2011) has solid-like properties (Gasser et al., 2008) and is composed of a fibrin mesh, traversed by a continuous network of interconnected canaliculi (Swedenborg and Eriksson, 2006), incorporated with blood cells, e.g., erythrocytes and neutrophils, aggregated platelets, blood proteins, and cellular debris. The ILT is thought to play an important role in the pathology and natural history of AAAs with a series of effects on the underlying aortic wall. Specifically it causes localized hypoxia, possibly leading to increased neovascularization, inflammation, and local wall weakening (Vorp et al., 2001). In addition, the changes to matrix-degrading protease expression (Kazi et al., 2005) and structural and cellular composition (Kazi et al., 2003) lead to a thinner wall compared to the aneurysm wall exposed to flowing blood (Swedenborg and Eriksson, 2006). The ILT has a significant structural impact on the biomechanics of AAAs and influences both the magnitude and the distribution of wall stress (Inzoli et al., 1993; di Martino et al., 1998; Wang et al., 2002; Polzer et al., 2011) and needs to be considered through a biomechanical rupture risk assessment (Gasser et al., 2010).
The hemostatic system maintains the integrity of the circulatory system in case of vascular damage. It maintains blood in a fluid state and responds to vessel injury by the rapid formation of a clot. The clot formation is the end result of a process initiated by the injury of a vessel wall and subsequent exposure of the subendothelium to blood flow. This triggers two interconnected processes: platelets (PLTs) aggregation and the CC. The first process involves platelets, anucleated cells that originate from bone marrow and circulate in blood as sentinels of vascular integrity. Platelets can be activated and adhere to the sites of exposed subendothelium. During this process the shape of the platelets modifies and chemicals are released in the blood stream so that new platelets are activated leading to aggregation with free or wall-bound ones.
The CC consists of a series of enzymatic reactions, in which a series of proenzymes (zymogens) is turned into their active enzyme form. The CC is triggered as soon as tissue factor (TF) binds to blood-borne factor VIIa. The series of activated reactions lead to the formation of thrombin, which in turn converts fibrinogen into fibrin (Gailani and Renné, 2007; Mackman et al., 2007). This process can lead to three outcomes: the impairment of clotting leading to bleeding disorders (like hemophilia), hemostasis, or hypercoagulability that in turn lead to thrombotic events such as the formation of an ILT. Platelet activation, aggregation, and adhesion can be regarded as the initiating response of thrombus formation that arrests hemorrhage in response to vascular injury and permits wound healing (Ruggeri and Mendolicchio, 2007), whereas actual hemostasis requires both platelets and the CC to occur.
It has long been known (Bluestein et al., 1996) that blood dynamics has a relevant influence on the clotting processes, consequently a meaningful ILT formation model should account for the biochemical and biomechanical interactions as it has been proposed by different studies in the literature (Hubbell and McIntire, 1986; Folie and McIntire, 1989; Kuharsky and Fogelson, 2001; Fogelson and Tania, 2005; Xu et al., 2008). However, simplistic geometrical, hemodynamical, and sometimes even biochemical, assumptions were applied. On the other hand, works such as Xu et al. (2008), despite the aforementioned limitations, present a multiscale approach able to integrate also the discrete blood elements and, therefore, are fundamental for a microscopic description of the coagulation process. The present work considers biochemical and biomechanical factors in the abdominal aorta through a macroscopic model of the CC and the Navier–Stokes equations, respectively. Specifically, the coupling is obtained by a set of convection-diffusion-reaction (CDR) equations added on top of the computed blood flow field. The proposed model predicts the distribution of chemicals in complex AAA blood flows and hence has the potential to improve our current understanding of ILT pathophysiology.
2. Materials and Methods
2.1. Geometry Representation and Discretization
An idealized 2D-axisymmetric AAA has been modeled to simulate a small fusiform aneurysm of 4.4 cm in (luminal) diameter, a size of clinical relevance sufficiently large to alter blood flow (Figure 1). We consider a thin-ILT-covered wall (Figure 3), a typical feature of AAA progression (Swedenborg, J. Personal Communication), with two regions of exposed subendothelium: Case (A) considers a large exposed area and Case (B) a focal one (Figure 1). The two cases aim at investigating the effect of varying size and position of the exposed subendothelium.
Figure 1. Axisymmetric fusiform abdominal aortic aneurysm (AAA) of 4.4 cm in (luminal) diameter with the dotted line denoting the axis of symmetry. Two cases of exposed subendothelium were introduced: Case (A) considers a large exposure and Case (B) a focal exposure, see also Section 2.2.3.
The domain has been discretized with an unstructured triangular mesh with higher density within concentration and velocity boundary layers. Details regarding the mesh are presented in Table 1 and in Figure 2 a region of the mesh is shown. The mesh was generated using the Advancing Front Method (Thompson et al., 1998) provided by COMSOL v3.5a (COMSOL AB).
Figure 2. Finite element mesh used to solve the fluid dynamical and chemical set of equations: (A) bulge region, (B) zoom of the near wall region.
A mesh convergence study verified that the computed results were mesh independent. For the fluid dynamical problem the velocity, pressure, and wall shear stress fields were sampled at different locations and times during the cardiac cycle (see Prakash and Ethier, 2001). For the CDR problem the concentration of thrombin, which is the final product of the modeled chemical reactions, was sampled.
2.2. Mathematical Models and Solution Procedure
2.2.1. Fluid dynamical model
The fluid dynamical field was rendered by the continuity and Navier–Stokes equations,
where u is the velocity vector, ρ the density, p the pressure, and τ = 2μD the deviatoric viscous stress tensor (Acheson, 2003). Boundary conditions, as previously used (Biasetti et al., 2010, 2011), were applied. At time t the inlet mass flow rate was prescribed and the inlet velocity profile u was calculated according to the no-slip boundary condition and the Navier–Stokes equations. At the outlet a prescribed pressure waveform with no-viscous stress boundary condition was applied. The prescribed waveform takes the form of p = pt(t), where pt(t) is a tabulated set of pressures in time, whilst the no-viscous stress condition prescribes that τn = 0. In the case of rigid wall and incompressible fluid the imposition of a pressure wave at the outlet is not strictly necessary since given the mass flow rate, the pressure differential (Δp) is calculated. This means that a constant pressure may also be applied. However, in order to get a meaningful pressure distribution (amplitude and phase) a physiological pressure wave need to be used. For a detailed explanation of different boundary conditions the reader is referred to Vignon-Clementel et al. (2006). The no-slip boundary condition was applied at the walls and at the symmetry axis an axial symmetry condition was set.
2.2.2. Constitutive modeling of blood
Blood has complex rheological properties involving shear-thinning, thixotropy, and viscoelasticity (Oka, 1981; McDonald, 2004). Its non-Newtonian behavior has recently been shown (Biasetti et al., 2011) to produce substantial differences in the local flow pattern compared to a Newtonian viscosity model; as such differences are mostly confined on a local level, a model dealing with local phenomena like platelets’ dynamics and chemicals’ concentrations should take the non-Newtonian behavior into proper consideration. We adopted the Carreau-Yasuda (CY) model for shear-thinning (Bird et al., 1987)
where denotes the scalar shear rate with D = (l+ lT)/2 being the symmetric part of the velocity gradient tensor l (Malvern, 1969). The antisymmetric part Ω = (l− lT)/2 of l, as it is used in Section 3, was also introduced. Note that at the lower and upper ends of the shear rate spectrum equation (3) describes blood as a Newtonian fluid of viscosities μ0 and μ∞, whereas the parameters λ, n, and a define the transition between these asymptotic conditions. Table 2 reports the Carreau-Yasuda model parameters applied in the present work, which were selected to fit experimental data of blood at 37.0°C (Leuprecht and Perktold, 2001; Abraham et al., 2005).
2.2.3. The coagulation cascade to thrombin formation
The CC is a series of subsequent reactions, where the activation of one clotting factor activates the ensuing one, finally leading to thrombin generation. The CC is triggered at sites of endothelial damage, i.e., when extravascular tissue factor (TF), which is present under the endothelial layer, binds with blood-borne factor VIIa (Figure 3A).
Figure 3. (A) Initiation of the CC at sites of subendothelial exposure, i.e., when extravascular tissue factor (TF) binds with blood-borne factor VIIa. (B) CC: subendothelial exposure triggers thrombin formation via the extrinsic (Tissue Factor, TF) pathway. Positive feedback (denoted by the dotted arrow) of the intrinsic pathway increases thrombin formation. The reaction from prothrombin to thrombin formally belongs to the common pathway, here indicated by the oval.
The CC can be divided into the extrinsic, the intrinsic (or contact), and the common pathways (Figure 3B). The extrinsic pathway is also referred as the tissue factor pathway and it is thought to be the primary activator of the CC (at least for cases where a foreign surface is not involved for which the intrinsic pathway is thought to start the CC; Basmadjian et al., 1997) while being also the essential one, meaning that without its reactions the CC cannot start (Hoffman and Monroe, 2007; Mackman et al., 2007).
We followed the CC model proposed in Jones and Mann (1994) and empirically supported by a companion work (Stram et al., 1994). It consists of 18 species and involves plasma-phase and surface-bound enzymes and zymogens (see Table 3) where the subscript “a” denotes the active complex. The model includes both plasma-phase and membrane-phase reactions. Note that for the present CC model, tissue factor (TF), and factor VIIa are not considered separately but they are already joined to form TF:VIIa.
A time-dependent version of the CC model has been implemented, where the 16 chemical reactions read:
Despite more recent models include also the stoichiometric (tissue factor pathway inhibitor (TFPI) and antithrombin) and dynamic (dynamic protein C (PC) system) inhibitory processes (see for example Hockin et al., 2002), AAA pathology motivates the choice of the CC model (4)–(15). In the vasculature TFPI is mainly bound to the endothelium (Österund, 2012), which is to a large extent not present in AAAs. However, TFPI is also present in blood at a very low concentration and may still play a role in inhibiting the clotting process. For the same reason, i.e., lack of endothelium, the dynamic protein C (PC) system has not been considered in the present model since it needs, in order to be activated, the integral membrane protein thrombomodulin which is expressed on the surface of endothelial cells. Antithrombin is a relatively inefficient inhibitor on its own and needs to bind to heparan sulfate proteoglycans expressed on the vascular endothelium or heparin in order to accelerate its function (Rosenberg and Bauer, 1994). Heparan sulfate proteoglycans are not available due to the lack of endothelium and the effect of heparin is debated (Marcum et al., 1986). These reasons led to the exclusion of this contribution.
To show the reaction kinetics behavior in a simple case, the chemical reactions in equations (4)–(15) were first computed for the case of a perfectly mixed physical environment with constant volume, i.e., in a batch reactor, and at a constant temperature of 37°C. Initially all chemical species are at physiological concentration, see Table 4.
Introducing a concentration vector c ([mol/m3])1, the governing equation of the CC model reads
where S denotes the stoichiometric matrix and r is the reaction rate vector ([mol/(m3 s)])2. Specifically, S is a 18 × 16 matrix that reads
where reversible reactions are considered through single reactions. The 16 columns in S represent the 16 chemical reactions in equations (4)–(15) that build up the model and c stores the concentration of the species in the following order: IX, TF:VIIa, IX:TF:VIIa, IXa, X, X:TF:VIIa, Xa, VIIIa:IXa, X:VIIIa:IXa, V, Va, VIII, VIIIa, IIa, II, Va:Xa, II:Va:Xa, mIIa. The rate constants, used in the CC model, are given in Table 5.
Table 5. Rate constants used in the CC model (see Jones and Mann, 1994).
2.2.4. Convection-diffusion-reaction model
In the case of the batch reactor, chemical concentrations can be considered as a function of time only, because of the so called well-mixed assumption. On the other hand, in blood flows, including flows in AAAs, the concentrations vary in space and time, leading therefore to a 3D time-dependent problem. This can be solved by coupling the fluid dynamical model (see Section 2.2.1) to the time-dependent CC model (see Section 2.2.3). Specifically, the coupled model defines a set of evolution equations for species’ concentrations ci, i = 1, …, 18 that takes the form of a CDR problem and reads
where Ri are the reaction terms (defined by the CC model, see Section 2.2.3), u the blood velocity vector (defined by the Navier–Stokes equations, see Section 2.2.1) and Di the isotropic diffusion coefficients. Values for Di were taken from Folie and McIntire (1989), Xu et al. (2008) and set equal, slightly modified (see Section 4), to 10−8 m2/s. Chemical species’ concentrations at the inlet are set equal to the values in Table 4 except for TF:VIIa which is set equal to zero. At sites of subendothelial exposure TF:VIIa is set equal to 1 × 10−6 mol/m3. The domain is initialized with the inlet values. The final CDR problem is clearly convection-dominated and requires numerical regularization as detailed in Section 2.2.5. Despite TF:VIIa, IX:TF:VIIa, and X:TF:VIIa are bound to the subendothelium, we model them also in the fluid domain by adding CDR equations for these three compounds. This choice is motivated by the fact that blood contains microparticles (MPs) that can express TF, allowing the formation of the TF:VIIa, IX:TF:VIIa, and X:TF:VIIa compounds (Owens and Mackman, 2011). These three compounds are not present in blood at the beginning of the simulation (see Section 3.2) but progressively appear from the site of exposed subendothelium, which is thought to model the gradual appearance of TF-bearing MPs from activated blood elements. Moreover, considering an intravascular source of TF is crucial in thrombosis, and in particular under conditions of low flow and impaired endothelium (Del Conde et al., 2005; Langer et al., 2008; Antoniak et al., 2009; Morel et al., 2009; Mann, 2011). Finally, it is noted that the present CC model is based on the assumption of an infinite supply of platelets, which in turn leads to an infinite binding site density. This assumption is satisfied in our computational model due to the continuous supply of fresh blood from upstream. For a discussion on the implications of a finite density of binding sites the reader is referred, for example, to Kuharsky and Fogelson (2001).
2.2.5. Finite element discretization and numerical algorithms
184.108.40.206. Time-dependent solution of the CC model. To solve the time-dependent CC model (see Section 2.2.3), the IDA solver for Differential-Algebraic Equation (DAE) systems (Hindmarsh et al., 2005), as implemented in COMSOL Reaction Engineering Lab v3.5a (COMSOL AB), was used.
220.127.116.11. Finite element discretization. The fluid dynamical problem was discretized with Lagrange P2P1 elements, Lagrange elements with quadratic interpolation for the velocity field and linear interpolation for the pressure field, whereas linear elements were used for the CDR problem (Zienkiewicz and Taylor, 2000). The Galerkin method was used to discretize the equations, which may become unstable for an element Peclet number (Pe) larger than one, i.e., Pe = |u|h/(2D) > 1 with h being the mesh element size. Since both problems, fluid dynamics and CDR, are convection-dominated (Pe > 1; Figure 4), stabilization techniques have been employed.
Figure 4. Peclet number contours for the CDR problem at peak systole (t = 0.25 s) and late diastole (t = 1.00 s). In both cases the Peclet number is much larger than one in the whole domain, thus requiring numerical stabilization.
The Navier–Stokes equations’ solution typically includes shear layers (embedded within the flow or at its boundaries) where strong gradients of the solution variables are present. The width of these layers can be smaller than the element size, and hence the solution variables cannot be resolved properly, leading to spurious oscillations (wiggles) in the solution. Galerkin least-squares (GLS) streamline diffusion combined with crosswind diffusion (Ck= 0.1) was used to stabilize the Navier–Stokes equations (Hauke, 2002; COMSOL, 2008).
The CDR problem is more challenging to stabilize. Wiggles, even of small amplitude, are not permissible in order to ensure non-negative species concentrations. Isotropic diffusion at the lowest possible level was introduced and δISO= 0.1 was set (COMSOL, 2008).
18.104.22.168. Time-advancing algorithm. The fluid dynamical problem and the CDR one are decoupled. Consequently, the Navier–Stokes equations were solved first and the solution (velocity vector u) then fed to the CDR set of equations. The Generalized-α method (Jansen et al., 2000) with the numerical damping parameter was used to solve the Navier–Stokes equations. The different rate constants in the chemical reactions in equations (4)–(15) define a CDR problem with temporal scales encompassing nine orders of magnitude requiring very small time steps for the non-regularized problem. The numerical parameter was found to be appropriate for canceling out the non-resolvable frequencies while retaining a good degree of accuracy in the computed solution.
The arising non-linear set of equations was solved using an affine invariant form of the damped Newton method (Deuflhard, 1974) together with the PARDISO direct solver (http://www.pardiso project.org/) available in COMSOL v3.5a (COMSOL AB).
All computations have been performed on a 64-bit PC equipped with 2 dual core processors Intel(R) Xeon(TM) CPU 3.40 GHz and 16 GB RAM, with Windows Server 2008 R2 Enterprise Operating System. The applied numerical parameters are summarized in Table 6.
Table 6. Numerical parameters used to regularize the fluid dynamical and convection-diffusion-reaction (CDR) problems.
2.3 Eduction Scheme: The λ2-method
The concept of coherent vortical structures or VSs enables a precise analysis of complex flow fields. Nonetheless, the lack of an accepted mathematical definition of what constitutes a coherent vortical structure did lead to issues in its identification. The most accepted and widely used method of eduction (the process of finding a VS in a flow field) is the λ2-method (Jeong and Hussain, 1995). The idea behind this technique is that a vortex is a flow region with predominance of vorticity surrounding a local pressure minimum. The method is based on a reduced version of the Navier–Stokes equation for incompressible flows which ignores contributions from unsteady straining and viscous effects:
This reduced version allows to analyze the pressure minima solely in relation to vortical motion; a vortex is identified as a connected region with two negative eigenvalues of D2 + Ω2. Since D2 + Ω2 is symmetric, its eigenvalues are real. Ranking the eigenvalues as λ1 ≥ λ2 ≥ λ3 defines the equivalent condition λ2 < 0 for the presence of a VS. Note that for visualization purposes, the condition λ2 = λ2tr is typically applied, where λ2tr denotes a threshold value defining the surface of the VS.
3.1. CC Model
The CC is triggered as soon as TF:VIIa comes into contact with factor IX or X, i.e., at sites of subendothelial exposure, and the following evolution of species concentrations is shown in Figure 5. Initially the production of thrombin is slow (time-lag phase) and increases then rapidly due to the activation of the intrinsic pathway (accelerated phase). This boost is commonly referred to be the effect of the positive thrombin feedback mechanism, that is indeed the intrinsic pathway. Finally, thrombin concentration reaches the plateau value of 1.4 × 10−3 mol/m3, i.e., equal to the initial concentration of prothrombin (plateau phase).
Figure 5. Evolution of species’ concentrations in the CC model. Thrombin production shows three distinct phases, i.e., time-lag, accelerated, and plateau phases. The concentrations of some chemical species are too low to be seen in the figure.
3.2. Fluid-Chemical Model
In this section the evolution of thrombin in relation to the blood flow field is studied. Specifically, the role of VSs was investigated keeping in mind that in the limit of inviscid fluid (in the present case this is therefore still valid although to some approximation) they represent material surfaces/volumes. Consequently, VSs generated at some position tend to entrain the matter present at its formation site, i.e., VSs generated at the wall capture chemicals at that position and transport them along their trajectory.
3.2.1. VSs’ dynamics and flow field
At the beginning of the systolic phase (t ≈ 0.15 s), the wall shear layer separates in the expansion region giving rise to a free vorticity layer [Figure 6 (0.2 s) and Figure 7]. This free shear layer quickly increases the bending of its profile through a process driven by the velocity difference between the two sides and augmented by self-induction (Kelvin–Helmholtz instability, see, e.g., Tritton, 1988). During this process it rolls up, generating a counter-clockwise VS, from now on indicated as VS1, which moves away from the wall [Figure 6 (0.3 s) and Figure 7 second row]. The velocity field induced by VS1 generates a negative velocity gradient in the direction parallel to the wall and pointing upward, thus accumulating the diffused vorticity in the wall region. In addition, the wall-normal induced velocity component promotes the detachment of the lump of vorticity thus generated at the wall. These two effects tend synergically to cause the formation of a secondary VS [VS2; Figure 6 (0.3 and 0.4 s) and Figure 8]. Once detached from the wall, VS2 moves upstream because of the velocity field induced by VS1 and the influence of the mean flow [Figure 6 (0.4 s) and Figure 7]. VS1, on the other hand, moves downstream because of self-induction. Note that, as the computation is axisymmetric, all VSs must be seen as the trace on the computational plane of a ring vortex. Therefore VS1 undergoes the standard self-induction process, i.e., every section of the vortex induces downstream (because of the sign of the vorticity) velocity on the opposite section, leading to the final result of a downstream movement of the whole vortex. While moving down the lumen, VS1 maintains a strong identity, although the intensity is slightly reduced via viscous diffusion. When VS1 approaches the distal bend, it induces the generation of a shear layer which rolls up in a fashion very similar to what explained above but backward (Figure 7). This shear layer causes a remarkable wall shear stress (WSS) peak illustrated in Figure 9.
Figure 6. Vorticity magnitude (s−1) (first row) and Vortical Structures (VSs), educed with the λ2-method (s−2), (second row) dynamics at three selected times t in the cardiac cycle. The sequence highlights the formation process of the VSs (VS1 and VS2) as a result of the vorticity sheet separation and roll-up (details in the text).
Figure 7. Birth and evolution of Vortical Structures (VSs) at different times t throughout the cardiac cycle. A vortex sheet formed at the proximal end of the AAA (t = 0.15–0.2 s) develops in a vortex (VS1) rotating counter-clockwise (t = 0.3 s). The vortex moves downward until it impinges on the distal contraction of the AAA. The genesis and motion of the second vortex (VS2) is also visible.
Figure 8. Mechanism of VS-induced wall-vorticity accumulation and Vortical Structures (VSs) detachment. The main VS (VS1), superimposed on a main flow directed from left to right, induces a velocity field in the near wall region which accumulates vorticity forming another VS (VS2) and, subsequently, promotes its detachment.
Figure 9. Spatially averaged wall shear stress (AWSS) in the distal contraction during the cardiac cycle. During Vortical Structures (VSs) impingement (between t = 0.15 s and t = 0.4 s) the averaged WSS increases more than four times with respect to the diastolic values.
WSS is highly affected by VSs, as demonstrated by Figure 9 and our previous work (Biasetti et al., 2011); in particular a correlation between WSS and near wall VSs exists. Figure 9 shows the space-averaged WSS (AWSS) over the distal vortex impinging region. AWSS hovers below 0.5 Pa throughout the cardiac cycle due to the larger passage area caused by the aneurysmatic bulge, except when VS1 impinges on the wall. This happens between t = 0.1 and t = 0.4 s, where the AWSS reaches values up to 2.1 Pa. Notice that this is a spatially averaged WSS: the WSS is observed to peak at around 3.5 Pa within the averaging region. The vortex dynamics described here is coherent with experimental observations presented in Salsac et al. (2006), thus giving confirmation of the validity of the present results. For a more detailed discussion of VSs’ dynamics in AAA and a comparison between Newtonian and non-Newtonian models on VSs’ behavior the reader is referred to Biasetti et al. (2011).
3.2.2. Chemical field
At the sites of subendothelial exposure TF:VIIa, present with the constant concentration of 1 × 10−6 mol/m3, comes into contact with factors IX and X and triggers the CC. Recall that in this case the concentration of TF:VIIa in the fluid domain is initially zero, compare with Table 4. Once triggered, the CC propagates inside the fluid and leads to thrombin production also within the lumen. After an initial transient period of about 15–20 cardiac cycles, depending on the initial subendothelial exposure area, thrombin concentration reaches a periodic solution. The combined effect of VSs shedding and recirculating region inside the aneurysmatic bulge convect thrombin to the distal AAA region in both cases of endothelial damage, i.e., case (A) and (B) (Figure 10). In the present case VS1 captures thrombin that is generated proximally and transports it distally due to the described VS-mediated transport phenomenon, see Section 3.2. The recirculation region, present during parts of the cardiac cycle, helps the produced chemicals to remain inside the aneurysmatic bulge. Noteworthy, VS1 moves downstream with an anti-clockwise rotation and creates a trail of thrombin that curls up in an anti-clockwise way (Figure 11). Fluid convection prevents thrombin over-accumulation in the distal area by washing chemicals downstream. Consequently, fluid flow might then be regarded as a “mechanical antagonist” of procoagulant accumulation.
Figure 10. Thrombin (IIa) distribution once the periodic state is reached in the two investigated cases. Considering a large [Case (A)] and a small [Case (B)] subendothelial exposure (see Section 2.1) both concentration patterns are strongly shifted toward the distal Abdominal Aortic Aneurysm (AAA) region. A larger area of high thrombin concentration is found in Case (A) due to the larger endothelial damage.
Figure 11. Vortical Structures (VSs) educed with the λ2-method (left) and thrombin (IIa) distribution (right) at time t = 0.4 s during the cardiac cycle. The thrombin distribution is shown using a logarithmic scale. As clearly shown, the moving VS is carrying thrombin and the characteristic counter-clockwise trail is due to the counter-clockwise rotation of the VS.
The link between biomechanical loads and biochemical signals is thought to be one of the key aspects of cardiovascular homeostasis and pathology development. Due to the complexity of the involved phenomena this link has not yet been fully elucidated. AAA is a complex pathology involving biological, mechanical, and fluid dynamical aspects and an integrated mechanochemical view may help understand its development and the role played by the ILT. Consequently, our previous fluid dynamical model (Biasetti et al., 2011) was coupled to the biochemistry allowing the study of the CC in complex flows. In earlier works (Biasetti et al., 2010, 2011) we explored the behavior of VSs and recirculating regions and we found indications that both promote favorable conditions for PLT activation, convection, and deposition at the wall in the distal AAA region. Related to this is the finding that thrombus thickness is largest in the distal portion of the aneurysmatic aorta (Figure 12). VSs capture chemicals (Figure 11), in their core and convey them until VSs’ burst (the breaking-up of the vortex), seen to happen at the distal contraction. The strong thrombin accumulation observed distally correlates with the thickness of the ILT layer, and suggests a possible correlation between fluid dynamics, biochemistry, and ILT growth. Thrombin converts soluble fibrinogen into insoluble strands of fibrin, which in turn form the clot in conjunction with platelets and other blood elements. Consequently, the accumulated thrombin in the distal AAA region might lead to an increased fibrin production enhancing ILT growth with respect to other regions. In particular, in whole blood activated with tissue factor, the production of a platelet-fibrin clot starts with thrombin concentration of approximately 2 × 10−6 mol/m3 (Mann, 2011) and during normal coagulation the concentration of thrombin can reach 5 × 10−4 mol/m3 (Wolberg and Campbell, 2008). In the cases discussed here, the thrombin concentration was observed to reach values as high as twice that reported above, i.e., ≈1 × 10−3 mol/m3.
Figure 12. ILT thickness measured on slices perpendicular to the centerline from 61 patients with small AAAs (<5.5 cm). Dots and bars denote the median and the lower and upper quartile respectively. The thickest ILT is located in the distal region of the aneurysm, at around 70% of its length starting from the renal arteries. Further details are given in Martufi et al. (submitted).
Thrombus formation in the venous and arterial systems differs markedly. In the venous system, the low flow rates and stasis allow the accumulation of activated coagulation factors leading to the formation of thrombin mostly without platelets contribution. On the other hand, in the normal arterial circulation, flow conditions prevent this accumulation (Brass et al., 2007) and platelets play a major role in arterial thrombus formation: they accelerate thrombin formation and provide a scaffold for fibrin accumulation. In AAAs, due to pathological flow conditions involving stasis and recirculation, the normal arterial clotting process might be altered and “switch” to a mechanism similar to the one observed in the venous circulation with an enhanced accumulation of procoagulant. In addition, the convective effect of VSs tends to favor the presence of activated platelets in regions of high chemical concentrations. The combination of these effects could cause the thicker ILT observed in the distal region of the lumen – this picture being consistent with the views put forward in Swedenborg (2008), Yoshimura et al. (2011).
Below we discuss the remaining open points and limitations of the present study that should be addressed in future research.
4.1. Chemical Model
Our model does not take into account the inhibitory processes for the reasons explained previously, but blood-borne anticoagulants from the proximal healthy arterial system might still play an important role. However to the best of the authors’ knowledge no indication has been reported in literature. Our model can then be enriched by considering this contribution too. Platelets recruited to the exposed subendothelium accumulate and cover it, blocking the source of TF:VIIa at the wall (Hathcock and Nemerson, 2005). Despite this, human thrombotic occlusion has been observed in vivo on a timescale of minutes (Ambrose et al., 1997) indicating that another source of procoagulant activity is available at the fluid-clot interface. Blood-borne TF has been indicated as one of the contributors (Hathcock and Nemerson, 2005). Our model does not take into account platelet accumulation at the wall allowing the continuous availability of TF:VIIa at the wall. Since the actual accretion of the ILT is not simulated in our model, the continuous availability of TF:VIIa at the wall substitutes, in first approximation, the blood-borne source of TF present at the fluid-clot interface.
To avoid numerical instabilities the diffusion coefficients used in the CDR model were increased from about 10−10 m2/s (Folie and McIntire, 1989; Xu et al., 2008) to 10−8 m2/s, which resulted in an over-diffusive solution. Since the underlying chemical field is strongly convection-dominated the influence of this numerical regularization is probably negligible, so that the general conclusions drawn here will still hold especially since lower diffusion coefficients are expected to result in even more segregation of thrombin in the near wall region.
4.2. ILT Growth
The present model does not consider the actual ILT growth which will modify the local fluid dynamics. This alteration is difficult to predict but we expect that the growing ILT will induce VSs’ burst location to move slightly more proximally, a feature seen in patient-specific geometries (Biasetti et al., 2011). The migration of the VSs’ burst location might also be a consequence of the lost axisymmetry of the fluid domain during ILT evolution – a feature that is constrained in the present model.
It should be noted that VSs’ burst happens distally both in the present case and in patient-specific geometries (Biasetti et al., 2011), leading to the postulation that this pattern of thrombin accumulation might be found also in patient-specific cases.
4.3. 2D-Axisymmetry Approximation
3D simulations of a similar flow (result not shown) reveal that in 3D the ring vortex undergoes the classical azimuthal bending instability (a waviness of its profile) which however does not modify the flow significantly. We therefore expect the 2D-axisymmetric model to provide the same results, in terms of chemical distributions, as the 3D one, with the advantage of saving computational time.
4.4. Parameter Space
Parameters like inlet flow rate, geometry of the AAA and initial concentration of TF:VIIa can influence the distribution of the chemicals involved. A parametric study would provide important insights regarding the sensitivity of the modeled system with respect to these parameters.
4.5. Stabilization Method
Isotropic diffusion, used to stabilize the CDR problem, is an inconsistent method and more sophisticated stabilization techniques are available (e.g., SUPG). Unfortunately, we observed that such techniques could not guarantee a non-negative concentration of the chemical species, probably due to the strong destabilizing effect of the reaction terms. For this reason, the choice of a less sophisticated but more robust method has been compulsory at the present stage.
4.6. Fluid-Structure Interaction
One limitation of the present work is the assumption of rigid wall. Despite the small motion of the aneurysm wall during the cardiac cycle (Länne et al., 1992; Long et al., 2004), and good agreement between Fluid-Structure Interaction (FSI) and rigid wall models (Wolters et al., 2005), the soft ILT tissue may deform significantly during the cardiac cycle influencing the fluid dynamics and hence the distribution of chemicals.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We are extremely thankful to Professor Jay D. Humphrey (Biomedical Engineering, Yale University, Connecticut, USA) for the fruitful exchange of ideas regarding platelets biochemistry and M.Sc. Fred “Nooj” Nugen (ICES, The University of Texas at Austin) for the helpful discussions regarding numerical methods. We would also like to thank M.Sc. Salvatore Nania [Department of Clinical Science, Intervention and Technology (CLINTEC), Karolinska Institutet, Stockholm, Sweden] for his help in the literature survey. This work has been supported by the Young Faculty Grant No. 2006-7568 provided by the Swedish Research Council, VINNOVA, and the Swedish Foundation for Strategic Research, which are gratefully acknowledged.
Ambrose, J. A., Almeida, O. D., Sharma, S. K., and Ratner, D. G. (1997). Angiographic evolution of intra-coronary thrombus and dissection following percutaneous transluminal coronary angioplasty (the Thrombolysis and Angioplasty in Unstable Angina [TAUSA] trial). Am. J. Cardiol. 79, 559–563.
Antoniak, S., Boltzen, U., Eisenreich, A., Stellbaum, C., Poller, W., Schultheiss, H. P., and Rauch, U. (2009). Regulation of cardiomyocyte full-length tissue-factor expression and microparticles release under inflammatory conditions in vitro. J. Thromb. Haemost. 7, 871–878.
Biasetti, J., Gasser, T. C., Auer, M., Hedin, U., and Labruto, F. (2010). Hemodynamics conditions of the normal aorta compared to fusiform and saccular abdominal aortic aneurysms with emphasize on thrombus formation. Ann. Biomed. Eng. 38, 380–390.
Biasetti, J., Hussain, F., and Gasser, T. C. (2011). Blood flow and coherent vortices in the normal and aneurysmatic aortas. A fluid dynamical approach to intra-luminal thrombus formation. J. R. Soc. Interface 8, 1449–1461.
Bluestein, D., Niu, L., Schoephoerster, R. T., and Dewanjee, M. K. (1996). Steady flow in an aneurysm model: correlation between fluid dynamics and blood platelet deposition. J. Biomed. Eng. 118, 280–286.
Choke, E., Cockerill, G., Wilson, W. R., Sayed, S., Dawson, J., Loftus, I., and Thompson, M. M. (2005). A review of biological factors implicated in abdominal aortic aneurysm rupture. Eur. J. Vasc. Endovasc. Surg. 30, 227–244.
Del Conde, I., Shrimpton, C. N., Thiagarajan, P., and Lopez, J. A. (2005). Tissue-factor-bearing microvesicles arise from lipid rafts and fuse with activated platelets to initiate coagulation. Blood 106, 1604–1611.
di Martino, E. S., Mantero, S., Inzoli, F., Melissano, G., Astore, D., Chiesa, R., and Fumero, R. (1998). Biomechanics of abdominal aortic aneurysm in the presence of endoluminal thrombus: experimental characterization and structural static computational analysis. Eur. J. Vasc. Endovasc. Surg. 15, 290–299.
Fleming, C., Whitlock, E. P., Beil, T., and Lederle, F. A. (2005). Review: screening for abdominal aortic aneurysm: a best-evidence systematic review for the U.S. preventive services task force. Ann. Intern. Med. 142, 203–211.
Folie, B. J., and McIntire, L. V. (1989). Mathematical analysis of mural thrombogenesis. Concentration profiles of platelet-activating agents and effects of viscous shear flow. Biophys. J. 56, 1121–1141.
Gasser, T. C., Auer, M., Labruto, F., Swedenborg, J., and Roy, J. (2010). Biomechanical rupture risk assessment of abdominal aortic aneurysms. Model complexity versus predictability of finite element simulations. Eur. J. Vasc. Endovasc. Surg. 40, 176–185.
Gasser, T. C., Görgülü, G., Folkesson, M., and Swedenborg, J. (2008). Failure properties of intraluminal thrombus in abdominal aortic aneurysm under static and pulsating mechanical loads. J. Vasc. Surg. 48, 179–188.
Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban, R., Shumaker, D. E., and Woodward, C. S. (2005). SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Trans. Math. Softw. 31, 363–396.
Jansen, K. E., Whiting, C. H., and Hulbert, G. M. (2000). A generalized-α method for integrating the filtered Navier-Stokes equations with a stabilized finite element method. Comput. Methods Appl. Mech. Eng. 190, 305–319.
Kazi, M., Thyberg, J., Religa, P., Roy, J., Eriksson, P., Hedin, U., and Swedenborg, J. (2003). Influence of intraluminal thrombus on structural and cellular composition of abdominal aortic aneurysm wall. J. Vasc. Surg. 38, 1283–1292.
Kazi, M., Zhu, C., Roy, J., Paulsson-Berne, G., Hamsten, A., Swedenborg, J., Hedin, U., and Eriksson, P. (2005). Difference in matrix-degrading protease expression and activity between thrombus-free and thrombus-covered wall of abdominal aortic aneurysm. Arterioscler. Thromb. Vasc. Biol. 25, 1341–1346.
Langer, F., Spath, B., Haubold, K., Holstein, K., Marx, G., Wierecky, J., Brummendorf, T. H., Dierlamm, J., Bokemeyer, C., and Eifrig, B. (2008). Tissue factor procoagulant activity of plasma microparticles in patients with cancer-associated disseminated intravascular coagulation. Ann. Hematol. 87, 451–457.
Länne, T., Sonesson, B., Bergqvist, D., Bengtsson, H., and Gustafsson, D. (1992). Diameter and compliance in the male human abdominal aorta: influence of age and aortic aneurysm. Eur. J. Vasc. Surg. 6, 178–184.
Long, A., Rouet, L., Bissery, A., Rossignol, P., Mouradian, D., and Sapoval, M. (2004). Compliance of abdominal aortic aneurysms: evaluation of tissue Doppler imaging. Ultrasound Med. Biol. 30, 1099–1108.
Michel, J. B., Martin-Ventura, J. L., Egido, J., Sakalihasan, N., Treska, V., Lindholt, J., Allaire, E., Thorsteinsdottir, U., Cockerill, G., and Swedenborg, J. (2011). Novel aspects of the pathogenesis of aneurysms of the abdominal aorta in humans. Cardiovasc. Res. 90, 18–27.
Morel, O., Pereira, B., Averous, G., Faure, A., Jesel, L., Germain, P., Grunebaum, L., Ohlmann, P., Freyssinet, J. M., Bareiss, P., and Toti, F. (2009). Increased levels of procoagulant tissue factor-bearing microparticles within the occluded coronary artery of patients with ST-segment elevation myocardial infarction: role of endothelial damage and leukocyte activation. Atherosclerosis 204, 636–641.
Polzer, S., Gasser, T. C., Swedenborg, J., and Bursa, J. (2011). Intra-luminal thrombus failure on the mechanical stress in the wall of abdominal aortic aneurysms. Eur. J. Vasc. Endovasc. Surg. 41, 467–473.
Rosenberg, R. D., and Bauer, K. A. (1994). “The heparin-antithrombin system: a natural anticoagulant mechanism,” in Hemostasis and Thrombosis: Basic Principles and Clinical Practice, 3rd Edn, eds R. W. Colman, J. Hirsch, V. J. Marder, and E. W. Salzman (Philadelphia, PA: Lippincott), 837–860.
Salsac, A. V., Sparks, S. R., Chomaz, J. M., and Lasheraz, J. C. (2006). Evolution of the wall shear stresses during the progressive enlargement of symmetric abdominal aortic aneurysms. J. Fluid Mech. 560, 19–51.
Swedenborg, J. (2008). “The role of the intraluminal thrombus in abdominal aortic aneurysms,” in Aortic Aneurysms, New Insights into an Old Problem, ed. N. Sakalihasan (Liége: Édition de l’Université de Liége), 343–346.
Vignon-Clementel, I. E., Figueroa, C. A., Jansen, K. E., and Taylor, C. A. (2006). Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Comput. Methods Appl. Mech. Eng. 195, 3776–3796.
Vorp, D. A., Lee, P. C., Wang, D. H., Makaroun, M. S., Nemoto, E. M., Ogawa, S., and Webster, M. W. (2001). Association of intraluminal thrombus in abdominal aortic aneurysm with local hypoxia and wall weakening. J. Vasc. Surg. 34, 291–299.
Wang, D. H., Makaroun, M. S., Webster, M. W., and Vorp, D. A. (2002). Effect of intraluminal thrombus on wall stress in patient-specific models of abdominal aortic aneurysm. J. Vasc. Surg. 36, 598–604.
Wolters, B. J. B. M., Rutten, M. C. M., Schurink, G. W. H., Kose, U., de Hart, J., and van de Vosse, F. N. (2005). A patient-specific computational model of fluid-structure interaction in abdominal aortic aneurysms. Med. Eng. Phys. 27, 871–883.
A simple example
Consider a series of four reactions involving five species
The reaction rate vector rT = [r1,r2,r3,r4] defines the reaction rates of every single reaction, i.e.
where ki,i = 1, …, 4 denote the reaction rate constant of the i-th reaction. Consequently, the system of reactions reads
where the 5 × 4 stoichiometric matrix S was introduced. The five rows and the four columns in S relate to the five species and four reactions, respectively.
Keywords: abdominal aortic aneurysm, coagulation cascade, computational fluid dynamics, convection-diffusion-reaction equations, intra-luminal thrombus, platelets, thrombin, vortical structures
Citation: Biasetti J, Spazzini PG, Swedenborg J and Gasser TC (2012) An integrated fluid-chemical model toward modeling the formation of intra-luminal thrombus in abdominal aortic aneurysms. Front. Physio. 3:266. doi: 10.3389/fphys.2012.00266
Received: 31 March 2012; Accepted: 25 June 2012;
Published online: 20 July 2012.
Edited by:Rajat Mittal, Johns Hopkins University, USA
Reviewed by:Haoxiang Luo, Vanderbilt University, USA
Daniel Goldman, The University of Western Ontario, Canada
Copyright: © 2012 Biasetti, Spazzini, Swedenborg and Gasser. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.
*Correspondence: Jacopo Biasetti, Department of Solid Mechanics, School of Engineering Sciences, KTH Royal Institute of Technology, Osquars Backe 1, SE-100 44 Stockholm, Sweden. e-mail: email@example.com