## Abstract

A coupled chemo-fluidic computational model for investigating flow-mediated thrombogenesis in infarcted left ventricles (LVs) is proposed. LV thrombus (LVT) formation after the acute myocardial infarction (AMI) may lead to thromboembolic events that are associated with high mortality and morbidity, and reliable stratification of LVT risk is the key to managing the treatment of AMI patients. There have been several studies emphasizing the importance of LV blood flow patterns on thrombus formation; however, given the complex interplay between ventricular flow dynamics and biochemistry of thrombogenesis, current understanding is mostly empirical. In the present model, blood flow in the LV is obtained by solving the incompressible Navier-Stokes equations, and this is coupled to the biochemical modeling of the coagulation cascade, platelet activation, and fibrinogen polymerization. The coupled model is used to examine the effect of ventricular flow patterns on thrombogenesis in modeled ventricles. It is expected that the method developed here will enable in-depth studies of thrombogenesis in patient-derived infarcted LV models.

- computational hemodynamics
- cardiac flow
- intraventricular flow
- blood clot
- coagulation cascade
- biochemical reaction

## NEW & NOTEWORTHY

A computational method that simultaneously resolves the left ventricular flow patterns as well as the biochemical coagulation cascade (CC) is proposed. The method enables a clear identification of the effect of flow patterns on left ventricular thrombogenesis and suggests a metric quantitatively correlated to thrombogenic risk.

one of the most feared complications for patients with acute myocardial infarction (AMI) is the occurrence of thromboembolic events associated with left ventricular (LV) thrombus (LVT) formation (11) and cardioembolic strokes, which account for 20–30% of all ischemic strokes (25, 27, 33). The identification of high-risk patients and the pharmacological prevention of LVT formation are the keys to preventing these embolic events. Thrombus formation results from the combination of disturbed blood flow patterns, endothelial injury, and hypercoagulability of the endothelium. After AMI, LV regional wall akinesia or dyskinesia changes the patterns of blood flow in the LV, and this can result in local stasis of blood flow (3). Prolonged ischemia also leads to subendocardial tissue injury with inflammatory changes and may lead to the development of a LV aneurysm (30, 34), which may cause further abnormalities in LV flow. Furthermore, patients with acute coronary syndrome show a hypercoagulable state with increased concentrations of prothrombin and fibrinopeptide (11). The above prerequisites in the LV with AMI lead to the formation of LVT composed of fibrin, platelets, and red blood cells.

Risk factors for LVT formation include infarct size, degree of apical dyssynchrony, and abnormalities in the blood flow pattern. In general, the incidence of LVT formation is higher with lower LV ejection fractions (EF) (7), which is generally accompanied by dilation of LV [increase of end-diastolic volume (EDV)] and/or the decrease in the stroke volume (SV) due to the infarction. Stratification of patients at risk for LVT formation is currently based primarily on a global assessment of ventricular function and image-based assessment of ventricular wall motion to identify LV regional wall akinesia and dyskinesia. Some previous studies have, however, shown that global parameters, such as EF, SV, and EDV, do not always correlate well with LVT formation (29, 40, 45), and existing approaches for stratifying LVT risk, therefore, lack sensitivity and specificity (44). The identification of high-risk patients and the pharmacological prevention of LVT formation are the key to preventing embolic events. Once an “at risk” patient is identified, treatment with heparin and then warfarin reduces the risk of LVT formation (43) and results in an 86% reduction in the rate of embolization (44). However effective, anticoagulation is not without risk. Most patients with an anterior ST elevation myocardial infarction (MI) in the current era undergo intracoronary stent implantation requiring “triple therapy” with warfarin, aspirin, and clopidogrel. Triple therapy significantly increases the fatal and nonfatal bleeding risk to 22.6% in the first 30 days and remains above 20% up to 90 days. Current risk stratification methods, therefore, result in the treatment of 10 patients to prevent LVT in 1 patient (22), while the risk of bleeding in the triple therapy era is 1 in 5. Clearly, more effective methods of risk stratification are needed.

Several studies have emphasized the importance of LV blood flow patterns on thrombus formation (3, 10, 23, 29, 40, 45), and it has even been suggested that the direct assessment of LV flow for the identification of abnormal flow patterns (10, 29, 45), or a metric closely related to the LV flow (40), could be a better marker for LVT formation. However, the correlation between the flow pattern and LVT risk reported in earlier studies is rather descriptive or qualitative. A quantitative metric related to the flow pattern in the LV and its correlation to LVT formation has been suggested (40), but the limitation of this prior study was that the flow patterns in LVT cases were assessed after thrombus formation, making it difficult to determine whether the metric represents conditions conducive to the genesis of LVT or the effect of the LVT on the LV flow.

The objective of the present study is to develop a computational method to investigate how flow patterns couple with the biochemistry of thrombogenesis in infarcted LVs to generate LVT. A coupled chemo-fluidic computational model is developed to simulate flow-mediated thrombogenesis in infarcted LV models derived from cardiac computerized tomography (CT) scan data. Blood flow in the LV is simulated by solving the incompressible Navier-Stokes equations using the immersed boundary method (37, 39), and this is coupled to a three-component biochemical model, which models the CC (4), fibrin formation (32), and platelet activation (28). The key feature of the current computational model is that the complex flow patterns, as well as the transport and accumulation of key cellular and biochemical factors associated with thrombogenesis, can be tracked simultaneously, enabling a clear identification of the effect of flow patterns on thrombogenesis. Based on these simulations, a metric that is quantitatively correlated to thrombogenic risk is also proposed, and this could contribute to the improvement of methods for LVT risk stratification.

## METHODS

#### Infarcted LV models.

Models of infarcted LVs that are typical of ST elevated, antero-apical infarcts are constructed for the simulation of blood flow and thrombus formation. First, a relatively simple baseline model of the LV is constructed using a three-dimensional (3D), high-resolution contrast CT scan of a normal human ventricle (37, 38). Given the prevalence of LV aneurysms in ST-elevated MIs and their association with thrombogenic risk, we employ three different infarcted LV geometries with a small (grade III), medium (grade IV), and large (grade V; Ref. 30) apical aneurysm, as shown in Fig. 1*A*. These aneurysms are generated “in silico” by artificially dilating the apical region of the baseline computational model. Since the EF is the other critical factor connected with thrombogenic risk, we have considered three different EFs (28, 35, and 37%) and the parameters for the six final cases are summarized in Table 1. These cases are designed to study the effects of EF, SV, LV volume, and aneurysm size on LV hemodynamics and thrombus formation. *Cases A–C* enable the study of the effects of SV and EF for the fixed end-systolic volume (ESV) with a large apical aneurysm, and *cases E* and *F* are for the small aneurysm. *Cases C*, *D*, and *F* enable the study of the effects of aneurysm size and SV for a fixed EF. Given the significant computational effort required for these simulations (each simulation requires about 30,000 central processing unit hours), a larger parameter study is quite difficult at the current time. The chosen set of cases, however, cover a fairly large range of primary LV parameters (SV: 27–62 ml; EF: 28–37%; and ESV: 70–105 ml) to enable a detailed study of LVT formation.

