Computational modeling of blood component transport related to coronary artery thrombosis in Kawasaki disease.

Coronary artery thrombosis is the major risk associated with Kawasaki disease (KD). Long-term management of KD patients with persistent aneurysms requires a thrombotic risk assessment and clinical decisions regarding the administration of anticoagulation therapy. Computational fluid dynamics has demonstrated that abnormal KD coronary artery hemodynamics can be associated with thrombosis. However, the underlying mechanisms of clot formation are not yet fully understood. Here we present a new model incorporating data from patient-specific simulated velocity fields to track platelet activation and accumulation. We use a system of Reaction-Advection-Diffusion equations solved with a stabilized finite element method to describe the evolution of non-activated platelets and activated platelet concentrations [AP], local concentrations of adenosine diphosphate (ADP) and poly-phosphate (PolyP). The activation of platelets is modeled as a function of shear-rate exposure and local concentration of agonists. We compared the distribution of activated platelets in a healthy coronary case and six cases with coronary artery aneurysms caused by KD, including three with confirmed thrombosis. Results show spatial correlation between regions of higher concentration of activated platelets and the reported location of the clot, suggesting predictive capabilities of this model towards identifying regions at high risk for thrombosis. Also, the concentration levels of ADP and PolyP in cases with confirmed thrombosis are higher than the reported critical values associated with platelet aggregation (ADP) and activation of the intrinsic coagulation pathway (PolyP). These findings suggest the potential initiation of a coagulation pathway even in the absence of an extrinsic factor. Finally, computational simulations show that in regions of flow stagnation, biochemical activation, as a result of local agonist concentration, is dominant. Identifying the leading factors to a pro-coagulant environment in each case-mechanical or biochemical-could help define improved strategies for thrombosis prevention tailored for each patient.

Introduction

Coronary artery thrombosis is the major risk associated with Kawasaki disease (KD) [1]. Long-term management of KD patients with persistent aneurysms requires thrombotic risk assessment and clinical decisions regarding the administration of anticoagulation therapy [1]. Typically, medication prescribed to these patients falls into two main categories of anticoagulants: vitamin K inhibitors such as Warfarin, acting on coagulation factors, and antiplatelet drugs such as Clopidogrel, acting on the P2Y receptor on the platelet surface related to Adenosine Diphosphate (ADP) mediated platelet activation. Also, oral Xa direct inhibitors or direct thrombin inhibitors may be used as an alternative anticoagulation therapy in KD [1]. These are currently not approved for pediatric use, although clinical trials are on-going [2]. Patient-specific computational modeling has provided new insight regarding thrombotic risk for KD patients [3-6], demonstrating that hemodynamic variables obtained from computational fluid dynamics (CFD) simulations have the potential to identify regions prone to thrombosis [3]. However, the underlying mechanisms of clot formation, which could help determine favorable strategies for thrombosis prevention, are not fully understood.

The coagulation cascade has been widely studied. In particular, in the setting of hemostasis. The two coagulation pathways, the intrinsic or contact pathway and the extrinsic pathway initiated by tissue factor (TF), converge to the common pathway leading to activation of factor X and generation of thrombin, fibrin, and ultimately a stable clot (Fig 1). Platelets play an essential role in thrombus formation, especially at early stages, and can be activated by thrombin and through interactions with collagen and von Willebrand factor (vWF). Upon activation, platelets undergo changes in their membrane that facilitate adhesion to the vessel wall and aggregation to other platelets, and release the content stored in their dense granules. Chemical agonists released from platelet dense granules include adenosine diphosphate (ADP), thromboxane A2 (TXA2), which are involved in platelet activation; and polyphosphate (PolyP) [7] reported in recent studies as a potential activator for the intrinsic coagulation pathway [8-11].

Fig 1

Simplified representation of the coagulation cascade.

Intrinsic and extrinsic pathways and proposed pathway for aneurysm induced thrombosis (increased platelet activation + PolyP). TF = Tissue Factor, AP = Activated platelet, RP = Resting platelet, PolyP = Polyphosphate.

