# FSGe: A fast and strongly-coupled 3D fluid-solid-growth interaction method

Martin R. Pfaller  Marcos Latorre  Erica L. Schwarz  Fannie M. Gerosa  Jason
M. Szafron  Jay D. Humphrey  Alison L. Marsden

###### Abstract

Equilibrated fluid-solid-growth (FSGe) is a fast, open source, three-
dimensional (3D) computational platform for simulating interactions between
instantaneous hemodynamics and long-term vessel wall adaptation through growth
and remodeling (G&R). Such models are crucial for capturing adaptations in
health and disease and following clinical interventions. In traditional G&R
models, this feedback is modeled through highly simplified fluid models,
neglecting local variations in blood pressure and wall shear stress (WSS).
FSGe overcomes these inherent limitations by strongly coupling the 3D Navier-
Stokes equations for blood flow with a 3D equilibrated constrained mixture
model (CMMe) for vascular tissue G&R. CMMe allows one to predict long-term
evolved mechanobiological equilibria from an original homeostatic state at a
computational cost equivalent to that of a standard hyperelastic material
model. In illustrative computational examples, we focus on the development of
a stable aortic aneurysm in a mouse model to highlight key differences in
growth patterns and fluid-solid feedback between FSGe and solid-only G&R
models. We show that FSGe is especially important in blood vessels with
asymmetric stimuli. Simulation results reveal greater local variation in
fluid-derived WSS than in intramural stress (IMS). Thus, differences between
FSGe and G&R models became more pronounced with the growing influence of WSS
relative to pressure. Future applications in highly localized disease
processes, such as for lesion formation in atherosclerosis, can now include
spatial and temporal variations of WSS.

††journal: Computer Methods in Applied Mechanics and Engineering

[inst1] organization=Department of Pediatrics – Cardiology, Stanford
Univeristy, city=Stanford, state=CA 94306, country=USA

[inst2] organization=Center for Research and Innovation in Bioengineering,
Universitat Politècnica de València, city=València, country=Spain

[inst3] organization=Department of Bioengineering, Stanford Univeristy,
city=Stanford, state=CA 94306, country=USA

[inst4] organization=Department of Biomedical Engineering, Yale Univeristy,
city=New Haven, state=CT 06520, country=USA

##  1 Introduction

The intricate interaction between hemodynamics and growth (change in mass) and
remodeling (change in microstructure) of a blood vessel plays a vital role in
development, homeostasis, and disease progression [1]. The constrained mixture
model (CMM) simulates the continuous deposition and degradation of tissue
constituents [2, 3], often based on the concept of mechanical homeostasis,
which is visualized in Figure 1. In this model of mechanobiological
homeostasis, changes in vascular composition and properties stem from
deviations in regulated variables from defined set points. These regulated
variables primarily include intramural stress (IMS), $\sigma_{I}$, and wall
shear stress (WSS), $\tau_{w}$, which have consistently been found to be
maintained at homeostatic values under physiologic conditions in vivo [4].

Hemodynamics $\mathcal{F}$ are a crucial determinant of the biomechanical
state of a blood vessel, characterized in part by its regulated variables.
Solid mechanics $\mathcal{S}$ within the vessel wall define blood vessel
geometry and microstructure. Changes in the solid mechanics thus trigger
changes in hemodynamics, $\mathcal{F}$, rendering the feedback loop a fluid-
solid-growth interaction (FSG) problem. Most prior work considered either
constant or reduced-order hemodynamics. The novelty of our work is a strongly
coupled feedback between hemodynamics $\mathcal{F}$ and solid mechanics
$\mathcal{S}$ in a fully three-dimensional (3D) framework based on the
equilibrated constrained mixture model (CMMe) [5]. In the following, we review
previous approaches in more detail and highlight the contributions of this
work.

Figure 1: Negative feedback characteristics of mechanical homeostasis,
including hemodynamics $\mathcal{F}$ and solid mechanics $\mathcal{S}$ .
Stimuli consist of deviations of intramural stress, $\Delta\sigma_{I}$, and
wall shear stress, $\Delta\tau_{w}$, from homeostatic set-point values. Tissue
changes are modulated by gain factors $K^{\alpha}_{\sigma_{I}}$ and
$K^{\alpha}_{\tau_{w}}$, which capture cell sensitivity to the particular
stimulus. G&R are influenced further by degradation rate $k^{\alpha}$,
homeostatic prestretch $G^{\alpha}$, and orientation angle $\alpha^{\alpha}$.
In this work, tissue constituents $\alpha$ are elastin $e$, four collagen
fiber families $c$, and smooth muscle $m$, which are visualized in Figure 2.
Figure adapted from [3].

Hemodynamics in vivo typically include complex flow phenomena such as flow
separation and recirculation induced by changes in cross-sectional area or
bifurcations [6]. Consequently, pressure and WSS typically show a high spatial
and temporal variation. Local hemodynamics can be quantified with
computational fluid dynamics (CFD) or fluid-structure interaction (FSI)
solutions [7]. However, the CMM commonly includes simplified assumptions for
axial pressure gradient $\Delta p$ and WSS $\tau_{w}$ based on
Hagen–Poiseuille law [5]

$\Delta p=\frac{8\mu lQ}{\pi r^{4}},\quad\tau_{w}=\frac{4\mu Q}{\pi r^{3}},$
(1)

with dynamic blood viscosity $\mu$, vessel length $l$ and radius $r$, and bulk
blood flow rate $Q$. These equations assume fully developed steady-state flow
through a straight cylindrical tube with axisymmetric flow and zero radial and
circumferential components, which is generally not upheld in in vivo settings.

Reference | Fluid | Solid | Stimuli | Coupling | Geometry
---|---|---|---|---|---
Figueroa et al. [8] |  3D Navier-Stokes | Membrane CMM | pulsatile mean | weak | idealized
Sheidaei et al. [9] |  3D Navier-Stokes | Membrane CMM | pulsatile mean | weak | in vivo
Watton et al. [10] |  3D Navier-Stokes |  3D hybrid | pulsatile mean | weak | in vivo
Grytsan et al. [11] |  3D Navier-Stokes | 3D hybrid | steady-state | weak | in vivo
Latorre et al. [12] | Control-volume |  3D CMMe | steady-state | strong | idealized
Schwarz et al. [13] |  3D Navier-Stokes |  3D CMM | steady-state | strong | in vivo
This work |  3D Navier-Stokes |  3D CMMe | steady-state | strong | idealized
Table 1: Overview of some fluid-solid-growth interaction (FSG) interaction
models.

The local interaction between hemodynamics and vessel growth and remodeling
(G&R) mechanisms is crucial in numerous diseases and their treatments, as, for
example, in tissue-engineered vascular grafts (TEVGs) [14, 15], vein graft
adaptation and failure following coronary artery bypass graft surgery [16,
17], and in congenital heart patients with pulmonary disease [18, 19]. For an
overview of G&R models, i.e., kinematic, constrained mixture model (CMM), and
hybrid, see [20]. A non-comprehensive overview of similar fluid-solid-growth
interaction models is given in Table 1. Figueroa et al. [8] developed a method
that couples three-dimensional (3D) hemodynamics to a membrane model of the
vessel wall. This theory is based on the theory of small on large [21],
superimposing "small" deformations during the cardiac cycle upon "large"
deformations due to G&R. The membrane wall formulation is related to the
coupled momentum method for fluid-structure interaction (FSI) in blood vessels
[22]. The long-term G&R deformation is coupled to pulsatile hemodynamics,
illustrated in an example of aneurysmal development. This approach was used in
a patient-specific aneurysm geometry in Sheidaei et al. [9]. The work of
Latorre et al. [12] couples 3D G&R with a simplified control-volume
formulation of fluid dynamics. Pressure and wall shear stress (WSS) are
derived from an axial node-based control-volume approach based on the Hagen-
Poiseuille law (1) with additional pressure loss terms. For a more detailed
discussion of reduced-order fluid dynamics, see Pfaller et al. [23]. Most
recently, Schwarz et al. [13] coupled 3D CMM and 3D hemodynamics. While this
constitutes the most detailed approach so far, it comes at high computational
cost due to the need to integrate heredity integrals point-wise in full CMM.

Most FSG approaches reviewed above, except for Latorre et al. [12] and Schwarz
et al. [13], employ weak coupling between fluid and solid G&R, also known as
loose or weak coupling. Here, the fluid and solid domains are solved
independently, with their solutions being exchanged only once in each load
step. This approach assumes that the interaction between the fluid and solid
is relatively weak, allowing for some degree of decoupling between the
domains. Loose coupling is typically used for problems where the fluid and
solid have a minor influence on each other or when computational efficiency is
prioritized over solution accuracy. On the other hand, strong coupling schemes
rely on an iterative procedure that ensures the convergence of fluid and solid
solutions.

This work introduces an equilibrated fluid-solid-growth (FSGe) method to
strongly couple 3D G&R based on the equilibrated constrained mixture model
(CMMe) with 3D hemodynamics. Using CMMe greatly reduces the computational cost
of evaluating the solid model and generally arrives at the desired model
prediction in much fewer steps. In the following, we introduce the governing
differential equations of the fluid and solid domains (Section 2). We then
review computational techniques employed to solve the coupled fluid-solid
problem (Section 3). In illustrative computational examples of an asymmetric
aneurysm, we compare key differences between a solid-only G&R model and our
FSGe formulation (Section 4). We close with a discussion of the results,
limitations, and future perspectives (Section 5).

##  2 Governing equations of fluid and solid domains

As outlined in Figure 1, the FSGe model relies on the solution of fluid (2.1)
and solid (2.2) mechanics, whose governing equations we briefly review in this
section. Additionally, for completeness we introduce both the full CMM (2.2.1)
and CMMe (2.2.2) G&R models. Both 3D fluid and solid were solved numerically
in our open-source multi-physics solver svFSIplus [24], the C++ version of
svFSI [25]; both solvers are released with the SimVascular project [26].

###  2.1 Fluid dynamics

We model blood flow in large vessels at high shear rates as a Newtonian fluid
in the quasi-static fluid domain $\Omega^{\text{f}}$ within time scale
$T^{\text{f}}$ with the incompressible Navier-Stokes equations,

$\displaystyle\nabla\cdot\bm{\sigma}^{\text{f}}(\bm{u},p)$
$\displaystyle=\rho\left[\frac{\partial\bm{u}}{\partial
t}+(\bm{u}\cdot\nabla)\bm{u}\right],\quad$
$\displaystyle\bm{x}\in\Omega^{\text{f}},\quad t\in[0,T^{\text{f}}],$ (2)
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0,\quad$
$\displaystyle\bm{x}\in\Omega^{\text{f}},\quad t\in[0,T^{\text{f}}],$ (3)

with density $\rho$, viscosity $\mu$, and fluid Cauchy stress tensor
$\bm{\sigma}^{\text{f}}=\mu(\nabla\bm{u}+\nabla\bm{u}^{\intercal})-p\bm{I}$,
where $\bm{u}$ and $p$ are fluid velocity and pressure, respectively. We
prescribe boundary conditions

$\displaystyle\bm{u}(\bm{x},t)$ $\displaystyle=\mathbf{0},\qquad$
$\displaystyle\bm{x}\in\Gamma_{\text{int}},\quad$ $\displaystyle
t\in[0,T^{\text{f}}],$ $\displaystyle\bm{u}(\bm{x},t)$
$\displaystyle=\bm{u}_{\text{in}}(\bm{x}),\qquad$
$\displaystyle\bm{x}\in\Gamma_{\text{f,in}},\quad$ $\displaystyle
t\in[0,T^{\text{f}}],$
$\displaystyle\bm{n}\cdot\bm{\sigma}^{\text{f}}(\bm{x},t)\cdot\bm{n}$
$\displaystyle=p_{\text{out}},\qquad$
$\displaystyle\bm{x}\in\Gamma_{\text{f,out}},\quad$ $\displaystyle
t\in[0,T^{\text{f}}],$

with a no-slip boundary condition at the fluid-solid interface
$\Gamma_{\text{int}}$, constant quadratic inflow profile $\bm{u}_{\text{in}}$
at the inlet $\Gamma_{\text{f,in}}$, and constant outlet pressure
$p_{\text{out}}$ at the outlet $\Gamma_{\text{f,out}}$ with normal $\bm{n}$.
The initial conditions are

$\displaystyle\bm{u}(\bm{x},t=0)$ $\displaystyle=\bm{u}_{0}(\bm{x}),\qquad$
$\displaystyle\bm{x}\in\Omega^{\text{f}},$ $\displaystyle p(\bm{x},t=0)$
$\displaystyle=p_{0}(\bm{x}),\qquad$
$\displaystyle\bm{x}\in\Omega^{\text{f}},$

with initial velocity field $\bm{u}_{0}$ and initial pressure field $p_{0}$.
See Esmaily Moghadam et al. [27, 28] for details on the P1-P1 variational
multiscale finite element solution of the Navier-Stokes equations.

###  2.2 Solid growth and remodeling

We follow the classic approach of finite strain theory to model the vessel
wall displacement, $\bm{d}=\bm{x}-\bm{X}$, with reference position $\bm{X}$
and current position $\bm{x}$. We calculate the deformation gradient $\bm{F}$,
the Jacobian $J$, and the right Cauchy–Green tensor $\bm{C}$ as