The heart rate is assumed to be 60 beats/min, a reasonable goal heart rate in a post-MI patient, and the flow rate, Q(*t*) (plotted in Fig. 1*B*) through the LV is based on published data (9). The wave forms for the early (E-) and atrial (A-) waves are based on a simple harmonic oscillator model, which has previously been verified with Doppler echocardiographic data (18, 19), and the E/A ratio is set to 1.2 (5). The flow rate wave forms normalized by the maximum flow rate through the mitral inlet (Q_{max}, see Table 1) are the same for all cases. Further details regarding the flow rate specification are given in appendix a. The flow into and out of the ventricle is driven by a prescribed expansion and contraction of the LV (37). The wall motions follow the descriptions given for the LV with aneurysm in Ref. 30. The motion of the apical wall is suppressed, representing the apical akinesia due to the aneurysm, and thus the volumetric expansion and contraction of the LV are mainly in the radial directions (*x* and *y*) rather than in the axial direction (*z*).

The ESV of the ventricle models with the large, medium, and small aneurysm are 105, 87, and 70 ml, respectively. The length of LV in the axial direction (*z*), *L*_{LV}, is 8.4 cm for the large-aneurysm model, 7.9 cm for the medium-aneurysm model, and 7.3 cm for the small-aneurysm model. The apical region of the LV wall inside the aneurysm is assumed to be damaged due to an acute infarction, as indicated in Fig. 1*A*, and the length of damaged wall region in the axial direction (*L*_{DM}) is ∼1.4 cm for all models.

The present LV models include a relatively simple model for the mitral valve (MV), which is shown to produce realistic LV flow patterns (38). The shape of the MV leaflets is based on the detailed anatomical measurements of Ranganathan et al. (35), and the motion of the MV leaflets during diastole is modeled by prescribing the opening angles (angles between the leaflet and the mitral annulus plane) for anterior (AL) and posterior leaflets (PL), ϕ_{AL} and ϕ_{PL}, respectively, as indicated in Fig. 1*A*. Further details of this motion are provided in appendix b.

#### Modeling of blood flow.

The incompressible Navier-Stokes equations are employed to simulate the dynamics of blood flow, since the blood in the ventricle can be assumed to be a Newtonian fluid. The equations are written as:
(1)where is a velocity vector, P is pressure, and ρ_{0} and ν_{0} are the density and kinematic viscosity of the blood, respectively. The equations are solved on the non-body-conformal Cartesian grid, and the moving boundary of the LV endocardium and MV leaflets is treated by a sharp-interface immersed boundary method (31). The current flow solver has been successfully applied to cardiac flows in a number of previous studies (8, 37–39, 47–49) and has been shown to generate flow features that are inline with in vivo observations (38). A recent validation study (46) for a model ventricle (15) with the same procedure also confirmed that our computational modeling procedure does produce the key features of the velocity and vorticity fields observed in ventricular flows. These previous simulations also provide guidance on the resolution required to model these flows.

In the present study, the endocardial wall and the MV leaflet surfaces are discretized with a total of 37,138 and 7,200 triangular elements, respectively, and a total of 128 × 128 × 256 (4.2 million) grid points are employed to resolve the flow (37). The time step size used is 1 × 10^{−4} s, which is more than sufficient to ensure high-temporal accuracy. The flow Reynolds number based on the average peak E-wave velocity and mitral annulus diameter ranges from 1,700 (SV = 27 ml) to 3,900 (SV = 62 ml), and the Womersley number based on the mitral annulus radius and heart rate is 14. These nondimensional parameters are in the range considered typical for an adult human heart (49).

#### Coagulation modeling.

For the modeling of blood coagulation and thrombus formation in the infarcted LV, we adopt thrombosis biochemical reaction models employed in previous studies (4, 28, 32) of vascular thrombosis. It is assumed here that thrombus formation in the LVs with AMI is initiated by the exposure of blood to tissue factor (TF) (42) emanating from the damaged ventricle wall. Once TF is exposed to the bloodstream, it combines with the active form of coagulation, factor VII, and this TF:VIIa complex initiates the CC (16). This extrinsic, TF pathway of CC is common in damaged vessels, and the biochemical reactions in this procedure are schematically shown in Fig. 2. To model the reactions between coagulation factors in Fig. 2, we employed the reaction model proposed by Biasetti et al. (4). This reaction model consists of a total of 16 reactions between 18 biochemical species and a detailed description of these reactions can be found in appendix c. Since those biochemical species are transported by the blood flow, we solve the convection-diffusion-reaction (CDR) equations to resolve the transport and evolution of the concentrations. This equation can be written as: (2)

where *C*_{i} is the concentration of *i*th species, *D* is the diffusion coefficient, and *R*_{i} is the reaction source term. The primary effect of flow on the biochemical reaction model is the transportation of species by the flow and the CDR equation (*Eq. 2*) is coupled with the hemodynamic simulation through the flow velocity vector . The diffusion coefficient *D* for the coagulation factors is set to 1.0e-8 m^{2}/s following Biasetti et al. (4). The detailed description of the reaction term *R*_{i} and the initial concentrations of the species are given in appendix c. At the wall, a zero diffusive flux boundary condition is applied except in the damaged region, where the TF is exposed and combined with factor VIIa. At the damaged wall, a Dirichlet boundary condition with a constant concentration (equal to 1e-4 mol/m^{3}) is applied for the TF:VIIa complex following Biasetti et al. (4).

As a result of the CC shown in Fig. 2, thrombin (factor IIa) is produced. Thrombin plays an important role in thrombus formation (1, 16). First, thrombin activates key coagulation factors in the CC (as shown in Fig. 2) and accelerates its own production. Second, thrombin also initiates the activation of platelets (1), which is another key cellular component of blood clots. Finally, thrombin converts the fibrinogen present in blood plasma into fibrin, which is a fibrous protein and a key component of thrombii (1). Therefore, in this study, we simulate the thrombin-based activation of platelets, as well as the polymerization of fibrinogen into fibrin.

For modeling the polymerization of the fibrinogen, we have employed the CDR equation model proposed by Neeves et al. (32). The evolution equations for the fibrinogen and fibrin concentrations are given by:
(3)
where *C*_{f} is the concentration of fibrinogen, *C*_{m} is the concentration of fibrin, and *k*_{t} is the conversion rate which is defined by:
(4)where *C*_{IIa} is the concentration of thrombin, and the other parameters are given (32) by *k*_{cat} = 84 s^{−1}, *K*_{m} = 7.2e-3 mol/m^{3}, and *k*_{p} = 8.2e2 mol·m^{−3}s^{−1}. Following Neeves et al. (32), the initial concentration of fibrinogen in the blood plasma is set to 9e-3 mol/m^{3}.