Computational studies aiming to model arterial thrombosis often rely on an arterial wall injury that exposes collagen, vWF, and other sub-endothelial components to the bloodstream, which facilitates platelet adhesion to the wall even in a non-activated state [12], and subsequent clot formation. Also, extravascular TF can come in contact with the bloodstream providing the trigger required to initiate the extrinsic coagulation cascade and generate thrombin. Several studies have modeled the extrinsic coagulation pathway and most of the chemical reactions involved using computational methods. Models by Fogelson and Kuharsky et al. [13-15] provide a comprehensive description of the coagulation cascade, including the interactions of coagulation factors, platelet deposition and aggregation; and flow-mediated transport to explore the role of hemodynamics on early coagulation events. Work by Leiderman et al. [16, 17] investigated the implications of blood components transport in thrombus growth. Later studies further demonstrated the essential role that transport physics plays in thrombosis [18, 19], in particular during clot formation and growth. Work by Chatterjee et al. [20] proposed a systems biology approach to model coagulation initiation that also considered the intrinsic coagulation pathway and the interplay between the extrinsic and intrinsic coagulation pathways. However, these models are not adequate to explain thrombosis in other settings where even in the absence of wall injury, blood clots may still form under complex flow conditions. Experimental data suggest the role of hemodynamics in thrombus initiation and growth is two-fold: shear forces can contribute to platelet activation, and flow stagnation facilitates platelet aggregation, driving thrombus formation and growth [21]. Also, under elevated shear conditions shear-induced platelet aggregation can lead to occlusive arterial thrombosis [12, 22]. This process is mediated by vWF and is independent of platelet activation. In addition, platelets can be activated by either mechanical (shear-induced platelet activation) or biochemical stimuli, which suggest different potential targets for thromboprophylaxis [15]. Data from both experimental [23, 24] and computational studies [25, 26] show that platelet activation depends not only on the level of shear rate exposure but also on the length of the exposure, suggesting that longer exposure to moderate shear levels could also induce platelets change to an activated state. Platelet activation potential and shear rate exposure have been investigated related to stenosis in coronary or carotid arteries, as a means to quantify the effect of shear rate on platelet activation [23, 25, 27, 28], and regarding high shear flows induced by medical devices [24, 29, 30]. Prior CFD studies have considered the effects of shear-induced platelet activation by using activation functions that rely on critical thresholds for shear stress, either low [28] or high [31]. Work by Di Achille et al. [32] used Lagrangian particle tracking to obtain information on the shear history of particles traveling through the domain. In this case, shear history was used to compute a platelet activation potential index, which was proposed as a predictor for areas prone to intraluminal thrombosis in abdominal aortic aneurysms. Other models, focused on the TF or extrinsic coagulation pathway considered platelet activation as a function of chemical agonists and thrombin concentrations and proposed a thrombus growth model, based on platelet-platelet and platelet-wall adhesion forces, at low and high shear rate levels [33].

Typically thrombosis models have been coupled with the flow, widely using a 2D approach, relying on idealized flow waveforms [16, 34–38] and considering small computational domains on the order of micrometers, localized in the proximity of a TF source. A few recent studies have presented 3D models to investigate low shear and high shear thrombosis yet the majority have considered idealized rather than patient-specific anatomies [22, 33, 39, 40]. Yadzani et al. [41] proposed a hemodynamic transport model to predict the location and extent of thrombus in aortic dissection in mice. Regarding the use of patient-specific anatomies in the context of thrombosis modeling, work by Di Achille et al. [32, 42] and Arzani et al. [43] proposed hemodynamic based metrics to identify regions prone to thrombosis in abdominal aortic aneurysms. However, these studies considered only a mechanical perspective of thrombus formation and did not consider the biochemical aspects related to thrombus initiation. Menichini et al. [28] presented a computational model to predict false lumen thrombosis in patient-specific aortic dissection anatomies, based on fluid shear rate, residence time, and platelet distribution. More recently, Li et al. [44] presented a predictive model for intraluminal thrombus formation in diabetic retinal microaneurysms, which combines a continuum and a particle model to simulate platelet aggregation on thrombogenic surfaces inside these microaneurysms. Our previous work in KD patients [3] also provided wall shear stress exposure and residence time as metrics for thrombotic risk stratification.