$\bm{F}=\frac{\partial\bm{x}}{\partial\bm{X}},\quad
J=\text{det}\,\bm{F},\quad\bm{C}=\bm{F}^{\text{T}}\bm{F}.$

The balance of linear momentum in the case of negligible inertial loads [29]
in the solid domain $\Omega^{\text{s}}$ yields the quasi-static problem

$\displaystyle\nabla\cdot\bm{\sigma}^{\text{s}}=\bm{0},\qquad$
$\displaystyle\bm{x}\in\Omega^{\text{s}},\quad$ $\displaystyle
t\in[0,T^{\text{s}}],$ (4)

with solid stress tensor $\bm{\sigma}^{\text{s}}$ provided constitutively by
the G&R material and its solid time scale $T^{\text{s}}$. We prescribe
boundary conditions

$\displaystyle d_{a}(\bm{x},t)$ $\displaystyle=0,\qquad$
$\displaystyle\bm{x}\in\Gamma_{\text{cap}},\quad$ $\displaystyle
t\in[0,T^{\text{s}}],$
$\displaystyle\bm{\sigma}^{\text{s}}(\bm{x},t)\cdot\bm{n}(\bm{x},t)$
$\displaystyle=k\,\bm{d}(\bm{x},t),\qquad$
$\displaystyle\bm{x}\in\Gamma_{\text{out}},\quad$ $\displaystyle
t\in[0,T^{\text{s}}],$

with zero Dirichlet boundary conditions in axial direction $z$ at the caps of
the blood vessel $\Gamma_{\text{cap}}$ and external tissue support with
outward normal $\bm{n}$ and stiffness $k$ [30]. Note that our external tissue
support does not include viscous damping since G&R is a quasi-static process.
The initial condition is

$\displaystyle\bm{d}(\bm{x},t=0)$ $\displaystyle=\bm{0},\qquad$
$\displaystyle\bm{x}\in\Omega^{\text{s}},$

We use a G&R constitutive law to relate deformation $\bm{F}$ to stresses
$\bm{\sigma}^{\text{f}}$. In the following, we briefly review the full CMM and
the CMMe. The CMM has been diversely applied, e.g., in altered flow and
pressure conditions [31], abdominal [32] and thoracic [33] aortic aneurysms,
vein grafts [34], or tissue-engineered vascular grafts [35]. However, due to
its computational complexity, there are only limited implementations in three
dimensions [36, 37]. In this work, we use the fast, rate-independent CMMe
method that directly predicts the long-term evolved solution that is in
mechanobiological equilibrium [38]. This material model had been previously
implemented in the open-source finite element solver FEBio [39], with
computational costs comparable to hyperelastic materials [5]. Prior work
relied on constant-pressure models [5, 33, 40] and control-volume hemodynamics
[12]. We briefly introduce the CMM in Section 2.2.1 and outline the CMMe in
Section 2.2.2 but refer interested readers to the original publications.

Figure 2: Blood vessel with length $l$, radius $r$, and thickness $h$. Its
constituents are elastin $e$ (isotropic), smooth muscle $m$ (red,
circumferential), and four collagen fiber families $c_{i}$ (circumferential,
axial, diagonal angle $\alpha_{0}$). Image created by Sebastian L. Fuchs and
licensed under the Creative Commons Attribution 4.0 International License.

####  2.2.1 Constrained Mixture Model

The CMM assumes that structurally significant constituents $\alpha$ with
individual reference configurations, material properties, mass balance
relations, and turnover rates are constrained to move together. The mass
production per unit volume at G&R time $\tau\in[0,s]$, where $s$ is the
current G&R time and $\tau$ is an intermediate time at which constituents are
deposited, is expressed as

$\displaystyle m^{\alpha}(\tau)=m^{\alpha}_{o}\,\Upsilon^{\alpha}(\tau),$

using an original rate $m^{\alpha}_{o}>0$ that is modulated by a stimulus
function, as, for example,

$\displaystyle\Upsilon^{\alpha}(\tau)=1+K^{\alpha}_{\sigma}\,\Delta\sigma_{I}(\tau)-K^{\alpha}_{\tau}\,\Delta\tau_{w}(\tau_{w}),$
(5)

with gain parameters $K^{\alpha}_{i}>0$ controlling the sensitivity to stress
deviations $\Delta\sigma_{I}$ and $\Delta\tau_{w}$, respectively. We define
these stimuli as

$\displaystyle\Delta\sigma_{I}(\tau)=\frac{\sigma_{I}(\tau)}{\sigma_{Io}}-1,\quad\Delta\tau_{w}(\tau)=\frac{\tau_{w}(\tau)}{\tau_{wo}}-1.$
(6)

governed by changes in intramural stress (IMS) $\sigma_{I}$ (induced by
pressure and axial force) and WSS $\tau_{W}$ (induced by blood flow) from
homeostatic set points $\sigma_{Io}$ and $\tau_{wo}$, respectively. The
constituent removal between deposition times $\tau$ and current time $s$ is
modeled by an exponential decay

$q^{\alpha}(s,\tau)=\exp\left(-\int_{\tau}^{s}k^{\alpha}(t)\,\mathrm{d}t\right),\quad
k^{\alpha}(t)=k^{\alpha}_{o}\,\left[1+K^{\alpha}\,\Delta\sigma^{2}(t)\right],$

where $k^{\alpha}$ is a rate parameter with homeostatic value
$k^{\alpha}_{o}$, and $K^{\alpha}>0$ is a gain parameter for IMS deviation. We
can express the mass density per current volume of each constituent
$\rho^{\alpha}_{R}$ by

$\rho^{\alpha}(s)=\int_{-\infty}^{s}m^{\alpha}(\tau)\,q^{\alpha}(s,\tau)\,\mathrm{d}\tau,$
(7)

where $m^{\alpha}$ is the mass density production rate. The constituent-
specific elastic strain energy density functions $\hat{W}^{\alpha}$ for
elastin $e$, collagen $c$, and smooth muscle $m$ are characterized by the neo-
Hookean and Fung-type exponential relations as

$\displaystyle\hat{W}^{e}=c^{e}(I^{e}_{1}-3),\quad\hat{W}^{m,c}=\frac{c_{1}^{m,c}}{4c_{2}^{m,c}}\left[\exp\left(c_{2}^{m,c}(I_{4}^{m,c}-1)^{2}\right)-1\right],$
(8)

where $c^{e},c_{1}^{m,c},c_{2}^{m,c}$ are material parameters and
$I^{e}_{1},I_{4}^{m,c}$ are constituent-specific measures of deformation that
are invariant to the coordinate system used. The constituents are visualized
in Figure 2. Collagen fibers are further classified into circumferentially,
axially, and diagonally oriented (at angle $\alpha_{0}$) populations with
fractions $\beta_{\theta}$, $\beta_{z}$, and $\beta_{d}$, respectively. The
constituent relations (8) are evaluated as
$\hat{W}^{\alpha}(\bm{F}^{\alpha}_{n(\tau)}(s))$ with the deformation gradient

$\displaystyle\bm{F}^{\alpha}_{n(\tau)}(s)=\bm{F}(s)\,\bm{F}^{-1}(\tau)\,\bm{G}^{\alpha}(\tau),$
(9)

where $\bm{G}^{\alpha}$ is the constituent-specific preferred homeostatic pre-
stretch, which can be reduced to the scalars $G^{m,c}$ for smooth muscle cells
and collagen fibers. The elastin prestretch tensor $\bm{G}^{e}$ is
characterized by the components $G^{e}_{\theta,r,z}$. We obtain the total
strain energy density

$\displaystyle W(s)=\sum_{\alpha}W^{\alpha}(s),\quad
W^{\alpha}(s)=\int_{-\infty}^{s}\frac{m^{\alpha}(\tau)\,q^{\alpha}(s,\tau)}{\rho(s)}\,\hat{W}^{\alpha}_{n(\tau)}(s)\,\mathrm{d}\tau.$

The constituent-specific mass fraction $\phi^{\alpha}$ at time $s$ is defined
as

$\displaystyle\phi^{\alpha}(s)=\frac{\rho^{\alpha}(s)}{\rho(s)}.$

####  2.2.2 Equilibrated Constrained Mixture Model

The CMMe directly predicts the long-term evolved state in mechanical and
biological equilibrium, i.e., at $s\to\infty$, where production again balances
removal. Note that CMMe can still be used to predict transient G&R phenomena
as long as the process is considered to be quasi-equilibrated, similar to
classical problems in mechanics of quasi-equilibrated problems. Toward this
end, it is assumed that the imposed perturbations to the homeostatic state are
slow relative to the adaptive G&R process. For this purpose, the G&R process
is broken up into several quasi-static load steps where perturbations are
applied in increments. This rate-independent approach does not require
computing hereditary integrals (7), thus drastically lowering the
computational demand. In CMMe, the stimulus function in Equation (5) reduces
to

$\displaystyle\Upsilon^{\alpha}_{h}(\Delta\sigma_{Ih},\Delta\tau_{wh})=1,$
(10)

depending on the equilibrated stimulus functions that can be simplified from
(6) as

$\displaystyle\Delta\sigma_{Ih}=\frac{\sigma_{Ih}}{\sigma_{Io}}-1,\quad\Delta\tau_{wh}=\frac{\tau_{wh}}{\tau_{wo}}-1,$
(11)

comparing the evolved equilibrium $h$ to the original homeostasis $o$, hence
eliminating the time-dependency on $\tau$. All values with index $o$ are
determined at a hyperelastic pre-loading stage ($t=0$), followed by the G&R
stage $h$ to determine the long-term evolved variables in mechanobiological
equilibrium ($t>0$). The equilibrated stimulus function yields

$\displaystyle\Delta\sigma_{Ih}=K_{\tau\sigma}\,\Delta\tau_{wh},\quad
K_{\tau\sigma}=\frac{K^{\alpha}_{\tau}}{K^{\alpha}_{\sigma}}\geq
0~{}\forall\alpha,\quad\eta=\frac{K^{m}}{K^{c}},$ (12)

with a single shear-to-intramural gain ratio $K_{\tau\sigma}$ for all
constituents $\alpha$ and smooth muscle-to-collagen turnover ratio $\eta$.
Note that for $K_{\tau\sigma}=0$ and $K_{\tau\sigma}\to\infty$ the blood
vessel perfectly maintains IMS $\sigma_{Ih}$ and WSS $\tau_{wh}$,
respectively. We obtain the solid Cauchy stress

$\displaystyle\bm{\sigma}^{s}=\bm{\sigma}^{x}-p_{h}\bm{I},\quad\bm{\sigma}^{x}=\sum_{\alpha}^{e,c,m}\phi_{h}^{\alpha}\hat{\bm{\sigma}}_{h}^{\alpha},$
(13)

where $\hat{\bm{\sigma}}^{\alpha}_{h}$ are the equilibrated constituent-
specific Cauchy stresses, $\bm{\sigma}^{x}$ is the mixture stress, and $p_{h}$
is the equilibrated Lagrange multiplier that enforces mechanobiological
equilibrium. We define the IMS invariant throughout this work as

$\displaystyle\sigma_{I}=\frac{\text{tr}\,\bm{\sigma}^{s}}{3},$ (14)

calculated from the solid stress tensor $\bm{\sigma}^{s}$. We can now express
the Lagrange multiplier, using (12), (13), and (14), as

$\displaystyle
p_{h}=\sigma^{x}_{Ih}-\sigma_{Io}\left[1+K_{\tau\sigma}\left(\frac{\tau_{wh}}{\tau_{wo}}-1\right)\right],\quad\sigma^{x}_{Ih}=\frac{\text{tr}\,\bm{\sigma}^{x}}{3}.$
(15)

The invariant $\sigma^{x}_{Ih}$ can be explicitly calculated from the mixture
stress $\bm{\sigma}^{x}$. The original homeostatic IMS $\sigma_{Io}$ is stored
locally from the pre-loading stage. Note that the tangential stiffness matrix
of the CMMe material model only has minor symmetries, no major symmetries.
This can pose challenges for iterative linear solvers and their
preconditioners. We thus use the distributed memory sparse direct linear
solver MUMPS [41] within the PETSc suite [42].

##  3 Fluid-solid-growth coupling

In classic FSI, the solid material behavior is purely determined by an elastic
response to boundary conditions and interface conditions at the fluid-solid
interface $\Gamma_{\text{int}}$. In contrast, the CMMe material behavior is
characterized by a long-term mechanobiological response to mechanical stimuli.
In this section, we first highlight key similarities and differences between
FSI and FSGe problems and then derive a robust coupling scheme. WSS is a
crucial stimulus in the models used in this work. Its local calculation is a
key difference between purely solid G&R models and FSGe:

$\displaystyle\text{\lx@glossaries@gls@link{acronym}{gr}{\leavevmode
G\&R}:}~{}\Delta\tau_{wh}$
$\displaystyle\approx\left(\frac{a_{h}}{a_{o}}\right)^{-3}-1\approx\left[\frac{r_{o}}{a_{o}}\,\lambda_{\theta}-\left(\frac{r_{o}}{a_{o}}-1\right)\,\lambda_{r}\right]^{-3}-1,\quad\lambda_{i}=\sqrt{\mathbf{C}_{h}:\mathbf{e}_{i}\otimes\mathbf{e}_{i}}.$
(16) $\displaystyle\text{\lx@glossaries@gls@link{acronym}{fsge}{\leavevmode
FSGe}:}~{}\Delta\tau_{wh}$
$\displaystyle=\frac{\|\bm{\tau}_{wh}\|}{\|\bm{\tau}_{wo}\|}-1.$ (17)