Platelets activated by thrombin may deposit on and attach to the damaged wall and also aggregate with other platelets to form a blood clot. In the present study, these procedures are modeled by using the evolution equations of the platelet number density proposed by Leiderman and Fogelson (28). The original equations of Ledierman and Fogelson are slightly modified, and we separate the platelets into three groups based on their status: *1*) inactivated, mobile platelets; *2*) activated, mobile platelets; and *3*) platelets bound to the wall or the other bound platelets. The evolution equations for the number densities of these platelet groups are written as:
(5)
(6)
(7)where *PT*_{m,u}, *PT*_{m,a}, and *PT*_{b} are the number densities of the inactivated, activated, and bound platelets, respectively; *PT*_{max} is the maximum platelet number density allowed in the specific volume and the reported value (28) is 6.67e7 mm^{−3}; and *D*_{P} is the diffusion coefficient which is set to 1.e-9 m^{2}/s. The term *A*_{IIa} is the activation rate by the thrombin and defined as (28):
(8)where *k*_{u} = 0.5 s^{−1} and *K*_{IIa} is set to 1.e-6 mol/m^{3}. The third and fourth terms on the right-hand side of *Eq. 6* are the sink terms associated with the adhesion of activated platelets to the damaged wall and the cohesion to the bound platelets, respectively, and these terms correspond to the source terms on the right-hand side of *Eq. 7*. In *Eqs. 6* and *7*, is a function indicating the proximity to the damaged wall and has a value of 1 adjacent (within 0.5 mm) to the damaged wall and 0 otherwise. The function *g*(*PT*_{b}) evaluates the number density of nearby bound platelets within one computational grid point by:
(9)where *i*, *j*, *k* are the grid indexes. The adhesion rate coefficient, *k*_{adh}, and the cohesion rate coefficient, *k*_{coh}, are both set to 5,000 s^{−1}. The number density of inactivated platelets in the blood is set to 2.5e5 mm^{−3} (28), and this value is used for the initial condition of *PT*_{m,u}, while the initial values for the other platelet groups are 0.

All of the evolution equations in the coagulation modeling (*Eqs. 2*, *3*, *5*, and *7*) are discretized with a second-order finite difference method and integrated in time using a four-stage Runge-Kutta method. These equations are solved along with the incompressible Navier-Stokes equation, *Eq. 1*, since they require the hemodynamic velocity field. appendix d describes our study of the sensitivity of the current CC model results to the initial concentration of key chemical factors, and the data indicate that the overall trends are quite insensitive to these factors.

#### Flow metrics for assessing LV hemodynamics.

A general notion regarding the effect of flow patterns on thrombus formation is that flow stasis (or very slow flow speed) in the LV will increase the chance of thrombus initiation (3, 29). This is because flow stasis allows sufficient time for the biochemical reactions between the coagulation factors and the damaged wall to occur. Thus, in the present study, we consider a number of characteristic quantities and metrics that can characterize the effect of ventricular hemodynamics on thrombus formation.

The formation as well as the propagation of vortexes in the ventricle has significant implications for ventricular hemodynamics. Healthy LVs generate vortexes that are by now well characterized (13, 15, 37), and the flow and dissipation induced by ventricular vortexes are expected to play a key role in modulating the conditions for thrombogenesis in infarcted ventricles. In the present study, 3D vortical structures are identified by the well-established Q-criterion (20). The effectiveness of vortex formation during the early (E-wave) and atrial (A-wave) fillings can also be characterized by the vortex formation number (17), which is defined by the length-to-diameter ratio of the ejected fluid column. For mitral inflow jet, this number can be calculated by:
(10)where VF is the vortex formation number, V_{filling} is the filling volume associated with E- or A-wave, *A*_{mitral} is the mitral inlet area, and *D*_{mitral} is the mean diameter of the mitral inlet. It is suggested that VF = 4 is optimal for the growth of vortex ring, and it has been reported that VF of the E-wave for normal LVs ranges from 3.3 to 5.5 (17). Another straightforward quantity that characterizes the flow dynamics in the LV is the time-averaged flow velocity and associated streamline pattern and vorticity field. In the present study, this time average is computed by averaging the blood flow velocity in the LV over a cardiac cycle. To assess the flow strength in the apical region near the damaged wall, the apical flow strength, *V*_{apical} is defined by the average velocity magnitude and calculated by:
(11)The volume integral is taken over the LV region, but the axial direction (*z*) is from the apex (*z*_{apex}) to the end of damaged wall region (*z*_{apex} + *L*_{DM}).

An additional metric that is considered in the current study is the wall shear stress (WSS) or wall shear rate (WSR). It is well established that the WSS or WSR may be responsible for the remodeling of the vascular endothelium, and WSR has been correlated with the location of vascular thrombus in the previous study (36). The WSR clearly indicates the strength of near-wall flow relative to the wall, since it is defined by the flow velocity gradient at the wall. Thus low WSR regions should be correlate with increased propensity for thrombogenesis (36). The time-averaged WSR is defined by:
(12)where *n* denotes the wall normal direction, and *T* is the duration of a cardiac cycle. Note that the WSR can easily be converted to WSS by multiplying by the viscosity of blood.

The final metric employed here is the residence time (RT), which is defined as the time for which the blood resides in a given region of interest (ROI). The RT can be used as more direct measurement of flow stasis. Stagnant or low-flow speed in the ROI should result in a high RT. It is noted that the RT can also be considered as a direct measure of the time available for the coagulation associated biochemical reactions in the ROI, and thus regions with high RT should also correlate with the thrombus location. In the present study, we employ the Eulerian method proposed by Esmaily-Moghadam et al. (14) for the calculation of RT. The ROI is set to be the region adjacent (within 0.5 mm) to the damaged wall, since the initiation of CC and the binding of activated platelets all occur in that region. We call this metric as the near (damaged) wall residence time (NWRT), and it is obtained by solving the following equation concurrently with the flow equations: (13)where τ is the NWRT and is the function, indicating the proximity to the damaged wall ROI and has the value of 1 in the ROI and 0 otherwise.

#### Correlation coefficient.

To quantify the spatial correlation between the hemodynamic and coagulation metrics, correlation coefficients are calculated between pairs of metrics over the damaged LV wall:
(14)where *r* is the correlation coefficient, φ is any hemodynamic or coagulation metric, and *A*_{s} is the area of the damaged wall. The surface integral is taken over the damaged wall (see Fig. 1*A*).

## RESULTS