One limitation of current thrombosis models is their reliance on assumptions about the thrombotic potential of the wall, for which there is very limited quantitative data. These models require that initial extrinsic factor, TF, or assumed collagen or vWF concentration is prescribed as the initiator for the coagulation cascade. Although these assumptions hold in the case of trauma or athero-thrombosis, there are other scenarios where no data on the thrombogenic potential of the arterial wall is available. Therefore, these models fail to accurately represent the process of clot formation in such cases. Coronary artery aneurysm thrombosis following KD is an example of the need to understand the mechanisms of thrombus initiation in the absence of an extrinsic factor. No obvious wall injury, which would enable a local release of TF, has been reported concerning thrombosis in KD aneurysms. Therefore, we hypothesized that the intrinsic or contact coagulation pathway may be responsible for thrombus formation in this setting. Our previous studies related to KD aneurysms, have related thrombosis risk with abnormal hemodynamics, suggesting that hemodynamics are likely the primary driver. We present a new model that incorporates information from patient-specific simulated velocity fields to track platelet activation and accumulation in patient-specific coronary artery aneurysms caused by KD. Increased concentration of activated platelets and agonists induces a pro-coagulant environment that may drive thrombus initiation and growth. This hypothesis is based on recent studies that have suggested that critical concentration of PolyP, released from platelets upon activation, can activate factor XII and initiate the intrinsic coagulation cascade [7-10] and data that relates critical concentrations of ADP with platelet activation and aggregation [45]. The main goal of this study is to investigate the potential for platelet activation from patient-specific hemodynamics and establish links between platelet activation—mechanical and biochemical—and initiation of the intrinsic coagulation pathway.

Methods

Ethics statement

This study was approved by the Institutional Review Board at Stanford University and The Hospital for Sick Children, University of Toronto. Written subject consent or assent and parental written consent were obtained as appropriate for the imaging and simulation studies.

Patient population

We considered six KD patients with aneurysms and one control patient with normal coronary arteries. The diseased group was drawn from our previous study on patient-specific hemodynamics and risk of thrombosis in KD [3]. This group included three patients classified as low risk for thrombosis according to all the hemodynamic metrics analyzed in [3], and three patients with confirmed thrombosis, featuring different aneurysms shapes, sizes, and locations.

Shear induced platelet activation: Cumulative shear rate

To quantify cumulative shear rate exposure in the 3D domain we used a Lagrangian approach, computing shear rate exposure along path lines calculated from a patient-specific velocity field. Shear rate is defined as: where E is the rate of strain tensor

Cumulative shear rate was defined as a linear function of shear rate exposure time

To initiate the calculation, massless particles representing platelets (Stokes number ≪ 1) were released at the most proximal section of the coronary artery model (ostium) and their trajectories were time-integrated using a fourth order Runge-Kutta scheme based on the velocity field obtained from a patient-specific blood flow simulation [3]. We obtained hemodynamic quantities, including instantaneous or cumulative shear rate, along these path lines. Considering a sufficient number of seeding points at the inlet, we were able to retrieve a full 3D reconstruction of the average cumulative shear rate. This 3D field represents the average cumulative shear rate that a platelet, at a specific location in the 3D model, would be exposed to, as it travels along the coronary artery (Fig 2). The number of seeds was determined to match the number of nodes from the finite element mesh at the inlet for each model, therefore the number of seeds was adjusted to the size of the domain for each case. The simulation was initiated with the release of the particles at the seed points (representing platelets entering the domain) and ended when one of the following was achieved: path line reached the distal outlet (platelet leaving the coronary model), stagnation (constant point coordinates of the path line for more than one cardiac cycle), the maximum number of cycles prescribed for the simulation was reached. For this study, we set the maximum number of cardiac cycles for the simulation to 30 cycles, based on previous work [3]. Path lines were calculated in parallel by distributing the number of seed points to multiple processors so that the integration of each path line was performed independently. The cumulative shear rate computed for each trajectory was combined to obtain an average 3D cumulative shear rate field in a post-processing step. Examples of the cumulative shear rate distributions for a normal right coronary artery and a right coronary artery with an aneurysm are presented in Fig 2. The cumulative shear rate in the normal coronary artery increases linearly as a function of distance from the ostium, as expected from the definition (Eq 3). On the other hand, in the aneurysm case, we observe a non-linear distribution, with regions of high cumulative shear rate in proximal and medial regions of the coronary artery. Visualization of the path lines (Fig 3) shows the contrast between high shear rates and high cumulative shear rates in different regions of the model.

Fig 2

Cumulative shear rate distribution in Kawasaki disease right coronary artery.

Normal coronary vs. coronary artery aneurysm. CumSR = Cumulative Shear rate.

Fig 3

Path lines of particles released at the beginning of the cardiac cycle in a left coronary branch (P4).