In the G&R model, local WSS is approximated with the Poiseuille flow-derived
WSS (1). The evolved local inner radius $a_{h}$ of the blood vessel is
approximated locally at each Gauß point of the finite elements. The
approximation uses the radial coordinate $r_{o}$ and inner radius $a_{o}$ in
original homeostasis as well as the evolved local circumferential stretch
$\lambda_{\theta}$ and radial stretch $\lambda_{r}$ with unit vectors
$\mathbf{e}_{i}$ [5]. In the FSGe model, we use the magnitude of the local WSS
vector $\bm{\tau}_{w}$ extracted at the fluid-solid interface from a steady-
state fluid solution. In response to local changes in WSS, the endothelium
alters its production of vasoactive molecules that rapidly diffuse into the
wall and cause vascular smooth muscle cell relaxation or contraction [43].
Instead of solving the diffusion equation, this work assumes that the WSS
stimulus at the interface $\Gamma_{\text{int}}$ effects each radial location
equally within the vessel wall. Thus, the FSGe model is coupled via

$\displaystyle\bm{\sigma}^{\text{s}}\vcentcolon=\bm{\sigma}^{\text{s}}(\bm{\tau}_{w})=\bm{\sigma}^{\text{s}}(\bm{\sigma}^{\text{f}}_{\Gamma_{\text{int}}}),\qquad\bm{x}\in\Omega^{\text{s}},$
(18)

to the solid stress tensor $\bm{\sigma}^{\text{s}}$. Modeling the effect of
vasoactive molecules, at everly location in the vessel the solid stress
$\bm{\sigma}^{\text{s}}$ depends on WSS $\bm{\tau}_{w}$, which, in turn,
depends on the fluid stress tensor
$\bm{\sigma}^{\text{f}}_{\Gamma_{\text{int}}}$ at the interface
$\Gamma_{\text{int}}$. This interaction differs from classic FSI problems,
where fluid and solid are only coupled at the interface. In FSGe, the fluid-
solid interface is coupled to the whole solid volume. Furthermore, in blood
flow, WSS $\tau_{w}$ is typically five orders of magnitude smaller than the
IMS $\sigma_{I}$. Thus, WSS has a negligible mechanical influence in vascular
FSI models, but, depending on the gain ratio $K_{\tau\sigma}$, a large
biological influence in the FSGe model. Also, solid G&R takes place on a time
scale $T^{\text{s}}\gg T^{\text{f}}$ of weeks or months, which is much larger
than the fluid dynamics time scale $T^{\text{f}}$. Thus, we can derive the
interface velocity from the interface displacements $\bm{d}$ as

$\displaystyle\bm{u}_{\Gamma_{\text{int}}}=\frac{\partial\bm{d}}{\partial
t^{\text{s}}}\approx\bm{0},\qquad\bm{x}\in\Gamma_{\text{int}},$ (19)

which differs from the typical FSI condition that non-zero fluid and solid
velocities must match at the interface. The interface condition (19) is
equivalent to a no-slip condition in a rigid-wall fluid simulation. The
traction interface condition

$\displaystyle\bm{\sigma}^{\text{s}}_{\Gamma_{\text{int}}}\cdot\bm{n}=\bm{\sigma}^{\text{f}}_{\Gamma_{\text{int}}}\cdot\bm{n},\qquad\bm{x}\in\Gamma_{\text{int}},$
(20)

is identical to most FSI problems, with interface normal vector $\mathbf{n}$.
While blood pressure is typically assumed to be constant in purely solid G&R
models, the FSGe couples the normal components of local solid and fluid stress
tensors, $\bm{\sigma}^{\text{s}}_{\Gamma_{\text{int}}}$ and
$\bm{\sigma}^{\text{f}}_{\Gamma_{\text{int}}}$, respectively, which spatially
vary according to the solution of the Navier-Stokes equations.

The choice of coupling scheme is crucial for the robustness and accuracy of
the FSGe model. For a comprehensive review of FSI coupling approaches, we
refer to Hou et al. [44]. We refer to Degroote [45] for an extensive
discussion of partitioned fluid-structure coupling approaches and reuse some
of the notations in this section. The FSGe model proposed in this work is
available open-source in the Python coupling code svFSGe [46]. In the
remainder of this section, we briefly review key differences between
monolithic (3.1) and partitioned (3.2) coupling schemes for their use in FSGe.
The monolithic approach simultaneously solves the space- and time-discretized
fluid and solid fields. In contrast, the partitioned approach solves fluid and
solid fields separately and can be implemented with either weak (3.2.1) or
strong (3.2.2) coupling.

###  3.1 Monolithic coupling

Let $\mathbf{f}$ and $\mathbf{s}$ denote the fluid and solid residual,
respectively. The vector $\mathbf{U}$ represents the discrete fluid variables
velocity and pressure; the vector $\mathbf{D}$ represents the discrete solid
displacements. The monolithic FSI problem is then given by

$\begin{cases}\begin{aligned} \mathbf{f}(\mathbf{U},\mathbf{D})&=\mathbf{0}\\\
\mathbf{s}(\mathbf{U},\mathbf{D})&=\mathbf{0}\end{aligned}\end{cases}.$ (21)

This system is commonly solved with the Newton-Raphson method. The update for
fluid and solid increments, $\Delta\mathbf{U}$ and $\Delta\mathbf{D}$,
respectively in iteration $k$ at a load step $t$ becomes

$\begin{bmatrix}\partial_{\mathbf{U}}\,\mathbf{f}&\partial_{\mathbf{D}}\,\mathbf{f}\\\
\partial_{\mathbf{U}}\,\mathbf{s}&\partial_{\mathbf{D}}\,\mathbf{s}\end{bmatrix}_{t,k}\begin{bmatrix}\Delta\mathbf{U}\\\
\Delta\mathbf{D}\end{bmatrix}_{t,k+1}=-\begin{bmatrix}\mathbf{f}\\\
\mathbf{s}\end{bmatrix}_{t,k},$ (22)

where $\partial_{\mathbf{U}}\,\mathbf{f}$ denotes the partial derivative of
$\mathbf{f}$ with respect to $\mathbf{U}$. All blocks in Equation (22) are
treated in the same mathematical framework and implemented in svFSIplus. For
FSI problems, the advantages of monolithic coupling can be faster convergence
and higher accuracy than the partitioned approach. However, it may require
more computational resources and a more intrusive software implementation.
Independent of these challenges, we found a monolithic coupling approach to be
intractable for an FSGe problem. Solving Equation (22) requires solving both
fluid and solid fields with the same time step size, typically
$T_{\text{FSI}}\sim 1\,\text{ms}$ or below in cardiovascular FSI problems
[47]. However, as discussed in Section 3, solid G&R takes place over a much
longer time scale than blood flow fluid dynamics. In fact, the CMMe returns
the mechanobiological equilibrated G&R response after an infinitely long
period $s\to\infty$ at every load step $t$. This limit renders solving
Equation 22 with a uniform fluid-solid time step size intractable, and we,
therefore, chose a partitioned scheme for FSGe. Note that we refer to $t$
throughout this work as a quasi-static G&R load step opposed to a time step in
a dynamic problem. All load steps represent a long-term equilibrium at
$s\to\infty$, typically for different insult severities.

###  3.2 Partitioned coupling

In partitioned FSI, fluid field $\mathcal{F}$ and solid field $\mathcal{S}$
solve for interface traction $\mathbf{u}$ and interface displacements
$\mathbf{d}$, respectively so that

$\displaystyle\mathbf{u}$ $\displaystyle=\mathcal{F}(\mathbf{d}),$ (23)
$\displaystyle\mathbf{d}$ $\displaystyle=\mathcal{S}(\mathbf{u}).$ (24)

In the fluid field $\mathcal{F}$, in addition to solving blood flow dynamics
(Section 2.1), we apply a deformation to the computational grid of
$\Omega^{f}$ to avoid excessive fluid mesh distortion. This is similar to the
arbitrary Lagrangian-Eulerian (ALE) formulation typically employed in FSI
problems where an arbitrary grid velocity is applied. However, due to (19),
the mesh deformation is quasi-static in FSGe. In the FSGe model, we apply an
interface mesh deformation $\mathbf{d}_{\text{int}}$ and solve a linear
elasticity problem to solve for the fluid mesh deformation

$\mathbf{K}\cdot\mathbf{d}_{f}=\mathbf{0},\quad\text{with~{}}\mathbf{d}_{f}=\mathbf{d}_{\text{int}}\text{~{}on~{}}\Gamma_{\text{int}},$
(25)

with a stiffness matrix $\mathbf{K}$. Inspired by Degroote [45], we outline
the steps within the fluid and solid solvers in Algorithms 1 and 2,
respectively. Partitioned FSI offers the advantage of reusing and combining
existing fluid and solid solvers with solution techniques targeted to the
respective domain. Note that while independent solvers could be used for
$\mathcal{F}$ and $\mathcal{S}$, we use our open-source multi-physics solver
svFSIplus for both [25]. Crucially, this allows us to natively split both time
scales arising in the FSGe problem. We thus focus on a partitioned coupling
approach for the remainder of this work.

1:Apply interface displacement $\mathbf{d}$

2:Deform the fluid domain (25)

3:Solve for fluid variables $\mathbf{U}$

4:Extract interface traction $\mathbf{u}$

Algorithm 1 $\mathcal{F}$: Steady-state fluid

1:Apply interface traction $\mathbf{u}$

2:Propagate WSS through solid domain

3:Solve for solid variables $\mathbf{D}$

4:Extract interface displacement $\mathbf{d}$

Algorithm 2 $\mathcal{S}$: Quasi-static solid

Note that in comparison to a classic FSI problem, in FSGe we additionally need
to radially propagate the effects of the WSS at the interface throughout the
thickness of the blood vessel in the solid domain (line 2 in Algorithm 2) as
discussed in Section 2.2. In a cylindrical reference domain with a structured
solid grid, this can be easily facilitated with a cylindrical coordinate
system. A radial coordinate can be defined in patient-specific geometries,
e.g., by solving a diffusion problem commonly used in cardiac geometries [48].

####  3.2.1 Weak coupling

The group of partitioned FSI coupling schemes can be further split into
implicit (or strong), semi-implicit, and explicit (or weak) schemes [49]. As
Section 1 outlines, most previous FSG approaches used weak coupling. Algorithm
3 illustrates the most basic case of a serial staggered weak coupling scheme.
Here, the fluid solution is updated once in every load step $t$ based on the
solid solution $\mathbf{d}_{t-1}$ from the previous load step. The updated
fluid solution is then input to update the solid solution $\mathbf{d}_{t}$.
For FSGe, this weak coupling method is unstable for $K_{\tau\sigma}>0$ and is
thus not further pursued in this work. Although weakly coupled FSI problems
can sometimes be stabilized by reducing the time step size, this strategy is
not possible in FSGe. While an applied stimulus can be spread over more load
steps, the time-independent solid CMMe model always yields a long-term evolved
equilibrium at $s\to\infty$.

1:$\mathbf{d}_{0}=\mathbf{0},~{}t=0$

2:do

3: $t++$

4: $\mathbf{d}_{t}=\mathcal{S}\circ\mathcal{F}(\mathbf{d}_{t-1})$

5:while $t\leq t_{\text{max}}$

Algorithm 3 Weak coupling scheme

####  3.2.2 Strong coupling

Strong coupling schemes enforce an equilibrium between fluid and solid fields
in every load step, which requires an additional iteration loop. Strong
coupling is essential to solve the FSGe problem when there is a strong
interaction between fluid and solid, which is the case for $K_{\tau\sigma}>0$,
that is, when the solid G&R material depends on local WSS. We control the
convergence towards an equilibrium in coupling iterations $k$ at every load
step $t$ with the displacement residual

$\mathbf{r}_{t,k}=\tilde{\mathbf{d}}_{t,k}-\mathbf{d}_{t,k},$ (26)

In the following, $\tilde{\mathbf{d}}_{t,k}$ denotes the output of a
subsequent fluid-solid solve

$\tilde{\mathbf{d}}_{t,k}=\mathcal{S}\circ\mathcal{F}(\mathbf{d}_{t,k}),$ (27)

in contrast to the input $\mathbf{d}_{t,k}$. Equilibrium between fluid and
solid is achieved when the relative displacement error norm

$\frac{\|\mathbf{r}_{t,k}\|_{2}}{\|\mathbf{d}_{t,k}\|_{2}}<\epsilon_{0},$ (28)

is smaller than a given tolerance $\epsilon_{0}$. The most basic form of
strong coupling is the Gauß-Seidel iteration scheme depicted in line 8 of
Algorithm 4. Here, fluid mesh deformation, fluid, and solid are executed
sequentially, and we use the output from the previous coupling iteration as
the input for the current iteration:

$\mathbf{d}_{t,k}=\tilde{\mathbf{d}}_{t,k-1}.$ (29)

