11affiliationtext: Department of Mathematics, University of Utah, Salt Lake
City, UT, USA22affiliationtext: Department of Mathematics, University of North
Carolina, Chapel Hill, NC, USA33affiliationtext: Advanced Medical Imaging Lab,
University of North Carolina Medical Center, Chapel Hill, NC,
USA44affiliationtext: University of North Carolina School of Medicine, Chapel
Hill, NC, USA55affiliationtext: Division of Cardiology, Department of
Medicine, University of North Carolina, Chapel Hill, NC, USA66affiliationtext:
Department of Biomedical Engineering, University of California Irvine, Irvine,
CA, USA77affiliationtext: Departments of Mathematics, Applied Physical
Sciences, and Biomedical Engineering, University of North Carolina, Chapel
Hill, NC, USA88affiliationtext: Carolina Center for Interdisciplinary Applied
Mathematics, University of North Carolina, Chapel Hill, NC,
USA99affiliationtext: Computational Medicine Program, University of North
Carolina, Chapel Hill, NC, USA1010affiliationtext: McAllister Heart Institute,
University of North Carolina, Chapel Hill, NC, USA1111affiliationtext:
Departments of Mathematics and Biomedical Engineering, University of Utah,
Salt Lake City, UT, USA**affiliationtext<EMAIL_ADDRESS>
# A Model of Fluid-Structure and Biochemical Interactions for Applications to
Subclinical Leaflet Thrombosis

Aaron Barrett  Jordan A. Brown  Margaret Anne Smith  Andrew Woodward  John P.
Vavalle  Arash Kheradvar  Boyce E. Griffith  Aaron L. Fogelson

###### Abstract

Subclinical leaflet thrombosis (SLT) is a potentially serious complication of
aortic valve replacement with a bioprosthetic valve in which blood clots form
on the replacement valve. SLT is associated with increased risk of transient
ischemic attacks and strokes and can progress to clinical leaflet thrombosis.
SLT following aortic valve replacement also may be related to subsequent
structural valve deterioration, which can impair the durability of the valve
replacement. Because of the difficulty in clinical imaging of SLT, models are
needed to determine the mechanisms of SLT and could eventually predict which
patients will develop SLT. To this end, we develop methods to simulate leaflet
thrombosis that combine fluid-structure interaction and a simplified
thrombosis model that allows for deposition along the moving leaflets.
Additionally, this model can be adapted to model deposition or absorption
along other moving boundaries. We present convergence results and quantify the
model’s ability to realize changes in valve opening and pressures. These new
approaches are an important advancement in our tools for modeling thrombosis
in which they incorporate both adhesion to the surface of the moving leaflets
and feedback to the fluid-structure interaction.

##  1 Introduction

Subclinical leaflet thrombosis is a potentially serious complication of
bioprosthetic aortic valve replacement and may occur following either surgical
or transcatheter aortic valve replacement. Although bioprosthetic heart valves
(BHVs) are remarkably less thrombogenic than mechanical heart valves (MHVs),
clinical valve thrombosis can occur as a life-threatening complication. Recent
studies [38, 39, 51] have suggested that the rate of subclinical leaflet
thrombosis (SLT) is as high as 13–38% [47]. SLT is associated with increased
risk of transient ischemic attacks and strokes, acute myocardial infarction,
and accelerated valve deterioration [45]. Further, if left untreated, SLT can
progress to clinical valve thrombosis. While a cardiac computed tomography
(CT) scan can detect SLT, predicting which patients will develop SLT is
currently not possible. Accordingly, there is a need for of computational
tools to model the fluid-structure and biochemical interactions that
predispose a particular patient to develop SLT.

Prior work to model leaflet thrombosis has focused on computational fluid
dynamics (CFD) simulations of blood flow through the valve. Plitman Mayo et
al. [44] performed CFD experiments of deployed transcatheter aortic valve
replacements (TAVRs) to determine areas of stagnated blood flow, suggesting
possible sites of thrombosis formation. Vahidkhah et al. [58] compared blood
residence times behind the coronary and non-coronary leaflets after a TAVR
procedure and determined similar residence times for all the leaflets. Kivi et
al. [29] performed two dimensional fluid-structure interaction (FSI)
simulations with leaflets of varying stiffness. A common finding in CFD and
FSI simulations is the presence of stagnant regions in the aortic sinus, in
which blood clots are thought to form. Hatoum et al. [23] combined a CFD model
of flow through patient specific geometry post-TAVR with a reduced order model
that predicted thrombus growth based on the wall shear stress and percent
stasis volume measurements. While they were able to determine a correlation
between circulation and amount of thrombosis, they concluded that finer flow
metrics or FSI analysis are needed to fully predict thrombosis.

Mathematical and computational models of thrombosis have also been developed,
but methods suitable for modeling thrombosis on dynamic flexible structures,
which is critical for describing leaflet thrombosis, are lacking. Fogelson et
al. [12, 36] developed a model of intravascular platelet deposition and
determined the sensitivity of thrombus formation due to various chemical and
platelet factors. Du and Fogelson [8] developed a multiphase model of platelet
aggregation in which the thrombus is modeled as a viscoelastic fluid. This
model can be seen as an extension of models by Fogelson and Guy [11] that were
created to study thrombus formation in a moving fluid. Models describing
flowing platelets and platelet deposition onto a stationary vessel wall have
been developed using a variety of multiscale modeling and computational
approaches [7, 59, 63, 55]. These models describe both fluid-phase transport
of platelets and the influence of platelet deposits on the hemodynamics
through and near the deposits. In these models, the platelets deposit over
stationary surfaces. However, to our knowledge, no thrombosis model has yet
been developed that allows for thrombus growth on a surface whose motion is
determined by solving an FSI problem, e.g., a heart valve leaflet.

There are several models that couple the advection and diffusion of chemical
species and their sources from immersed boundaries [50, 49, 25]. Typically,
these models use sources that are then spread from the immersed boundary to
the surrounding fluid using the regularized delta function. Restricting
species from diffusing across the interface remains a challenge. While many
different methods have been proposed to restrict diffusion and enforce Robin
boundary conditions across a moving interface [56, 27, 60, 54, 24], there are
far fewer that have tested the method in the context of an immersed boundary
model. Chen and Lai [6] used a diffuse domain approach to model absorption of
surfactants on a surface. Their approach is based on the methods introduced by
Li et al. [35] who demonstrated that this method enforces the boundary
condition at first order accuracy. In the methods used herein, we enforce a
boundary condition without smoothing the interface, leading to second order
accuracy up to and including the boundary [5].

The present study introduces new numerical methods to simulate the deposition
of material onto thin moving leaflets. The leaflet and fluid motion are
determined through an FSI calculation and the material deposition feeds back
onto the FSI calculation by modifying the leaflet’s mechanical properties.
While we refer to the fluid as blood and the deposited material as platelets,
the current paper deals only with a prototype of a situation that would arise
in modeling leaflet thrombosis. In a complete model of leaflet thrombosis, the
deposited material would consist of platelets, fibrin, and potentially other
inflammatory cells. The model would include mechanisms for activating
platelets through contact with the leaflet surface, exposure to high shear
stress, or encounter of soluble platelet agonists [14, 8, 9]. It would also
include treatment of coagulation biochemistry [37, 36, 34, 31] coupled with
fibrin polymerization [15, 13]. The current work is a major step towards
simulating the dynamics of such a model.

##  2 Continuous Equations

We consider an FSI model of thrombus formation on the aortic valve leaflets.
The valve geometry is created via slicing a three-dimensional reconstruction
of a typical trileaflet aortic valve, as will be discussed further in section
3.1. In this simplified model, fluid phase platelets can bind to the leaflet
surface while the surface-bound platelets stiffen the leaflets and can also
dissociate back into the fluid.

###  2.1 Fluid-Structure Interaction

The fluid-structure system is modeled using the immersed finite element/finite
difference method [21]. In this approach, a fixed computational domain
$\Omega$ is partitioned into a time-dependent fluid subdomain
$\Omega_{t}^{\text{f}}$ and a time-dependent solid subdomain
$\Omega_{t}^{\text{s}}$, so that
$\Omega=\Omega_{t}^{\text{f}}\cup\Omega_{t}^{\text{s}}$. The fluid domain is
further subdivided into the lumen $\Omega_{t}^{\text{f}^{-}}$ (i.e., the space
occupied by the blood, in which platelets are free to advect and diffuse) and
the space outside the aortic root $\Omega_{t}^{\text{f}^{+}}$, with
$\Omega_{t}^{\text{f}}=\Omega_{t}^{\text{f}^{-}}\cap\Omega_{t}^{\text{f}^{+}}$;
see Figure 1. We denote Eulerian physical coordinates with $\mathbf{x}$. The
solid domain is tracked using Lagrangian material coordinates $\mathbf{X}$,
and the mapping between the reference and current coordinates is
$\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$. The motion of the
fluid-structure system is described by

$\displaystyle\rho\mathopen{}\left(\frac{\partial\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}}{\partial
t}+\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}\cdot\nabla\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}\right)\mathclose{}$
$\displaystyle=-\nabla
p\mathopen{}\left(\mathbf{x},t\right)\mathclose{}+\mu\nabla^{2}\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}+\mathbf{f}\mathopen{}\left(\mathbf{x},t\right)\mathclose{},$
(1a)
$\displaystyle\nabla\cdot\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$
$\displaystyle=0,$ (1b)
$\displaystyle\mathbf{f}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$
$\displaystyle=\int_{\Omega_{0}^{\text{s}}}\mathbf{F}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}\,\delta\mathopen{}\left(\mathbf{x}-\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}\right)\mathclose{}\,d\mathbf{X},$
(1c) $\displaystyle\frac{\partial\bm{\chi}}{\partial
t}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$
$\displaystyle=\int_{\Omega}\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}\,\delta\mathopen{}\left(\mathbf{x}-\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}\right)\mathclose{}\,d\mathbf{x}=\mathbf{u}\mathopen{}\left(\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{},t\right)\mathclose{},$
(1d)