a) Shear Rate (SR); b) Cumulative shear rate (CumSR) along path lines; c) Cumulative shear rate 3D distribution (CumSR), computed as a time average for the cardiac cycle.

To describe platelet activation related to high cumulative shear rate, we used an activation function, based on the cumulative shear rate levels observed in the healthy coronary model, and reported literature values for platelet activation indices found in regions prone to thrombosis [23, 25]. The activation function was defined as: where and is the rate for platelet activation due to cumulative shear rate.

We chose four time points equally distributed along the cardiac cycle for particles release to compute path lines (0, 0.25T, 0.5T, 0.75T, where T is the period of the cardiac cycle). Although we observed differences in cumulative shear rate distribution when comparing results from particles released at the beginning of the cardiac cycle to other time points, we did not observe substantial differences when increasing the number of time points up to ten.

Reaction-Advection-Diffusion equations

We use a continuum approach to model the transport of the main blood components and chemical agonists involved in thrombus initiation based on the models presented by Ananad [31] and Leiderman and Fogelson [16]. We solve the Reaction-Advection-Diffusion (RAD) equations using a stabilized Finite Element Method suitable for advection dominated problems [46], to compute the temporal evolution of the concentrations of non-activated platelets [RP] and activated platelets [AP], and local concentrations of Adenosine Diphosphate [ADP] and PolyPhosphate [PolyP]. To perform this simulation, sub-models of the right and left coronary branches were considered independently. Activation of platelets is modeled considering two different contributions to the reaction term: mechanical activation, a function of cumulative shear rate, and biochemical activation, a function of the local concentration of activated platelets and agonists such as ADP. We define the advection-diffusion operator as: where [A] represents the concentration of the i chemical specie considered and u is the patient-specific pre-computed fluid velocity. The velocity fields were computed using Simvascular [47] as described in our previous work [3]. Also, the velocity fields obtained from KD simulations have been previously validated against flow fields measured experimentally in vitro using phase-contrast MRI [48]. We express the RAD equations as where represents the different contributions to production or depletion of A.

The definition of the RAD equation, for each of the chemical species considered in our model is presented below:

Resting Platelets: where A(SR) is the activation function defined in Eq 4.

Active Platelets:

ADP:

Polyphosphate (PolyP): where [PolyP] and [ADP] account for ADP and PolyP released from dense granules upon platelet activation.

Boundary conditions and initial conditions

A Dirichlet boundary condition was imposed at the inlet of the domain (proximal segment of the coronary artery) to match the circulating concentration of platelets and other chemical species considered in our model. A no flux Neumann boundary condition was prescribed at the walls. The 3D mesh was initialized to the average circulating values of resting platelets, activated platelets, ADP, and PolyP according to literature data (Table 1). No abnormal platelet count or hematocrit is usually reported in patients with KD, therefore we assumed average population values for platelet count (2.5 1011 platelets/l) with 5% of activated platelets circulating in the bloodstream [49]. Chemical reaction rates were set to literature values, [14], s−1[31]. Platelet activation via thrombin is neglected () since we are considering thrombus initiation before thrombin production. The diffusivity constant for platelets considering the average shear rate levels observed in our hemodynamic simulations was on the order of 10−7. This value was estimated by combining Brownian diffusion and enhanced-shear diffusion using the definition given by Anand et al. [31]. Enhanced-shear diffusion was also considered for ADP and PolyP. Diffusion coefficients and initial conditions for the simulation are summarized in Table 1.

Table 1

Initial conditions and diffusion coefficient values for Reaction-Advection-Diffusion (RAD) equations.

Species	Initial concentration [nM]	Diffusion Coefficient [cm2/s]	Reference
RP	10	10−7	[14], [31]
AP	0.5	10−7	[49], [31]
ADP	250	25.7⋅10−7	[50]
PolyP	0	10−7	[8, 51, 52, 31]

RP = resting platelets, AP = activated platelets, ADP = adenosine diphosphate, PolyP = polyphosphate

Biochemical platelet activation

One of the primary mediators of biochemical platelet activation is the release of chemical agonists from platelet dense granules upon activation. We assumed the following values for the contents of platelet dense granules based on the literature [7]: [ADP] = 3 mol/1017 plt and [PolyP] = 0.743 mol/1017 plt. Assuming at least 75% percent of the platelets dense granules content is released, we can estimate the corresponding concentrations of released [ADP] and [PolyP]:

Although reaction constant rates are given in terms of [AP] and [RP] we included an activation function that limits the activity of this reaction based on a minimum ADP concentration, [ADP] > 0.2 μM, which is reported in the literature as the threshold value for reversible changes in platelets [45]. The activation function is defined as a ramp function between [ADP] = 0.2 μM (reversible activation) and [ADP] = 1 μM (irreversible activation).

Post-processing and quantities of interest

We performed a qualitative comparison of the distribution of [AP] in a normal coronary tree, and low and high thrombosis risk cases, and quantified the concentration of resting and activated platelets within the aneurysms. In addition, we quantified the concentration of agonists such as ADP and PolyP to determine if their concentration reached the critical values associated with irreversible changes in platelets and platelet aggregation in the case of ADP, and activation of factor XII (intrinsic coagulation pathway) in the case of PolyP. Critical values were set to [ADP] = 1μM [PolyP] = 2μgr / ml according to experimental values reported in the literature [8, 45]. To quantify the potential for thrombus initiation and the pro-coagulant environment induced by the aneurysm, we considered the following quantities of interest, based on concentration fields obtained from the RAD simulations:

Volume fraction of the domain where the concentration of activated platelets is greater than or equal to a specified level of platelet activation. We define AP as the volume fraction of the aneurysm where platelet activation is greater or equal to x% activation, with full platelet activation at [AP] ≥ 10nM.

Volume fraction of the aneurysm where [ADP] and [PolyP] are above the specified critical thresholds.

For patients with more than one consecutive aneurysms with no normal coronary segment in between, consecutive aneurysms were considered in aggregate.

Results

High concentration of activated platelets correlates with regions prone to thrombosis

Qualitative comparison of [AP] spatial distribution between coronary arteries with aneurysms and the control case showed that high [AP] correlates with regions of confirmed thrombosis (Figs 4 and 5). In the case of patient P2, we were able to perform a direct validation of our predictions since the thrombus was still visible on the image data (Cardiac MRI). Volumetric segmentation of the clot overlapped with the area of the highest [AP] (Fig 4). The concentration of activated platelets in most of the aneurysms classified as high risk of thrombosis reached values substantially higher than in the normal coronary case, where [AP] remained close to the initial value for the simulation, which assumed a 5% platelet activation. Coronary arteries with aneurysms classified as low risk for thrombosis, presented with [AP] generally lower than the cases classified as high risk of thrombosis and confirmed clot formation. In patients with confirmed thrombosis, levels of [AP] were equal or greater than the initial concentration of non-activated platelets (10 nM/l), at the specific locations where thrombosis was reported, which implies full platelet activation. The volume fraction of the aneurysm where [AP] > 10 nM, representing full platelet activation, ranged from 10–100% for the high thrombotic risk patients. On the low-risk cases, this volume fraction was smaller, ranging from 0–10% (Fig 6). Volume fraction where [AP] >10M for the normal coronary case was 0%.

Fig 4

Qualitative comparison of [AP] spatial distribution for Kawasaki disease patients classified as high thrombotic risk.

Arrow indicates region of confirmed thrombosis. Yellow shaded area represents the thrombus segmentation from image data. [AP] = Activated platelet concentration.

Fig 5

Qualitative comparison of [AP] spatial distribution for Kawasaki disease patients classified as low thrombotic risk and a normal coronary artery control.

[AP] = Activated platelet concentration.

Fig 6

Pro-coagulant environment quantification based on concentration of activated platelets and agonists.

a) Changes in fraction of aneurysm volume where [AP] >5, 8, 10 nM; b) fraction of aneurysm volume where [AP] >10 nM. Pro-coagulant environment quantification based on ADP and PolyP critical thresholds. c) Fraction of aneurysm volume where [ADP] > [ADP]; d) fraction of aneurysm volume where [PolyP] > [PolyP]. Th = confirmed thrombosis, NonTh = non confirmed thrombosis, [AP] = activated platelet concentration.

ADP and PolyP concentration reaches critical thresholds in areas of confirmed thrombosis

Quantification of aneurysm volume fraction where ADP and PolyP were greater or equal to [ADP] and [PolyP] respectively suggested a more pro-coagulant environment in those aneurysms with confirmed thrombosis (Fig 6). Except for patient P7, pro-coagulant volume fractions measured were almost zero for non-thrombosis branches. Note that for patient P3 the aneurysm in series with the aneurysms where the clot was reported also exhibits larger volume fraction of pro-coagulant environment measured by [ADP] and [PolyP].