Gauß-Seidel schemes typically converge slowly or might not converge at all
[45]. A straightforward approach to stabilize the Gauß-Seidel scheme is to
introduce static relaxation

$\mathbf{d}_{t,k}=\omega\,\tilde{\mathbf{d}}_{t,k-1}+(1-\omega)\,\mathbf{d}_{t,k-1},\quad
0<\omega\leq 1,$ (30)

with a static relaxation parameter $\omega$, which reduces to the Gauß-Seidel
scheme (29) for $\omega=1$. While usually improving numerical stability for
$0<\omega<1$, this coupling scheme typically converges slowly and highly
depends on the choice of the relaxation parameter $\omega$. A much more
efficient scheme is obtained by dynamically choosing $\omega_{k}$ in each
iteration $k$ through Aitken relaxation [50, 51]:

$\omega_{k}=-\omega_{k-1}\frac{\mathbf{r}_{k-1}\cdot\Delta\mathbf{r}_{k}}{\|\Delta\mathbf{r}_{k}\|_{2}},\quad\text{with}~{}\Delta\mathbf{r}_{k}=\mathbf{r}_{k}-\mathbf{r}_{k-1}.$
(31)

While easy to implement, Aitken relaxation can use only information from one
previous coupling iteration and is thus limited in its potential to introduce
stabilization to the coupling scheme. The interface quasi-Newton with inverse
Jacobian from a least-squares model (IQN-ILS) can use information from up to
$q$ previous coupling iterations [52, 45]. It is based on applying the Newton-
Raphson method to solve the equation $\mathbf{r}_{t}=\mathbf{0}$, which yields

$\mathbf{d}_{t,k}=\tilde{\mathbf{d}}_{t,k-1}-\left(\left.\frac{\partial\mathbf{r}}{\partial\mathbf{d}}\right|_{t,k-1}\right)^{-1}\mathbf{r}_{k-1}.$
(32)

However, the Jacobian matrix $\partial\mathbf{r}/\partial\mathbf{d}$ is not
known explicitly, and an approximation with finite differences would require
several additional solver evaluations. In IQN-ILS, the Jacobian is
approximated from a linear combination of known previous displacements,
computing the input displacements as

$\displaystyle\mathbf{d}_{t,k}=\tilde{\mathbf{d}}_{t,k-1}+\mathbf{W}\mathbf{c},\quad\text{with}~{}\mathbf{W}=[\Delta\tilde{\mathbf{d}}_{k-1},\Delta\tilde{\mathbf{d}}_{k-2},\dots],\quad\Delta\tilde{\mathbf{d}}_{k-1}=\tilde{\mathbf{d}}_{k}-\tilde{\mathbf{d}}_{k-1},$
(33)

with coefficient vector $\mathbf{c}$ which we obtain by solving the least
squares problem

$\displaystyle\min_{\mathbf{c}}\|\mathbf{V}\mathbf{c}+\mathbf{r}_{k}\|_{2},\quad\text{with}~{}\mathbf{V}=[\Delta\mathbf{r}_{k-1},\Delta\mathbf{r}_{k-2},\dots].$
(34)

This minimization problem is solved using the QR-decomposition
$\mathbf{Q}\mathbf{R}=\mathbf{V}$. The matrices $\mathbf{V}$ and $\mathbf{W}$
are updated in every coupling iteration $k$ and truncated to contain a maximum
of $q$ entries. We reuse information from previous load steps in $\mathbf{V}$
and $\mathbf{W}$ but never calculate increments across load steps.
Furthermore, filtering ensures that linearly dependent columns $i_{\epsilon}$
are excluded from $\mathbf{V}$ and $\mathbf{W}$ when
$\mathbf{R}_{ii}<\epsilon_{QR}$. The IQN-ILS approach is summarized in
Algorithm 6.

1:$\mathbf{V}=\emptyset,~{}\mathbf{W}=\emptyset,~{}t=0$

2:do

3: $k=0$

4: $\mathbf{d}_{t,0}$ from Algorithm 5

5: do

6: $k++$

7: if $t=0$ then

8: $\mathbf{d}_{t,k}=\tilde{\mathbf{d}}_{t,k-1}$

9: else if $(k=1)$ or ($t=1$ and $k\leq 5$) then

10:
$\mathbf{d}_{t,k}=\omega\,\tilde{\mathbf{d}}_{t,k-1}+(1-\omega)\,\mathbf{d}_{t,k-1}$

11: else

12: $\mathbf{d}_{t,k}$ from Algorithm 6

13: end if

14: $\tilde{\mathbf{d}}_{t,k}=\mathcal{S}\circ\mathcal{F}(\mathbf{d}_{t,k})$

15: $\mathbf{r}_{t,k}=\tilde{\mathbf{d}}_{t,k}-\mathbf{d}_{t,k}$

16: while $\|\mathbf{r}_{t,k}\|_{2}\,/\,\|\mathbf{d}_{t,k}\|_{2}>\epsilon_{0}$

17: $t++$

18:while $t\leq t_{\text{max}}$

Algorithm 4 Strongly coupled FSGe

1:if $t=0$ then

2: $\mathbf{d}_{t,\,0}=\mathbf{d}_{0}$

3:else

4: $\mathbf{d}_{t,\,0}=2\mathbf{d}_{1}-\mathbf{d}_{0}$

5:end if

Algorithm 5 Polynomial predictor

1:$\mathbf{V}~{}=(\mathbf{r}_{t,k}\,-\mathbf{r}_{t,k-1},\mathbf{V})$

2:$\mathbf{W}=(\tilde{\mathbf{d}}_{t,k}-\tilde{\mathbf{d}}_{t,k-1},\mathbf{W})$

3:$\mathbf{V}=\mathbf{V}\setminus\mathbf{V}_{i>q},~{}\mathbf{W}=\mathbf{W}\setminus\mathbf{W}_{i>q}$

4:do

5: $\mathbf{Q}\mathbf{R}=\mathbf{V}$

6: $i_{\epsilon}=\arg_{i}\left(|\mathbf{R}_{ii}|<\epsilon_{\text{QR}}\right)$

7:
$\mathbf{V}=\mathbf{V}\setminus\mathbf{V}_{i_{\epsilon}},~{}\mathbf{W}=\mathbf{W}\setminus\mathbf{W}_{i_{\epsilon}}$

8:while $i_{\epsilon}\neq\emptyset$

9:$\mathbf{c}=\mathbf{R}^{-1}\,\mathbf{b}$ with
$\mathbf{b}=-\mathbf{R}^{-\text{T}}\,(\mathbf{V}\mathbf{r}_{t,k})$

10:$\mathbf{d}_{t,k}=\tilde{\mathbf{d}}_{t,k-1}+\mathbf{W}\mathbf{c}$

Algorithm 6 IQN-ILS

The strongly coupled FSGe method is outlined in Algorithm 4. We use a Gauß-
Seidel scheme exclusively at $t=0$ for the pre-loading stage in CMMe [5] in
line 8. Here, the displacements are negligibly small, and there is only weak
coupling as pre-loading does not depend on local WSS. The G&R stage in CMMe at
$t>0$ requires the more robust IQN-ILS coupling approach in line 12. Note that
the performance of the algorithm also crucially depends on an accurate
prediction of displacements $\mathbf{d}_{t,0}$ at the beginning of each load
step. If a previous load step is available, we use the linear predictor in
Algorithm 5. Further, note that at the beginning of each coupling iteration
$k=1$, we obtain input displacements $\mathbf{d}_{t,k}$ from a static
relaxation step (30). Similarly, we use static relaxation for the first five
coupling iterations at $t=1$, i.e., the first G&R load step, to collect an
initial number of iterations for IQN-ILS.

##  4 Illustrative computational examples

This section compares our equilibrated fluid-solid-growth (FSGe) model to a
purely solid growth and remodeling (G&R) model. As an illustrative example, we
focus on the growth of an aortic aneurysm in a mouse. Aortic aneurysms have
been studied in prior work using the equilibrated constrained mixture model
(CMMe) [5, 33] and a prior fluid-solid-growth coupling method [8]. Given the
inherently large deformation in aortic aneurysms, they locally alter blood
flow [53] and are thus a good example to demonstrate differences between solid
G&R and FSGe and a good starting step for more complex in vivo flow patterns
in future research.

Solid G&R parameters [5] |  |
---|---|---
Constituents | $e,m,c$ | elastin, smooth muscle, collagen
Inner radius, thickness, length | $a_{o},h_{o},l_{o}$ | 0.647 mm, 0.04 mm, 15 mm
Mass fractions | $\phi^{e}_{o},\phi^{m}_{o},\phi^{c}_{o}$ | 0.34, 0.33, 0.33
Collagen orientation fractions | $\beta_{\theta},\beta_{z},\beta_{d}$ | 0.056, 0.067, 0.877
Diagonal collagen orientation | $\alpha_{0}$ | 29.9°
Elastin stiffness | $c^{e}_{o}$ | 89.71 kPa
Smooth muscle parameters | $c^{m}_{1},c^{m}_{2}$ | 261.4 kPa, 0.24
Collagen parameters | $c^{c}_{1},c^{c}_{2}$ | 234.9 kPa, 4.08
Deposition stretches | $G^{e}_{\theta},G^{e}_{z},G^{e}_{r},G^{m},G^{c}$ |  1.90, 1.62, $1/(G^{e}_{\theta}G^{e}_{z})$, 1.20, 1.25
Smooth muscle-to-collagen turnover ratio | $\eta$ | 1.0
Shear-to-intramural gain ratio | $K_{\tau\sigma,o}$ | {0.0, 0.2, 0.4, 0.6, 0.8, 1.0}
External tissue support stiffness | $k$ | 2.0 kPa/mm
Aneurysm parameters |  |
Circumferential extent | $\theta_{od}$ | $0.55\,\pi$
Circumferential decay | $v_{\theta}$ | 6
Axial extent | $z_{od}$ | $l_{o}/4$
Axial decay | $v_{z}$ | 2
Maximum elastin degradation | $\varphi^{e}_{hm}$ | 0.7
Fluid parameters |  |
Maximum inflow velocity | $u_{\text{in}}$ | 1000 mm/s
Outlet pressure [5] | $p_{\text{out}}$ | 104.9 mmHg
Dynamic viscosity | $\mu$ |  $4.0\cdot 10^{-6}$ kg/mm/s
Density | $\rho$ |  $1.06\cdot 10^{-6}$ kg/mm3
Numerical parameters |  |
Number of elements | $n_{\theta}\times n_{r}\times n_{z}$ | $64\times 1\times 40$
Number of G&R load steps [5] | $t_{\text{max}}$ | 10
Coupling tolerance | $\epsilon_{0}$ | $10^{-3}$
Static relaxation parameter | $\omega$ | 0.1
IQN-ILS filtering tolerance | $\epsilon_{\text{IQN-ILS}}$ | $10^{-1}$
IQN-ILS number of old iterations | $q$ | 20
Table 2: Parameters used in all computational examples in this work.

For an overview of all parameters, see Table 2, where the solid G&R parameters
are taken from Latorre and Humphrey [5]. All parameters in the G&R model are
identical to the FSGe model, with the FSGe having additional fluid and
coupling parameters. The inflow velocity $u_{\text{in}}$ was chosen as an
upper bound for systolic velocities observed in mouse aortas [54]. Aneurysmal
growth is associated with elastin degradation, modeled here by decreasing
elastin stiffness [55], i.e., $c^{e}_{h}\to 0$, and loss of mechanosensation,
i.e., $K_{\tau\sigma,h}\to 0$ [56]. We apply a spatial insult profile to
trigger asymmetric aneurysmal growth,

$\displaystyle c^{e}_{h}(\theta,z,t)$
$\displaystyle=c^{e}_{o}\,[1-\varphi^{e}_{hm}\,f(\theta,z,t)],$ (35)
$\displaystyle K_{\tau\sigma,h}(\theta,z,t)$
$\displaystyle=K_{\tau\sigma,o}\,[1-f(\theta,z,t)],$ (36)

reducing long-term evolved elastin stiffness $c^{e}_{h}$ and gain ratio
$K_{\tau\sigma,h}$ with the spatial and temporal factor $f$ and maximum
elastin degradation factor $\varphi^{e}_{hm}$. In our computational results,
we sampled aneurysmal growth for different original homeostatic gain ratios
$K_{\tau\sigma,o}\in[0,1]$ as it proved crucial for the G&R vs. FSGe
comparison. The circumferential, axial, and temporal factors are

$\displaystyle f(\theta,z,t)=f_{\theta}\,f_{z}\,f_{t},\quad
f_{\theta}=\exp\left(-\left|\frac{\theta_{o}-\pi}{\theta_{od}}\right|^{v_{\theta}}\right),\quad
f_{z}=\exp\left(-\left|\frac{z_{o}}{z_{od}}-\frac{1}{2}\right|^{v_{z}}\right),\quad
f_{t}=\frac{\text{tanh}\,(2t/t_{\text{max}})}{\text{tanh}\,(2)},\quad t>0,$
(37)