The coupled chemo-fluidic simulations of thrombogenesis in each of the six infarcted LV models are carried out for a total of six cardiac cycles. While the flow reaches a repeatable state by the third cycle, the additional cycles allow us to obtain reliable statistics on the growth rate of the various biochemical components of the CC. These simulations are conducted on a multiprocessors cluster and simulation of one cardiac cycle takes about 6,000 central processing unit hours. We note here that the computational expense of these simulation prohibits significant longer time integration, and the current models are best suited for examining the very early trends in LVT formation.

#### Flow patterns.

The time evolution of the vortex structures for the six infarcted cases is plotted in Fig. 3. For *case A*, the SV is 62 ml and the E-wave vortex formation number (*Eq. 10*) is 4.1. Consequently, the formation of the vortex ring below the MV is clearly visible at *t* = 0.12 s (around the peak of E-wave). The vortex ring generated by the E-wave separates from the valve lip and propagates down to the apex. Due partly to the interaction with the LV wall, the vortex ring breaks down into smaller vortex structures in the apical aneurysm (*t* = 0.36 s) and then slowly dissipates due to viscous effects (*t* = 0.48 s). The A-wave jet generates a similar but weaker vortex ring that is visible around *t* = 0.6 s.

*Case B* has the same ESV and aneurysm size as *case A*, but the SV is slightly reduced to 50 ml, and this decreases the E-wave vortex formation number to 3.3. The vortex ring structure is still visible at the peak of the E-wave (*t* = 0.12 s), but it is not as strong as *case A*. The overall flow structure is similar to *case A*, but the vortical strength is slightly weaker than *case A* due to the decreased SV.

*Case C* also has the same ESV and aneurysm size as *cases A* and *B*, but the SV is further reduced to 41 ml, and the E-wave vortex formation number is 2.7. The vortex ring structure is, therefore, not clearly visible at the peak of the E-wave (*t* = 0.12 s), but, by the end of the E-wave (*t* = 0.24 s), the vortex ring does form, and it propagates downward to the apex. The vortex ring, however, starts to break down in the middle of ventricle, before it reaches the apex (*t* = 0.36 s), and the resulting small-scale vortex structures dissipate rapidly there (*t* = 0.48 s). Overall, the strength of circulatory flow is clearly weaker than the previous cases. The A-wave generates a vortex ring structure at the MV leaflet tips, but it is not clearly detached from the tip until the end of the A-wave (*t* = 0.6 s).

*Case D* has a medium-sized apical aneurysm, and thus the ESV is smaller than in *cases A–C*. The SV is reduced to 34 ml to maintain the same EF as in *case C*. The VF number for the E-wave is 2.25, and, therefore, the vortex ring structure is barely visible at the peak of E-wave (*t* = 0.12 s). The vortical flow of the E-wave jet is observed at the end of the E-wave (*t* = 0.24 s), but, similar to *case C*, it starts to break down in the middle of ventricle (*t* = 0.36 s). Although *case D* has a smaller aneurysm than *case C*, the vortical flow does not reach the apex before it dissipates.

The ESV as well as the aneurysm size for *case E* are smaller than the above cases, but the SV (and therefore the vortex formation number) is the same as for *case C*. Therefore, the overall vortex ring formation and propagation patterns for *case E* are almost identical to those of *case C*, as one can see at *t* = 0.12 and 0.24 s. Similarly, the vortex ring starts to break at around *t* = 0.36 s, but, since *case E* has the smaller LV size, the vortex structures reach in to the apical aneurysm before dissipating (*t* = 0.48 s). The A-wave vortex ring structure is again very similar to that for *case C* (0.6 s).

*Case F* has the same ESV and aneurysm size as *case E*, but the SV is further decreased to 27 ml. The vortex formation number for the E-wave is only 1.78, and thus the vortex ring structure is not clearly visible during the E-wave due to lower early diastolic filling. A ring-like vortex structure is observed at the end of the E-wave (*t* = 0.24 s) at around the middle of the ventricle, but the structure soon breaks (0.36 s) and dissipates (0.48 s) with relatively little penetration into the ventricle. The vortex structure for the A-wave is formed following the MV leaflets, but it does not develop into a vortex ring (0.6 s).

The time-averaged streamline pattern and the apical flow strength are shown in Fig. 4. Figure 4*A* shows the velocity magnitude averaged over a cardiac cycle for the *cases A–F*. The averaged velocity is high along the line of mitral jet propagation and also near the outflow tract toward the aorta, while, due to the overall swirling flow pattern, the velocity magnitude at the center of the LV is low. The apical flow strength, i.e., the velocity magnitude inside the apical aneurysm, is clearly different for the cases considered. To assess the overall strength of the flow in the apical region, the apical flow strength, *V*_{apical} is calculated by *Eq. 11* over the region near the damaged wall (under the dashed line in Fig. 4*A*) and presented in Fig. 4*A* for all cases. With the highest SV, *case A* shows the strongest apical flow (*V*_{apical} = 11.3 cm/s), and we observe the mitral jet smoothly connecting the mitral inflow and the flow toward the outflow tract. For *case B*, which has a reduced SV, the *V*_{apical} is decreased to 7.3 cm/s. For all other cases, *C*–*F*, however, *V*_{apical} is as low as ∼4 cm/s, and the localized zone of nearly stagnant flow is observed along the septal wall inside the apical aneurysm. For *case F*, this low speed zone at the septal-apical wall extends to the middle of the LV.

The streamline patterns and vorticity for the time-averaged flow are shown in Fig. 4*B*. For all of the cases, an overall clockwise circulating (negative vorticity) flow pattern is observed. The mitral inflow is deflected to the lateral wall and then redirected toward the outflow tract inside the apical aneurysm. Complete flow stasis in the aneurysm is not observed for any of the case, and thus the apical blood pool is being “washed out” by the flow for all the cases. However, as evident in Fig. 4*A*, the overall flow speed in the aneurysm is quite different between the cases, and this will result in different degree of washout inside the aneurysm. For the low apical flow strength, blood constituents will reside in the proximity of the damaged LV wall inside the aneurysm for a longer time, and this will increase the accumulation and concentration of the procoagulant biomolecules. As the septal-apical wall region inside the aneurysm is found to have low-velocity magnitude for the *cases C–F*, this region may be subjected to a higher propensity for coagulation and thrombogenesis.

#### Thrombogenesis.

The simulation of the CC and the fibrin and platelet activation and deposition are performed along with the hemodynamic simulations for a total of six cardiac cycles. On the damaged LV wall inside the aneurysm indicated in Fig. 1*A*, the TF:VIIa complex initiates the reactions, and the coagulation factors accumulate as the reactions continue. At the same time, these biochemical species are transported by flow, primarily along the directions indicated by the streamlines. As indicated in Fig. 4, the blood volume inside the aneurysm is washed out by the flow and finally ejected through the outflow tract. If the flow speed near the damaged wall is high, the local blood volume is washed out rapidly, and this limits the accumulation of the coagulation factors and the deposition of activated platelets. On the other hand, if the flow speed is low, the blood constituents reside longer near the damaged wall, resulting in higher concentrations of procoagulant factors. Therefore, the concentration and distribution of the coagulation factors are expected to depend strongly on the flow field in the vicinity of the damaged LV wall.