in which $\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$ and
$p\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$ are the Eulerian velocity
and pressure, respectively,
$\mathbf{f}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$ is the Eulerian
force density, $\mathbf{F}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$ is
the Lagrangian force density, which is determined in a manner specfied below,
and $\delta\mathopen{}\left(\mathbf{x}\right)\mathclose{}$ is the Dirac delta
function. The fluid density $\rho$ and viscosity $\mu$ are assumed to be
constant. Equations 1a and 1b are the well known Navier-Stokes equations and
hold across the entire computational domain $\Omega$. Equations 1c and 1d
couple the Lagrangian and Eulerian variables. The integral in equation 1c is
over the reference configuration of the solid subdomain while that in equation
1d is over the entire computational domain.

For the current study, the aortic walls are treated as approximately rigid
while the valve leaflets are elastic and deformable. For rigid structures, we
use a penalty formulation that is intended to tether the structure in place
using the force

$\mathbf{F}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}=\kappa\mathopen{}\left(\mathbf{X}-\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}\right)\mathclose{},$
(2)

in which $\kappa$ is the stiffness parameter [32]. In practice, we choose
$\kappa$ to be the largest stable value permitted by the numerical scheme so
that the structure’s motion is minimized.

Figure 1: The domain is decomposed into a solid subdomain
$\Omega_{t}^{\text{s}}$ (denoted in the black and tan curves) and a fluid
subdomain $\Omega_{t}^{\text{f}}$, which itself is partitioned into interior
subregion $\Omega_{t}^{\text{f}^{-}}$, which corresponds to the lumen, and an
exterior subregion $\Omega_{t}^{\text{f}^{+}}$, which corresponds to the space
outside the vessel. The lumen boundary
$\Gamma_{t}^{\text{f}^{-}}=\partial\Omega_{t}^{\text{f}^{-}}$ is composed of
regions where platelet binding can occur (shown in blue) and no penetration
conditions for the platelets (shown in red). The inlet to the domain is the
left ventricle outflow tract (LVOT) and the outlet is the ascending aorta.
Boundary conditions at the inlet are determined using a time-dependent
elastance based model of the heart, including the left ventricle (LV), mitral
valve (MV), and the left atrium (LA). The outlet boundary conditions are
determined using a three element Windkessel model [33, 22].

For elastic structures which in this study are the valve leaflets, the
response is given by the first Piola-Kirchhoff stress $\mathbb{P}$, which is
determined from a strain-energy function
$\Psi\mathopen{}\left(\mathbb{F}\right)\mathclose{}$ via
$\mathbb{P}=\frac{\partial\Psi}{\partial\mathbb{F}},$ in which
$\mathbb{F}=\frac{\partial\bm{\chi}}{\partial\mathbf{X}}$ is the deformation
gradient tensor. Following the immersed finte element/difference approach of
Vadala-Roth et al. [57], we split the strain energy functional into deviatoric
and dilational parts,
$\Psi\mathopen{}\left(\mathbb{F}\right)\mathclose{}=W\mathopen{}\left(\bar{\mathbb{F}}\right)\mathclose{}+U\mathopen{}\left(J\right)\mathclose{},$
in which $J=\text{det}~{}\mathbb{F}$ is the Jacobian of the deformation tensor
and $\bar{\mathbb{F}}=J^{-1/3}\mathbb{F}$. In what follows, we choose the
dilational part of the energy to be

$U\mathopen{}\left(J\right)\mathclose{}=\frac{\kappa_{\text{stab}}}{2}\mathopen{}\left(\log
J\right)\mathclose{}^{2},$ (3)

in which $\kappa_{\text{stab}}$ is the numerical bulk modulus. The Lagrangian
force density is then computed by requiring

$\int_{\Omega_{0}^{\text{s}}}\mathbf{F}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}\cdot\mathbf{V}\mathopen{}\left(\mathbf{X}\right)\mathclose{}\,d\mathbf{X}=-\int_{\Omega_{0}^{\text{s}}}\mathbb{P}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}:\nabla_{\mathbf{X}}\mathbf{V}\mathopen{}\left(\mathbf{X}\right)\mathclose{}\,d\mathbf{X},$
(4)

for all smooth test functions
$\mathbf{V}\mathopen{}\left(\mathbf{X}\right)\mathclose{}$.

The leaflets are modeled as a hyperelastic material that follows an
exponential neo-Hookean model [18, 41]. The deviatoric strain energy
functional for this model is given by

$W\mathopen{}\left(\bar{\mathbb{F}}\right)\mathclose{}=C_{10}\mathopen{}\left(e^{C_{01}\mathopen{}\left(\bar{I}_{1}-3\right)\mathclose{}}-1\right)\mathclose{},$
(5)

in which $C_{10}$ and $C_{01}$ are material constants and $\bar{I}_{1}$ is the
first deviatoric strain invariant of the modified right Cauchy-Green tensor,
$\bar{I}_{1}=\text{tr}\mathopen{}\left(\bar{\mathbb{C}}\right)\mathclose{}$,
which is defined in terms of the modified deformation tensor,
$\bar{\mathbb{C}}=\bar{\mathbb{F}}^{\text{T}}\bar{\mathbb{F}}$. The material
parameter $C_{10}$ is set to be a function of the bound platelet
concentration, as described in section 2.3.

###  2.2 Boundary Conditions

We use reduced order models to determine pressure-flow relationships in the
ascending aorta and left ventricular outflow tract (LVOT); see Figure 1. These
reduced order models are connected to the FSI model through boundary
conditions imposed on the fluid [33, 22]. We use a coupling scheme in which
the net flow rate through each of the boundary surfaces serves as an input to
the corresponding reduced order model. In turn, each reduced order model
determines a pressure that is prescribed on the corresponding boundary
surface. The net flow rate through the LVOT boundary is the integral of the
vertical component of the velocity over the portion of the bottom boundary of
the computation domain between the vesesl walls. The net flow rate for the
aortic boundary is defined in a similar way. On the remainder of the
computational domain’s boundary, we use zero velocity boundary conditions.

We use a three-element Windkessel model [53] for the downstream reduced order
model that models the aortic outflow,

$\displaystyle C\frac{dP_{\text{Wk}}}{dt}$
$\displaystyle=Q_{\text{Ao}}-\frac{P_{\text{Wk}}}{R_{\text{p}}},$ (6)
$\displaystyle P_{\text{Ao}}$
$\displaystyle=P_{\text{Wk}}+Q_{\text{Ao}}R_{\text{c}},$ (7)

in which $C$ is the compliance, $R_{\text{c}}$ is the characteristic
resistance, $R_{\text{p}}$ is the peripheral resistance, $P_{\text{Wk}}$ is
the Windkessel pressure, $Q_{\text{Ao}}$ is the computed volumetric flow rate
at the outlet of the ascending aorta, and $P_{\text{Ao}}$ is pressure at the
outlet of the ascending aorta which is then prescribed as a boundary condition
for the fluid.

For the upstream model to model the inflow from the heart, we employ a time-
dependent elastance-based left heart model [53],

$\displaystyle\frac{d\mathopen{}\left(C_{\text{LA}}P_{\text{LA}}\right)\mathclose{}}{dt}$
$\displaystyle=Q_{\text{vein}}-Q_{\text{MV}},$ (8)
$\displaystyle\frac{d\mathopen{}\left(C_{\text{LV}}P_{\text{LV}}\right)\mathclose{}}{dt}$
$\displaystyle=Q_{\text{MV}}-Q_{\text{LVOT}},$ (9) $\displaystyle
P_{\text{LVOT}}$ $\displaystyle=P_{\text{LV}}-Q_{\text{LVOT}}R_{\text{LVOT}},$
(10) $\displaystyle Q_{\text{MV}}$
$\displaystyle=\left\\{\begin{array}[]{ll}0,&\mbox{ if }\quad
P_{\text{LA}}\leq P_{\text{LV}},\\\
\frac{P_{\text{LA}}-P_{\text{LV}}}{R_{\text{MV}}},&\mbox{ if }\quad
P_{\text{LA}}>P_{\text{LV}},\end{array}\right.$ (13)

in which $C_{\text{LA}}$ and $C_{\text{LV}}$ are the time-dependent
compliances of the left atrium and left ventricle, respectively.
$R_{\text{LVOT}}$ and $R_{\text{MV}}$ are the resistances of the LVOT and
mitral valve, the latter of which is modeled as a diode. $P_{\text{LA}}$,
$P_{\text{LV}}$, and $P_{\text{LVOT}}$ are the left atrial, left ventricular,
and the LVOT pressures, and $Q_{\text{vein}}$, $Q_{\text{MV}}$, and
$Q_{\text{LVOT}}$ are the volumetric flow rates of the pulmonary veins, mitral
valve, and LVOT. In this model, $Q_{\text{vein}}$ is prescribed as a constant
inflow rate into the left atrium. $Q_{\text{LVOT}}$ is the computed flow rate
at the inlet of the computational domain. $P_{\text{LVOT}}$ is then prescribed
as a boundary condition for the momentum equation 1a. We determine the time-
dependent compliances $C(t)$ from specified elastance functions $E(t)$ via
$C(t)=1/E(t)$. We use the “two-Hill” elastance waveform given by Mynard et
al.[43],

$\displaystyle
E\mathopen{}\left(t\right)\mathclose{}=\mathopen{}\left(E_{\text{max}}-E_{\text{min}}\right)\mathclose{}\alpha\mathopen{}\left(t\right)\mathclose{}+E_{\text{min}},$
(14)
$\displaystyle\alpha\mathopen{}\left(t\right)\mathclose{}=\frac{k\frac{g_{1}}{1+g_{1}}\frac{1}{1+g_{2}}}{\max\mathopen{}\left(\frac{g_{1}}{1+g_{1}},\frac{1}{1+g_{2}}\right)\mathclose{}},$
(15) $\displaystyle g_{i}=\left(\frac{t}{\tau_{i}}\right)^{m_{i}}.$ (16)

We use the elastance parameters in $E\mathopen{}\left(t\right)\mathclose{}$
for the left atrium from Mynard et al. [43]. The remaining parameters are fit
to experimental measurements of human aortic pressures $P_{\text{Ao}}$ and
aortic flow rates $Q_{\text{Ao}}$ from Murgo et al. [42] by taking the
experimental measurements of $Q_{\text{Ao}}$ as input to the Windkessel model,
and comparing the resulting values of $P_{\text{Ao}}$ to its experimental
values. We calculate the best-fit parameters to data from Murgo et al. [42]
for a “Type A” beat for the upstream model. The downstream model is fit using
the corresponding downstream data from Murgo et al. [42]. The fits were
created using MATLAB’s `fmincon`, a nonlinear optimization tool.

###  2.3 Mass Deposition Model

We couple the FSI model to a mass deposition model that includes a fluid-phase
$c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$ concentration
measured per unit volume and a surface-bound
$C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$ concentration
field measured per unit reference area. Although this model does not include
the cellular and biochemical interactions describing thrombosis, it does
include fields which we view as platelet populations, and the conversion of
fluid-phase platelets in $c_{\text{f}}$ to surface-bound platelets in
$C_{\text{b}}$ as platelet adhesion. The fluid-phase species diffuses and
advects with the local fluid velocity in the interior fluid domain
$\Omega_{t}^{\text{f}^{-}}$ and can be converted into the surface-bound
species along the boundary
$\Gamma_{t}\subset\Gamma_{t}^{\text{f}^{-}}=\partial\Omega_{t}^{\text{f}^{-}}$.
In the results in section 4, $\Gamma_{t}$ is the downstream side of one or
both of the leaflets. The surface-bound species moves with the structure and
can dissociate to become the fluid-phase species. The model equations are

$\displaystyle\frac{\partial
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}}{\partial
t}+\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}\cdot\nabla
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}=$ $\displaystyle
D\nabla^{2}c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{},$
$\displaystyle\mathbf{x}\in\Omega_{t}^{\text{f}^{-}},$ (17a)
$\displaystyle\frac{\partial
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}}{\partial\mathbf{n}}=$
$\displaystyle 0,$
$\displaystyle\mathbf{x}\in\Gamma_{t}^{\text{f}^{-}}\setminus\Gamma_{t},$
(17b) $\displaystyle-D\frac{\partial
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}}{\partial\mathbf{n}}=$
$\displaystyle
k_{\text{on}}\mathopen{}\left(C_{\text{b}}^{\text{max}}-C_{\text{b}}\mathopen{}\left(\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{},t\right)\mathclose{}\right)\mathclose{}J_{\text{s}}c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$
$\displaystyle-
k_{\text{off}}C_{\text{b}}\mathopen{}\left(\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{},t\right)\mathclose{}J_{\text{s}},$
$\displaystyle\mathbf{x}\in\Gamma_{t},$ (17c) $\displaystyle\frac{\partial
C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}}{\partial t}=$
$\displaystyle
k_{\text{on}}\mathopen{}\left(C_{\text{b}}^{\text{max}}-C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}\right)\mathclose{}c_{\text{f}}\mathopen{}\left(\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{},t\right)\mathclose{}$
$\displaystyle-
k_{\text{off}}C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{},$
$\displaystyle\mathbf{X}\in\Gamma_{0},$ (17d)