Biochemical platelet activation is dominant in regions of flow stagnation

We compared the contributions of mechanical and biochemical mechanisms to platelet activation. We observed that in regions of flow stagnation, the dominant reaction term represented the contribution due to local [AP]. The analysis of the temporal evolution of these contributions showed that, even if at the beginning of the simulation, mechanical activation, based on cumulative shear rate, was the main driver for platelet activation, as [AP] increases, the biochemical activation term becomes dominant (Fig 7). The contribution of the cumulative shear rate in the cases we analyzed, even if dominant during the first or second cardiac cycles, was small. This suggests that the contribution of shear rate to platelet activation in KD aneurysms is almost negligible, and platelet activation is instead mediated by the local increase of chemical agonist concentrations in the aneurysm as a result of flow stagnation (Fig 8).

Fig 7

Temporal evolution of the contribution of mechanical and biochemical factors to platelet activation.

Contribution to platelet activation as defined in Eq 9. Volume fraction where SR is dominant is negligible in regions of flow stagnation. Data corresponding to P2 right coronary artery. BioChem = Platelet activation due to concentration of agonists, SR = platelet activation based on cumulative shear rate activation function (Eq 4).

Fig 8

Qualitative comparison of platelet activation due to mechanical and biochemical contributions as a function of time.

Thrombus 3D segmentation in yellow.

Pro-coagulant environment quantification shows correlation with hemodynamic variables

We observed a qualitative correlation between areas of low shear stress and high residence time with regions of increased pro-coagulant activity, measured by [AP], [ADP] and [PolyP]. Quantitative correlation of AP10 with hemodynamic variables showed high correlation with residence time (r2 = 0.97). On the other hand, we observed low correlation of AP10 (volume fraction) with low wall shear stress exposure, A (surface area fraction, TAWSS < 1 dyne/cm2, r2 = 0.23).

Discussion

This paper is the first computational study to consider the biochemical aspects of thrombosis for a patient-specific quantification of aneurysm pro-coagulant environment in KD. The model is also novel in the sense that it investigates the conditions for thrombus initiation without relying on an extrinsic trigger to the coagulation cascade. This model can be extended to investigate other cardiovascular problems, such as aortic dissection, that present with stagnation-driven thrombus as a result of complex flow patterns and without a focal injury site. Also, the timeline for clot formation may be distinct from other coagulation and thrombosis scenarios where the clot forms in seconds or hours, usually as a result of an extrinsic component such as TF being released into the bloodstream, or the exposure of sub-endothelial collagen and vWF. The model presented here relies solely on the anatomy of the patient, at a given point in time, and the physics and biochemistry that govern blood flow and platelet activation to predict where and if the clot is likely to form.

Previous studies have demonstrated that hemodynamics, in particular, wall shear stress and residence time correlate with regions prone to thrombosis [3, 5]. Here we provided added value to residence time calculation since we not only identify areas of potential accumulation of platelets and agonists, but also we quantify the local concentration of chemical agonists and compare with critical thresholds related to platelet aggregation in the case of ADP, and factor XII activation in the case of PolyP.

Our results show qualitatively high correlation between regions of low WSS and high residence time and regions with pro-coagulant activity. However, surface metrics (area exposed to low WSS, A) are not well correlated with volume metrics (volume fraction where AP concentration is high, AP10). The fraction of surface area exposed to low WSS (A) is calculated with a threshold of TAWSS < 1 dynes/cm2. This suggests that wall quantities such as low WSS exposure (A) may not be enough to quantify the abnormal hemodynamics and pro-coagulant environment in this case. Even if the fraction of the wall exposed to pathological WSS is smaller, there could still be areas of high recirculation and flow stagnation within the aneurysms. The focus of many studies on wall quantities in relation to thrombosis is typically motivated by the assumed thrombogenic potential of the aneurysm wall and the hypothesis that platelets concentrate near the wall. The latter is referred to as platelet margination, affected by the Fåhraeus-Lindqvist effect [53]. However, limitations arise when applying this assumption in large arteries in the cardiovascular system, since the scales at which this effect has been usually observed are on the order of μm which is several orders of magnitude below the typical size of coronary aneurysms. Although in some cases platelet margination has been observed at the mm scale, the development of near the wall excess of platelets required ideal flow conditions and lengths well beyond the typical coronary artery length [54].