The growth of the volume-averaged concentrations of thrombin, fibrin, and the number density of activated and bound platelets are plotted in Fig. 5. The values are normalized by the initial concentration of their preactivation form, respectively, i.e., prothrombin for the thrombin, fibrinogen for the fibrin, and inactivate platelet for the activated (and bound) platelets. Since it is observed that the coagulation factors grow exponentially in time, we take a logarithm to plot the time histories. Regression analysis indicates that the growth of these quantities is best fit to the exponential function: (15)

where *C** is the normalized concentrations, and *K* and *B* are constants. The line plots in Fig. 5 are the results of curve fitting to *Eq. 15*, and *R*^{2} for all of the regression fits are >0.9.

It is interesting to note that, if we take a time derivative of *Eq. 15*, the growth of normalized coagulation factor concentrations can be expressed as a system of two evolution equations with a single rate coefficient:
(16)where ϕ* can be considered as a reaction potential, which is abundant initially and decreases as reaction continues. *Equation 16* can be adopted as a model of a reduced-order reaction for the volume-averaged coagulation factors, and the constant *K* (unit of s^{−1}) represents the overall reaction rate coefficient. *Equations 15* and *16* can also be used to extrapolate the coagulation results over a longer elapsed time. The reaction rate *K* for the thrombin, fibrin, and platelets are listed in Table 2 for all cases. Apparently, different LV conditions and flow patterns for the cases result in a different value of model reaction rate, *K*. The fibrin and activated platelets are produced directly by the thrombin, and thus the changes in the model reaction rate of fibrin and platelet for different cases follow the same trend as thrombin. A higher *K* indicates a faster growth of the corresponding coagulation factor, and thus the value of *K* in Table 2 may be used to determine the overall thrombosis risk.

While the overall growth of the coagulation factors in the LV volume is quantified in Fig. 5, the localized accumulation of those coagulation factors is also important for the onset of thrombus formation. Therefore, the spatial distribution of thrombin concentration averaged over the last cardiac cycle is examined in Fig. 6 for all cases.

Figure 6 shows that thrombin is mainly concentrated near the septal-apical wall inside the aneurysm for all of the cases. *Cases A–C* show a similar distribution of thrombin with high concentrations at the septal wall inside the large aneurysm, although the maximum concentration for *case C* is about five orders of magnitude higher than for *case A*. For *case D*, higher concentration of thrombin in the septal-apical region is spread over a larger volume; the volume-averaged concentration is consequently of the same order of magnitude as *case C*. For *case E*, the distribution of thrombin is highly localized near the apex due to the smaller aneurysm size, and the maximum thrombin concentration is about three orders of magnitude higher than for *case A*. For *case F*, thrombin is distributed over a large area, which is caused by an extended region of low-flow velocity, as shown in Fig. 4. The local maximum thrombin concentration for *case F* is about an order of magnitude smaller than for *case C* and *D*, although the overall growth rate of the volume-averaged concentration is comparable to those cases (see Table 2).

#### Quantitative correlation between the flow and coagulation metrics.

As described in the above sections, the accumulation of coagulation factors seems to be related to the apical flow strength and the resulting washout. In this section, we examine and quantify the correlation between the flow and coagulation patterns. The apical flow strength near the damaged wall is quantified using two hemodynamic metrics: WSR (*Eq. 12*) and NWRT (*Eq. 13*). These quantities are calculated as described in the methods section. WSR is proportional to the magnitude of the flow velocity relative to the adjacent wall motion; it, therefore, measures the flow strength near the wall. The NWRT measures how long blood volume resides near the wall, so it can be used for the evaluation of the degree of flow stasis as well as washout quality near the wall. A lower apical flow speed will result in a higher (1/WSR) and a higher NWRT, which will correlate with greater accumulation of coagulation factors. Note that both 1/WSR and NWRT are expressed in units of time (i.e., seconds). The WSR is a metric that has to be determined on the surface, while the NWRT can be evaluated both on the surface and in the volume. To investigate the correlation between these hemodynamic metrics and the coagulation results, we focus on the damaged LV wall region inside the aneurysm. For the coagulation metrics, we pick the thrombin concentration (which is the primary contributor to the thrombus formation) and the bound platelet number density (the direct constituent of the thrombus). The bound platelets are the activated platelets that are bound to the damaged wall and can, therefore, be evaluated easily on the infarcted LV wall.

The distributions of hemodynamic metrics (1/WSR and NWRT) and the coagulation results (thrombin concentration and bound platelets number density) averaged over the last cardiac cycle on the apical region are plotted in Fig. 7 for all of the cases using a cardiac polar plot. One can see that two hemodynamic metrics, 1/WSR and NWRT, have very similar distributions. Similarly, the distributions of the thrombin and the bound platelets are almost identical, since the platelets are directly activated by thrombin in the present model. The locations of the maximum values for both hemodynamic and coagulation metrics are found to be shifted toward the septal direction, as noted earlier, but the plots also show a bias toward the anterior direction. While there are differences in the distributions of the hemodynamic and coagulation metrics, the peak locations of the hemodynamic metrics match relatively well with the peak locations of the coagulation factors, especially for the cases with a higher accumulation of coagulation factors (*cases B–F*). The spatial correlations between the hemodynamic and coagulation metrics are quantified by the correlation coefficients (*Eq. 14*), and the calculated correlation coefficients are listed in Table 3.

The spatial correlation between 1/WSR and coagulation results are moderate for the cases with high local accumulation (*C* and *D*) (*r* = 0.53–0.68), but low for other cases (*r* = 0.05–0.46). In contrast, the correlations between NWRT and the coagulation metrics are moderate to high for all of the cases (*r* = 0.55–0.79). This is because, while the WSR is only proportional to the local relative velocity magnitude, the NWRT is directly related to the washout quality and the time allowed for reaction.