in which $D$ is the diffusion coefficient, $k_{\text{on}}$ and
$k_{\text{off}}$ are the reaction rates for adhesion and dissociation,
respectively, $C_{\text{b}}^{\text{max}}$ is the carrying capacity of
$C_{\text{b}}$ per unit undeformed area along the boundary $\Gamma_{0}$, and
$J_{\text{s}}=\frac{dA}{da}$ is the surface Jacobian, which is the ratio of
reference and current areas. The first term on the right hand side of equation
17c gives the rate of binding of fluid-phase platelets with concentration
$c_{\text{f}}$ to the valve leaflet where the available binding sites have
surface density
$\mathopen{}\left(C_{\text{b}}^{\text{max}}-C_{\text{b}}\mathopen{}\left(\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{},t\right)\mathclose{}\right)\mathclose{}J_{\text{s}}$
with respect to the current leaflet configuration. The second term gives the
rate at which absorbed platelets with surface density
$C_{\text{b}}\mathopen{}\left(\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{},t\right)\mathclose{}J_{\text{s}}$
detach from the leaflet.

To model the effect of thrombosis over the valve leaflets, we set the
stiffness coefficient of the leaflets $C_{10}$ to be a function of the surface
concentration $C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$.
Because $C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$ is
defined only on the surface of the leaflet, we use a harmonic interpolation
procedure to extend the surface concentration into the interior of the
leaflet, where the Lagrangian forces are calculated. Specifically, we solve
Laplace’s equation

$\displaystyle\nabla^{2}C_{\text{b}}^{\text{in}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$
$\displaystyle=0,\quad\mathbf{X}\in\Omega^{\text{leaf}}_{0},$ (18)
$\displaystyle
C_{\text{b}}^{\text{in}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$
$\displaystyle=\left\\{\begin{array}[]{ll}C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{},&\mbox{
if }\mathbf{X}\in\Gamma_{0},\\\ 0,&\mbox{ otherwise},\end{array}\right.$ (21)

in which $\Omega^{\text{leaf}}_{0}$ is the leaflet domain in the initial
configuration. Having found
$C_{\text{b}}^{\text{in}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$, we
then set the stiffness of the leaflet to be

$C_{10}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}=C_{10}^{\text{base}}\mathopen{}\left(\frac{\beta+1}{2}-\frac{\beta-1}{2}\cos\mathopen{}\left(\frac{\pi
C_{\text{b}}^{\text{in}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}}{C_{\text{b}}^{\text{max}}}\right)\mathclose{}\right)\mathclose{},$
(22)

in which $C_{10}^{\text{base}}$ is the stiffness with no accumulation and
$\beta C_{10}^{\text{base}}$ is the maximum stiffness.

The parameters of the mass deposition model are chosen so that the reactions
occur on a similar time scale as the fluid-structure interactions. These
values are several orders of magnitude larger than those used in a similar
thrombosis model as described previously [34, 11]. Use of physiologically
relevant reaction rates would require performing simulations over thousands of
computational cycles, which is currently not feasible. We are actively working
on a temporal multiscale method to meet this challenge and allow use of
realistic reaction rates.

##  3 Computational Models and Numerical Methods

The model is implemented in IBAMR, which provides implementations of the
immersed boundary method and several of its extensions along with support for
adaptive mesh refinement [19]. IBAMR utilizes libMesh for the finite element
representation of the structural deformations [28] and PETSc for linear
solvers [2, 3, 4]. Support for structured adaptive mesh refinement is provided
by SAMRAI [26]. While the model can be naturally extended to three spatial
dimensions, we describe the numerical implementation and results in two
spatial dimensions.

###  3.1 Imaged Model and Mesh Generation

Our two-dimensional aortic root geometry is informed by a three-dimensional
patient-specific aortic root model based on pre-procedural computed tomography
(CT) image data of a female patient preparing for aortic valve replacement at
UNC Medical Center. The images used in this study were obtained under a
protocol approved by the UNC Institutional Review Board (study number
18-0202). The CT scan was performed using a Siemens SOMATOM Definition CT
Scanner with an image resolution of $512\times 512\times 226$ and a voxel size
of
$$0.441\text{\,}\frac{\mathrm{mm}}{}$\times$0.441\text{\,}\frac{\mathrm{mm}}{}$\times$0.6\text{\,}\frac{\mathrm{mm}}{}$$.
The patient images are segmented by a semi-automated method in ITK-SNAP [62],
which implements an active contour model that minimizes an energy functional
of voxel intensities [61]. The aortic root measures
$28\text{\,}\frac{\mathrm{mm}}{}$ in diameter and
$7.68\text{\,}\frac{\mathrm{cm}}{}$ in length, and the thickness of the aortic
wall is $1.0\text{\,}\frac{\mathrm{mm}}{}$. The inflow boundary of the model
is truncated at the LVOT, and the outflow boundary of the model is truncated
downstream of the aortic valve before the first arterial bifurcation.
Artificial circular extensions are added at both boundaries using SOLIDWORKS
(Dassault Systèmes SOLIDWORKS Corporation, Waltham, MA, USA) to simplify the
application of boundary conditions to the computational model. The radius of
the vessel at both truncations is $21\text{\,}\frac{\mathrm{mm}}{}$. Idealized
aortic valve replacement leaflets with a thickness of
$0.7\text{\,}\frac{\mathrm{mm}}{}$ are created based on the measurements from
Sahasakul et al. [48] and trimmed to fit within the reconstructed aortic root
in SOLIDWORKS. To derive our two-dimensional aortic root geometry from the
three-dimensional model, we extract a slice through the diameter of the aorta
using Coreform Cubit (Computational Simulation Software, LLC, American Fork,
UT, USA), which is a software application based on CUBIT from Sandia National
Laboratory. We then use Cubit to smooth the angles in both the aortic root and
leaflet surfaces and to generate structural meshes consisting of triangular
elements.

###  3.2 Fluid-Structure Interaction

The fluid equations (1a) and (1b) are solved using a second-order Cartesian
staggered-grid finite difference method. The nonlinear term is approximated
using a piecewise parabolic method [46]. The resulting saddle point system is
solved using GMRES with a projection method as a preconditioner [20].

The solid subdomain $\Omega_{t}^{\text{s}}$ is discretized using
$\mathcal{C}^{0}$ finite elements. A triangulation $\mathcal{T}_{h}$ of the
structure is constructed. The size of each element in the triangulation is
chosen so that there is approximately one node per Cartesian grid cell. On
$\mathcal{T}_{h}$, we define Lagrangian basis functions
$\left\\{\phi_{l}\mathopen{}\left(\mathbf{X}\right)\mathclose{}\right\\}_{l=1}^{m}$,
in which $m$ is the total number of nodes in the triangulation. We approximate
the structural deformation and force using the basis functions via

$\displaystyle\bm{\chi}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}=\sum_{l=1}^{m}\bm{\chi}_{l}\mathopen{}\left(t\right)\mathclose{}\,\phi_{l}\mathopen{}\left(\mathbf{X}\right)\mathclose{},$
(23)
$\displaystyle\mathbf{F}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}=\sum_{l=1}^{m}\mathbf{F}_{l}\mathopen{}\left(t\right)\mathclose{}\,\phi_{l}\mathopen{}\left(\mathbf{X}\right)\mathclose{}.$
(24)

Coupling between the fluid and structure is mediated using regularized delta
functions in equations 1c and 1d. Recently, Lee and Griffith [32] suggested
using delta functions with smaller support for structures in shear driven
regimes. Therefore, in this work, we use the three-point $B$-spline kernel for
the flexible valve leaflets, and a two-point piecewise linear kernel for the
nearly rigid walls of the aortic root.

###  3.3 Mass Deposition Model

The fluid phase concentration field is approximated using a hybrid semi-
Lagrangian cut-cell method [5]. For brevity, we omit the details and only
highlight the changes of the discretization. To summarize, we split equation
17 into an advection step,

$\frac{\partial
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}}{\partial
t}+\mathbf{u}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}\cdot\nabla
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}=0,\quad\mathbf{x}\in\Omega_{t}^{\text{f}^{-}},$
(25)