Importantly, the margination effect has been mostly observed under ideal or healthy flow conditions, in which the flow is primarily unidirectional. Therefore it may not be straightforward to extend this behavior to situations of complex flow. Some studies have suggested shear stress levels may impact platelet margination. In particular, low wall shear stress may limit the margination effect by slowing down margination and increasing the margination length [55-59]. Therefore, based on these data, we did not expect significant margination effect in this study due to the low shear rates observed in KD aneurysms (100 s−1) from our hemodynamic simulations and the typical length of coronary arteries. In addition, geometry changes such as expansions or constrictions have been reported to significantly effect particle margination [60]. For these reasons, here we considered the full 3D fluid domain in our calculations of platelet activation and platelet accumulation. We observed areas of high platelet concentrations growing into the lumen of the artery in those aneurysms exhibiting exceptionally complex flow patterns. This could explain the low correlation between A and the pro-coagulant volume fraction AP10 and demonstrates the added value of combining wall-surface and volumetric metrics. Since, the aneurysmal wall is abnormal in KD patients, especially after remodeling has occurred [61, 62], it is possible that the probability of platelet adhesion to the wall is higher than in normal coronary arteries. Histology and Optical Coherence Tomography studies have reported intimal hyperplasia and irreversible damage to the aneurysm wall [61, 63], however, the relationship between this damage and subsequent thrombus formation remains incompletely understood. Experimental animal models that replicate aneurysm formation in vivo could be useful to guide more realistic modeling of KD thrombus progression in future work.

Since we considered a linear approximation to compute the cumulative shear rate, we expected that it would increase linearly with the distance along the vessel. This holds for the normal coronary case, however, the diseased anatomy shows non-linear concentrations of higher cumulative shear rate. This is a manifestation of the complex flow patterns created in the pathological anatomy. Of note, in the diseased coronary arteries (Figs 3 and 4) we observe areas of high cumulative shear rate in both proximal and distal regions of the coronary artery. This demonstrates that the cumulative shear rate encodes information of high shear rate magnitudes as well as high exposure times, to moderate or even normal shear rates.

Our simulation results show that the contribution of shear rate to platelet activation is almost negligible for the cases considered in this study. However, calculated cumulative shear rate distributions suggest that stenosis proximal to the aneurysm may substantially increase the contribution of shear rate to platelet activation. The biochemical activation term is dominant in areas of high flow stagnation, where the accumulation of circulating pro-coagulants results in concentrations that may be sufficient to induce platelet aggregation or initiate the intrinsic coagulation cascade via factor XII activation [8]. In the absence of flow stagnation, the accumulation of activated platelets, and the increase of [ADP] and [PolyP] in the coronary models analyzed was negligible. We observed that although higher cumulative shear rate correlated with thrombus location to some extent (Fig 2), such values of cumulative shear rate were not enough to result in substantial platelet activation according to the shear-induced platelet activation function (Eq 5). To assess sensitivity of our results to changes in the activation function definition (Eq 5), we also considered a lower threshold for cumulative shear rate activation (SR > 250) and our results did not differ substantially. Completely removing the mechanical contribution to platelet activation also did not substantially change the result. This was expected since the dominant contribution in the cases analyzed here was due to agonist-induced platelet activation.

Shear stress levels in localized areas of high activated platelet accumulation were typically < 1 dyne/cm2. According to experimental assessments of thrombus deformation and growth at different shear rates [37], this mechanical environment would facilitate clot growth and would induce small deformations. Future studies that include platelet dynamics similar to Lu et al. [38] would provide more insight into the potential risk of growth to total occlusion.

In our computational experiments, we did not find substantial differences between [AP] distribution from patient-specific velocity fields obtained from rigid wall or fluid-solid interaction simulations. However, we observed that as the flow becomes more complex, especially in the presence of large flow stagnation areas, wall deformation may have a more significant effect on near-wall transport. This should be considered in future studies, in particular, if large wall deformations are expected, such as in aortic dissection.