where the insult is ramped up over the load steps $t_{\text{max}}$. Recall
that load step $t=0$ corresponds to the CMMe pre-loading stage, whereas $t>0$
corresponds to G&R. We modulate the temporal factor $f_{t}$ with the
hyperbolic tangent to yield approximately equal displacement increments in all
load steps; this enables using a linear predictor that initializes each load
step with a solution close to the converged one. The spatial insult profile
$f_{\theta}\,f_{z}$ is visualized in Figure 3 (left). We introduce an
aneurysmal region locally at the top of the blood vessel, whereas the ends and
the bottom maintain healthy parameters.

Figure 3: Mesh (left) with fluid domain $\mathcal{F}$ (yellow) and solid
domain $\mathcal{S}$, colored by the spatial insult profile
$f_{\theta}\,f_{z}$ (white healthy, purple initiating insult profile). Plot
locations (right), top, right, bottom, left, with a cylindrical coordinate
system $(r,\theta,z)$.

We generated the cylindrical structured mesh in Figure 3 (left) with
hexahedral linear displacement-based finite elements that are refined at the
location of the aneurysm in the middle. Further, the mesh has a conforming
interface between fluid and solid domains. The global $(\theta,r,z)$
coordinate system and four different plot locations within the blood vessel
are shown in Figure 3 (right). All results in this work are extracted at the
fluid-solid interface, i.e., at the endothelium. We show plots over the axial
length of the vessel and point-wise results at the apex of the aneurysm.

###  4.1 Computational performance

Figure 4 shows the cumulative number of FSGe coupling iterations in each load
step for varying gain ratios $K_{\tau\sigma}\in[0,1]$ required to reach the
coupling tolerance $\epsilon_{0}$. Load step $t=0$ denotes the CMMe pre-
loading step, which takes only three iterations to converge in all simulations
since the displacements $\mathbf{d}$ are close to zero during pre-loading. Our
IQN-ILS coupling scheme from Algorithm 6 relies on solutions from previous
coupling iterations, which are not yet available at the beginning of load step
$t=1$. In addition, since there is only one prior converged load step, we
cannot use extrapolation for the predictor. Thus, load step $t=1$ requires
more coupling iterations than other load steps in all simulations.
Importantly, there is a clear trend that more coupling iterations are required
with higher $K_{\tau\sigma,o}$. The average number of coupling iterations for
load steps 2 to 10 ranges from 2.0 for $K_{\tau\sigma,o}=0.0$ to 7.6 for
$K_{\tau\sigma,o}=1.0$. For the longest simulation $K_{\tau\sigma,o}=1.0$ with
84 coupling iterations, the total wall times on an Apple M2 Max CPU were 2.7
min (one core) for linear-elastic mesh deformation and 52 min (six cores) for
fluid dynamics in the fluid domain $\mathcal{F}$ and 1.7 min (one core) for
G&R with CMMe in the solid domain $\mathcal{S}$. This renders the steady-state
fluid simulation the computational bottleneck of FSGe.

Figure 4: Cumulative number of FSGe coupling iterations per load step for
different gain ratios $K_{\tau\sigma,o}$.

###  4.2 Solid mechanics

This section compares the solid mechanics between solid-only G&R simulations
and coupled FSGe simulations, both using CMMe as the solid model. Figure 5
shows the displacement fields for G&R (left column) and FSGe (right column)
for increasing original homeostatic gain ratios $K_{\tau\sigma,o}\in[0,1]$
from top to bottom. The red color corresponds to the evolved gain ratio
$K_{\tau\sigma,h}$ after applying the aneurysm insult. Note that the evolved
blood vessel retains its healthy original homeostatic parameters outside the
top part of the aneurysm where the aneurysm insult was applied. Generally, a
smaller original homeostatic gain ratio $K_{\tau\sigma,o}$ leads to larger
aneurysm growth. Comparing the top row with uniform $K_{\tau\sigma,o}=0.0$,
there is little difference between G&R and FSGe results. However, the
aneurysms are smaller with FSGe than with G&R for $K_{\tau\sigma,o}>0$. The
differences become more pronounced for increasing gain ratio
$K_{\tau\sigma,o}\to 1.0$, especially in the bottom aneurysmal region with
intact elastin stiffness and mechanosensation. Interestingly, the radial
deformation in this region is close to zero for $K_{\tau\sigma,o}=0.2$ and
becomes negative, i.e., inward, for $K_{\tau\sigma,o}>0.2$.

G&R $K_{\tau\sigma,o}=0.0$ FSGe

G&R $K_{\tau\sigma,o}=0.2$ FSGe

G&R $K_{\tau\sigma,o}=0.4$ FSGe

G&R $K_{\tau\sigma,o}=0.6$ FSGe

G&R $K_{\tau\sigma,o}=0.8$ FSGe

G&R $K_{\tau\sigma,o}=1.0$ FSGe

Flow direction $\Longrightarrow$

Figure 5: Long-term evolved solid domains in G&R (left) vs. FSGe (right) for
different original gain ratios $K_{\tau\sigma,o}\in[0,1]$ (from top to
bottom). The geometries are colored by the evolved gain ratio
$K_{\tau\sigma,h}$.

Figure 6 shows the local thickness over the axial length, $z$, of the blood
vessel for $K_{\tau\sigma,o}\in\\{0.0,1.0\\}$, i.e., for the top and bottom
rows in Figure 5. Again, G&R and FSGe models were very similar for
$K_{\tau\sigma,o}=0.0$ (top row). While all G&R results are symmetric, we
observed minor thickening of about 3% downstream at $z=11$ compared to
upstream at $z=4$ for the FSGe model. With $K_{\tau\sigma,o}=1.0$, there was a
pronounced axial asymmetry in the FSGe model (bottom right) compared to the
G&R model (bottom left). Comparing the same locations, we observed 12%
downstream thinning at the top of the aneurysm (blue) and 24% thickening at
the bottom (red). Most importantly, the thickness at the bottom of the
aneurysm (red) reaches an evolved equilibrial thickness of more than three
times the original homeostatic value.

Flow direction $\Longrightarrow$

Figure 6: Thickness according to locations in Fig. 3: top (blue), bottom
(red), sides (orange/green).

The local displacements $d_{\\{\theta,r,z\\}}$ in the cylindrical coordinate
system are shown in Figure 7. The comparison is shown for the highest sampled
gain ratio $K_{\tau\sigma,o}=1.0$ (bottom row of Figure 5), which showed the
most pronounced differences between the two models. Due to the symmetry of the
aneurysm, there were no circumferential displacements (top row) at the top and
bottom (blue/red dotted) of the aneurysm in both models. However, at the sides
of the aneurysm (orange, green), there was double the circumferential
displacement, i.e. rotation, towards the bottom with FSGe than with G&R. The
radial displacements (middle row) at the bottom (red) and the sides
(orange/green dotted) of the aneurysm became negative in the FSGe model,
indicating inward growth in these regions. Given the symmetry of the G&R
model, axial displacements (bottom row) are symmetric around $z=7.5$ mm but
asymmetric in the FSGe model.

Flow direction $\Longrightarrow$

Figure 7: Displacements $d_{\theta},d_{r},d_{z}$ (top to bottom) at the
endothelium according to locations in Fig. 3: top (blue), right (orange),
bottom (red), left (green) for varying original homeostatic gain ratios
$K_{\tau\sigma,o}$.

The local mechanobiology in both G&R and FSGe for $K_{\tau\sigma,o}=1.0$ is
visualized in Figure 8. The prescribed insult profile (top row) of
mechanosensation loss towards the top of the aneurysm is identical in both
models. However, the same prescribed gain ratios
$K_{\tau\sigma,h}=\Delta\sigma_{I}/\Delta\tau_{w}$ are achieved with different
levels of intramural stimuli $\Delta\sigma_{I}$ (middle row) and wall shear
stress (WSS) stimuli $\Delta\tau_{w}$ (bottom row). Since $K_{\tau\sigma,h}=0$
at the top peak of the aneurysm (blue curves at $z=7.5$), it follows for both
models that $\Delta\sigma_{I}=0$ at this location. Furthermore, recall that we
maintain intact mechanosensation at the bottom of the aneurysm (red curves at
$z=7.5$). While in the G&R model, the mechanobiologically equilibrated
solution yields $K_{\tau\sigma,h}=1.0\approx(-0.1)/(-0.1)$, the FSGe model
yields $K_{\tau\sigma,h}=1.0\approx(-0.45)/(-0.45)$. This means that in the
G&R model’s evolved equilibrium, both intramural stimulus $\Delta\sigma_{I}$
and WSS stimulus $\Delta\tau_{w}$ are reduced by 10% compared to the original
homeostatic state. However, in the FSGe model, these were reduced by 45%, thus
arriving at a significantly different long-term evolved equilibrium that’s
characterized by the observed differences in local thickness (Figure 6) and
displacements (Figure 7).

Flow direction $\Longrightarrow$

Figure 8: Prescribed evolved equilibrium gain ratio
$K_{\tau\sigma,h}=\Delta\sigma_{I}/\Delta\tau_{w}$ (top) and resulting
intramural $\Delta\sigma_{I}$ (middle) and WSS stimuli $\Delta\tau_{w}$
(bottom) in G&R (left) and FSGe (right) models. Locations according to Figure
3: top (blue), bottom (red), sides (orange/green).

We visualize the evolved local collagen mass $\phi^{c}_{h}J_{h}$ in Figure 9
at the aneurysm mid-section over the sampled original gain ratios
$K_{\tau\sigma,o}$. Overall, collagen is deposited most at the top of the
aneurysm (blue) and, to a lesser extent, at the healthy bottom (red). Notably,
with increasing gain ratio $K_{\tau\sigma,o}\to 1$ in the FSGe model (right),
collagen is degraded at the bottom (red) below its initial mass fraction of
$\phi^{c}_{o}=0.33$ (black). The degradation of collagen at the bottom (red)
of the aneurysm is associated with inward growth around this region. This is
visible in Figure 10, where radial displacements at the sides (green/orange)
and bottom (red) become negative for $K_{\tau\sigma}\geq 0.4$. These
mechanisms are much less evident in the G&R model, where there is little
dependence on the original gain ratio $K_{\tau\sigma,o}$ and no collagen
degradation.

Figure 9: Evolved referential collagen mass fraction (mass of collagen per
unit reference volume) $\phi^{c}_{h}J_{h}$ at the endothelium at $z=7.5$ mm
for varying original gain ratios $K_{\tau\sigma,o}$. The black line indicates
the reference collagen mass of $\phi^{c}_{o}=0.33$. Locations according to
Fig. 3: top (blue), bottom (red), sides (orange/green).
Figure 10: Radial displacement $d_{r}$ at the endothelium at $z=7.5$ mm for
varying original gain ratios $K_{\tau\sigma,o}$. Displacements $d_{r}<0$
(black line) indicate an inward growth of the aneurysm. Locations according to
Fig. 3: top (blue), bottom (red), sides (orange/green).

###  4.3 Fluid dynamics

We compared solid mechanics in G&R and FSGe models in Section 4.2. This
section compares the fluid dynamics in both models. Figure 11 shows the WSS
stimulus $\Delta\tau_{w}$ for G&R (left column) and FSGe (right column)
simulations. The plot locations are identical to Figure 5. Whereas the WSS
stimulus in the FSGe model stems from solving the incompressible Navier-Stokes
equations (3), it is approximated from the Poiseuille solution (1) in the G&R
model. Purple regions with a value of $\Delta\tau_{w}\approx 0$ correspond to
regions that maintain WSS during the development of the aneurysm. As expected,
a higher gain ratio $K_{\tau\sigma,o}\to 1$ (from top to bottom) leads to more
regions with $\Delta\tau_{w}\to 0$ due to the increasing weight of WSS in
finding the mechanobiological equilibrium. In G&R, $\Delta\tau_{w}$ is
strictly coupled to the vessel’s local change in endothelial radius, leading
to strong circumferential variations in the apex of the aneurysm. In FSGe,
$\Delta\tau_{w}$ is circumferentially more uniform around the deformed
centerline of the blood vessel.

G&R $K_{\tau\sigma,o}=0.0$ FSGe

G&R $K_{\tau\sigma,o}=0.2$ FSGe

G&R $K_{\tau\sigma,o}=0.4$ FSGe

G&R $K_{\tau\sigma,o}=0.6$ FSGe

G&R $K_{\tau\sigma,o}=0.8$ FSGe

G&R $K_{\tau\sigma,o}=1.0$ FSGe

Flow direction $\Longrightarrow$

Figure 11: WSS stimulus $\Delta\tau_{w}$ in G&R (left) vs. FSGe (right) for
different original gain ratios $K_{\tau\sigma,o}\in[0,1]$ (from top to
bottom).

Figure 12 shows velocity (left) and pressure (right) for the vessel’s
longitudinal midsection for different gain ratios $K_{\tau\sigma,o}$ in the
FSGe. Note that, due to the absence of local fluid dynamics in the G&R model,
we cannot compare the models. The geometries correspond to the right column in
Figure 5. Recalling the fluid boundary conditions from Section 4, it is
evident that all flows have identical (parabolic) inflow conditions on the
very left and identical outflow pressures of 104.9 mmHg on the very right.
While there is little variation in pressure fields for the sampled gain ratios
$K_{\tau\sigma,o}$, there are some notable differences in the velocity field.
There is a small recirculation zone visible at the top apex of the aneurysm.
Due to the incompressibility of blood, for $K_{\tau\sigma,o}=0.0$ (top), the
doubling in radius inside the aneurysm drastically decreases local velocity.
This effect is less pronounced as the maximum radius decreases with increasing
gain ratio $K_{\tau\sigma,o}$. Interestingly, as the gain ratio approaches
$K_{\tau\sigma,o}\to 1$, the velocity is deflected upward into the top part
upstream of the aneurysm.