and a diffusion step,

$\frac{\partial
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}}{\partial
t}=D\nabla^{2}c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{},\quad\mathbf{x}\in\Omega_{t}^{f^{-}},$
(26)

along with the boundary conditions equations 17b and 17c and the surface
concentration equation 17d. During the diffusion step, the domain
$\Omega_{t}^{\text{f}^{-}}$ is assumed to be fixed. The advective step is
treated with a semi-Lagrangian method using polyharmonic splines to
reconstruct the function $c_{\text{f}}$. The diffusion step is treated with a
cut-cell finite volume method. The surface concentration
$C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$ is solved for
by extrapolating the fluid-phase field
$c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$ to the boundary
and approximating the ODE in equation 17d.

####  3.3.1 Diffusion

To approximate the diffusion step in equation 26, we employ a cut-cell finite
volume method in which the domain $\Omega_{t}^{f^{-}}$ is considered fixed for
the duration of this step. Integrating equation 26 over a grid cell
$\mathbf{c}_{i,j}$ that is entirely or partially interior to
$\Omega_{t}^{f^{-}}$ and dividing by that cell’s volume, we obtain

$\frac{1}{\left|\mathbf{c}_{i,j}\cap\Omega_{t}^{f^{-}}\right|}\int_{\mathbf{c}_{i,j}\cap\Omega_{t}^{f^{-}}}\frac{\partial
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}}{\partial
t}\text{d}\mathbf{x}=\frac{1}{\left|\mathbf{c}_{i,j}\cap\Omega_{t}^{f^{-}}\right|}\int_{\mathbf{c}_{i,j}\cap\Omega_{t}^{f^{-}}}D\Delta
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}\text{d}\mathbf{x}.$
(27)

We define $C_{\text{f},i,j}$ as the cell average of
$c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$ in the cell
$\mathbf{c}_{i,j}\cap\Omega_{t}^{f^{-}}$. Replacing the cell average in the
left hand side of equation 27 and employing the divergence theorem on the
right hand side, we obtain

$\frac{\text{d}C_{\text{f},i,j}}{\text{d}t}=\frac{1}{\left|\mathbf{c}_{i,j}\cap\Omega_{t}^{f^{-}}\right|}\int_{\partial\mathopen{}\left(\mathbf{c}_{i,j}\cap\Omega_{t}^{f^{-}}\right)\mathclose{}}D\frac{\partial
c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}}{\partial\mathbf{n}}\cdot\text{d}\mathbf{A}.$
(28)

The integral in equation 28 consists of two parts, an integral over the
boundary $\Gamma_{t}^{\text{f}^{-}}$ that is interior to cell
$\mathbf{c}_{i,j}$ and an integral over the portion of the boundary of the
cell $\mathbf{c}_{i,j}$ that is interior to $\Omega_{t}^{\text{f}^{-}}$. The
first type consists of an integral over the physical boundary and using the
provided boundary conditions in equations 17b and 17c, can be computed using
techniques described in the next section. The second integral is discretized
using second order finite differences. This discretization requires the
computation of the cut cell volume
$\left|\mathbf{c}_{i,j}\cap\Omega_{t}^{f^{-}}\right|$, which is described in
section 3.3.4.

####  3.3.2 Surface Reactions

Along part of the surface $\Gamma_{t}^{\text{f}^{-}}$, we allow for binding of
the fluid-phase species to the boundary and for unbinding of the surface-bound
species into the fluid, as described by equations 17c and 17d. We extract a
boundary mesh represented by $C^{0}$ elements from the volumetric leaflet mesh
as described in section 3.2. We maintain a representation of both the surface
concentration $C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$
per unit reference area and the fluid concentration
$c_{\text{f}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}$ per unit volume
restricted to the boundary. These values are represented using Lagrangian
basis functions
$\left\\{\psi_{l}\mathopen{}\left(\mathbf{X}\right)\mathclose{}\right\\}_{l=1}^{n_{\text{bd}}}$
in which $n_{\text{bd}}$ is the number of nodes of the boundary mesh. We note
these are the same basis functions used for the structural deformation, but
restricted to the surface. The concentrations along the boundary are
accordingly

$\displaystyle
C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}=\sum_{l=1}^{n_{\text{bd}}}C_{\text{b}}^{l}\mathopen{}\left(t\right)\mathclose{}\psi_{l}\mathopen{}\left(\mathbf{X}\right)\mathclose{},$
(29) $\displaystyle
c_{\text{f}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}=\sum_{l=1}^{n_{\text{bd}}}c_{\text{f}}^{l}\mathopen{}\left(t\right)\mathclose{}\psi_{l}\mathopen{}\left(\mathbf{X}\right)\mathclose{}.$
(30)

The values $c_{\text{f}}^{l}$ are found by using a radial basis function
interpolant as described in section 3.3.3 to extrapolate the values of
$c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$ to the surface
nodes. The nodal values $C_{\text{b}}^{l}$ are found by solving the ODE in
equation 17d using a two stage Runge Kutta method.

This finite element representation allows for easy evaluations of the flux
defined in equation 17c from the boundary to the fluid. To evaluate this flux,
we require the value of the Jacobian
$J_{\text{s}}=\frac{d\mathbf{a}}{d\mathbf{A}}$ that converts areas in the
reference configuration to areas in the current configuration. Because we are
using a $C^{0}$ representation of the surface, the Jacobian is discontinuous
at nodes. To obtain a continuous representation, we project $J_{\text{s}}$
onto the finite element basis [30]. In practice, this amounts to computing the
Jacobian at quadrature points along the surface.

####  3.3.3 Reconstructions

Both the semi-Lagrangian step and the surface reactions involve reconstructing
$c_{\text{f}}\mathopen{}\left(\mathbf{x},t\right)\mathclose{}$ at various
points $\hat{\mathbf{x}}$ in the computational domain. The details of the
reconstruction procedure depend on where the reconstruction is being performed
within this domain. Away from the boundary, we use the four closest grid
points to $\hat{\mathbf{x}}$ to form a bilinear interpolant. If there are too
few points to form the bilinear interpolant (e.g., near cut-cells), we use a
radial basis function (RBF) interpolant [10, 52]. The RBF interpolant is
constructed via a polyharmonic spline

$q\mathopen{}\left(\mathbf{x}\right)\mathclose{}=\sum_{j=1}^{k}\lambda_{j}\lVert\mathbf{x}-\mathbf{x}_{j}\rVert^{m}+\sum_{j=1}^{s}\beta_{j}p_{j}\mathopen{}\left(\mathbf{x}\right)\mathclose{},$
(31)

in which $m$ is an odd integer and
$p_{j}\mathopen{}\left(\mathbf{x}\right)\mathclose{}$ form a set of $s$
polynomial basis functions. The total number of points in the stencil $k$ is
chosen so that $k=2m+1$. The points $\mathbf{x}_{j}$ are the $k$ closest
points to the location $\hat{\mathbf{x}}$. We find the coefficients
$\lambda_{j}$ and $\beta_{j}$ by requiring

$\displaystyle q\mathopen{}\left(\mathbf{x}_{j}\right)\mathclose{}$
$\displaystyle=f_{j}$ $\displaystyle\mbox{ for }j=1,\ldots,k,$ (32a)
$\displaystyle\sum_{i=1}^{s}\lambda_{i}p_{i}\mathopen{}\left(\mathbf{x}_{j}\right)\mathclose{}$
$\displaystyle=0$ $\displaystyle\mbox{ for }j=1,\ldots,k.$ (32b)

Equation 32a are the interpolation conditions, and equation 32b are the
orthogonality conditions to ensure a unique interpolant. This results in a
linear system for the coefficients, which is solved using a QR algorithm. In
our computations, we set the integer $m=3$ and use up to quadratic
polynomials.

####  3.3.4 Cut Cell Geometries

In the cut-cell finite volume discretization of equation 26, we require the
computation of the geometries of cells cut by the boundary $\Gamma^{f^{-}}$.
We denote the node of a Carteisan grid cell by
$\mathbf{x}_{i+\frac{1}{2},j+\frac{1}{2}}$. To find cut cell volumes, we first
calculate the signed distance function to the surface at each node
$\mathbf{x}_{i+\frac{1}{2},j+\frac{1}{2}}$. To do this, we first find
intersections of the $C^{0}$ representation of the immersed structure with the
background Eulerian grid, and for each element of the immersed structure, we
calculate outward facing normals. We note that this requires a consistent
traversal (e.g., counter-clockwise) of the structure to ensure a consistent
facing normal. Then, for each cell node
$\mathbf{x}_{i+\frac{1}{2},j+\frac{1}{2}}$ of the background grid, we find the
projection of the point onto each element and compute its distance from
$\mathbf{x}_{i+\frac{1}{2},j+\frac{1}{2}}$. If multiple minimal distance
projections exist, we use the angle weighted average of the projections [1].
The sign of the distance is computed using the previously computed structure
normal. Once we have the signed distances at each cell node, we can compute
partial cell volumes. Following Min and Gibou [40], we compute cell volumes by
decomposing the cell into simplices, for which analytic formulas for the
volume exist.