To assess the proportionality between the hemodynamic and coagulation metrics, the maximum and averaged hemodynamic and coagulation metrics over the damaged wall are listed in Table 4 for all of the cases. Overall, the hemodynamic and coagulation metrics are positively correlated; the higher 1/WSR and NWRT values correspond to a higher thrombin concentration and bound platelet number density. As expected, *case A*, which has the lowest thrombin and bound platelet accumulation, shows the lowest 1/WSR and NWRT as well. The peak NWRT for *case A* is just 0.4 (s), which means that the blood volume stays near the damaged wall no longer than one-half of a cardiac cycle. This confirms the higher degree of apical washout and explains the low accumulation of coagulation factors. The hemodynamic metrics for *case C* are about three to five times higher than for *case A*. The peak NWRT is about 2 (s), which indicates a much degraded apical washout. Interestingly, the coagulation metrics for *case C* are about five orders of magnitude higher than for *case A*. *Case D* shows almost identical NWRT and coagulation metrics with *case C*, while those two cases have different LV geometries with different aneurysm sizes. *Cases B* and *E* also display similar behavior, but the peak thrombin concentration is about two orders of magnitude lower than for *cases C* and *D*. Interestingly, *case F* has a similar mean NWRT as that for *cases C* and *D*, but the peak NWRT is about one-half, and thus the peak thrombin concentration is about one order of magnitude lower. These results show that, although higher 1/WSR and NWRT values are indicative of a greater accumulation of coagulation factors, the actual relation between the hemodynamic and coagulation metrics is not linear. The reason for this is that the produced thrombin itself also activates the key coagulation factors (factors V and VIII) for further thrombin production, as shown in Fig. 2, and, once these factors are abundant, the production rate of the thrombin, as well as the species activated by thrombin, is increased significantly by those factors.

To assess the relationship between the allowed local reaction time and the hemodynamic metrics (1/WSR and NWRT), the relationship between the peak value of the coagulation factor concentrations and the peak value of hemodynamic metrics is subjected to regression analysis. If the two hemodynamic metrics are related to the time allowed for reactions to occur, the relationships should correspond to the expression of *Eq. 15*, which describes the function for the temporal growth of the coagulation factors. The fitted curves for the 1/WSR and NWRT vs. the peak thrombin concentration are plotted in Fig. 8, *A* and *B*. Both plots show that the present data fit reasonably well to *Eq. 15*, and this suggests that NWRT or 1/WSR indeed represents the reaction time allowed in the ROI. The goodness of fit is better for the NWRT (*R*^{2} = 0.91) than 1/WSR (*R*^{2} = 0.8). We note that the peak thrombin concentration does not correlate well with LV parameters such as SV, EF, and LV size (ESV and EDV), which is reflective of the fact that SV and EF can only provide a gross and, therefore, imprecise measure of LV hemodynamics, especially vis-à-vis the degree of flow strength in the apical region, which is key to LV thrombogenesis.

Given the fact that 1/WSR and NWRT both correlate with the local reaction time, the average value of these metrics should have positive correlations with the overall reaction rate presented in Table 2. The relationship between these variables is plotted in Fig. 8, C and D, and they show a reasonably clear linear relation, especially for NWRT (*R*^{2} = 0.95). This also suggests that these hemodynamic metrics may be used to estimate the overall LVT risk.

## DISCUSSION

Several studies have emphasized the role of flow pattern on LVT formation (10, 29, 40, 45). These studies have suggested that weak flow strength or the flow stasis in the apical region results in thrombus formation. The correlations between the flow pattern and the thrombus formation presented in those studies were somewhat qualitative, since it was difficult to measure simultaneously the details of the flow patterns and LVT risk in such in vivo settings. In the present study, we have performed high-fidelity computational simulations with simple but representative models of infarcted ventricles. The computational model includes both hemodynamics and biochemical reactions for the coagulation, and the results, therefore, allow us to investigate the detailed correlation between the flow patterns and the thrombus formation.

The results indicate that the accumulation of coagulation factors is primarily due to flow stasis and the associated RT of blood near the infarct site in the apical region, and this is in line with the observations from previous studies (3, 29). Weak flow near the infarction location enables the blood constituents to reside longer in that region, allowing sufficient time for the generation of active coagulation factors. The coagulation factors generated are also not transported away from the infarct location due to weak flow, and the concentration increases exponentially in time (see Fig. 5), since the produced coagulation factors continuously activate other factors. This increased concentration eventually leads to activated fibrin and platelets, which would tend to form a thrombus in that region. Thus the local thrombus risk directly depends on the local washout quality associated with the flow velocity. We have observed that the coagulation factors are mainly concentrated near the septal-apical wall inside the aneurysm for all of the present cases. This is because, although the TF:VIIa is available all over the damaged LV wall inside the aneurysm, the produced biochemical species are transported from the lateral to the septal direction by the overall clockwise swirling flow pattern, as shown in Fig. 4. Also, the septal-apical region generally has a low-flow velocity magnitude, and thus a longer RT is allowed here for reactions to occur.

Regression analysis indicates that the complex multistep reaction kinetics are well represented by a double-exponential curve, suggesting that low-dimensional models of the reaction kinetics could eventually be employed for representing the CC. The key factor in this lower dimensional model is the model reaction coefficient *K* (Table 2) that represents overall growth rate of the coagulation factors. The model reaction coefficient *K* for *case A* is much smaller than for other cases, although it has the large aneurysm. This is because of the effective washout of the apical blood volume near the damaged wall by the fast apical flow, as examined in Fig. 4*A*. *Case C* has the same aneurysm size as *case A*, but SV and EF are reduced, and thus the apical flow strength is also decreased. As a result, the apical washout quality is degraded, and this increases the growth rate. *Case D* has a smaller aneurysm than *case C*, but the same EF and the apical flow strength for *case D* is the lowest, resulting in the highest reaction rate. *Case E* has the same EF as *case A*, and the aneurysm size is smaller than that in *case A*, but the model reaction coefficient is noticeably higher than in *case A* due to the decreased SV. *Case F* has the smallest SV and EF among the cases. Although the aneurysm size for *case F* is the smallest, the reaction rate *K* is as high as that in *cases C* and *D*.

It has been reported that the concentration of key coagulation factors could serve as biomarkers for the LVT risk (11). Based on the quantitative coagulation results obtained from the coupled hemodynamic and biochemical simulations, the relative thrombus risk for the cases can, therefore, be determined based on the results shown in Table 2 and Figs. 5 and 6. First, the risk for *case A* would be very low due to the extremely small concentration of thrombus factors and the relatively low overall growth rate. *Cases C*, *D*, and *F* may have a similar thrombus risk based on the comparable overall growth rates. However, *case C* would have a much higher risk for the localized thrombus formation based on the highest maximum concentration value. The risk for *cases B* and *E* is judged to be moderate. These preliminary evaluations of the relative thrombus risk based on simulated coagulation factor concentrations clearly show that thrombus risk does not simply depend on global parameters such as SV, EF, and aneurysm size. Instead, thrombus risk is driven by apical washout quality, which may not have a simple relationship to these global parameters.