FSGe, $K_{\tau\sigma,o}=0.0$

FSGe, $K_{\tau\sigma,o}=0.2$

FSGe, $K_{\tau\sigma,o}=0.4$

FSGe, $K_{\tau\sigma,o}=0.6$

FSGe, $K_{\tau\sigma,o}=0.8$

FSGe, $K_{\tau\sigma,o}=1.0$

Flow direction $\Longrightarrow$

Figure 12: Fluid domains in FSGe models with velocity and flow profiles (left)
and pressure (right) with vessel wall (black) for original homeostatic gain
ratios $K_{\tau\sigma,o}\in[0,1]$ (from top to bottom).

Figure 13 shows fluid velocity $\bm{u}$ (left) and pressure $p$ (right) along
the deformed centerline of the vessel for all gain ratios in Figure12. Again,
prescribed boundary conditions are evident for the inlet velocity at $z=0$ mm
and outlet pressure at $z=15$ mm. The centerline velocity also coincides with
the region of maximum velocity. The drop in velocity inside the aneurysm and
the downstream acceleration distal to the aneurysm are visible. We also
observe a linear drop in pressure along $z$ upstream and downstream of the
aneurysm and a pressure plateau in the aneurysmal region. While there’s a more
than 60% fluctuation in centerline velocity throughout all $K_{\tau\sigma,o}$
samples, the pressure changes only about 3.5% compared to the overall mean
pressure.

Figure 13: Velocity magnitude (left) and pressure (right) along the deformed
vessel centerline from $K_{\tau\sigma,o}=0.0$ (light) to
$K_{\tau\sigma,o}=1.0$ (dark).

##  5 Discussion

We presented a novel equilibrated fluid-solid-growth (FSGe) method, where we
strongly couple the three-dimensional (3D) Navier-Stokes equations and the 3D
equilibrated constrained mixture model (CMMe). We demonstrated crucial
differences in the solution compared to a solid-only growth and remodeling
(G&R) model in illustrative computational examples of an asymmetric aortic
aneurysm. We identified the original homeostatic gain ratio $K_{\tau\sigma,o}$
as a crucial factor determining the degree of interaction between fluid and
solid growth in the long-term evolution. In the remainder of this section, we
discuss the results from Section 4, limitations, and future perspectives. We
first review computational aspects of the FSGe method proposed in this work:
its equilibrated nature (5.1), its strong coupling (5.2), stability of CMMe
(5.3), and original homeostasis (5.4). We further review mechanobiological
aspects of FSGe, namely wall shear stress (WSS) influence (5.5) and the
directionality and pulsatility of blood flow (5.6).

###  5.1 Equilibrated vs. transient growth and remodeling

Latorre and Humphrey [5] introduced the CMMe as a fast, rate-independent
version of the constrained mixture model (CMM). It directly yields the long-
term equilibrated G&R response, i.e., at $s\to\infty$. The computational cost
of CMMe is comparable to a solid mechanics FEM solution with a hyperelastic
material. Together with our implicit coupling scheme, this makes the FSGe
highly efficient. On the other hand, it may be critical to model transient G&R
in cases where the insult is applied slowly compared to the G&R response. For
example, tissue-engineered vascular grafts (TEVGs) in Fontan patients can
develop a severe stenosis, which sometimes spontaneously resolves [14, 15].
Here, the time-resolved development of G&R is crucial to studying the
mechanisms behind the stenosis formation and resolution [13], which can be
modeled with FSGe [12]. However, we can only study phenomena that will yield a
long-term mechanobiological equilibrium. Yet, many clinically interesting
applications involve G&R processes that do not necessarily yield an evolved
equilibrial state [57]. For example, theories have been developed to quantify
the stability or instability of G&R processes such as aneurysm formation [58,
59]. Thus, transient G&R models like the one developed by Schwarz et al. [13]
will be necessary to study fluid-solid-growth interaction (FSG) in unstable
cases.

###  5.2 Strong vs. weak fluid-solid-growth coupling

We proposed a strongly coupled partitioned fluid-structure interaction model
where the solid G&R model is given by the CMMe. A fully implicit monolithic
coupling scheme is intractable in this case due to the discrepancy in time
scales, $T_{f}\sim 1\,\text{ms}$ for the fluid and $T_{s}\to\infty$ for the
solid CMMe. Our coupling algorithm converged in only two coupling iterations
for $K_{\tau\sigma,o}=0.0$, allowing for a weakly coupled partitioned scheme.
However, the number of coupling iterations rose with $K_{\tau\sigma,o}\to 1$,
highlighting the strong interaction between fluid and solid. We experimented
with static relaxation for FSGe coupling. However, this required a very low
relaxation parameter $\omega$ to stabilize the coupling scheme, which led to
excessive coupling iterations. We achieved considerable speedups by
introducing dynamic Aitken relaxation [50] but could not find a solution for
high $K_{\tau\sigma,o}\to 1$. Thus, IQN-ILS remains our coupling algorithm of
choice [45]. More speedup could be achieved by taking into account the
iterations internal to each solver [60].

###  5.3 Computational stability of the equilibrated constrained mixture model

In addition to potential instabilities arising from the fluid-structure
coupling in FSGe (Section 5.2), we discovered instabilities in the underlying
solid CMMe model for $K_{\tau\sigma,o}\to 1$, leading to oscillations in
displacements, Jacobian, and internal G&R variables. These are triggered by
the growth of the aneurysm and can be interpreted as a mechanical instability
problem. There may exist other biological equilibria of $\Delta\tau_{w}$ and
$\Delta\sigma_{I}$ that satisfy a certain gain ratio $K_{\tau\sigma,o}>0$.
Only the special case $K_{\tau\sigma,o}=0=\Delta\sigma_{I}$ has a unique
biological solution for $\Delta\sigma_{I}$. A mechanobiologically equilibrated
solution is obtained by simultaneously solving for a mechanical equilibrium
(4). We believe these oscillations originate from evaluating the
mechanobiological equilibrium (12) point-wise on each Gauß-point, without
having any continuity conditions for internal G&R variables. Thus, it is
possible that Gauß-points within the same element may converge to different
mechanobiological equilibria. We currently prevent these oscillations by
choosing only one element across the vessel thickness. A more versatile
solution could be using an internal G&R variable, e.g., collagen mass fraction
$\phi^{c}_{h}$, as an additional unknown in addition to the displacements and
choosing a lower polynomial order for its shape function. Element formulations
could be adopted from (quasi)-incompressibility formulations, e.g., mixed-
field approaches such as Q1P0 [61] or MINI elements [62]. Alternatively, a
functional could be added to the solid weak form preventing oscillations in
collagen mass fraction. This idea is derived from edge stabilization in
computational fluid dynamics (CFD), where a least-squares term penalizes the
jump in the gradient of the discrete solution over element boundaries [63].

###  5.4 Definition of original homeostasis

In computational biomechanics, the solid reference configuration, often
obtained from in vivo medical imaging, is typically not stress-free due to
fluid or other loads. This requires either finding a stress-free reference
configuration, which might not exist [64], or local residual stresses so that
the initial configuration is in equilibrium with fluid loads. Many strategies
have been developed to address this issue [65, 66, 67]. In G&R, we need to
satisfy a biological equilibrium, i.e., original homeostasis $o$, in addition
to the aforementioned mechanical equilibrium. In this work, we achieve a
mechanobiological equilibrated original homeostasis by solving a non-linear
equation for a thin-walled cylinder for the initial in vivo blood pressure
$P_{o}$. Here, the inner and outer radii are fixed, as well as the homeostatic
parameters in Table 2. This approach works well in a given ideal cylinder with
an unknown uniform in vivo blood pressure that can be freely chosen in the
model. However, it cannot be generalized to in vivo blood vessel geometries
that were imaged under measured blood pressure. Given the limited prior work
with 3D CMM, there are few references for pre-loading in vivo G&R models. The
algorithm proposed by Wu et al. [68] solves a minimization problem by locally
varying fiber alignment. Laubrie et al. [69] prescribed a transmurally varying
prestretch gradient for constituents to establish an initial homeostatic
state. Gebauer et al. [70] introduced an iterative algorithm to prestress a
fixed patient-specific cardiac geometry under fixed ventricular pressures. The
algorithm solves for the local isotropic elastin prestretch until the model is
in mechanobiological equilibrium up to a specified displacement tolerance. The
choice of elastin prestretch was motivated by the fact that elastin is the
only constituent that is assumed to not turn over in the G&R time scale
considered, given its long half-life [71]. In the future, this approach could
be used to define the original homeostasis in FSGe problems by either locally
solving for the elastin prestretch or other quantities, such as vessel
thickness.

###  5.5 Influence of local wall shear stress

Previous FSG studies did not investigate differences across G&R models or
found only minor differences [13]. We observed in Section 4.2 that both G&R
and FSGe models delivered similar results for $K_{\tau\sigma,o}\to 0$ despite
the underlying assumptions in the purely solid G&R model. However, for high
gain ratios $K_{\tau\sigma,o}\to 1$, the different long-term evolved
equilibria between G&R and FSGe models became more pronounced. These
differences were especially prevalent in the bottom healthy region of the
blood vessel with intact elastin stiffness and mechanosensation. Compared to
G&R, the FSGe model grew inward with a three-fold increase in vessel wall
thickness. Furthermore, in Section 4.3, we observed that the inward growth of
the bottom region of the blood vessel is associated with an upward deflection
of blood flow toward the top of the aneurysm.

As mentioned previously, earlier versions of G&R relied on the Poiseuille and
thin-walled cylinder assumptions to calculate mechanical stimuli. Assuming
Poiseuille flow (1) and a thin-walled cylinder for the blood vessel, analytic
relationships can be derived for WSS $\tau_{w}$ and wall stress $\sigma$
components, respectively,

$\displaystyle\tau_{w}\sim Q\cdot\frac{1}{r^{3}},\quad\sigma_{\theta}\sim
P\cdot\frac{r}{h},\quad\sigma_{r}\sim-P,\quad\sigma_{z}\sim
F\cdot\frac{1}{rh},$ (38)

with radius $r$, wall thickness $h$, flow rate $Q$, pressure $P$, and axial
force $F$. We specifically developed our FSGe coupling method to not rely on
Poiseuille flow and thin-walled assumptions but rather the numerical solution
of the 3D Navier-Stokes equations and 3D finite strain elasticity,
respectively. Yet, we studied aneurysmal development under steady-state flow
with a low Reynolds number of $\text{Re}\approx 170$ which resulted in a flow
profile similar to parabolic even in aneurysmal regions of the blood vessel.
Thus, comparing the fluid and solid field solutions derived from our FSGe
model to these simplifying assumptions gives us insights into the feedback
loop of biomechanically regulated variables (Figure 1) and helps us explain
the differences in mechanobiological equilibria. In our examples, the flow
rate $Q$ was constant, and there was little change in local pressure $P$
throughout the fluid domain. We can thus derive the following relationships
for the change in WSS and intramural stimuli

$\displaystyle\Delta\tau_{w}$
$\displaystyle\sim\frac{r_{o}^{3}}{r_{h}^{3}}-1,$ (39)
$\displaystyle\Delta\sigma_{I}$
$\displaystyle\sim\frac{h_{o}}{r_{o}}\frac{r_{h}}{h_{h}}-1.$ (40)

The differences in aneurysmal growth were associated with these different
definitions of stimuli and, therefore, different values of mechanobiological
equilibria in the two models. Our illustrative computational examples
highlight the importance of incorporating local fluid-derived WSS in the FSGe
model. In the FSGe model, the fluid couples the deformation in the aneurysmal
top part of the blood vessel with the healthy bottom part. In the purely solid
G&R model, the WSS in these regions is evaluated independently from local
deformations. The stimuli $\Delta\tau_{w}<0$ and $\Delta\sigma_{I}<0$ were
downregulated to a much larger extent in FSGe than in G&R at the bottom of the
aneurysm. In the G&R model, relying on relationship (39) resulted in an
increase in evolved radius $r_{h}$ and thickness $h_{h}$. In the FSGe model,
the fluid-derived WSS, $\Delta\tau_{w}$ decreased with decreasing $r_{h}$
which is contrary to (39). This was possible due to the Navier-Stokes-derived
WSS (Figure 11, right) and the upward deflection of the blood flow inside the
aneurysm (Figure 12). In both models, this change was accompanied by an
increase in thickness $h_{h}$ according to (40).

###  5.6 Pulsatility and directionality of blood flow