###  3.4 Time Stepping

In summary, the steps to advance the solution from time $t^{n}$ to time
$t^{n+1}$ are:

  1. 1.

Compute the cut cell geometries.

  2. 2.

Perform a half step of the diffusion solve, evolving both the fluid-phase and
the structure-bound concentration fields.

  3. 3.

Solve the Navier-Stokes equations and update the position of the immersed
structure.

  4. 4.

Update the cut cell geometries using the new position of the immersed
structure.

  5. 5.

Perform a full step of the semi-Lagrangian method, using the velocities from
the Navier-Stokes solve.

  6. 6.

Perform a half step of the diffusion step, evolving the fluid-phase and
surface-bound concentrations.

Our use of an explicit time stepping scheme for several of these steps limits
our time step size to resolve the fastest time scale. In this case, the
fastest time scale is that of the leaflet elasticity. We determine an
empirical scaling relationship between the time step size and the stiffness of
the leaflet that maintains numerical stability under increasing leaflet
stiffness. Specifically, we choose the time step such that

$\Delta
t=\frac{C_{\text{ts}}}{\sqrt{\text{max}\mathopen{}\left(C_{10}\right)\mathclose{}}},$
(33)

in which $C_{\text{ts}}$ is chosen to be as large as possible.

##  4 Results

Table 1 provides the values of all relevant physical and numerical parameters.
At $t=0$, the initial fluid phase concentration $c_{\text{f}}$ is set to be 1
throughout the domain $\Omega_{t}^{\text{f}^{-}}$. During the first two
cycles, the binding and unbinding coefficient, $k_{\text{on}}$ and
$k_{\text{off}}$, are set to zero; afterward they are reset to their non-zero
values. We emphasize that the binding and unbinding coefficients and the
diffusion coefficient are artificially increased by several orders of
magnitude compared to other clotting models [34, 11] to ensure that sufficient
binding can occur within the duration of the simulation.

Table 1: Values of the parameters used in the simulation. Structure parameters | Deposition and Fluid parameters
---|---
$\kappa$ | $20.17\text{\,}\frac{\mathrm{GPa}}{{\mathrm{cm}}^{2}}$ | $k_{\text{on}}$ | $0.033\,21\text{\,}\frac{{\mathrm{cm}}^{3}}{\mathrm{s}\text{\,}\mathrm{platelet}}$
$k_{\text{stab}}$ | $58.4\text{\,}\frac{\mathrm{MPa}}{}$ | $k_{\text{off}}$ | $0.01\text{\,}\frac{}{\mathrm{s}}$
$C_{01}$ | $3.25\text{\,}\frac{}{}$ | $D$ | $0.1\text{\,}\frac{{\mathrm{cm}}^{2}}{\mathrm{s}}$
$C_{10}^{\text{min}}$ | $2.264\text{\,}\frac{\mathrm{MPa}}{}$ | $C_{\text{b}}^{\text{max}}$ | $1.41\text{\times}{10}^{7}\text{\,}\frac{\mathrm{platelet}}{{\mathrm{cm}}^{2}}$
$\beta$ |  varies, between $1-600$ | $c_{\text{f}}^{\text{max}}$ | $1.5\text{\times}{10}^{5}\text{\,}\frac{\mathrm{platelet}}{{\mathrm{cm}}^{3}}$
|  | $\rho$ | $1\text{\,}\frac{\mathrm{g}}{{\mathrm{cm}}^{3}}$
|  | $\mu$ | $0.035\text{\,}\frac{\mathrm{g}}{\mathrm{cm}\text{\,}\mathrm{s}}$
Boundary model parameters
$R_{\text{p}}$ | $0.9046\text{\,}\frac{\mathrm{mmHg}\text{\,}\mathrm{s}}{\mathrm{mL}}$ | $\tau_{1,\text{LA}}$ | $0.097\,89\text{\,}\frac{\mathrm{s}}{}$
$C$ | $1.950\text{\,}\frac{\mathrm{mL}}{\mathrm{mmHg}}$ | $\tau_{1,\text{LV}}$ | $0.0887\text{\,}\frac{\mathrm{s}}{}$
$R_{\text{c}}$ | $0.042\text{\,}\frac{\mathrm{mmHg}\text{\,}\mathrm{s}}{\mathrm{mL}}$ | $m_{1,\text{LA}}$ | $1.32\text{\,}\frac{}{}$
$Q_{\text{vein}}$ | $6.2\text{\,}\frac{\mathrm{L}}{\mathrm{min}}$ | $m_{1,\text{LV}}$ | $2.404\text{\,}\frac{}{}$
$Q_{\text{LVOT}}$ | $0.015\text{\,}\frac{\mathrm{mmHg}\text{\,}\mathrm{s}}{\mathrm{mL}}$ | $\tau_{2,\text{LA}}$ | $0.1602\text{\,}\frac{\mathrm{s}}{}$
$R_{\text{MV}}$ | $0.005\text{\,}\frac{\mathrm{mmHg}\text{\,}\mathrm{s}}{\mathrm{mL}}$ | $\tau_{2,\text{LV}}$ | $0.4461\text{\,}\frac{\mathrm{s}}{}$
$E_{\text{max,\text{LA}}}$ | $0.17\text{\,}\frac{\mathrm{mmHg}}{\mathrm{mL}}$ | $m_{2,\text{LA}}$ | $13.1\text{\,}\frac{}{}$
$E_{\text{min,\text{LA}}}$ | $0.08\text{\,}\frac{\mathrm{mmHg}}{\mathrm{mL}}$ | $m_{2,\text{LV}}$ | $20.952\text{\,}\frac{}{}$
$E_{\text{min,\text{LV}}}$ | $0.0265\text{\,}\frac{\mathrm{mmHg}}{\mathrm{mL}}$ | $E_{\text{max,\text{LV}}}$ | $0.16\text{\,}\frac{\mathrm{mmHg}}{\mathrm{mL}}$

###  4.1 Convergence Study

The flow regime during peak systole is turbulent, with the largest Reynolds
number being approximately 5000. Because of the chaotic nature of the
simulation, convergence of the numerical method is not well defined. Small
changes in the simulations (e.g. grid spacing and time step size) can lead to
large changes in the flow velocities. Further, the fluid-phase and surface
concentrations and hence the stiffness of the leaflets is directly affected by
the turbulent flow. Therefore, to assess the accuracy of the simulation, we
compare the average fluid velocity near peak systole across grid sizes. We
modify the model to use a parabolic velocity profile which corresponds to
three-quarters systole. We generate three different resolutions of meshes with
maximum element edge lengths of 0.32, 0.24, and
$0.18\text{\,}\frac{\mathrm{mm}}{}$, which correspond to 2, 3, and 4 elements
across the width of the leaflets, respectively. The background Cartesian grid
is refined such that there is approximately one structural mesh node per grid
cell. The modified model is then run without accumulation. Figure 2 shows the
average fluid velocity from time $t=0.5$ to a final time $T=3.5$. We observe
consistent values across all grid resolutions tested. While convergence is not
clear, we expect the average flow velocity to show convergence in the limit as
$T\rightarrow\infty$. In the full model, we do observe grid independence of
the surface concentration field. Figure 3 shows the total bound concentration
$\int_{\Gamma_{0}}C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}d\mathbf{X}$
for all three grid resolutions. For the results presented below, we use the
coarse mesh, consisting of two elements across the leaflet.

(a) |  (b) |  (c)
---|---|---
(d) |  (e)
Figure 2: To assess convergence, we perform simulations at approximately
three-quarters systole, and we compute the average velocity from time $0.5$ to
time $T=3.5$. Panel (a) shows the average fluid velocity magnitude across the
time interval. Panels (b)-(e) show slices of the average fluid velocity
magnitude for three different grid resolutions. The coarsest grid is shown in
light blue, the medium grid is shown in blue, and the finest grid is shown in
black. We observe consistent values across all grid resolutions tested. Figure
3: The total accumulation
$\int_{\Gamma_{0}}C_{\text{b}}\mathopen{}\left(\mathbf{X},t\right)\mathclose{}d\mathbf{X}$
over time for three different mesh sizes. Although we do not expect pointwise
convergence of the fluids phase concentration and velocity fields for the
turbulent flow regime considered in this study, we do observe grid
independence of the bound concentration field.

###  4.2 Leaflet Deposition

Figure 4 shows fluid-phase concentrations $c_{\text{f}}$ and velocity
magnitude snapshots from the last cycle of a simulation with deposition only
on the right leaflet. At the times of these plots, the right leaflet is
substantially stiffer than the left, and the predominant flow through the
valve is shifted towards the left leaflet. Figure 5 shows fluid-phase
concentrations $c_{\text{f}}$ and velocity magnitude snapshots with deposition
on both leaflets. Here, both leaflets become stiff and open less over time. We
observe higher velocity magnitudes when deposition occurs on both leaflets as
opposed to a single leaflet.

---
(a) |  (b) |  (c) |  (d)
|  |  |

Figure 4: For a simulation in which deposition happens only on the right
leaflet, snapshots during (a) middle of diastole, (b) peak systole, (c), end
systole, and (d) middle of diastole the last cycle showing (top) the fluid-
phase concentration $c_{\text{f}}$ and the surface concentration
$C_{\text{b}}J_{\text{s}}$ and (bottom) the fluid velocity magnitude. This
simulation uses $\beta=600$. Notice that the right leaflet is considerably
stiffer than the left one and that the fluid concentration $c_{\text{f}}$ is
depleted in and downstream of the aortic sinus.

---
(a) |  (b) |  (c) |  (d)
|  |  |

Figure 5: For a simulation in which deposition happens on both leaflets,
snapshots during (a) middle of diastole, (b) peak systole, (c), end systole,
and (d) middle of diastole in the last cycle showing (top) the fluid-phase
concentration $c_{\text{f}}$ and the surface concentration
$C_{\text{b}}J_{\text{s}}$ and (bottom) the fluid velocity magnitude. This
simulation uses $\beta=600$. Notice that both leaflets stiffen in this
simulation and that the peak velocity magnitudes are larger than those in
Figure 4.