The overall flow pattern resolved by the simulations can be characterized by the big and slow circulation, which is centered near the middle of the ventricle. This pattern is very similar to the one observed by Beppu et al. (3) for ventricles with apical akinesia using contrast echocardiography. However, the flow strength in the apical aneurysm is very different from case to case, and depends significantly on the strength of the mitral jet during diastole and its propagation to the apex. This difference in apical flow strength leads to significant changes in the production of coagulation factors. The present simulation results for various SV and EF show that the strength of the mitral jet is primarily proportional to the SV, but, even for the same SV, the local washout quality at the apex can differ based on the LV and aneurysm size. The time evolution of the vortex structures (Fig. 3) shows that the strength of the vortex ring and its propagation speed into the apex depends primarily on the SV. A higher SV results in a faster penetration of the vortex into the ventricular cavity. If the initial vortex ring propagates further down toward the apex, it will increase the flow strength in the apical aneurysm. However, the flow strength in the apical region is in fact determined by the combination of flow characteristics and the LV geometry. For example, if the LV is significantly dilated, a higher SV is required to make the vortex ring propagate all the way down to the apex. The results of *cases C* and *E* clearly demonstrate this. These cases have the same SV and vortex formation number that result in almost identical vortex ring dynamics, but, since *case C* has a larger LV size due to the larger aneurysm, the vortexes do not reach the apex before they dissipate.

Although EF is the SV normalized by EDV and can represent the global washout quality (12, 37), it is not an effective single metric to represent the “local” washout, as one can see from the results of *cases A* and *E* and *cases C* and *F*. These observations explain why the stratification of LVT risk based on the conventional echo-Doppler parameters, such as SV, EF, etc., is difficult (29, 40, 45). The local flow strength based on the velocity magnitude may better represent the washout quality and thrombus risk, as suggested in the previous study (29), but defining the location and size of the region for the consistent velocity measurement may be difficult. Moreover, it is also difficult to measure the flow velocity very near the wall. For vascular thrombus, Rayz et al. (36) showed a remarkable correlation between the region of high blood RT and thrombus location. This supports the idea that the accumulation of thrombus factors depends on the local washout quality. In the present study, two hemodynamic metrics, 1/WSR and NWRT, are used to quantify the local washout in the infarcted LVs. With the direct simulation of the CC, we were able to extract quantitative correlations between the hemodynamic metrics and the coagulation factors. Based on the correlation with the accumulated coagulation factors, we have found that NWRT is a very good metric for local washout responsible for coagulation.

The quantitative spatial correlation between the NWRT and the coagulation results found in the present computational study suggests that the NWRT could be used for the estimation of possible thrombus location. The NWRT also has a positive correlation with the peak concentration of thrombin, but the actual relation is not linear. The relation between the NWRT and thrombin concentration is well fit to the model of temporal growth of coagulation factors, *Eq. 15*, and this suggested that the NWRT actually represent the local reaction time allowed in ROI. It is also shown that the mean NWRT is linearly correlated with the overall reaction rate coefficient for the averaged coagulation factors. Based on this, the hemodynamic metric, NWRT could be used for the stratification of overall LVT risk. With advanced medical imaging technique, such as contrast echo and CT, and especially echo-particle image velocimetry (24), it may be possible to develop methods to estimate NWRT in vivo in the near future.

A number of study limitations are worth noting. First, the coagulation model is limited to the extrinsic pathway, and the intrinsic (41) and inhibition pathways (2, 6) are not considered. Second, the chemical reaction models employed in this study are based on previous studies, but there could be patient-specific variations in the reaction coefficients, as well as the initial concentration of procoagulant biomolecules. All of the above could potentially change the absolute magnitude of the coagulation factor concentrations obtained from the simulations but, as shown in appendix d, would likely not change the comparative correlation between the various cases simulated here. For instance, inclusion of a fibrinolysis pathway (6) could be analogous to reducing the forward reaction rate of fibrin production and would reduce the final concentration of fibrin predicted by the simulations. However, given the strong correlation that we have found between active coagulation factor concentration and NWRT, which is unrelated to biochemistry, we do not expect that this change in effective reaction rate would change the comparative correlation between the cases. Given that the primary objective of the current study is the investigation of correlations between the LV flow dynamics and the coagulation patterns, rather than the prediction of the exact concentrations of the coagulation factor, the precise choice of model coefficients might, therefore, not be critical here.

It is also noted that, while the present simulations are carried out for six cardiac cycles due to the excessive computational time for the coupled hemodynamic and biochemical systems, the formation of an actual thrombus occurs over a much longer time. Thus the present study only provides a measure of the very initial stages of the thrombus formation process, and the current one-way coupled approach, where the effect of thrombus on the flow is not considered, may also be valid only for the early stage of thrombus growth. In this regard, however, the low-dimensional model in *Eq. 16* or *Eq. D5* could eventually be used to extrapolate the growth of the thrombus to a longer time. Our correlation study is focused on the near wall region adjacent to the damaged site, and thus the identified correlation may not apply to a large, volumetric thrombus occurring in the later stage. In the present study, we considered three sizes of apical aneurysms, but the effect of the location of the infarcts has not been investigated. Finally, the LV geometries employed are highly simplified and do not account for the natural variability that exists in healthy and diseased ventricles. In this regard, we are currently conducting the study with patient-specific LV geometries obtained via an ongoing clinical patient study.

## GRANTS

This research is supported by the National Science Foundation (NSF) through Grants IOS-1124804, IIS-1344772, and CBET-1511200. This work used the Extreme Science and Engineering Discovery Environment, which is supported by NSF Grant TG-CTS100002. This work was also made possible by support from the Johns Hopkins Medicine Discovery Fund.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

J.H.S., R.T.G., and R.M. conception and design of research; J.H.S. performed experiments; J.H.S. and T.A. analyzed data; J.H.S., T.A., R.T.G., and R.M. interpreted results of experiments; J.H.S. prepared figures; J.H.S. and R.M. drafted manuscript; J.H.S., R.T.G., and R.M. edited and revised manuscript; J.H.S., T.A., R.T.G., and R.M. approved final version of manuscript.

## APPENDIX A: FLOW RATE SPECIFICATION

Temporal parameters for the flow rate specification shown in Fig. 1*B* are based on the data by Chung et al. (9) for the heart rate of 60 beats/min, and given by E-wave acceleration time (AT) = 0.12, deceleration time (DT) = 0.136, diastasis duration (DIA) = 0.205, A-wave duration (A_{dur}) = 0.139, and systole duration (SYS) = 0.3 (s). We assumed that the peak systole time is at one-third of SYS. The flow rate waveforms for E-wave (*Q*_{E}), A-wave (*Q*_{A}), and systole (*Q*_{S}) were obtained from the harmonic oscillator model (18, 19) and given by *Eqs. A1*, *A2*, and *A3*, respectively:
(A1)
(A2)
(A3)where α_{E} and α_{S} are defined by:
(A4)where SYS_{p} = SYS/3, and ω_{A} is obtained from the solution of *Q*_{A} (A_{dur}) = 0. Amplitude coefficients, *A*_{E}, *A*_{A}, and *A*_{S} are related with the SV and E/A ratio (EAR) by:
(A5)
(A6)
(A7)