Our FSGe model considered steady-state flow conditions as an input for the
intramural and WSS G&R stimuli, $\Delta\sigma_{I}$ and $\Delta\tau_{w}$,
respectively. This is equivalent to some previous FSG methods in Table 1,
whereas others used time-averaged G&R stimuli from a pulsatile flow
simulation. Pulsatility is an important determinant of the blood vessel wall
geometry, structure, and properties [72, 73]. For example, Eberth et al. [74]
found that increased wall thickness in hypertension correlated more strongly
with pulse than with mean or systolic pressure. Extending our FSGe approach to
incorporate pulsatility would require solving an fluid-structure interaction
(FSI) problem with an elastic wall in each coupling iteration. We could then
extract the instantaneous elastic properties of the CMMe with a fixed G&R
state from the evolved equilibrial state, as shown by Latorre and Humphrey
[5]. In the aneurysmal region, we currently introduce dysfunctional
mechanosensing at the endothelium by regionally setting $K_{\tau\sigma,o}\to
0$ in the top part of the aneurysm [56] as a model input. In the future,
mechanosensing could be a function of pulsatility and thus as part of the
model prediction.

Flow direction $\Longrightarrow$

Figure 14: Recirculation zone in an FSGe simulation with
$K_{\tau\sigma,o}=0.0$ (top left in Figure 12), with backflow velocity vectors
for $u_{z}<0$ at the top of the aneurysm.

Another crucial input would be the variability of WSS over a cardiac cycle.
Figure 14 shows a recirculation zone with backflow against the predominant
flow direction at the top of the aneurysm. This flow separation only occurs in
the FSGe model with $K_{\tau\sigma,o}=0.0$ (top left in Figure 12) due to its
nearly two-fold increase in radius. Note, however, that the maximum backflow
velocity is only 1.6% of the maximum inlet velocity. This effectively reverses
the direction of the local WSS vector $\bm{\tau}_{w}$ acting on the vessel
wall. However, in solid G&R models, $\tau_{w}$ can only be estimated as a
scalar quantity to inform the WSS stimulus, typically derived from Poiseuille
flow (1). It is well established that not only the magnitude but also the
pulsatility and direction of WSS are constantly sensed by the endothelium [75,
76]. For example, the directionality of WSS plays an important role in the
development of atherosclerosis [77]. Mohamied et al. [78] proposed a new
metric,

$\text{transWSS}=\frac{1}{T}\int_{0}^{T}\frac{\bm{\tau}_{w}\cdot\left(\bm{n}\times\bm{\tau}_{\text{mean}}\right)}{\|\bm{\tau}_{\text{mean}}\|}\,\text{d}t,\quad\bm{\tau}_{\text{mean}}=\frac{1}{T}\int_{0}^{T}\bm{\tau}_{w}\,\text{d}t,$
(41)

to quantify the change of direction of WSS over the course of a cardiac cycle
$T$, with vessel normal vector $\bm{n}$ and mean WSS vector
$\bm{\tau}_{\text{mean}}$. This metric was shown to correlate better with
atherosclerotic lesion prevalence than other metrics like time-averaged or
oscillatory WSS index. Our FSGe approach allows locally quantifying
directionality measures like transWSS. In the future, incorporating local and
directional WSS will enable the testing of hypotheses in FSGe-mediated
aneurysm growth and atherosclerotic lesion formation using FSGe. Since current
G&R models were designed with the static and scalar nature of reduced fluid
dynamics in mind, we expect future biological insights from FSGe to arise from
the advancements outlined in this work.

##  6 Acknowledgements

This work was supported by NIH Grants K99HL161313, R01HL139796, R01HL159954,
the Stanford Maternal and Child Health Research Institute, and the Additional
Ventures Foundation Cures Collaborative. We further thank Dr. Richard
Schussnig for his help with IQN-ILS and Drs. Alice Cortinovis and Luca
Pegolotti for their help with linear solvers. In addition, we thank Drs. Yuri
Bazilevs, Fabian Bräu, David Li, Matteo Salvador, and Vijay Vedula for
stimulating discussions on the stability of the CMMe model.

## References

  * [1] S. Loerakker, J. D. Humphrey, Computer model-driven design in cardiovascular regenerative medicine, Annals of Biomedical Engineering 51 (1) (2022) 45–57.  doi:10.1007/s10439-022-03037-5.
  * [2] J. D. Humphrey, K. R. Rajagopal, A constrained mixture model for growth and remodeling of soft tissues, Mathematical Models and Methods in Applied Sciences 12 (03) (2002) 407–430.  doi:10.1142/s0218202502001714.
  * [3] J. D. Humphrey, Constrained mixture models of soft tissue growth and remodeling – twenty years after, Journal of Elasticity 145 (1-2) (2021) 49–75.  doi:10.1007/s10659-020-09809-1.
  * [4] D. Dajnowiec, B. Langille, Arterial adaptations to chronic changes in haemodynamic function: coupling vasomotor tone to structural remodelling, Clinical Science 113 (1) (2007) 15–23.  doi:10.1042/cs20060337.
  * [5] M. Latorre, J. D. Humphrey, Fast, rate-independent, finite element implementation of a 3d constrained mixture model of soft tissue growth and remodeling, Computer Methods in Applied Mechanics and Engineering 368 (2020) 113156\.  doi:10.1016/j.cma.2020.113156.
  * [6] T. W. Secomb, Hemodynamics (mar 2016).  doi:10.1002/cphy.c150038.
  * [7] E. L. Schwarz, L. Pegolotti, M. R. Pfaller, A. L. Marsden, Beyond cfd: Emerging methodologies for predictive simulation in cardiovascular health and disease, Biophysics Reviews 4 (1) (2023) 011301.  arXiv:https://doi.org/10.1063/5.0109400, doi:10.1063/5.0109400.