To assess the opening of the valve, we project the leaflets onto the valve
ring, as shown in Figure 6. The opening area is then normalized by the area of
the fully open valve. Figure 6 shows the normalized open valve area over each
cycle for deposition on both leaflets, as we increase the maximum stiffness
factor $\beta$. For lower maximal stiffnesses, we observe similar normalized
open valve areas compared to a simulation with no deposition. For larger
maximum stiffness, we observe a smaller normalized open valve area as more
deposition occurs. Figure 6 compares the normalized open valve area for
deposition on both leaflets versus on only the right leaflet for the same
maximum stiffness. When deposition occurs on only the right leaflet, the
normalized open valve area still decreases compared to no accumulation, but
does not realize as dramatic reductions as when deposition on both leaflets is
allowed. The left leaflet, which has a constant stiffness over time,
compensates and opens more as the right leaflet stiffens.

Figure 7 shows the maximum and minimum accumulations
$C_{\text{b}}J_{\text{s}}$ across the leaflet. Because the diffusion
coefficient is large compared to the reaction rates, there is always
sufficient fluid-phase platelets to bind to the leaflets, and accordingly, a
consistent rate of binding to the surface and an increasing surface
concentration as the simulation progresses. By the end of the simulation,
bound platelets occupy approximately 23% of the carrying capacity. The minimum
concentration periodically jumps while the maximum concentration is
monotonically increasing. The periodic jumps are due to the physical location
of the minimum and maximum which affects the amount of platelets per current
area. The minimum occurs near the position where the leaflet attaches to the
aortic wall. This location sees the largest changes in area as the valve
opens. The maximum concentration is found on the tips of the leaflet, which
move through regions of high fluid-phase concentration. The tip of the
leaflets deform less than the rest of the leaflet upon opening and closing of
the valve, leading to a steadily increasing surface bound concentration field.
While the fluid-phase concentration $c_{\text{f}}$ is not completely depleted,
we do observe reductions in the concentration than initially.

Figure 8 shows the velocity magnitude during the final cycle near peak
systole. We observe a vortex in the sinus region that grows in strength as we
increase the maximum stiffness. This vortex is not present when accumulation
occurs exclusively on the right leaflet, but is present when there is no
accumulation.

(a) (b) (c)

Figure 6: The computed normalized open valve area over time. Panel (a) depicts
the computation of the normalized open valve area. The open area of the valve
is projected onto the valve ring, and then normalized by the area of the fully
opened valve. Panel (b) depicts the normalized open valve area over time for
accumulation on both leaflets as $\beta$ increases during the entire
simulation (top) and during only the last cycle (bottom). Notice that as the
leaflets get stiffer, the normalized open valve area decreases. Panel (c)
depicts the normalized open valve area over time for accumulation on the right
leaflet or both leaflets during all the cycles (top) and during only the last
cycle (bottom). Notice that if accumulation occurs on both leaflets, the
normalized open valve area decreases more than if accumulation occurs on a
single leaflet.

(a)

(c)

(b)

(d)

Figure 7: (a) The surface concentration $C_{\text{b}}J_{\text{s}}$ along the
leafets at the end of the simulation with $\beta=600$. (b) The minimum and
maximum surface concentration $C_{\text{b}}J_{\text{s}}$ on the leaflets.
Panel (c) highlights the accumulation over the last three cycles for
accumulation on both leaflets for three values of $\beta$. In panel (d), the
accumulation is shown for accumulation only on the right leaflet versus both
leaflets with $\beta=600$. The vertical lines denote the beginning of systole.
There is a consistent accumulation of material on the leaflets. The jumps in
the minimum concentration are due to the the change in $J_{\text{s}}$, which
changes most where the leaflets attach to the aortic wall.

(a)

(b) (c)

(d)

(e)

Figure 8: The magnitude of the velocity field at peak systole during the last
cycle for (a) no accumulation, (b) $\beta=100$ with accumulation on both
leaflets, (c) $\beta=300$ with accumulation on both leaflets, (d) $\beta=600$
with accumulation on both leaflets, and (e) $\beta=600$ with accumulation on
only the right leaflet. Notice the presence of a vortex in the right sinus
that is absent if accumulation occurs only on the right leaflet. The choice of
colorbar is intended to highlight vortex formation in the sinus region. We
observe peak flow rates through the valve of
$170\text{\,}\frac{\mathrm{cm}}{\mathrm{s}}$ for $\beta=1$ and
$275\text{\,}\frac{\mathrm{cm}}{\mathrm{s}}$ for $\beta=600$.

###  4.3 Pressures and Flow Rates

Here we quantify the valve’s resistance to the flow at different maximum
stiffnesses. We measure the pressures at locations just upstream and
downstream of the valve. We use a Gaussian filter to smooth the curves in both
space and time, yielding the results shown in Figures 9 and 10. We observe
marginal increases in the pressures upstream of the valve of about
$14\text{\,}\frac{\mathrm{mmHg}}{}$; however, there are sharp decreases of
about $515\text{\,}\frac{\mathrm{mmHg}}{}$ in the aortic pressures downstream
of the leaflets as we increase stiffness. A similar trend is observed with
deposition on only the right leaflet, although the differences are not as
pronounced.

(a)

(b)

Figure 9: The pressures just upstream and downstream of the valve. Panel (a)
compares the pressures for no accumulation and accumulation on both leaflets
with $\beta=600$. Panel (b) compares pressures for accumulation on only the
right leaflet and both leaflets, with $\beta=600$. While the upstream pressure
increases mildly, the downstream pressures decrease by between $5$ and
$10\text{\,}\frac{\mathrm{mmHg}}{}$.

(a)

(b)

Figure 10: The pressures during the last three cycles for accumulation on both leaflets with $\beta=1$, 100, 300, and 600. Note that $\beta=1$ corresponds to a baseline where the leaflet stiffness is constant. The aortic side pressures are shown in panel (a). The left ventricle side pressures are shown in panel (b). While the upstream pressure increases by less than $5\text{\,}\frac{\mathrm{mmHg}}{}$, the downstream pressures decrease by between $5$ and $15\text{\,}\frac{\mathrm{mmHg}}{}$. Table 2: Pressures just upstream of the leaflets (LVOT) and downstream of the leaflets (aorta) during peak systole. While we observe increases of $12\text{\,}\frac{\mathrm{mmHg}}{}$ in the pressure upstream of the valve, there are greater decreases of $1015\text{\,}\frac{\mathrm{mmHg}}{}$ in the pressure downstream of the valve. | Base | Both Leaflets | Right Leaflet
---|---|---|---
|  | $\beta=100$ | $\beta=300$ | $\beta=600$ | $\beta=600$
LVOT | $117\text{\,}\frac{\mathrm{mmHg}}{}$ | $117\text{\,}\frac{\mathrm{mmHg}}{}$ | $118\text{\,}\frac{\mathrm{mmHg}}{}$ | $119\text{\,}\frac{\mathrm{mmHg}}{}$ | $117\text{\,}\frac{\mathrm{mmHg}}{}$
Aorta | $111\text{\,}\frac{\mathrm{mmHg}}{}$ | $105\text{\,}\frac{\mathrm{mmHg}}{}$ | $101\text{\,}\frac{\mathrm{mmHg}}{}$ | $94\text{\,}\frac{\mathrm{mmHg}}{}$ | $104\text{\,}\frac{\mathrm{mmHg}}{}$

We additionally compute the effective orifice area (EOA). The EOA
$A_{\text{AV}}$ is computed using conservation of mass in equation 1b by
assuming the relation
$V_{\text{LVOT}}A_{\text{LVOT}}=V_{\text{AV}}A_{\text{AV}},$ in which
$V_{\text{LVOT}}$ is the average time integral flow rate through the left
ventricle outflow tract during each cycle when the valve is open,
$A_{\text{LVOT}}$ is the area of the left ventricle outflow tract, and
$V_{\text{AV}}$ is the average time integral flow rate through the aortic
valve. $V_{\text{LVOT}}$ is computed from the boundary condition model
described in section 2.2. To compute $V_{\text{AV}}$, we interpolate the
velocity to the midpoint between the two points on each leaflet that are
closest during systole. $V_{\text{AV}}$ is then computed as the time integral
of the component of this interpolated velocity normal to the valve ring.
Figure 11 plots the EOA for each cycle. We observe a general decrease in EOA
as more the total surface concentration $C_{\text{b}}$ increases. This
indicates that the fluid velocity increases to compensate for the stiffening
of the valve. The EOA decreases more when accumulation occurs on both leaflets
when compared to accumulation on only the right leaflet. Because the left
leaflet remains at a constant stiffness, the left leaflet opens more to
compensate for the stiffening of the right leaflet. This causes the jet to
shift towards the left leaflet, as shown in Figure 4.

(a) (b)

Figure 11: The effective orifice area (EOA) for each cycle. Panel (a) shows
the EOA with $\beta=1$, 100, 300, and 600 for accumulation on both leaflets.
Note that $\beta=1$ corresponds to a baseline where the leaflet stiffness is
constant. Panel (b) shows the EOA for accumulation on no leaflets (base),
accumulation on the right leaflet (right), and accumulation on both leaflets
(both). We observe a general decrease of EOA over time as the accumulation
increases. There is a greater decrease of EOA if accumulation occurs on both
leaflets than just the right leaflet.

##  5 Conclusions

This study presents new numerical methods incorporating both deposition and
fluid-structure interaction to simulate leaflet thrombosis. The simplified
thrombosis model serves as a stepping stone to demonstrate the capabilities of
our simulation approach that includes concentration fields describing fluid-
phase platelets and structure-bound platelets. Platelets can deposit onto the
leaflet surface, and bound platelets can dissociate into the fluid. In our
model, the stiffness of the leaflet is a function of the bound platelet
concentration. We have shown that our model is capable of realizing drops in
pressure and decreases in effective orifice area, without fully occluding the
aortic valve. The results also show that the stiffness of the valve can lead
to a variety of flow features in the sinus of Valsalva region. These flow
features affect the amount of material that is locally present to deposit over
the leaflets.