## Appendix B: MITRAL VALVE OPENING ANGLE SPECIFICATION

Temporal functions of the opening angle of MV leaflets, ϕ_{i}, where *i* = AL or PL, are specified by:
(B1)where ϕ_{i,0} is the initial angle just after the separation of two leaflets from the coaptation, ϕ_{i,f} is the angle at the fully opened position, ϕ_{i,m} is the angle during diastasis, and *t*_{EP} = 0.12, *t*_{EE} = 0.256, *t*_{AS} = 0.461, *t*_{AP} = 0.539, and *t*_{ED} = 0.65 (s) correspond to the times of peak E-wave, end E-wave, beginning A-wave, peak A-wave, and end diastole, respectively. In the present study, we set ϕ_{AL,0} = 24.78, ϕ_{PL,0} = 70.6, ϕ_{AL,f} = ϕ_{PL,f} = 90, ϕ_{AL,m} = 45, and ϕ_{PL,m} = 70.6 (degrees). These angles are based on the measurements from the particular echocardiogram. The measured angles are slightly modified to be fitted to the present computational model of the LV.

## APPENDIX C: BIOCHEMICAL MODELING OF COAGULATION CASCADE

The TF pathway of the CC shown in Fig. 2 schematically can be described in detail with total 16 biochemical reaction equations for 18 species (4):
(C1)where TF is the tissue factor; V, VIII, IX, and X are the coagulation factors; and Va, VIIa, VIIIa, IXa, and Xa are the active forms of the coagulation factors. II and IIa are the prothrombin and active thrombin, respectively. Note that the complexes TF:VIIa, VIIIa:IXa, Va:Xa, IX:TF:VIIa, X:VIIIa:IXa, and X:TF:VIIa are treated as separate species. Eighteen species are ordered as listed in Table C1, and the reaction coefficients *k* are given in Table C2. The detailed description for the reaction coefficient can be found in Biasetti et al. (4) and Jones and Mann (21). The rate coefficients are determined by Jones and Mann (21) based on the experimental study of Lawson et al. (26).

The reaction source term in the CDR equation (*Eq. 2*) for *i*th species *R*_{i} is given by:
(C2)where *r*_{j} is the reaction rate of the *j*th reaction, and *S*_{ij} is the element of stoichiometric matrix *S*, which relates the species and reactions. If the *i*th species presents on the left-hand side of the *j*th reaction, *S*_{ij} = −1, and if it is on the right side, *S*_{ij} = 1. Otherwise, *S*_{ij} = 0. The reaction rate for *j*th reaction *r*_{j} is calculated as:
(C3)where *k*_{f,j} and *k*_{b,j} are the reaction coefficients for the forward and backward reactions, respectively; *C*_{a} and *C*_{b} are the concentrations of the species on the left-hand side of the reaction equation; and *C*_{c} and *C*_{d} are for the species on the right-hand side. If there is only one species on the left- or right-hand side, the multiplication of two concentrations is replaced by the concentration of the only species. For example, for the first reaction in *Eq. C1*, the reaction rate is given by *r*_{1} = *k*_{6}[*C*_{IX}][*C*_{TF:VIIa}] − *k*_{16}[*C*_{IX:TF:VIIa}].

## APPENDIX D: SENSITIVITY TEST OF COAGULATION CASCADE MODEL

The sensitivity of the current CC model to the initial concentration of the coagulation factors in the bloodstream and the available TF:VIIa concentration on the damaged wall is examined here. Given the large number of parameters to be tested, we employ for this analysis, a low-dimensional model, instead of performing full 3D simulations.

The one-dimensional (1D) model is based on a small control volume very near the damaged wall, as shown in Fig. D1. Based on the CDR equation, *Eq. 2*, the evolution equation for the volume-averaged concentration of species in the control volume can be written as:
(D1)

where *C̄*_{i} is the volume averaged concentration over the control volume; *C*_{i,in}*U*_{in} and *C*_{i,out}*U*_{out} are the concentration fluxes into and out of the control volume, respectively; Δ*s* is the length of control volume; *R*_{i} is the reaction source term; and *F*_{i} is the diffusion term. Since the control volume is small and very near the wall, we apply the following assumptions:
(D2)where *U* is the streamwise velocity magnitude, and *C*_{i,0} is the free stream concentration. The reaction term *R*_{i} is the same as given by *Eq. C2*. Since the diffusion coefficients for the biochemical species are very small, the diffusion term is neglected in the present simplified model, except the net diffusive flux of TF:VIIa from the damaged wall, which is modeled as:
(D3)where *C*_{2,wall} is the concentration of TF:VIIa at the lumen wall, and Δ*n* is the height of the control volume and is set to 0.5 mm, which is the same value used in the 3D simulations.

The blood will stay in the small control volume for a duration of Δ*s*/*U*; thus we can also assume that:
(D4)

With the above assumptions, *Eq. D1* can be reduced to:
(D5)

*Equation D5* is a 1D model of CC near the damaged wall, and the effect of local washout is incorporated via the NWRT in the model. Since the reaction mechanism remains the same with the reaction term, *R*_{i}, we can perform the sensitivity test to the reaction parameters using this 1D model with much reduced computational cost. The overall validity of the current 1D model is also checked in the later of this section.

To test the sensitivity on thrombin production, we vary the initial concentrations of factor IX, X, V, VIII, VIIIa, and II, and the wall concentration of TF:VIIa by taking 1/10, 1/3, 1, 3, and 10 times the baseline values (see Table C1). The 1D model, *Eq. D5* is solved with NWRT = 1 s for six cardiac cycles (*t* = 6 s), and the results are plotted in Fig. D2*A*. The data show that the produced thrombin concentration is not very sensitive to the change of initial concentration of these factors, except TF:VIIa. The result is understandable, since TF:VIIa is the main initiator of CC, and it activates several other coagulation factors. The second-most sensitive parameter is the factor X concentration, but a two orders-of-magnitude variation in this parameter results only in about a threefold change in the thrombin concentration.

Given the sensitivity of thrombin concentration to the wall concentration of TF:VIIa, we have examined the correlation between the NWRT and the thrombin concentration with our 1D model for three different values of TF:VIIa wall concentration: 10, 100, and 1,000 μmol/m^{3}. The 1D model is solved with NWRT = 0.2–2.5 s for these cases, and the results are shown in Fig. D2*B*. First we note that the 1D model result with the baseline value (100 μmol/m^{3}) compares well with the results from the 3D simulations presented in Fig. 8, and this confirms the general validity of the 1D model. Second, the 1D simulations for the other values of TF:VIIa show an overall qualitative behavior that is very similar to the baseline case. This strongly suggests that, while quantitative values of thrombogenic products might depend on wall concentration of TF:VIIa, the overall trends are quite robust to changes in these values.

- Copyright © 2016 the American Physiological Society