The low velocities and low shear rates typically observed in KD aneurysms suggest diffusion transport may have an important role. However, we did not find significant differences in the total volume of full platelet activation (AP10) comparing advection only vs. advection and diffusion solutions from our simulations. This was also confirmed by calculating the Peclet number. Considering coronary artery diameter as the characteristic length scale, Peclet number was >>1 for all cases, demonstrating that our models, at a macroscopic scale, are advection dominated. To investigate the ratio of advection to diffusion, we re-computed the Peclet number using the size of platelets as the characteristic length, approximating flow moving through an aggregate of platelets where the size of the opening is approximately the platelet size. Using this characteristic dimension, representative of the flow at the microscale, there are small regions near the aneurysm wall where diffusion may dominate. To capture diffusion of chemical agonists appropriately, an adaptive meshing method, which refines the mesh as platelets begin aggregating, could be integrated into the solver. In this case, the already low velocities observed in some regions of the aneurysms may be reduced substantially, increasing the diffusion contribution to the transport of agonists at this location.

We performed a pro-coagulant environment quantification as a post-processing step of patient-specific flow simulations. It is important to note that high concentration of activated platelets contributes to thrombus initiation, however this criterion alone is likely not a sufficient condition for thrombus development. For this reason, we also quantify other agonists related to thrombosis (ADP, PolyP). Future extensions of this model will consider the full multi-scale process of clot formation to assess how concentration of activated platelets and agonists directly affect this process. Because of the short time scale considered, it is reasonable to assume that the effect of the fluid-thrombus interaction at this point is very small. Therefore this data provides an initial condition for a future multi-scale thrombus growth simulation identifying regions likely to develop thrombosis due to induced pro-coagulant environments. Our analysis showed that a simulation time of 10 cardiac cycles was sufficient to capture activated platelet and agonist concentration distributions from a given patient-specific velocity field. Future studies should consider a coupled approach that solves the fluid and the RAD equations simultaneously and that accounts for the changes in the flow field as the clot grows. Multi-phase models such as those proposed by Xu et al. [35] and Zheng et al. [40] could potentially be adapted for patient-specific simulations. Experimental data describing clot composition and material properties would be key to increase the realism of these models and would allow validation of the simulated clot formation process under flow.

Our model simulations did not include the full intrinsic coagulation cascade, but the estimation of the critical conditions that may trigger the intrinsic coagulation pathway ([PolyP] > [PolyP]). Empirical PolyP biochemical data including activation rate constant for factor XII activation needs to be incorporated to achieve a patient-specific intrinsic coagulation cascade model from patient-specific anatomy to clot formation. In addition, future studies aiming to simulate thrombus growth and concentrated on the microscale region near the wall would need to consider a more precise diffusivity value for PolyP, which is specific for PolyP chains released upon platelet activation.

Conclusion

This computational modeling study illustrates possible mechanisms for thrombus initiation under low shear conditions, such as those induced by aneurysms. In the absence of extremely high shear rates, prolonged exposure to moderate shear rates has the potential to activate platelets. In the normal coronary case we observe that even if platelets are activated due to high cumulative shear rates, the activated platelets and chemical agonists released due to platelet activation are typically washed out of the domain after one cardiac cycle. That is, platelets do not accumulate anywhere in the coronary model. On the contrary, in the case of aneurysms, we observe a significant accumulation of both activated platelets and chemical agonists due to flow stagnation. Furthermore, the local concentration of agonists in some aneurysms reaches levels comparable to those reported as critical thresholds to induce platelet activation, amplifying the effect of shear rate activation and creating an environment prone to thrombus initiation, potentially though the intrinsic coagulation pathway. Current guidelines rely only on the diameter of the aneurysm as the criterion for initiating systemic anticoagulation. Our results suggest that it is precisely for those patients, such as patient P2, who are near the cutoff Z-score value that alternative metrics, based on hemodynamic and biochemical simulations, such as pro-coagulant environment volume quantification, can add value to thrombotic risk prediction. In addition, we showed that coronary artery aneurysm thrombosis may depend to a great degree on platelet activation, suggesting some type of antiplatelet therapy as the most efficient antithrombotic therapy.

This work provides a new framework to investigate thrombus initiation from a patient-specific perspective. The model presented here could help identify regions at higher risk of thrombosis based on patient-specific hemodynamics, even in the absence of a focal lesion, as well as the best strategies for thrombosis prevention by analyzing the most significant contributions to platelet activation in each case. Results presented here would motivate prospective studies of larger patient cohorts. Future work should include information on patient-specific coagulation profiles to increase patient-specificity and clinical relevance. This approach could also enable the addition of anticoagulation drugs into the model to evaluate their efficiency and dose management, and design treatment strategies tailored for each patient.