Extensions of this model to three dimensions require an efficient method for
solving the advection-diffusion equation in complex, time-evolving domains.
The method utilized here requires the computation of cut-cell volumes and
intersections, which remain challenging in three spatial dimensions. Recent
approaches to this class of problems include mesh-free RBF-FD methods [52] and
volume penalization methods [54]. The implementation of a more physiological
model of thrombosis remains important future work. A primary roadblock is the
disparate time scales present in thrombosis. While the heart beats on the
order of seconds, blood clots can form in hours to days. The use of
conditionally stable time stepping limits the numerical methods to time steps
that resolve the fastest timescale, which in this model is that of the fluid-
structure interaction. Recent work in multiscale time stepping algorithms [16,
17] could enable extensions of our modeling framework to enable such long-time
simulations. Further, with multiscale time stepping algorithms, this model
could be extended to study the affect of saturation of the bound concentration
field. While platelet deposition is important and the beginning step of
thrombus formation, a significant portion of the clot may be from coagulation
and fibrin mesh formation. However, a complete model of thrombosis will
require a computational model in which the blood clot on the moving valve
leaflets grows into the fluid [11, 34]. The development of such a model that
incorporates FSI is ongoing.

The new approaches described herein should be considered an important
steppingstone for thrombosis models in many different contexts. This model is
the first of its kind to incorporate both adhesion of a surface concentration
to the surface of the leaflets and feedback into the fluid-structure
interaction. Further, this model can be adapted to model deposition or
absorption along other moving boundaries, such as for particulate flow in the
lungs or drug absorption in the gut.

## Acknowledgements

We acknowledge funding from the NIH (Awards U01HL143336 and R01HL157631) and
NSF (OAC 1652541 and OAC 1931516).

