https://orcid.org/0009-0004-7664-8484
Figures
Abstract
Clot formation is a crucial process that prevents bleeding, but can lead to severe disorders when imbalanced. This process is regulated by the coagulation cascade, a biochemical network that controls the enzyme thrombin, which converts soluble fibrinogen into the fibrin fibers that constitute clots. Coagulation cascade models are typically complex and involve dozens of partial differential equations (PDEs) representing various chemical species’ transport, reaction kinetics, and diffusion. Solving these PDE systems computationally is challenging, due to their large size and multi-scale nature. We propose a multi-fidelity strategy to increase the efficiency of coagulation cascade simulations. Leveraging the slower dynamics of molecular diffusion, we transform the governing PDEs into ordinary differential equations (ODEs) representing the evolution of species concentrations versus blood residence time. We then Taylor-expand the ODE solution around the zero-diffusivity limit to obtain spatiotemporal maps of species concentrations in terms of the statistical moments of residence time, , and provide the governing PDEs for
. This strategy replaces a high-fidelity system of N PDEs representing the coagulation cascade of N chemical species by N ODEs and p PDEs governing the residence time statistical moments. The multi-fidelity order (p) allows balancing accuracy and computational cost providing a speedup of over N/p compared to high-fidelity models. Moreover, this cost becomes independent of the number of chemical species in the large computational meshes typical of the arterial and cardiac chamber simulations. Using a coagulation network with N = 9 and an idealized aneurysm geometry with a pulsatile flow as a benchmark, we demonstrate favorable accuracy for low-order models of p = 1 and p = 2. The thrombin concentration in these models departs from the high-fidelity solution by under 20% (p = 1) and 2% (p = 2) after 20 cardiac cycles. These multi-fidelity models could enable new coagulation analyses in complex flow scenarios and extensive reaction networks. Furthermore, it could be generalized to advance our understanding of other reacting systems affected by flow.
Author summary
The coagulation cascade is an intricate biochemical process that prevents excessive bleeding while maintaining vascular integrity. Modeling this process involves dozens of interdependent chemical reactions with disparate kinetics. Moreover, the reacting species are transported by the flow, leading to complex spatio-temporal dynamics. Consequently, the computational cost of modeling the coagulation cascade in flowing blood prohibits realistic simulations. To overcome this challenge, we introduce a new multi-fidelity approach that exploits the slow diffusion of chemical species to decouple simulating their flow transport and chemical reaction. This approach achieves a significant reduction in computational requirements while maintaining the accuracy of our simulations. We anticipate this new multi-fidelity approach will make coagulation cascade simulations in physiologically relevant scenarios accessible to many researchers. This technique will also enable studies necessitating multiple parametric runs, such as sensitivity or uncertainty quantification analyses. This advancement is poised to benefit medical professionals and researchers, opening new horizons in our understanding of coagulation processes and the effects of blood thinners.
Citation: Guerrero-Hurtado M, Garcia-Villalba M, Gonzalo A, Martinez-Legazpi P, Kahn AM, McVeigh E, et al. (2023) Efficient multi-fidelity computation of blood coagulation under flow. PLoS Comput Biol 19(10): e1011583. https://doi.org/10.1371/journal.pcbi.1011583
Editor: Alison Marsden, Stanford University, UNITED STATES
Received: June 1, 2023; Accepted: October 9, 2023; Published: October 27, 2023
Copyright: © 2023 Guerrero-Hurtado et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper, its Supporting information files, and on Zenodo (DDOI:10.5281/zenodo.8344615).
Funding: MGH, MGV and OF have been partially supported by the Spanish Research Agency and the European Regional Development Fund, under grant number PID2019-107279RB-I00. MGH, MGV, PML, JB and OF have been partially supported by the Comunidad de Madrid and the European Regional Development Fund, under grant number Y2018/BIO-4858 PREFI-CM, and by the Instituto de Salud Carlos III and the European Regional Development Fund, under grant numbers PI15/02211-ISBITAMI and DTS/1900063-ISBIFLOW. AG, EMcV, AK and JCdA have been partially supported by the US National Institutes of Health, under grant 1R01HL160024. JCdA has been partially supported by the US National Insitutes of Health, under grant number 1R01HL158667. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Blood coagulation, or clotting, is a highly regulated mechanism vital in sealing wounded blood vessels to prevent bleeding. Abnormal or excessive clotting can result in serious medical conditions, such as stroke and deep vein thrombosis. Therefore, blood coagulation has been the subject of extensive research, and understanding its mechanisms is crucial to diagnose and manage numerous diseases.
The initiation of blood coagulation is an enzymatic cascade that amplifies thrombin concentration in plasma, activating fibrin polymerization to form a clot [1, 2]. This cascade can be started via extrinsic and intrinsic pathways. The extrinsic pathway is triggered by a vessel-injury-mediated release of coagulation factor VII and tissue factor (TF) into the bloodstream. The intrinsic pathway is auto-initiated, i.e., it does not require exposure to an extravascular tissue factor, and begins with the activation of plasma factors XII, XI, IX, and VIII. Both pathways eventually activate factor X, converging into the common pathway that amplifies thrombin, which in turn converts fibrinogen into fibrin filaments and activates factor XIII, which cross-links the fibrin mesh [3, 4]. The initiation, propagation, and inhibition of this complex process involve a network of over 80 known biochemical reactions [5].
Since thrombosis is a ubiquitous complication of cardiovascular diseases and device implantation [6, 7], there is an abundance of computational models considering coagulation in diverse physiological and anatomical settings. These models may be coupled to computational fluid dynamics (CFD) solvers to capture the interactions between the abnormal blood flow, hypercoagulability, and vessel injury characteristic of thrombus formation. However, the problem is very complicated, and state-of-the-art studies including all interactions are limited. Intra-luminal thrombogenesis in an abdominal aortic aneurysm has been studied using an 18-equation model based on tissue factor activation [8, 9], providing an integrated mechanochemical picture of the process [10]. The same model was used later to study thrombogenesis in an infarcted left ventricle [11]. More complex models have also been used in the literature, either increasing the number of reaction equations involved in the model [12], or including the effects of tissue factor, platelet activation, and clot porosity on thrombus growth [13]. Other models have analyzed the contribution of intrinsic and extrinsic pathways at different timescales, highlighting the multi-stage character of the coagulation process [14].
Eulerian-Lagrangian approximations can be used to integrate coagulation cascade reaction-advection-diffusion equations with platelet activation and deposition. Researchers have used this approach to describe the evolution of non-activated and activated platelet concentrations [15], using the simulated velocity fields to track platelet activation and accumulation. The equations were solved using a stabilized finite element method. The accumulation model [16] accounts for various factors, including plasma-phase and membrane-phase reactions, coagulation inhibitors, and the presence of activated and unactivated platelets. Other groups have used similar strategies, like coupling a calibrated platelet aggregation model, which accounts for adhesion forces between platelet-platelet and platelet-wall at low and high shear rate levels, with an extrinsic coagulation cascade initiation model [17]. In this case, the coagulation cascade was based on a model using 23 chemical species [18].
The large number of coupled partial differential equations (PDEs) representing the reaction, advection, and diffusion of the species involved in the coagulation cascade creates stringent requirements for numerical simulation. This problem is aggravated by the high numerical cost of solving each PDE, owing to the disparate timescales associated with the flow, reaction kinetics, and diffusion [19–21]. The dimensional parameters involved in these processes for mid-size arteries are: flow velocity, Uc ∼ 10cm/s; vessel diameter, Lc ∼ 1cm; the cardiac cycle’s period, tc ≈ 1s; timescales of enzymatic reaction kinetics, ranging from a few seconds to hundreds of seconds (tr ∼ 102s) [9, 20, 22]; and mass diffusivity coefficients for species in blood, Di ∼ 10−6cm2/s [10, 23, 24]. The relatively slow reaction times in these systems require running simulations over many cardiac cycles to reach convergence, i.e., tr ≫ tc. Moreover, the slow diffusion of reactive species creates extremely thin layers in their concentration fields since the Peclet number, Pe = UcLc/Di = td/ta ≈ 107, which measures the ratio between convective and diffusive transport, is very large. For reference, the Schmidt number, which measures the ratio between the viscous diffusion acting in the Navier-Stokes equations and the mass diffusivity acting in the coagulation system’s equations, is Sc = ν/Di ∼ 104 where ν ≈ 4 × 10−2cm2/s is the kinematic viscosity of blood. Consequently, the spatial discretization of the coagulation system’s equations requires much finer computational grids than the ones used to solve the Navier-Stokes equations. This problem is well described in CFD literature and it is common to many other reactive and non-reactive cardiovascular transport problems [25].
Previous simulations of the coagulation cascade under flow have proposed concessions to reduce computational cost. First, the multi-scale nature of reaction kinetics have been simplified by assuming that fast-reacting species are in equilibrium, leading to reduced models with fewer chemical species [23], or by using phenomenological models [26]. The latter have been used to investigate hypercoagulability in the left heart by considering fibrin production from fibrinogen and thrombin [27]. Second, in most if not all studies, the Peclet number has been decreased explicitly by prescribing unphysically high values for Di [10, 11], or implicitly by using a diffusive numerical discretization (e.g., upwinding first-order finite differences). In non-reactive transport problems, this concession has been justified on the basis of accounting for additional noise sources and as long as the effective Pe remains ≫ 1 [28]. However, its adequacy is more questionable in reacting problems like the coagulation cascade, where the reaction timescale is intermediate between the transport and diffusive timescales. Another compromise adopted by some authors is to shorten the time integration to a few cardiac cycles, focusing only on the initial phases of thrombin activation [11]. Finally, many studies have considered idealized two-dimensional geometries to save computational cost [10, 13, 20, 21, 27]. In summary, there is an unmet need for computationally efficient strategies to model the coagulation cascade under flow.
We introduce a multi-fidelity modeling approach to significantly reduce the computational cost of coagulation cascade simulations in flowing blood and make this cost independent of the number of chemical species. The multi-fidelity approach transforms the reaction-advection-diffusion equations for species concentrations (ui, i = 1, …, N) into a system of ODEs by using blood residence time () as the independent variable. The resulting model requires integrating only one PDE for
and N ODEs for the coagulatory species, whose concentration fields can be mapped as
. The transformation is exact for zero diffusivity. For small, finite diffusivity, the model can be Taylor-expanded in terms of the residence time statistical moments, i.e.,
, to derive a family of customizable, multi-fidelity models that offer a balance between cost and accuracy. We compare 1st-order (MuFi-1: 1 PDE, N ODEs) and 2nd-order (MuFi-2: 2 PDEs, N ODEs) multi-fidelity models with the high-fidelity (HiFi: N-PDEs) model for pulsatile flow through an aneurysm-like geometry, using a 9-species coagulation system [23] as a benchmark. Overall, the MuFi-1 and HiFi models show good agreement up to t ≈ 10tc cardiac cycles, while this agreement is improved and extended to longer times (t ≈ 20tc) for MuFi-2. The proposed family of multi-fidelity coagulation models could benefit researchers in the field by enabling them to simulate and analyze complex blood coagulation phenomena more quickly and accurately, thus advancing our understanding of the underlying mechanisms and informing clinical practice.
Methods
This section presents a multi-fidelity model to reduce the cost of simulating coagulation networks of N species in flowing blood. To facilitate the model’s presentation, we first review the standard high-fidelity model (HiFi: N-PDE) and introduce its first-order (MuFi-1: 1-PDE, N-ODE) approximation, which neglects diffusion. We then introduce the second-order approximation (MuFi-2: 2-PDE, N-ODE) accounting for small, finite diffusion, define our benchmark flow problem and coagulation reaction system, and describe the numerical discretization methods.
High-fidelity and first-order multi-fidelity coagulation models
We consider blood as a continuum in space and time (x, t) flowing with velocity v(x, t), and model its coagulation by the system of reaction-advection-diffusion equations
(1)
where D/Dt denotes material derivative, the subindex i indicates each of the N species involved in coagulation system, ui(x, t) its concentration field, Ri(u1, u2, …, uN) its reaction rate from chemical kinetics, and Di its diffusivity coefficient. This system of N PDEs is denoted the high-fidelity (HiFi) model. Assuming that v(x, t) is known, this HiFi model can be solved with some appropriate initial and boundary conditions for ui. For simplicity, we consider uniform initial conditions ui(x, 0) = ui,0, Dirichlet boundary conditions at the domain flow inlets (ui = ui,0), and homogeneous Neumann boundary conditions (∂ui/∂n = 0) at solid surfaces and flow outlets.
The Eq (1) can be written in non-dimensional form using the flow velocity scale Uc and vessel length scale Lc,
(2)
where τ = tUc/Lc is a dimensionless time variable,
is a dimensionless reaction rate, the Damköhler number Da = Lc/(trUc) measures the relative importance of reaction kinetics (rate
) and convective terms, and the Péclet number Pe = UcLc/Di measures the relative importance of convection over diffusion. Using typical values corresponding to mid-size arteries and the reaction rate and diffusivity of coagulation cascade species (i.e. Uc ∼ 10 cm/s, Lc ∼ 1 cm, tr ∼ 102 s, Di ∼ 10−6 cm2/s) yields Da ∼ 10−3 and 1/Pe ∼ 10−7, so both terms on the right-hand side of Eq (2) are small. However, the reaction term is the only forcing in the equation and cannot be neglected, while the diffusive term is even smaller and can be neglected. Doing so simplifies Eq (1) to
(3)
where gi approximates ui in the limit of zero molecular diffusivity. The simplified model has no mixing [29], making the reaction rate within each fluid element independent of the species concentrations in surrounding elements, and exclusively dependent on its age. Consequently, it should be possible to write an ODE system for the coagulation system of each fluid element in a Lagrangian frame that follows the element as it moves with the flow. To avoid the complication of tracking Lagrangian trajectories, we follow previous works [30–33] and resort to the PDE governing the residence time
(4)
which can be used to calculate the age of fluid elements in the region of interest. Applying the chain rule on Eq (3) and taking into account Eq (4), we obtain
(5)
which simplifies the HiFi PDE system Eq (1) into the ODE system
(6)
We note that this ODE system is the same for all fluid elements, as it has no explicit dependence on x, and each fluid element’s pathline information is implicitly encoded by the spatial dependence of
. This model involves solving a system of N ODEs (i.e., Eq (6) for i = 1 … N) to obtain
, solving one PDE to calculate
, and mapping
. Because it combines solving ODEs and PDEs, we consider this approximate model a multi-fidelity (MuFi) model. When N is large, the MuFi model can be significantly cheaper to run than the HiFi model, which involves solving N PDEs.
Higher-order multi-fidelity approximations
In the previous section, we used overline notation for residence time to emphasize that as defined in Eq (4) is the ensemble average age of all molecules within a fluid element, defined as the first integral moment of its probability density function, fT:
(7)
Considering residence time as a stochastic variable is important because, as discussed in the Introduction, the numerical solutions to transport PDEs like Eq (4) explicitly include unphysically large diffusivities or employ discretization methods that implicitly introduce numerical diffusivity. Numerical dissipation helps control instabilities and spurious oscillations but it is a source of numerical error [34]. This error makes
satisfy an equivalent differential equation (EDE) instead of the theoretical PDE given by Eq (4). The specific form of the EDE depends on the details of the temporal integration scheme and the approximation used to discretize the spatial derivatives. For commonly used, first-order, dissipative methods, the EDE for
takes the form
(8)
where the coefficient Dn ∼ cΔx/2 represents the diffusivity of the numerical method, where c is the flow characteristic velocity and Δx the mesh spatial resolution [34].
In this effective scenario, diffusivity creates uncertainty in the residence time. This phenomenon can be shown using Itô’s differentiation [35] to derive the EDE for the second-order moment of the residence time, , as
(9)
where
(10)
Then, it is possible to express this equation in terms of the residence time variance,
, as
(11)
The interested reader can find the derivation of Eqs (8) and (11) in S1 Appendix. Of note, when Dn = 0 and fT is a Dirac delta function, Eqs (8) and (9) yield
and zero residence time variance, i.e.,
. But when Dn ≠ 0, the non-negative forcing term in Eq (11) causes
to increase unless residence time is constant. Note also that for higher order numerical methods, the differential operator multiplying Dn in the EDEs (8) and (9) will involve higher order derivatives, with a dispersive or dissipative character depending on the order of the temporal and spatial discretizations. This will result in more complex evolution equations for
, without changing the fact that numerical diffusion results in an increase on
in regions with strong gradients of
.
The growth of causes errors in the MuFi model derived in the previous section as time increases. To illustrate these errors and derive higher-order corrections, it is convenient to express the concentration of chemical species as
(12)
where Ui is the concentration’s statistical variable,
is its probability density function, and fT(T;x, t) is the probability density function of the residence time. The substitution Ui = gi(T) is warranted because, in the absence of diffusion, σT is zero and fT is a Dirac delta, yielding
, consistent with the definition of gi in Eq (3).
Assuming next that gi(T) is second-order differentiable, we can Taylor-expand it around and integrate the expansion to obtain
(13)
where the primes denote derivatives. This result demonstrates that the MuFi model
error grows with the residence time variance for non-linear reaction systems (i.e., those with g″ ≠ 0). But, more important, it provides a high-order correction that is also an inexpensive MuFi model as it only requires solving one additional PDE.
In summary, we introduce two multi-fidelity models that balance cost with accuracy:
- First-order (MuFi-1):
- Second-order (MuFi-2):
This procedure can be extended to higher order approximations by retaining additional terms in the Taylor expansion of gi around , and solving additional PDEs for higher-order moments of the residence time. For example, the MuFi model of third order (involving N-ODEs and three PDEs) is derived in S2 Appendix. Finally, we note that the PDEs to be solved in the MuFi-2 model are the true PDEs (i.e., Eqs (4) and (14)), not the EDEs (Eqs (8) and (9)). If the PDEs explicitly include diffusivity, then a diffusive term should be added. If not, the resulting
will capture the effect of the discretization’s numerical diffusivity, Dn, regardless of the order of the numerical method employed to solve Eqs (4) and (14). We emphasize that explicit knowledge of Dn is not required, which is advantageous since this diffusivity may be constant or vary in space and time depending on the numerical discretization.
Computational cost estimates
We estimate the computational cost of running a HiFi model in D dimensions to be φNnD+1 floating point operations (FLOPs), where N is the number of species, φ ∼ O(102) is a parameter that depends on the numerical discretization scheme and the function evaluations needed to calculate the reaction rates, and n is the number of elements in the spatial mesh along each direction. Note that the exponent D + 1 reflects that the temporal resolution is linked to the spatial resolution by the CFL condition to guarantee an accurate temporal integration. As a consequence, the number of time steps required to reach a finite integration time is proportional to n. On the other hand, a MuFi-p model is estimated to require βNn + θpnD+1 FLOPs, where the first term is the cost of integrating the N ODEs for the species reaction kinetics and the second term is the cost of solving the p PDEs governing the statistical moments of residence time. Like in the HiFi model, the parameters β and θ depend on the numerical implementation and the right-hand-side terms in the governing equations. It is worth noting that θ < φ because the forcing terms in the residence time equations (see e.g., Eqs 4 and 14) are simpler than the reaction rate terms in the equations governing species concentration (see S4 Appendix).
Based on these estimates and assuming that n is a large number as in, e.g., spatially resolved simulations of flow through arteries or the cardiac chambers, we make two remarks. First, we note that the cost of MuFi models becomes effectively independent of the number of species, N, since D ≥ 1. Second, we note that MuFi models achieve a speedup ∼ (φ/θ)(N/p) > N/p. Furthermore, MuFi models significantly reduce the memory allocation necessary to run simulations, which could lead to additional speedups in parallel implementations by avoiding the overheads associated with message passing and loss of cache memory coherence.
Test case: Coagulation cascade in an idealized aneurysm
To compare the multi-fidelity models MuFi-1 and MuFi-2 with the high-fidelity model based on the PDE system (1), we considered a simplified coagulation cascade model under pulsatile flow through an idealized two-dimensional geometry (Fig 1). This flow geometry broadly resembles a cerebral aneurysm or the left atrial appendage, two cardiovascular sites associated with thrombosis [36–40]. The parent vessel is modeled as a straight tube of diameter H, and the aneurysm is modeled as a circular cavity of radius 0.75H. The center of the cavity is located such that the aneurysm neck size is H. The corners at the aneurysm neck are smoothed with a radius of curvature equal to 0.067H, to avoid sharp corners in the geometry. The pulsatile flow was driven by imposing a two-dimensional Womersley flow as the inflow boundary condition (see S3 Appendix). The Reynolds and Womersley numbers are Re = UcH/ν = 500 and , respectively, where Uc is the maximum velocity. These values of Re and α are representative of intracranial saccular aneurysms [41]. Fig 2 shows the time history of the mass flow rate and the inlet velocity profiles at two time instants, corresponding to the minimum (t/tc = 0.5) and maximum (t/tc = 1) flow rates through the vessel. While the waveform of Q(t) does not include all the temporal complexity of physiological waveforms, it does include the acceleration and deceleration phases needed to drive the flow in the cavity, as described in section Results.
Idealized aneurysm geometry.
A: Time evolution of the mass flow rate through the vessel. B and C: Velocity profile at the inlet at t = 0.5tc and t = tc.
To model the coagulation cascade, we chose a 9-species system considering prothrombin (II), thrombin (IIa), fibrin (Ia), PCa, and factors XIa, IXa, Xa, VIIIa, and Va [23]. The corresponding source terms in Eq (1) are detailed in S4 Appendix. The equations describe the activation of prothrombin (II) into thrombin (IIa) by factors Va and Xa. In turn, thrombin activates the production of factors Va, VIIIa, XIa, and the system’s main inhibitor, PCa. The reaction rates for the model are adopted from previous works [23, 42] and are reported in Table A in S4 Appendix.
Fig 3 illustrates the evolution of the 9 species in stagnant blood (i.e., v = 0) for uniform initial concentrations (see Table 1), chosen within physiologically plausible ranges to ensure substantial thrombin growth within 20 cardiac cycles. This timescale aligns with the peak residence time values observed in the left atrial appendage (LAA) [43, 44]. These conditions lead to an accumulation of thrombin and factors Xa and VIIIa over the initial 10–17 cardiac cycles, followed by a rapid decrease in thrombin concentration over the subsequent 15 cycles.
Obtained by solving Eq (6) for the 9-species coagulation model.
Numerical methods
We used the in-house code TUCAN [45, 46] to solve the Navier-Stokes equations for Newtonian, incompressible flow in the configuration described in the previous section. The numerical discretization used second-order finite differences on a Cartesian staggered grid. The temporal integration was performed with a three-stage, low-storage, semi-implicit Runge-Kutta scheme. The no-slip boundary condition at the vessel walls was modeled by the immersed boundary method [47]. After discarding initial transient effects, the velocity field (v(x, t)) computed by TUCAN was sampled at constant time intervals (tsamp = tcycle/35) and stored to be linearly interpolated for integrating Eqs (1), (4) and (14). To assess the convergence of the velocity field, we performed a grid refinement study employing four resolutions, Δx/H = 1/38, Δx/H = 1/75, Δx/H = 1/150 and Δx/H = 1/300. Each simulation was run with a constant time step Δt that ensured the Courant number to be CFL = max(|u(x)|)Δt/Δx ≈ 0.1. The relative error was defined with respect to the case with Δx/H = 1/300,
(15)
where
is the averaged absolute value of the vorticity in the cavity computed with a spatial resolution Δx = H/k, and Ωcav is the volume of the cavity. Table 2 displays the values of the relative error for each resolution at three time points (t/tc = 0, t/tc = 0.33, and t/tc = 0.67). The case with resolution Δx/H = 1/150 was selected for the present study because it yielded a relative error consistently below 1% throughout the cardiac cycle.
A third-order weighted essentially non-oscillatory (WENO) scheme [48] was implemented to integrate the advection-reaction-diffusion PDEs (Eqs (1), (4) and (14)) as in previous works [43, 44]. This scheme locally adjusts numerical diffusivity to damp convective fluxes perpendicular to sharp scalar fronts, preventing spurious oscillations while at the same time keeping the overall numerical diffusivity low. The systems of PDEs (1), (4) and (14), and the system of ODEs (6) were integrated in time using an explicit, low-storage, 3-stage Runge Kutta scheme. The grid resolution employed to solve these PDEs was the same used for generating the velocity fields, Δx/H = 1/150, as in previous works [10, 11]. The interested reader can find the corresponding grid resolution study in the S5 Appendix. In the system of PDEs, uniform initial conditions were used for all variables, ui(x, 0) = ui,0, , while for the ODE system the corresponding initial conditions were imposed, ui(0) = ui,0.
Results
Flow patterns and residence time
Pulsatile flow in the parent channel has two distinct phases coinciding with the acceleration and deceleration of the inflow profile prescribed at the inlet (Fig 2). The deceleration phase occurs for 0 ≲ t/tc ≲ 0.5, whereas the acceleration phase comprises the rest of the cycle. Fig 4 shows instantaneous vorticity fields (panels A–C) and streamlines (panels D–F) at three different instants of the cycle. At the onset of deceleration (t = 0, first column in Fig 4), the velocity is maximum in the parent channel and a counter-clockwise vortex is the dominant pattern inside the cavity. As deceleration proceeds (t = 0.33tc, center column in Fig 4), the counter-clockwise vortex is partially sucked into the parent channel, creating a thin jet that drives fluid into the cavity in the downstream neck region while fluid slowly exits the cavity along the rest of the neck region. Finally, during acceleration (t = 0.67tc, last column in Fig 4), a vortex pair appears inside the cavity and pulls fluid from the parent channel near the upstream neck region, ejecting fluid to the channel near the downstream neck region.
A-C: Vorticity and D-E: instantaneous streamlines colored by velocity magnitude, both normalized by its maximum value across the cardiac cycle. G-I: Residence time during the 21th simulation cycle. Each variable is plotted at three different phases of the cardiac cycle, as indicated on top of each panel.
This alternating transport of fluid into and out of the cavity repeats every cycle, generating a layered structure in residence time values reminiscent of the growth rings in a tree trunk. Fig 4G, 4H and 4I highlight the developed layer pattern in the 21th cycle of the simulation, demonstrating that the WENO scheme effectively reproduces sharp gradients over long periods of time without introducing excessive numerical diffusivity. However, albeit small, the WENO scheme does introduce a non-zero Dn and, consequently, this scheme leads to σT > 0. This behavior is apparent in Fig 5, which displays snapshots of
and σT over the period 10tc ≤ t ≤ 21tc. Both variables grow with time, while keeping approximately the same spatial organization.
A-C: spatial distribution of . D-F: spatial distribution of σT/tc. Each variable is plotted at three different cycles. A,D: 11th cycle. B,E: 16th cycle. C,F: 21th cycle.
Fig 6 presents the temporal evolution of and σT at three points along the horizontal diameter of the cavity (i.e., the crosses in Fig 5A). In all three points, the residence time exhibits an initial phase of linear growth,
and σT ≈ 0, as in our previous work [43]. At points that do not receive “fresh” fluid from the parent channel (e.g., the wall point in Fig 6A), this phase should last indefinitely. In our simulations, a small departure from
and σT = 0 becomes noticeable for t ≳ 15tc due to the WENO scheme’s numerical diffusivity. Points exchanging fluid with the parent channel experience a different behavior characterized by two principal features. First,
rises and falls every cardiac cycle as pockets of stagnant and fresh fluid move back and forth over the point of interest. Second, the envelope of
saturates to a maximum value
indicating the time needed for fluid exchange with the parent channel to wash out the local blood pool. At these points, σT follows a similar behavior since fresh blood from the parent channel has not only low
but also low σT. We note that, for long enough simulation times, σT can grow at certain points to be comparable in magnitude to
. This result has implications for multi-fidelity modeling of the coagulation cascade but also for the uncertainty of
itself.
A,C,E: Temporal evolution of (
). B,D,F: Temporal evolution of σT/tc (
). Three locations are considered, indicated with × in Fig 5A: A,B at (x/H, y/H) = (1.78, 1.19); C,D at (x/H, y/H) = (2.5, 1.19); and E,F at (x/H, y/H) = (3.15, 1.19).
Multi-fidelity modeling of the coagulation cascade
The initiation of the coagulation cascade is characterized by the rapid growth of thrombin (IIa), peaking at about 15 cycles (see Fig 3). Fig 7 depicts the spatio-temporal structure of uIIa over the 16th simulation cycle, as obtained by the MuFi-1 (panels A–C) and MuFi-2 (panels D–F) models, as well as by the reference HiFi model (panels G–I). The three models capture the rise of uIIa inside the cavity and the formation of a layered structure similar to that of the residence time. There is a trend for MuFi-1 to underestimate uIIa, which is for the most part corrected by MuFi-2.
A-C: MuFi-1. D-F: MuFi-2. G-I: HiFi reference model. Three phases within the 16th cycle are plotted for each case, as indicated on the top row. A,D,G: t/tc = 15. B,E,H: t/tc = 15.34. C,F,I: t/tc = 15.67.
To study each model’s behavior in more detail, Fig 8 shows the thrombin concentration vs. time at the same three points considered in Fig 6, representative of scenarios where ,
, or
after running the simulation for a large number of cycles, t ≫ tc. For reference, the figure also shows the thrombin concentration obtained from the HiFi model and the ODE system corresponding to no flow and no diffusion, i.e.,
and σT = 0 (Eq 6).
Each line correspons to a different model: MuFi-1 (━), MuFi-2 () and HiFi model (
). For reference, the solution of the 9-ODE system Eq (6) is also included (
). Three locations are considered, indicated with × in Fig 5A: A,B at (x/H, y/H) = (1.78, 1.19); C,D at (x/H, y/H) = (2.5, 1.19); and E,F at (x/H, y/H) = (3.15, 1.19).
Similar to the residence time, uIIa experiences peaks and valleys every cycle due to the periodic fluid exchange between the cavity and the parent channel. In addition, the peak value of uIIa increases from one cycle to the next as the coagulation cascade progresses in the fluid trapped in the cavity. The no-flow ODE model (dashed lines in Fig 8) fails to capture the oscillations of uIIa and severely overestimates the growth of its peak values. The MuFi-1 model captures the oscillatory nature of uIIa, but it begins to underpredict the growth of its peak values after a number of cardiac cycles that varies from point to point. At the first sampled point, where , the MuFi-1 model is almost exact for
and remains fairly accurate for
(Fig 8A). At the second and third sampled points, where
, the MuFi-1 model significantly departs from the HiFi model beyond t ∼ 10tc (Fig 8C and 8E). In contrast, the MuFi-2 model, which incorporates both
and σT, remains in reasonable agreement with the HiFi model much longer than the MuFi-1 model, producing fairly accurate results for simulation times well over 10 cycles at the three sampled points (Fig 8B, 8D and 8F). By t = 20tc, the MuFi-1 model underestimates the peak uIIa by 4.8%, 21.5%, and 15% while the MuFi-2 model does so by 2.7%, 1%, and 4.1%, respectively, at the first, second, and third sampled points.
Next, we compare the spatio-temporal behavior of the high-fidelity and multi-fidelity models by representing in Fig 9 the concentrations of thrombin, prothrombin, and factor Xa (i.e., uIIa, uII, and uXa) along the horizontal cavity diameter depicted in Fig 7A, with y/H = 1.19. Spatial concentration profiles are plotted at time-points t/tc = 10, 15, 20. Factor Xa displays the best agreement between HiFi and MuFi models, followed by the prothrombin (II) and thrombin (IIa). Of note, the MuFi-2 model replicates the HiFi behavior for all three species, capturing their complex spatial oscillations. In comparison, the accuracy of the MuFi-1 model deteriorates faster with simulation time, although this model still retains the qualitative spatial dependence of species concentration.
Data is show at y/H = 1.19 for MuFi-1 (━), MuFi-2 () and HiFi model (
). A,B,C: thrombin, uIIa. D,E,F: prothrombin, uII. G,H,I: factor Xa, uXa. Three different times are plotted for each species, as indicated on the top row.
Multi-fidelity model error analysis
While interesting to understand the performance of the multi-fidelity models, Figs 8 and 9 include examples of extreme behavior that do not reflect these models’ overall accuracy. To systematically quantify the discrepancies between the high-fidelity and multi-fidelity models, we compute the relative errors
(16)
where i stands for species and p for MuFi order. Fig 10 shows the thrombin relative errors of the MuFi-1 and MuFi-2 models,
and
, for three instants along the 16th simulation cycle. These errors are negligible in the parent channel but reach appreciable values inside the cavity, where they exhibit a layered pattern with strong gradients, similar to
and σT. Also,
is significantly higher than
, consonant with the data shown in Figs 7–9.
The relative error is defined in Eq (16). A-C: MuFi-1. D-F: MuFi-2. Three phases within the 16th cycle are plotted for each case, as indicated on the top row. A,D: t/tc = 15. B,E: t/tc = 15.34. C,F: t/tc = 15.67.
To characterize the dependence of the MuFi errors on the residence time and its variance inside the cavity, we divide the plane in bins, ensemble-average
inside each bin, and plot the resulting error maps in Fig 11 together with the corresponding normalized bin counts [i.e., the probability density function
]. The error maps obtained for different values of t/tc indicate that the MuFi error increases with σT while it is much less sensitive to
. The MuFi-2 model particularly outperforms the MuFi-1 model in areas of large σT, since MuFi-1 assumes σT = 0. Similar observations can be made for the error maps for factor Xa and prothrombin (II) concentrations, which are provided in S6 Appendix.
A-C: Relative error in thrombin concentration in MuFi-1, , as a function of residence time and its standard deviation. D-F: Same, but for MuFi-2,
. G-I: Joint probability density function of residence time (
) and its standard deviation (σT). Data for all panels is compiled inside cavity during three different cycles, as indicated on the top row. A,D,G: 10th cycle. B,E,H: 15th cycle. C,F,I: 20th cycle.
Inspection of the probability density functions shows that the majority of the points inside the cavity are circumscribed to two regions in the plane. The region near the point
corresponds to locations within the cavity that receive periodic inflows of fresh flow from the parent channel. This periodic infusion sustains a low level of
throughout each cycle, since fresh flow has small values of
and
. The region near
corresponds to the stagnant areas of recirculating flow inside the cavity (see Fig 5). The modal errors in this region (black crosses in Fig 11G–11I) are 0.065, 0.141, and 0.187 for MuFi-1 and 0.004, 0.002, and 0.002 for MuFi-2 at t/tc = 9.5, 14.5, and 19.5, respectively. These stagnant flow areas displace upwards and towards the right in the
plane as the simulation time advances (see bottom row of Fig 11). Consequently, we expect σT and
to grow with t/tc inside the cavity and the MuFi-2 model to outperform the MuFi-1 model.
To test this hypothesis and quantify the overall error of MuFi models, we calculate the spatial average of over the cavity region,
(17)
and plot it vs. simulation time in Fig 12 for all 9 species. The results show differences among species but, overall, the MuFi-1 error
grows sharply with time for
, then more gradually for longer times, even plateauing or tapering off in some cases. For most species, the MuFi-2 error
is significantly lower than
over extended periods of time or both errors are very low. For example, after 20 integration cycles, MuFi-2 is about an order of magnitude more accurate than MuFi-1 for thrombin and pro-thrombin, i.e.,
vs.
, and
vs.
, whereas both models yield similar low errors for Factor XIa, i.e.,
. The species showing the highest relative errors is PCa, reaching
and
at t = 20tc. However, both the MuFi-1 error at t = 10tc and the MuFi-2 error at t = 20tc were well under 0.1 for all other species.
Each line correspons to a different model: MuFi-1 solid, MuFi-2 dashed. Dashed-dot lines correspond to and
. A: Thrombin (IIa) B: Prohrombin (II) C: Factor Xa D: Factor IXa E: Factor Xa F: Activated protein C (PCa) G: Factor VIIIa H: Factor Va I: Fibrin (Ia).
To assess the dependence of the MuFi errors on coagulation model complexity, we repeated the numerical experiments described above using a 3-species system considering thrombin, factor XIa, and the inhibitor PCa [23]. This simplified system can be derived from the 9-equation system under the assumption that the faster chemical equations are in equilibrium. Overall, the performance of MuFi models was similar on the 3-equation and 9-equation systems (see S7 Appendix). For thrombin, the MuFi-1 and MuFi-2 models yielded respectively and
at t = 20tc in the 3-equation system, comparable to the 9-equation values
and
.
Computational cost
The computational cost of solving 20 cycles of the 9-equation coagulation system using the HiFi approach was approximately 66.67 hours, using the MATLAB codes provided in [49] on an Intel(R) Xeon(R) Silver 4208 CPU @ 2.10GHz. On the same hardware, the MuFi-1 approach required around 4.24 hours, while the MuFi-2 approach took approximately 7.36 hours. This results in a speedup of 15.7 for MuFi-1, and 9.1 for MuFi-1, compared to the HiFi model. When the speedup estimated in Computational cost estimates is taken into account (i.e., speedup ≈ (φ/θ)(N/p)), these numbers imply that the overhead associated to computing the reaction terms Ri in HiFi (instead of the forcing terms for and
) is φ/θ ≈ 2.
Discussion
The computational modeling of the coagulation cascade in flowing blood poses significant challenges due to the multi-scale nature of the process and the large number of involved chemical species. High-fidelity (HiFi) continuum mechanics models of coagulation lead to dozens of reaction-advection-diffusion partial differential equations (PDEs) [50]. The reaction kinetics in these equations are much slower than the cardiac cycle, requiring long simulation times to cover the coagulation process [9, 20, 22]. Moreover, the diffusion of chemical species is orders of magnitude slower than their convection and reaction kinetics, causing sharp scalar fronts that require ultra-high-resolution spatial meshes [25]. Despite the vast advances in computational software and hardware achieved in past decades, these joint requirements still impede the fast simulations of chemo-fluidic coagulation models in arterial or intracardiac domains, hindering the adoption of these models in clinical decision-making. Modelers usually resort to surrogate metrics associated with thrombogenesis, derived from blood residence time or wall shear stress [33, 51–55]. On the other hand, chemo-fluidic models of thrombosis are rare [10, 11] and often suffer from excessive diffusivity (numerical or explicit), short simulation times, and/or simplified coagulation models with few species.
We introduce a family of tailorable multi-fidelity (MuFi) models to effectively decouple the computational cost of simulating the coagulation cascade under flow from the number of reacting species, N. The MuFi models are designed to approximate the HiFi model in the limit of vanishing molecular diffusivity. In this limit, the N-PDEs of the HiFi model can be transformed into ordinary differential equations (ODEs) by changing variables between simulation time and blood residence time. Consequently, the HiFi model is replaced by N ODEs representing the reaction kinetics and p PDEs representing the ensemble mean residence time within each fluid particle, , and p − 1 higher-order statistical moments (namely,
) [56]. We provide a procedure to incrementally derive MuFi models of arbitrary order starting from the first-order model obtained with p = 1, corresponding to zero diffusivity. In the presence of small diffusivity, natural or numerical, higher-order models can be derived by Taylor-expanding the HiFi model around the zero diffusivity limit.
We assess the performance of the MuFi-1 and MuFi-2 models obtained for p = 1, 2 in a well-characterized, simplified coagulation system representing nine species: prothrombin (II), thrombin (IIa), fibrin (Ia), PCa, and factors XIa, IXa, Xa, VIIIa and Va [23]. This coagulation system is evaluated in a pulsatile flow through a two-dimensional geometry consisting of a parent channel driven by a Womersley inflow profile with a laterally protruding cavity where blood becomes stagnant. The non-dimensional parameters governing the flow, i.e., the Reynolds (Re = 500) and Womersley (α = 10) numbers, are representative of a saccular aneurysm in an intermediate-size artery of the adult human circulatory system [57], although the instantaneous flow rate driving the flow does not include all the temporal complexity of physiological waveforms. For the residence time and its variance, the Peclet number is nominally set to be infinitely high by not including mass diffusivity terms in the corresponding PDE equations.
This configuration creates a cyclic influx of fresh fluid from the parent channel to the cavity and efflux of stagnant fluid from the cavity to the channel, together with a swaying fluid motion inside the cavity. Flow transport leads to a complex residence time pattern with thin layers separated by strong gradients, which intensifies as grows with simulation time. Our WENO scheme resolve the thinly layered
pattern for long simulation times. However, albeit low, the WENO scheme has a non-zero numerical diffusivity [48], and thus, the effective Pe in our HiFi simulations is finite, affecting the solution to the transport equations for
and
. Consequently, not only
but also its standard deviation, σT, increases with simulation time.
Overall, the MuFi-1 model compares well with the HiFi model, producing spatially averaged errors for thrombin inside the cavity that remain below 10% for up to 10 simulation cycles even if this model was derived under the assumption of Pe → ∞ and σT = 0. Nevertheless, MuFi-1 starts to underestimate the concentration of all species beyond that point. By t = 20tc, the spatially averaged errors for thrombin inside the cavity are as high as ≈20%. The MuFi-2 model yields significantly lower errors than the MuFi-1 model, with spatially averaged errors for thrombin that remain below 2% for up to t ≈ 20tc. These errors are comparable to those associated to solving the HiFi model’s PDEs (see S5 Appendix). The accuracy of the MuFi models is worst for PCa, the species with the fastest growth rate in the considered coagulation model. Statistical analysis shows that the MuFi errors depend almost exclusively on σT and that this dependence is steeper for MuFi-1 than for MuFi-2.
Because σT increases with numerical diffusivity, Dn, the performance of MuFi implementations is expected to depend on the numerical scheme and the spatio-temporal resolution used to discretize the model PDEs. While an exhaustive analysis of the numerical diffusivity of numerical discretizations of advection-reaction-diffusion problems is beyond the scope of this study, we have outlined a general methodology to formulate MuFi models of arbitrary order, derive each numerical scheme’s EDE for σT, and outlined a step-by-step process to compute σT without the need to derive its EDE. With these tools, it should be straightforward to tailor MuFi models to the peculiarities of each flow and numerical solver. These step-by-step procedures should also apply to models that compute residence time using non-classical numerical representations of the governing PDEs such as, e.g., neural networks.
The MuFi models (N ODEs, p PDEs) are much more efficient than the HiFi model (N PDEs) because solving ODEs is significantly cheaper than solving PDEs. In physiologically relevant computational meshes containing hundreds of elements in each direction, MuFi models achieve a speedup > N/p. We anticipate this speedup to exceed an order of magnitude, considering the favorable accuracy achieved by low-order MuFi models (p = 1, 2) and the dozens of species involved in realistic coagulation cascade models. An alternate interpretation is that the computational cost of MuFi models becomes independent of the number of species N. Additionally, one can run any number of different MuFi models inexpensively for a specified flow baseline. The costly part of MuFi models, i.e., solving the p-PDEs representing residence time and its higher-order moments, does not need to be redone unless the anatomy or inflow/outflow conditions change. Thus, MuFi models are highly efficient for sensitivity analyses, uncertainty quantification, kinetics model comparisons, and any type of study requiring multiple evaluations of the coagulation cascade model.
In the coagulation cascade, the diffusive and chemical timescales are, overall, significantly different, justifying the idea of MuFi models. However, the non-linearity of the reaction rate terms makes it difficult to rule out that diffusion became dominant for some species in some regions of state space, especially as the number of species increases. While this situation could make it necessary to increase the MuFi order as N increases, deteriorating the MuFi speedup, we have seen no indication of it being significant. Specifically, MuFi models reproduced the HiFi results for thrombin with similar accuracy on a classic 9-equation model [23], and its 3-equation simplified version. Investigation of MuFi models with N > 9 should clarify this matter further.
For simplicity and to establish proof of concept of the MuFi strategy, this study focuses on the intrinsic coagulation cascade, ignoring the extrinsic cascade initiated by the release of procoagulatory or inhibitory species from the vessel walls. While both the intrinsic and extrinsic coagulation pathways are often activated in cardiovascular conditions associated with thrombosis [58–60], focusing on the intrinsic pathway is not uncommon in the literature, motivated by the paucity of information about the wall’s prothrombotic potential in many clinically relevant scenarios [15]. Generating MuFi models for the extrinsic pathway would involve additional PDEs representing the residence time of wall-released chemical species or the time spent by blood near damaged wall regions, e.g., the near-wall residence time proposed by others [11]. As long as the number of additional PDEs remains lower than the total number of species, N, the resulting MuFi models would still save computational time. Furthermore, applications requiring multiple HiFi runs, e.g., uncertainty quantification studies, simulations pharmacological treatments, etc., can still be accelerated by MuFi models, even if this required solving a significant number of residence-time-like PDEs. A notable exception would be, however, applications requiring two-way coupling between flow and coagulation.
Finally, note that residence time and its higher order moments can not only be obtained in silico using CFD analysis [43, 44, 51], but also in vitro using experimental techniques like particle image velocimetry [61–63]. In vitro experiments offer opportunities to validate HiFi and MuFi approaches by, e.g., tracking the trajectories of single particle tracers or the evolution of dye injections. Moreover, residence time can also be estimated in vivo using medical imaging modalities like phase-contrast magnetic resonance imaging, Doppler ultrasound, and perfusion computed tomography [64–67]. Residence time maps obtained in vivo could be fed to the ODE component of the MuFi models, providing a computationally tractable method to calculate the spatio-temporal evolution of the coagulation species non-invasively in the clinical setting.
Conclusion
We present a novel multi-fidelity approach to mitigate the computational burden of simulating the coagulation cascade under flow. Using residence time and its statistical moments, the multi-fidelity (MuFi) modeling achieves a favorable trade-off between computational cost and accuracy. The multi-fidelity approach allows for the integration of various data streams, including CFD analysis, in vitro experiments, and non-invasive imaging in vivo techniques, enabling a comprehensive understanding of the spatial and temporal progression of coagulation species and offering promise for clinical translation. Finally, we note that the same multi-fidelity strategy developed here can be adopted to increase the efficiency of simulating other systems biology processes influenced by blood flow.
Supporting information
S2 Appendix. Multi-fidelity model of third order.
https://doi.org/10.1371/journal.pcbi.1011583.s002
(PDF)
S3 Appendix. Womersley inflow boundary condition.
https://doi.org/10.1371/journal.pcbi.1011583.s003
(PDF)
S5 Appendix. Grid resolution study for the HiFi model.
https://doi.org/10.1371/journal.pcbi.1011583.s005
(PDF)
S6 Appendix. Error maps for prothrombin and factor Xa.
https://doi.org/10.1371/journal.pcbi.1011583.s006
(PDF)
S7 Appendix. Averaged relative errors for N = 3.
https://doi.org/10.1371/journal.pcbi.1011583.s007
(PDF)
References
- 1. Macfarlane RG. An enzyme cascade in the blood clotting mechanism, and its function as a biochemical amplifier. Nature. 1964;202(4931):498–499. pmid:14167839
- 2. Davie EW, Ratnoff OD. Waterfall sequence for intrinsic blood clotting. Science. 1964;145(3638):1310–1312. pmid:14173416
- 3. Green D. Coagulation cascade. Hemodial Int. 2006;10(S2):S2–S4. pmid:17022746
- 4.
Robbins SL, Cotran RS. Pathologic basis of disease. Saunders; 1979.
- 5. Runyon MK, Johnson-Kerner BL, Ismagilov RF. Minimal functional model of hemostasis in a biomimetic microfluidic system. Angew Chem. 2004;116(12):1557–1562.
- 6. Wendelboe AM, Raskob GE. Global burden of thrombosis: epidemiologic aspects. Circ Res. 2016;118(9):1340–1347. pmid:27126645
- 7. Raskob GE, Angchaisuksiri P, Blanco AN, Buller H, Gallus A, Hunt BJ, et al. Thrombosis: a major contributor to global disease burden. Arterioscler Thromb Vasc Biol. 2014;34(11):2363–2371. pmid:25304324
- 8. Lawson JH, Kalafatis M, Stram S, Mann KG. A model for the tissue factor pathway to thrombin. I. An empirical study. J Biol Chem. 1994;269(37):23357–23366. pmid:8083241
- 9. Jones KC, Mann KG. A model for the tissue factor pathway to thrombin. II. A mathematical simulation. J Biol Chem. 1994;269(37):23367–23373. pmid:8083242
- 10. Biasetti J, Spazzini PG, Swedenborg J, Gasser TC. An integrated fluid-chemical model toward modeling the formation of intra-luminal thrombus in abdominal aortic aneurysms. Front Physiol. 2012;3:266. pmid:22934022
- 11. Seo JH, Abd T, George RT, Mittal R. A coupled chemo-fluidic computational model for thrombogenesis in infarcted left ventricles. Amer J Physiol-Heart Circul Physiol. 2016;310(11):H1567–H1582. pmid:27016582
- 12. Hockin MF, Jones KC, Everse SJ, Mann KG. A model for the stoichiometric regulation of blood coagulation. J Biol Chem. 2002;277(21):18322–18333. pmid:11893748
- 13. Leiderman K, Fogelson A. Grow with the flow: a spatial–temporal model of platelet deposition and blood coagulation under flow. Math Med Biol. 2011;28(1):47–84. pmid:20439306
- 14. Panteleev MA, Ovanesov MV, Kireev DA, Shibeko AM, Sinauridze EI, Ananyeva NM, et al. Spatial propagation and localization of blood coagulation are regulated by intrinsic and protein C pathways, respectively. Biophys J. 2006;90(5):1489–1500. pmid:16326897
- 15. Grande Gutiérrez N, Alber M, Kahn AM, Burns JC, Mathew M, McCrindle BW, et al. Computational modeling of blood component transport related to coronary artery thrombosis in Kawasaki disease. PLoS Comput Biol. 2021;17(9):e1009331. pmid:34491991
- 16. Kuharsky AL, Fogelson AL. Surface-mediated control of blood coagulation: the role of binding site densities and platelet deposition. Biophys J. 2001;80(3):1050–1074. pmid:11222273
- 17. Yazdani A, Li H, Humphrey JD, Karniadakis GE. A general shear-dependent model for thrombus formation. PLoS Comput Biol. 2017;13(1):e1005291. pmid:28095402
- 18. Anand M, Rajagopal K, Rajagopal KR. A model for the formation, growth, and lysis of clots in quiescent plasma. A comparison between the effects of antithrombin III deficiency and protein C deficiency. J Theor Biol. 2008;253(4):725–738. pmid:18539301
- 19. Wootton DM, Ku DN. Fluid mechanics of vascular systems, diseases, and thrombosis. Annu Rev Biomed Eng. 1999;1(1):299–329. pmid:11701491
- 20. Lobanov AI, Starozhilova TK. The effect of convective flows on blood coagulation processes. Pathophysiol Haemos Thromb. 2005;34(2-3):121–134. pmid:16432313
- 21. Ermakova EA, Panteleev MA, Shnol EE. Blood coagulation and propagation of autowaves in flow. Pathophysiol haemost thromb. 2005;34(2-3):135–142. pmid:16432314
- 22. Ataullakhanov FI, Krasotkina YV, Sarbash VI, Volkova RI, Sinauridse EI, Kondratovich AY. Spatio-temporal dynamics of blood coagulation and pattern formation: an experimental study. Int J Bifurcat Chaos. 2002;12(09):1969–1983.
- 23. Zarnitsina VI, Ataullakhanov F, Lobanov AI, Morozova OL. Dynamics of spatially nonuniform patterning in the model of blood coagulation. Chaos. 2001;11(1):57–70. pmid:12779441
- 24. Ratto N, Tokarev A, Chelle P, Tardy-Poncet B, Volpert V. Clustering of thrombin generation test data using a reduced mathematical model of blood coagulation. Acta Biotheor. 2020;68(1):21–43. pmid:31853681
- 25. Kaazempur-Mofrad MR, Ethier CR. Mass transport in an anatomically realistic human right coronary artery. Ann Biomed Eng. 2001;29:121–127. pmid:11284666
- 26. Ataullakhanov FI, Zarnitsyna VI, Kondratovich AY, Lobanova ES, Sarbash VI. A new class of stopping self-sustained waves: a factor determining the spatial dynamics of blood coagulation. Phys-Usp+. 2002;45(6):619.
- 27.
Qureshi A, Balmus M, Nechipurenko D, Ataullakhanov F, Williams S, Lip G, et al. Left Atrial Appendage Morphology Impacts Thrombus Formation Risks in Multi-Physics Atrial Models. In: 2021 Computing in Cardiology (CinC). vol. 48. IEEE; 2021. p. 1–4.
- 28. Ford MD, Stuhne GR, Nikolov HN, Habets DF, Lownie SP, Holdsworth DW, et al. Virtual angiography for visualization and validation of computational models of aneurysm hemodynamics. IEEE Trans Med Imaging. 2005;24(12):1586–1592. pmid:16350918
- 29. Villermaux E. Mixing versus stirring. Annu Rev Fluid Mech. 2019;51:245–273.
- 30.
Józsa J, Krámer T. Modelling residence time as advection-diffusion with zero-order reaction kinetics. In: Proceedings of the Hydrodynamics 2000 Conference, International Association of Hydraulic Engineering and Research. Citeseer; 2000. p. 23–27.
- 31. Mangual JO, Domenichini F, Pedrizzetti G. Describing the highly three dimensional right ventricle flow. Ann Biomed Eng. 2012;40:1790–1801. pmid:22396043
- 32. Esmaily-Moghadam M, Hsia TY, Marsden AL. A non-discrete method for computation of residence time in fluid mechanics simulations. Phys Fluids. 2013;25(11). pmid:24046509
- 33. Rossini L, Martinez-Legazpi P, Vu V, Fernandez-Friera L, Perez del Villar C, Rodriguez-Lopez S, et al. A clinical method for mapping and quantifying blood stasis in the left ventricle. J Biomech. 2016;49(11):2152–2161. pmid:26680013
- 34.
Hirsch C. Numerical computation of internal and external flows. (2nd edition). Elsevier; 2007.
- 35.
Itô K, Henry P Jr, et al. Diffusion processes and their sample paths: Reprint of the 1974 edition. Springer Science & Business Media; 1996.
- 36. Ngoepe MN, Frangi AF, Byrne JV, Ventikos Y. Thrombosis in cerebral aneurysms and the computational modeling thereof: a review. Front Physiol. 2018;9:306. pmid:29670533
- 37. Cohen JE, Yitshayek E, Gomori JM, Grigoriadis S, Raphaeli G, Spektor S, et al. Spontaneous thrombosis of cerebral aneurysms presenting with ischemic stroke. J Neurol Sci. 2007;254(1-2):95–98. pmid:17258773
- 38. Vanninen R, Manninen H, Ronkainen A. Broad-based intracranial aneurysms: thrombosis induced by stent placement. Am J Neuroradiol. 2003;24(2):263–266. pmid:12591645
- 39. Al-Saady NM, Obel OA, Camm AJ. Left atrial appendage: structure, function, and role in thromboembolism. Heart. 1999;82(5):547–554. pmid:10525506
- 40. Goldman ME, Pearce LA, Hart RG, Zabalgoitia M, Asinger RW, Safford R, et al. Pathophysiologic correlates of thromboembolism in nonvalvular atrial fibrillation: I. Reduced flow velocity in the left atrial appendage (The Stroke Prevention in Atrial Fibrillation [SPAF-III] study). J Am Soc Echocardiog. 1999;12(12):1080–1087.
- 41. Ku DN. Blood flow in arteries. Ann Rev Fluid Mech. 1997;29(1):399–434.
- 42. Zarnitsina VI, Pokhilko AV, Ataullakhanov FI. A mathematical model for the spatio-temporal dynamics of intrinsic pathway of blood coagulation. I. The model description. Thromb Res. 1996;84(4):225–236. pmid:8948047
- 43. Garcia-Villalba M, Rossini L, Gonzalo A, Vigneault D, Martinez-Legazpi P, Durán E, et al. Demonstration of patient-specific simulations to assess left atrial appendage thrombogenesis risk. Frontiers Physiol. 2021;12:596596. pmid:33716763
- 44. Gonzalo A, García-Villalba M, Rossini L, Durán E, Vigneault D, Martínez-Legazpi P, et al. Non-Newtonian blood rheology impacts left atrial stasis in patient-specific simulations. Int J Numer Method Biomed Eng. 2022;38(6):e3597. pmid:35344280
- 45. Moriche M, Flores O, Garcia-Villalba M. On the aerodynamic forces on heaving and pitching airfoils at low Reynolds number. J Fluid Mech. 2017;828:395–423.
- 46.
Flores O, Rossini L, Gonzalo A, Vigneault D, Bermejo J, Kahn AM, et al. Evaluation of blood stasis in the left atrium using patient-specific direct numerical simulations. In: ERCOFTAC Workshop Direct and Large Eddy Simulation. Springer; 2019. p. 485–490.
- 47. Uhlmann M. An immersed boundary method with direct forcing for the simulation of particulate flows. J Comput Phys. 2005;209(2):448–476.
- 48. Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. J Comput Phy. 1996;126(1):202–228.
- 49.
Guerrero-Hurtado M, Flores O. MultiFidelity models for coagulation cascade (in MATLAB); 2023. Available from: https://doi.org/10.5281/zenodo.8344615.
- 50. Leiderman K, Fogelson A. An overview of mathematical modeling of thrombus formation under flow. Thromb Res. 2014;133:S12–S14. pmid:24759131
- 51. Rayz VL, Boussel L, Ge L, Leach JR, Martin AJ, Lawton MT, et al. Flow residence time and regions of intraluminal thrombus deposition in intracranial aneurysms. Ann Biomed Eng. 2010;38(10):3058–3069. pmid:20499185
- 52. Di Achille P, Tellides G, Figueroa CA, Humphrey JD. A haemodynamic predictor of intraluminal thrombus formation in abdominal aortic aneurysms. Proc Math Phys Eng Sci. 2014;470(2172):20140163.
- 53. Kelsey LJ, Powell JT, Norman PE, Miller K, Doyle BJ. A comparison of hemodynamic metrics and intraluminal thrombus burden in a common iliac artery aneurysm. Int J Numer Method Biomed Eng. 2017;33(5):e2821.
- 54. O’Rourke MJ, McCullough JP, Kelly S. An investigation of the relationship between hemodynamics and thrombus deposition within patient-specific models of abdominal aortic aneurysm. Proc Inst Mech Eng H. 2012;226(7):548–564. pmid:22913102
- 55. Arzani A, Suh GY, Dalman RL, Shadden SC. A longitudinal comparison of hemodynamics and intraluminal thrombus deposition in abdominal aortic aneurysms. Am J Physiol Heart Circ Physiol. 2014;307(12):H1786–H1795. pmid:25326533
- 56. del Alamo JC, Rossini L, Kahn A, Bermejo J, Martínez-Legazpi P, Alvarez RY. Mapping and quantifying blood stasis and thrombus risk in the heart; 2020.
- 57.
Vlachopoulos C, O’Rourke M, Nichols WW. McDonald’s blood flow in arteries: theoretical, experimental and clinical principles. CRC press; 2011.
- 58. Ding WY, Gupta D, Lip GYH. Atrial fibrillation and the prothrombotic state: revisiting Virchow’s triad in 2020. Heart. 2020;106(19):1463–1468. pmid:32675218
- 59. Camaj A, Fuster V, Giustino G, Bienstock SW, Sternheim D, Mehran R, et al. Left ventricular thrombus following acute myocardial infarction: JACC state-of-the-art review. J Am Coll Cardiol. 2022;79(10):1010–1022. pmid:35272796
- 60. Cameron SJ, Russell HM, Owens AP III. Antithrombotic therapy in abdominal aortic aneurysm: beneficial or detrimental? Am J Hematol. 2018;132(25):2619–2628. pmid:30228233
- 61. Mareels G, Kaminsky R, Eloot S, Verdonck PR. Particle image velocimetry–validated, computational fluid dynamics–based design to reduce shear stress and residence time in central venous hemodialysis catheters. Asaio J. 2007;53(4):438–446. pmid:17667228
- 62. Falahatpisheh A, Kheradvar A. High-speed particle image velocimetry to assess cardiac fluid dynamics in vitro: From performance to validation. Eur J Mech B/Fluids. 2012;35:2–8.
- 63. Tomaszewski M, Sybilski K, Baranowski P, Małachowski J. Experimental and numerical flow analysis through arteries with stent using particle image velocimetry and computational fluid dynamics method. Biocybern Biomed Eng. 2020;40(2):740–751.
- 64. Steingoetter A, Weishaupt D, Kunz P, Mäder K, Lengsfeld H, Thumshirn M, et al. Magnetic resonance imaging for the in vivo evaluation of gastric-retentive tablets. Pharm Res. 2003;20:2001–2007. pmid:14725366
- 65. Li Y, Amili O, Moen S, Van de Moortele P, Grande A, Jagadeesan B, et al. Flow residence time in intracranial aneurysms evaluated by in vitro 4D flow MRI. J Biomech. 2022;141:111211. pmid:35780698
- 66. Hendabadi S, Bermejo J, Benito Y, Yotti R, Fernández-Avilés F, del Álamo JC, et al. Topology of blood transport in the human left ventricle by novel processing of Doppler echocardiography. Ann Biomed Eng. 2013;41:2603–2616. pmid:23817765
- 67. Maidu B, Gonzalo A, Bargellini C, Rossini L, Vigneault D, Martinez-Legazpi P, et al. Inferring left atrial appendage (LAA) hemodynamics from 4D CT contrast dynamics by physics informed neural networks (PINNs); 2022. Bulletin of the American Physical Society.