URL https://doi.org/10.1063/5.0109400

  * [8] C. A. Figueroa, S. Baek, C. A. Taylor, J. D. Humphrey, A computational framework for fluid–solid-growth modeling in cardiovascular simulations, Computer Methods in Applied Mechanics and Engineering 198 (45-46) (2009) 3583–3602.  doi:10.1016/j.cma.2008.09.013.
  * [9] A. Sheidaei, S. Hunley, S. Zeinali-Davarani, L. Raguin, S. Baek, Simulation of abdominal aortic aneurysm growth with updating hemodynamic loads using a realistic geometry, Medical Engineering & Physics 33 (1) (2011) 80–88.  doi:10.1016/j.medengphy.2010.09.012.
  * [10] P. N. Watton, N. B. Raberger, G. A. Holzapfel, Y. Ventikos, Coupling the hemodynamic environment to the evolution of cerebral aneurysms: Computational framework and numerical examples, Journal of Biomechanical Engineering 131 (10) (Sep. 2009).  doi:10.1115/1.3192141.
  * [11] A. Grytsan, P. N. Watton, G. A. Holzapfel, A thick-walled fluid–solid-growth model of abdominal aortic aneurysm evolution: Application to a patient-specific geometry, Journal of Biomechanical Engineering 137 (3) (mar 2015).  doi:10.1115/1.4029279.
  * [12] M. Latorre, J. M. Szafron, A. B. Ramachandra, J. D. Humphrey, In vivo development of tissue engineered vascular grafts: a fluid-solid-growth model, Biomechanics and Modeling in Mechanobiology 21 (3) (2022) 827–848.  doi:10.1007/s10237-022-01562-9.
  * [13] E. L. Schwarz, M. R. Pfaller, J. M. Szafron, M. Latorre, S. E. Lindsey, C. K. Breuer, J. D. Humphrey, A. L. Marsden, A fluid–solid-growth solver for cardiovascular modeling, Computer Methods in Applied Mechanics and Engineering 417 (2023) 116312.  doi:10.1016/j.cma.2023.116312.
  * [14] E. L. Schwarz, J. M. Kelly, K. M. Blum, K. N. Hor, A. R. Yates, J. C. Zbinden, A. Verma, S. E. Lindsey, A. B. Ramachandra, J. M. Szafron, J. D. Humphrey, T. Shin’oka, A. L. Marsden, C. K. Breuer, Hemodynamic performance of tissue-engineered vascular grafts in fontan patients, npj Regenerative Medicine 6 (1) (jul 2021).  doi:10.1038/s41536-021-00148-w.
  * [15] K. M. Blum, J. C. Zbinden, A. B. Ramachandra, S. E. Lindsey, J. M. Szafron, J. W. Reinhardt, M. Heitkemper, C. A. Best, G. J. M. Mirhaidari, Y.-C. Chang, A. Ulziibayar, J. Kelly, K. V. Shah, J. D. Drews, J. Zakko, S. Miyamoto, Y. Matsuzaki, R. Iwaki, H. Ahmad, R. Daulton, D. Musgrave, M. G. Wiet, E. Heuer, E. Lawson, E. Schwarz, M. R. McDermott, R. Krishnamurthy, R. Krishnamurthy, K. Hor, A. K. Armstrong, B. A. Boe, D. P. Berman, A. J. Trask, J. D. Humphrey, A. L. Marsden, T. Shinoka, C. K. Breuer, Tissue engineered vascular grafts transform into autologous neovessels capable of native function and growth, Communications Medicine 2 (1) (jan 2022).  doi:10.1038/s43856-021-00063-7.
  * [16] A. B. Ramachandra, A. M. Kahn, A. L. Marsden, Patient-specific simulations reveal significant differences in mechanical stimuli in venous and arterial coronary grafts, Journal of Cardiovascular Translational Research 9 (4) (2016) 279–290.  doi:10.1007/s12265-016-9706-0.
  * [17] M. O. Khan, J. S. Tran, H. Zhu, J. Boyd, R. R. S. Packard, R. P. Karlsberg, A. M. Kahn, A. L. Marsden, Low wall shear stress is associated with saphenous vein graft stenosis in patients with coronary artery bypass grafting, Journal of Cardiovascular Translational Research 14 (4) (2020) 770–781.  doi:10.1007/s12265-020-09982-7.
  * [18] I. S. Lan, W. Yang, J. A. Feinstein, J. Kreutzer, R. T. Collins, M. Ma, G. T. Adamson, A. L. Marsden, Virtual transcatheter interventions for peripheral pulmonary artery stenosis in williams and alagille syndromes, Journal of the American Heart Association 11 (6) (mar 2022).  doi:10.1161/jaha.121.023532.
  * [19] J. M. Szafron, W. Yang, J. A. Feinstein, M. Rabinovitch, A. L. Marsden, A computational growth and remodeling framework for adaptive and maladaptive pulmonary arterial hemodynamics, Biomechanics and Modeling in Mechanobiology (sep 2023).  doi:10.1007/s10237-023-01744-z.
  * [20] D. Ambrosi, M. Ben Amar, C. J. Cyron, A. DeSimone, A. Goriely, J. D. Humphrey, E. Kuhl, Growth and remodelling of living tissues: perspectives, challenges and opportunities, Journal of The Royal Society Interface 16 (157) (2019) 20190233\.  doi:10.1098/rsif.2019.0233.
  * [21] S. Baek, R. Gleason, K. Rajagopal, J. Humphrey, Theory of small on large: Potential utility in computations of fluid–solid interactions in arteries, Computer Methods in Applied Mechanics and Engineering 196 (31-32) (2007) 3070–3078.  doi:10.1016/j.cma.2006.06.018.
  * [22] C. A. Figueroa, I. E. Vignon-Clementel, K. E. Jansen, T. J. Hughes, C. A. Taylor, A coupled momentum method for modeling blood flow in three-dimensional deformable arteries, Computer Methods in Applied Mechanics and Engineering 195 (41-43) (2006) 5685–5706.  doi:10.1016/j.cma.2005.11.011.
  * [23] M. R. Pfaller, J. Pham, A. Verma, L. Pegolotti, N. M. Wilson, D. W. Parker, W. Yang, A. L. Marsden, Automated generation of 0D and 1D reduced-order models of patient-specific blood flow, International Journal for Numerical Methods in Biomedical Engineering 38 (10) (aug 2022).  doi:10.1002/cnm.3639.
  * [24] svFSIplus (https://github.com/SimVascular/svFSIplus) (2024).
URL https://github.com/SimVascular/svFSIplus

  * [25] C. Zhu, V. Vedula, D. Parker, N. Wilson, S. Shadden, A. Marsden, svFSI: A multiphysics package for integrated cardiac modeling, Journal of Open Source Software 7 (78) (2022) 4118.  doi:10.21105/joss.04118.
  * [26] A. Updegrove, N. M. Wilson, J. Merkow, H. Lan, A. L. Marsden, S. C. Shadden, SimVascular: An open source pipeline for cardiovascular simulation, Annals of Biomedical Engineering 45 (3) (2016) 525–541.  doi:10.1007/s10439-016-1762-8.
  * [27] M. E. Moghadam, I. E. Vignon-Clementel, R. Figliola, A. L. Marsden, A modular numerical method for implicit 0d/3d coupling in cardiovascular finite element simulations, Journal of Computational Physics 244 (2013) 63–79.  doi:10.1016/j.jcp.2012.07.035.
  * [28] M. Esmaily-Moghadam, Y. Bazilevs, A. L. Marsden, A new preconditioning technique for implicitly coupled multidomain simulations with applications to hemodynamics, Computational Mechanics 52 (5) (2013) 1141–1152.  doi:10.1007/s00466-013-0868-1.
  * [29] J. D. Humphrey, S. Na, Elastodynamics and arterial wall stress, Annals of Biomedical Engineering 30 (4) (2002) 509–523.  doi:10.1114/1.1467676.
  * [30] P. Moireau, N. Xiao, M. Astorino, C. A. Figueroa, D. Chapelle, C. A. Taylor, J.-F. Gerbeau, External tissue support and fluid–structure simulation in blood flows, Biomechanics and Modeling in Mechanobiology 11 (1–2) (2011) 1–18.  doi:10.1007/s10237-011-0289-z.
  * [31] A. Valentín, L. Cardamone, S. Baek, J. Humphrey, Complementary vasoactivity and matrix remodelling in arterial adaptations to altered flow and pressure, Journal of The Royal Society Interface 6 (32) (2008) 293–306.  doi:10.1098/rsif.2008.0254.
  * [32] J. S. Wilson, S. Baek, J. D. Humphrey, Importance of initial aortic properties on the evolving regional anisotropy, stiffness and wall thickness of human abdominal aortic aneurysms, Journal of The Royal Society Interface 9 (74) (2012) 2047–2058.  doi:10.1098/rsif.2012.0097.
  * [33] M. Latorre, J. D. Humphrey, Numerical knockouts–in silico assessment of factors predisposing to thoracic aortic aneurysms, PLOS Computational Biology 16 (10) (2020) e1008273.  doi:10.1371/journal.pcbi.1008273.
  * [34] A. B. Ramachandra, M. Latorre, J. M. Szafron, A. L. Marsden, J. D. Humphrey, Vascular adaptation in the presence of external support - a modeling study, Journal of the Mechanical Behavior of Biomedical Materials 110 (2020) 103943.  doi:10.1016/j.jmbbm.2020.103943.
  * [35] J. M. Szafron, R. Khosravi, J. Reinhardt, C. A. Best, M. R. Bersi, T. Yi, C. K. Breuer, J. D. Humphrey, Immuno-driven and mechano-mediated neotissue formation in tissue engineered vascular grafts, Annals of Biomedical Engineering 46 (11) (2018) 1938–1950.  doi:10.1007/s10439-018-2086-7.
  * [36] A. Valentín, J. Humphrey, G. Holzapfel, A finite element-based constrained mixture implementation for arterial growth, remodeling, and adaptation: Theory and numerical verification, International Journal for Numerical Methods in Biomedical Engineering 29 (8) (2013) 822–849.  doi:10.1002/cnm.2555.
  * [37] N. Horvat, L. Virag, G. A. Holzapfel, J. Sorić, I. Karšaj, A finite element implementation of a growth and remodeling model for soft biological tissues: Verification and application to abdominal aortic aneurysms, Computer Methods in Applied Mechanics and Engineering 352 (2019) 586–605.  doi:10.1016/j.cma.2019.04.041.
  * [38] M. Latorre, J. D. Humphrey, A mechanobiologically equilibrated constrained mixture model for growth and remodeling of soft tissues, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 98 (12) (2018) 2048–2071.  doi:10.1002/zamm.201700302.
  * [39] S. A. Maas, B. J. Ellis, G. A. Ateshian, J. A. Weiss, FEBio: Finite elements for biomechanics, Journal of Biomechanical Engineering 134 (1) (jan 2012).  doi:10.1115/1.4005694.
  * [40] S. Goswami, D. S. Li, B. V. Rego, M. Latorre, J. D. Humphrey, G. E. Karniadakis, Neural operator learning of heterogeneous mechanobiological insults contributing to aortic aneurysms, Journal of The Royal Society Interface 19 (193) (aug 2022).  doi:10.1098/rsif.2022.0410.
  * [41] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications 23 (1) (2001) 15–41.  doi:10.1137/s0895479899358194.
  * [42] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, et al., Petsc users manual, Tech. Rep. No. ANL-95/11 - Revision 2.3.2 (2019).
  * [43] J. D. Humphrey, S. L. O’Rourke, An Introduction to Biomechanics, Springer New York, 2015.  doi:10.1007/978-1-4939-2623-7.
  * [44] G. Hou, J. Wang, A. Layton, Numerical methods for fluid-structure interaction — a review, Communications in Computational Physics 12 (2) (2012) 337–377.  doi:10.4208/cicp.291210.290411s.
  * [45] J. Degroote, Partitioned simulation of fluid-structure interaction, Archives of Computational Methods in Engineering 20 (3) (2013) 185–238.  doi:10.1007/s11831-013-9085-5.
  * [46] svFSGe (https://github.com/StanfordCBCL/svFSGe) (2024).
URL https://github.com/StanfordCBCL/svFSGe

  * [47] R. Schussnig, D. R. Pacheco, M. Kaltenbacher, T.-P. Fries, Semi-implicit fluid–structure interaction in biomedical applications, Computer Methods in Applied Mechanics and Engineering 400 (2022) 115489.  doi:10.1016/j.cma.2022.115489.
  * [48] P. C. Africa, R. Piersanti, M. Fedele, L. Dede’, A. Quarteroni, lifex-fiber: an open tool for myofibers generation in cardiac computational models, BMC Bioinformatics 24 (1) (apr 2023).  doi:10.1186/s12859-023-05260-w.
  * [49] M. A. Fernández, Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit, SeMA Journal 55 (1) (2011) 59–108.  doi:10.1007/bf03322593.
  * [50] U. Küttler, W. A. Wall, Fixed-point fluid–structure interaction solvers with dynamic relaxation, Computational Mechanics 43 (1) (2008) 61–72.  doi:10.1007/s00466-008-0255-5.
  * [51] U. Küttler, Effiziente lösungsverfahren für fluid-struktur-interaktions-probleme, Dissertation, Technische Universität München, München (2009).
  * [52] J. Degroote, K.-J. Bathe, J. Vierendeels, Performance of a new partitioned procedure versus a monolithic procedure in fluid–structure interaction, Computers & Structures 87 (11-12) (2009) 793–801.  doi:10.1016/j.compstruc.2008.11.013.
  * [53] A. S. Les, S. C. Shadden, C. A. Figueroa, J. M. Park, M. M. Tedesco, R. J. Herfkens, R. L. Dalman, C. A. Taylor, Quantification of hemodynamics in abdominal aortic aneurysms during rest and exercise using magnetic resonance imaging and computational fluid dynamics, Annals of Biomedical Engineering 38 (4) (2010) 1288–1313.  doi:10.1007/s10439-010-9949-x.
  * [54] A. Feintuch, P. Ruengsakulrach, A. Lin, J. Zhang, Y.-Q. Zhou, J. Bishop, L. Davidson, D. Courtman, F. S. Foster, D. A. Steinman, R. M. Henkelman, C. R. Ethier, Hemodynamics in the mouse aortic arch as assessed by mri, ultrasound, and numerical modeling, American Journal of Physiology-Heart and Circulatory Physiology 292 (2) (2007) H884–H892.  doi:10.1152/ajpheart.00796.2006.
  * [55] P. N. Watton, N. A. Hill, M. Heil, A mathematical model for the growth of the abdominal aortic aneurysm, Biomechanics and Modeling in Mechanobiology 3 (2) (2004) 98–113.  doi:10.1007/s10237-004-0052-9.
  * [56] J. D. Humphrey, D. M. Milewicz, G. Tellides, M. A. Schwartz, Dysfunctional mechanosensing in aneurysms, Science 344 (6183) (2014) 477–479.  doi:10.1126/science.1253026.
  * [57] A. B. Ramachandra, S. Sankaran, J. D. Humphrey, A. L. Marsden, Computational simulation of the adaptive capacity of vein grafts in response to increased pressure, Journal of Biomechanical Engineering 137 (3) (Mar. 2015).  doi:10.1115/1.4029021.
  * [58] C. Cyron, J. Humphrey, Vascular homeostasis and the concept of mechanobiological stability, International Journal of Engineering Science 85 (2014) 203–223.  doi:10.1016/j.ijengsci.2014.08.003.
  * [59] M. Latorre, J. D. Humphrey, Mechanobiological stability of biological soft tissues, Journal of the Mechanics and Physics of Solids 125 (2019) 298–325.  doi:10.1016/j.jmps.2018.12.013.
  * [60] T. Spenke, N. Delaissé, J. Degroote, N. Hosters, On the number of subproblem iterations per coupling step in partitioned fluid-structure interaction simulations, arXiv (Mar. 2023).  arXiv:2303.08513, doi:10.48550/ARXIV.2303.08513.
  * [61] O. Gültekin, H. Dal, G. A. Holzapfel, On the quasi-incompressible finite element analysis of anisotropic hyperelastic materials, Computational Mechanics 63 (3) (2018) 443–453.  doi:10.1007/s00466-018-1602-9.
  * [62] E. Karabelas, G. Haase, G. Plank, C. M. Augustin, Versatile stabilized finite element formulations for nearly and fully incompressible solid mechanics, Computational Mechanics 65 (1) (2019) 193–215.  doi:10.1007/s00466-019-01760-w.
  * [63] E. Burman, P. Hansbo, The edge stabilization method for finite elements in cfd, in: Numerical Mathematics and Advanced Applications, Springer Berlin Heidelberg, 2004, pp. 196–203.  doi:10.1007/978-3-642-18775-9_17.
  * [64] M. W. Gee, C. Förster, W. A. Wall, A computational strategy for prestressing patient-specific biomechanical problems under finite deformation, International Journal for Numerical Methods in Biomedical Engineering 26 (1) (2009) 52–72.  doi:10.1002/cnm.1236.
  * [65] M. K. Rausch, M. Genet, J. D. Humphrey, An augmented iterative method for identifying a stress-free reference configuration in image-based biomechanical modeling, Journal of Biomechanics 58 (2017) 227–231.  doi:10.1016/j.jbiomech.2017.04.021.
  * [66] M. Genet, A relaxed growth modeling framework for controlling growth-induced residual stresses, Clinical Biomechanics 70 (2019) 270–277.  doi:10.1016/j.clinbiomech.2019.08.015.
  * [67] N. Barnafi, F. Regazzoni, D. Riccobelli, Reconstructing relaxed configurations in elastic bodies: Mathematical formulations and numerical methods for cardiac modeling, Computer Methods in Applied Mechanics and Engineering 423 (2024) 116845.  doi:10.1016/j.cma.2024.116845.
  * [68] J. Wu, C. M. Augustin, S. C. Shadden, Reconstructing vascular homeostasis by growth-based prestretch and optimal fiber deposition, Journal of the Mechanical Behavior of Biomedical Materials (2020) 104161doi:10.1016/j.jmbbm.2020.104161.
  * [69] J. D. Laubrie, S. J. Mousavi, S. Avril, About prestretch in homogenized constrained mixture models simulating growth and remodeling in patient-specific aortic geometries, Biomechanics and Modeling in Mechanobiology 21 (2) (2022) 455–469.  doi:10.1007/s10237-021-01544-3.
  * [70] A. M. Gebauer, M. R. Pfaller, F. A. Braeu, C. J. Cyron, W. A. Wall, A homogenized constrained mixture model of cardiac growth and remodeling: analyzing mechanobiological stability and reversal, Biomechanics and Modeling in Mechanobiology (jul 2023).  doi:10.1007/s10237-023-01747-w.
  * [71] A. J. Cocciolone, J. Z. Hawes, M. C. Staiculescu, E. O. Johnson, M. Murshed, J. E. Wagenseil, Elastin, arterial mechanics, and cardiovascular disease, American Journal of Physiology-Heart and Circulatory Physiology 315 (2) (2018) H189–H205.  doi:10.1152/ajpheart.00087.2018.
  * [72] J. F. Eberth, N. Popovic, V. C. Gresham, E. Wilson, J. D. Humphrey, Time course of carotid artery growth and remodeling in response to altered pulsatility, American Journal of Physiology-Heart and Circulatory Physiology 299 (6) (2010) H1875–H1883.  doi:10.1152/ajpheart.00872.2009.
  * [73] M. Kalenik, M. Chalecki, P. Wichowski, Real values of local resistance coefficients during water flow through welded polypropylene t-junctions, Water 12 (3) (2020) 895.  doi:10.3390/w12030895.
  * [74] J. F. Eberth, V. C. Gresham, A. K. Reddy, N. Popovic, E. Wilson, J. D. Humphrey, Importance of pulsatility in hypertensive carotid artery growth and remodeling, Journal of Hypertension 27 (10) (2009) 2010–2021.  doi:10.1097/hjh.0b013e32832e8dc8.
  * [75] P. F. Davies, Hemodynamic shear stress and the endothelium in cardiovascular pathophysiology, Nature Clinical Practice Cardiovascular Medicine 6 (1) (2008) 16–26.  doi:10.1038/ncpcardio1397.
  * [76] C. Wang, B. M. Baker, C. S. Chen, M. A. Schwartz, Endothelial cell sensing of flow direction, Arteriosclerosis, Thrombosis, and Vascular Biology 33 (9) (2013) 2130–2136.  doi:10.1161/atvbaha.113.301826.
  * [77] N. Baeyens, M. A. Schwartz, Biomechanics of vascular mechanosensation and remodeling, Molecular Biology of the Cell 27 (1) (2016) 7–11.  doi:10.1091/mbc.e14-11-1522.
  * [78] Y. Mohamied, E. M. Rowland, E. L. Bailey, S. J. Sherwin, M. A. Schwartz, P. D. Weinberg, Change of direction in the biomechanics of atherosclerosis, Annals of Biomedical Engineering 43 (1) (2014) 16–25.  doi:10.1007/s10439-014-1095-4.