## References

  * [1] J. Bærentzen and Henrik Aanæs “Signed distance computation using the angle weighted pseudonormal” In _IEEE Transactions on Visualization and Computer Graphics_ 11.3, 2005, pp. 243–253 DOI: 10.1109/TVCG.2005.49
  * [2] S Balay et al. “PETSc Users Manual” ISBN: ANL-95/11 - Revision 3.0.0 Publication Title: ReVision Volume: 2, 2010, pp. 1–211 URL: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.123.6845&rep=rep1&type=pdf
  * [3] Satish Balay et al. “PETSc Web page”, 2019 URL: https://www.mcs.anl.gov/petsc
  * [4] Satish Balay, William D. Gropp, Lois Curfman McInnes and Barry F. Smith “Efficient Management of Parallelism in Object-Oriented Numerical Software Libraries” In _Modern Software Tools for Scientific Computing_ Birkhäuser Press, 1997, pp. 163–202 DOI: 10.1007/978-1-4612-1986-6˙8
  * [5] Aaron Barrett, Aaron L. Fogelson and Boyce E. Griffith “A hybrid semi-Lagrangian cut cell method for advection-diffusion problems with Robin boundary conditions in moving domains” In _Journal of Computational Physics_ 449, 2022, pp. 110805 DOI: 10.1016/j.jcp.2021.110805
  * [6] Kuan-Yu Chen and Ming-Chih Lai “A conservative scheme for solving coupled surface-bulk convection–diffusion equations with an application to interfacial flows with soluble surfactant” In _Journal of Computational Physics_ 257, 2014, pp. 1–18 DOI: 10.1016/j.jcp.2013.10.003
  * [7] Scott L. Diamond, Jeremy Purvis, Manash Chatterjee and Matthew H. Flamm “Systems biology of platelet-vessel wall interactions” In _Frontiers in Physiology_ 4, 2013 DOI: 10.3389/fphys.2013.00229
  * [8] Jian Du and Aaron L. Fogelson “A Two-phase mixture model of platelet aggregation” In _Mathematical Medicine and Biology_ 35.2, 2018, pp. 225–256 DOI: 10.1093/imammb/dqx001
  * [9] Jian Du et al. “Clot Permeability, Agonist Transport, and Platelet Binding Kinetics in Arterial Thrombosis” In _Biophysical Journal_ 119.10, 2020, pp. 2102–2115 DOI: 10.1016/j.bpj.2020.08.041
  * [10] Natasha Flyer, Bengt Fornberg, Victor Bayona and Gregory A. Barnett “On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy” In _Journal of Computational Physics_ 321, 2016, pp. 21–38 DOI: 10.1016/j.jcp.2016.05.026
  * [11] Aaron L. Fogelson and Robert D. Guy “Immersed-boundary-type models of intravascular platelet aggregation” In _Computer Methods in Applied Mechanics and Engineering_ 197.25-28, 2008, pp. 2087–2104 DOI: 10.1016/j.cma.2007.06.030
  * [12] Aaron L. Fogelson, Yasmeen H. Hussain and Karin Leiderman “Blood Clot Formation under Flow: The Importance of Factor XI Depends Strongly on Platelet Count” Publisher: Cell Press In _Biophysical Journal_ 102.1, 2012, pp. 10–18 DOI: 10.1016/J.BPJ.2011.10.048
  * [13] Aaron L. Fogelson and James P. Keener “Toward an understanding of fibrin branching structure” In _Physical Review E - Statistical, Nonlinear, and Soft Matter Physics_ 81.5, 2010, pp. 1–9 DOI: 10.1103/PhysRevE.81.051922
  * [14] Aaron L. Fogelson and Keith B. Neeves “Fluid mechanics of blood clot formation” In _Annual Review of Fluid Mechanics_ 47.1, 2015, pp. 377–403 DOI: 10.1146/annurev-fluid-010814-014513
  * [15] Aaron L. Fogelson, Anna C. Nelson, Cheryl Zapata-Allegro and James P. Keener “Development of Fibrin Branch Structure Before and After Gelation” In _SIAM Journal on Applied Mathematics_ 82.1, 2022, pp. 267–293 DOI: 10.1137/21M1401024
  * [16] S. Frei and T. Richter “Efficient Approximation of Flow Problems With Multiple Scales in Time” Publisher: Society for Industrial and Applied Mathematics _eprint: 1903.12234 In _Multiscale Modeling & Simulation_ 18.2, 2020, pp. 942–969 DOI: 10.1137/19M1258396
  * [17] S. Frei, T. Richter and T. Wick “Long-term simulation of large deformation, mechano-chemical fluid-structure interactions in ALE and fully Eulerian coordinates” Publisher: Academic Press In _Journal of Computational Physics_ 321, 2016, pp. 874–891 DOI: 10.1016/J.JCP.2016.06.015
  * [18] T. Gasser, Ray W Ogden and Gerhard A Holzapfel “Hyperelastic modelling of arterial layers with distributed collagen fibre orientations” In _Journal of The Royal Society Interface_ 3.6, 2006, pp. 15–35 DOI: 10.1098/rsif.2005.0073
  * [19] Boyce E Griffith “IBAMR: an adaptive and distributed-memory parallel implementation of the immersed boundary method”, 2014 URL: https://github.com/IBAMR/IBAMR
  * [20] Boyce E. Griffith “An accurate and efficient method for the incompressible Navier-Stokes equations using the projection method as a preconditioner” ISBN: 0021-9991 Publisher: Elsevier Inc. In _Journal of Computational Physics_ 228.20, 2009, pp. 7565–7595 DOI: 10.1016/j.jcp.2009.07.001
  * [21] Boyce E. Griffith and Xiaoyu Luo “Hybrid finite difference/finite element immersed boundary method” Publisher: John Wiley & Sons, Ltd _eprint: 1612.05916 In _International Journal for Numerical Methods in Biomedical Engineering_ 33.12, 2017, pp. e2888 DOI: 10.1002/cnm.2888
  * [22] Boyce E. Griffith, Xiaoyu Luo, David M. McQueen and Charles S. Peskin “Simulating the Fluid Dynamics of Natural and Prosthetic Heart Valve Using the Immersed Boundary Method” In _International Journal of Applied Mechanics_ 01.01, 2009, pp. 137–177 DOI: 10.1142/S1758825109000113
  * [23] Hoda Hatoum et al. “Predictive Model for Thrombus Formation After Transcatheter Valve Replacement” In _Cardiovascular Engineering and Technology_ 12.6, 2021, pp. 576–588 DOI: 10.1007/s13239-021-00596-x
  * [24] Ásdís Helgadóttir, Yen Ting Ng, Chohong Min and Frédéric Gibou “Imposing mixed Dirichlet–Neumann–Robin boundary conditions in a level-set framework” In _Computers & Fluids_ 121, 2015, pp. 68–80 DOI: 10.1016/j.compfluid.2015.08.007
  * [25] Matthew M. Hopkins and Lisa J. Fauci “A computational model of the collective fluid dynamics of motile micro-organisms” In _Journal of Fluid Mechanics_ 455, 2002, pp. 149–174 DOI: 10.1017/S0022112001007339
  * [26] Richard D. Hornung and Scott R. Kohn “Managing application complexity in the SAMRAI object-oriented framework” In _Concurrency and Computation: Practice and Experience_ 14.5, 2002, pp. 347–368 DOI: 10.1002/cpe.652
  * [27] Huaxiong Huang, Kazuyasu Sugiyama and Shu Takagi “An immersed boundary method for restricted diffusion with permeable interfaces” In _Journal of Computational Physics_ 228.15, 2009, pp. 5317–5322 DOI: 10.1016/j.jcp.2009.04.040
  * [28] Benjamin S. Kirk, John W. Peterson, Roy H. Stogner and Graham F. Carey “libMesh : a C++ library for parallel adaptive mesh refinement/coarsening simulations” In _Engineering with Computers_ 22.3-4, 2006, pp. 237–254 DOI: 10.1007/s00366-006-0049-3
  * [29] Araz R. Kivi et al. “Fluid structure interaction modelling of aortic valve stenosis: Effects of valve calcification on coronary artery flow and aortic root hemodynamics” In _Computer Methods and Programs in Biomedicine_ 196, 2020, pp. 105647 DOI: 10.1016/j.cmpb.2020.105647
  * [30] Ebrahim M. Kolahdouz, Amneet Pal Singh Bhalla, Brent A. Craven and Boyce E. Griffith “An immersed interface method for discrete surfaces” In _Journal of Computational Physics_ 400, 2020, pp. 108854 DOI: 10.1016/j.jcp.2019.07.052
  * [31] Andrew L. Kuharsky and Aaron L. Fogelson “Surface-mediated control of blood coagulation: The role of binding site densities and platelet deposition” In _Biophysical Journal_ 80.3, 2001, pp. 1050–1074 DOI: 10.1016/S0006-3495(01)76085-7
  * [32] Jae H Lee and Boyce E Griffith “On the Lagrangian-Eulerian coupling in the immersed finite element/difference method” In _Journal of Computational Physics_ 457, 2022, pp. 111042 DOI: https://doi.org/10.1016/j.jcp.2022.111042
  * [33] Jae H. Lee et al. “Bioprosthetic aortic valve diameter and thickness are directly related to leaflet fluttering: Results from a combined experimental and computational modeling study” Publisher: Elsevier BV In _JTCVS Open_ 6, 2021, pp. 60–81 DOI: 10.1016/J.XJON.2020.09.002/ATTACHMENT/1012F996-5894-4BCC-86B2-6F05A1B63BAD/MMC6
  * [34] Karin Leiderman and Aaron L. Fogelson “Grow with the flow: A spatial-temporal model of platelet deposition and blood coagulation under flow” In _Mathematical Medicine and Biology_ 28.1, 2011, pp. 47–84 DOI: 10.1093/imammb/dqq005
  * [35] X. Li, J. Lowengrub, A. Ratz and A. Voigt “Solving PDEs in ComplexGeometries: A Diffuse Domain Approach” In _Communications in mathematical sciences_ 7.1, 2009, pp. 81–107 URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3097555/
  * [36] Kathryn G. Link et al. “A local and global sensitivity analysis of a mathematical model of coagulation and platelet deposition under flow” In _PLOS ONE_ 13.7, 2018, pp. e0200917 DOI: 10.1371/journal.pone.0200917
  * [37] Kathryn G. Link et al. “A mathematical model of coagulation under flow identifies factor V as a modifier of thrombin generation in hemophilia A” In _Journal of Thrombosis and Haemostasis_ 18.2, 2020, pp. 306–317 DOI: 10.1111/jth.14653
  * [38] Raj R. Makkar and Tarun Chakravarty “Transcatheter Aortic Valve Thrombosis” In _JACC: Cardiovascular Interventions_ 10.7, 2017, pp. 698–700 DOI: 10.1016/j.jcin.2017.02.041
  * [39] Raj R. Makkar et al. “Possible Subclinical Leaflet Thrombosis in Bioprosthetic Aortic Valves” In _New England Journal of Medicine_ 373.21, 2015, pp. 2015–2024 DOI: 10.1056/NEJMoa1509233
  * [40] Chohong Min and Frédéric Gibou “Geometric integration over irregular domains with application to level-set methods” In _Journal of Computational Physics_ 226.2, 2007, pp. 1432–1443 DOI: 10.1016/j.jcp.2007.05.032
  * [41] Kyle Murdock, Caitlin Martin and Wei Sun “Characterization of mechanical properties of pericardium tissue using planar biaxial tension and flexural deformation” In _Journal of the Mechanical Behavior of Biomedical Materials_ 77, 2018, pp. 148–156 DOI: 10.1016/j.jmbbm.2017.08.039
  * [42] J P Murgo, N Westerhof, J P Giolma and S A Altobelli “Aortic input impedance in normal man: relationship to pressure wave forms.” In _Circulation_ 62.1, 1980, pp. 105–116 DOI: 10.1161/01.CIR.62.1.105
  * [43] J.. Mynard, M.. Davidson, D.. Penny and J.. Smolich “A simple, versatile valve model for use in lumped parameter and one-dimensional cardiovascular models” In _International Journal for Numerical Methods in Biomedical Engineering_ 28.6-7, 2012, pp. 626–641 DOI: 10.1002/cnm.1466
  * [44] Romina Plitman Mayo et al. “Numerical models for assessing the risk of leaflet thrombosis post-transcatheter aortic valve-in-valve implantation” In _Royal Society Open Science_ 7.12, 2020, pp. 201838 DOI: 10.1098/rsos.201838
  * [45] Ravi Ramana et al. “Calcification and Thrombosis as Mediators of Bioprosthetic Valve Deterioration” In _Structural Heart_ 3.2, 2019, pp. 106–109 DOI: 10.1080/24748706.2018.1562265
  * [46] William J. Rider, Jeffrey A. Greenough and James R. Kamm “Accurate monotonicity- and extrema-preserving methods through adaptive nonlinear hybridizations” In _Journal of Computational Physics_ 225.2, 2007, pp. 1827–1848 DOI: 10.1016/j.jcp.2007.02.023
  * [47] Liesbeth Rosseel, Ole De Backer and Lars Søndergaard “Clinical Valve Thrombosis and Subclinical Leaflet Thrombosis Following Transcatheter Aortic Valve Replacement: Is There a Need for a Patient-Tailored Antithrombotic Therapy?” In _Frontiers in Cardiovascular Medicine_ 6, 2019 DOI: 10.3389/fcvm.2019.00044
  * [48] Yongyuth Sahasakul, William D. Edwards, James M. Naessens and A.Jamil Tajik “Age-related changes in aortic and mitral valve thickness: Implications for two-dimensional echocardiography based on an autopsy study of 200 normal human hearts” In _The American Journal of Cardiology_ 62.7, 1988, pp. 424–430 DOI: 10.1016/0002-9149(88)90971-X
  * [49] Matea Santiago, Nicholas A. Battista, Laura A. Miller and Shilpa Khatri “Passive concentration dynamics incorporated into the library IB2d, a two-dimensional implementation of the immersed boundary method” Publisher: IOP Publishing In _Bioinspiration &amp$\mathsemicolon$ Biomimetics_ 17.3, 2022, pp. 036003 DOI: 10.1088/1748-3190/ac4afa
  * [50] Matea Santiago, Kevin A. Mitchell and Shilpa Khatri “Numerical method for modeling photosynthesis of algae on pulsing soft corals” In _Physical Review Fluids_ 7.3, 2022, pp. 033102 DOI: 10.1103/PhysRevFluids.7.033102
  * [51] Gerhard Schymik et al. “How to Adapt the Implantation Technique for the New SAPIEN 3 Transcatheter Heart Valve Design” In _Journal of Interventional Cardiology_ 28.1, 2015, pp. 82–89 DOI: 10.1111/joic.12165
  * [52] Varun Shankar and Aaron L. Fogelson “Hyperviscosity-based stabilization for radial basis function-finite difference (RBF-FD) discretizations of advection–diffusion equations” _eprint: 1806.03798 In _Journal of Computational Physics_ 372, 2018, pp. 616–639 DOI: 10.1016/j.jcp.2018.06.036
  * [53] Nikos Stergiopulos, Berend E. Westerhof and Nico Westerhof “Total arterial inertance as the fourth element of the windkessel model” In _American Journal of Physiology-Heart and Circulatory Physiology_ 276.1, 1999, pp. H81–H88 DOI: 10.1152/ajpheart.1999.276.1.H81
  * [54] Ramakrishnan Thirumalaisamy, Neelesh A. Patankar and Amneet Pal Singh Bhalla “Handling Neumann and Robin boundary conditions in a fictitious domain volume penalization framework” Publisher: Academic Press _eprint: 2101.02806 In _Journal of Computational Physics_ 448, 2022, pp. 110726 DOI: 10.1016/J.JCP.2021.110726
  * [55] A. Tosenberger et al. “Modelling of platelet–fibrin clot formation in flow with a DPD–PDE method” In _Journal of Mathematical Biology_ 72.3, 2016, pp. 649–681 DOI: 10.1007/s00285-015-0891-2
  * [56] John D. Towers “A source term method for Poisson problems on irregular domains” In _Journal of Computational Physics_ 361, 2018, pp. 424–441 DOI: 10.1016/j.jcp.2018.01.038
  * [57] Ben Vadala-Roth et al. “Stabilization approaches for the hyperelastic immersed boundary method for problems of large-deformation incompressible elasticity” _eprint: 1811.06620 In _Computer Methods in Applied Mechanics and Engineering_ 365, 2020, pp. 112978 DOI: 10.1016/j.cma.2020.112978
  * [58] Koohyar Vahidkhah et al. “Blood Stasis on Transcatheter Valve Leaflets and Implications for Valve-in-Valve Leaflet Thrombosis” Publisher: Elsevier In _The Annals of Thoracic Surgery_ 104.3, 2017, pp. 751–759 DOI: 10.1016/J.ATHORACSUR.2017.02.052
  * [59] Ziheng Wu, Zhiliang Xu, Oleg Kim and Mark Alber “Three-dimensional multi-scale model of deformable platelets adhesion to vessel wall in blood flow” In _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 372.2021, 2014, pp. 20130380 DOI: 10.1098/rsta.2013.0380
  * [60] Sheng Xu and Z. Wang “An immersed interface method for simulating the interaction of a fluid with moving boundaries” In _Journal of Computational Physics_ 216.2, 2006, pp. 454–493 DOI: 10.1016/j.jcp.2005.12.016
  * [61] Paul A Yushkevich et al. “User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability” In _NeuroImage_ 31.3, 2006, pp. 1116–1128 DOI: https://doi.org/10.1016/j.neuroimage.2006.01.015
  * [62] Paul A. Yushkevich et al. “User-Guided 3D Active Contour Segmentation of Anatomical Structures: Significantly Improved Efficiency and Reliability” In _Neuroimage_ 31.3, 2006, pp. 1116–1128
  * [63] Peng Zhang et al. “Multiscale Particle-Based Modeling of Flowing Platelets in Blood Plasma Using Dissipative Particle Dynamics and Coarse Grained Molecular Dynamics” In _Cellular and Molecular Bioengineering_ 7.4, 2014, pp. 552–574 DOI: 10.1007/s12195-014-0356-5