#  Neural operator learning of heterogeneous mechanobiological insults
contributing to aortic aneurysms

Somdatta Goswami111These authors contributed equally to this work. David S.
Li222These authors contributed equally to this work. Bruno V. Rego  Marcos
Latorre
Jay D. Humphrey  George Em Karniadakis  Division of Applied Mathematics, Brown
University, Providence, RI, USA  Department of Biomedical Engineering, Yale
University, New Haven, CT, USA  Center for Research and Innovation in
Bioengineering, Valencia Polytechnic University, Valencia, Spain  School of
Engineering, Brown University, Providence, RI, USA

###### Abstract

Thoracic aortic aneurysm (TAA) is a localized dilatation of the aorta
resulting from compromised wall composition, structure, and function, which
can lead to life-threatening dissection or rupture. Several genetic mutations
and predisposing factors that contribute to TAA have been studied in mouse
models to characterize specific changes in aortic microstructure and material
properties that result from a wide range of mechanobiological insults. By
contrast, assessments of TAA progression _in vivo_ are largely limited to
measurements of aneurysm size and growth rate. It has been shown, however,
that aortic geometry alone is not sufficient to predict the patient-specific
progression of TAA but computational modeling of the evolving biomechanics of
the aorta could predict future geometry and properties from initiating
insults. In this work, we present an integrated framework to train a deep
operator network (DeepONet)-based surrogate model to identify contributing
factors for TAA by using synthetic finite element-based datasets of aortic
growth and remodeling (G&R) resulting from prescribed mechanobiological
insults. For training data, we investigate multiple types of TAA risk factors
and spatial distributions within a computationally efficient constrained
mixture model to generate axial–azimuthal maps of aortic dilatation and
distensibility. The trained network is then capable of predicting the initial
distribution and extent of the insult from a given set of dilatation and
distensibility information, which in turn can be used to determine subsequent
aortic geometry and mechanical properties. Two DeepONet frameworks are
proposed, one trained on sparse information and one on full-field grayscale
images, to gain insight into a preferred neural operator-based approach.
Performance of the surrogate models is evaluated through multiple simulations
carried out on insult distributions varying from fusiform (analytically
defined) to complex (randomly generated). We show that this integrated
continuous learning modeling approach can predict the patient-specific
mechanobiological insult profile associated with any given dilatation and
distensibility map with a high accuracy, particularly when based on full-field
images. Our findings demonstrate the feasibility of applying DeepONet to
support transfer learning of patient-specific inputs (e.g., age, hypertension,
diabetes, Marfan syndrome) to predict TAA progression with clinical images of
the aorta.

###### keywords:

Operator-based neural network , deep learning , growth and remodeling ,
thoracic aortic aneurysm

##  1 Introduction

Thoracic aortic aneurysms (TAAs) are localized dilatations of the aorta that
associate with a higher risk of life-threatening aortic dissection or rupture;
they can initiate from a variety of biomechanical and genetic factors, often
developing over several years [1, 2, 3, 4, 5, 6]. Treatment of TAAs may
involve surgical replacement with a synthetic graft or repair via the
placement of an endovascular stent [7, 8]. Determination of the optimal
approach and time of intervention depends on myriad factors, including the
size, growth rate and location of the aneurysm, as well as the individual’s
genetic history [9]; ultimately, to facilitate systematic improvement of the
patient’s prognosis and therapeutic design, there is a pressing need to
understand the complex roles these aspects play in the development of
aneurysms [10].

Assessment of aortic health in the clinic is largely limited to _in vivo_
anatomical information and hemodynamic measurements, with few biomarkers
available. On the other hand, murine models of aortic aneurysm have provided
valuable insight into the structural and biomechanical properties of the
normal and diseased thoracic aorta via _in vitro_ experimentation of excised
tissue specimens [11, 12] that complement _in vivo_ studies. Several genetic
mutations leading to compromised aortic structure and function have been
identified as critical predisposing factors driving TAA formation, including
those affecting extracellular matrix integrity, smooth muscle contractile
dysfunction, and aberrant intracellular signaling [13, 14, 2, 15]. Although
much has been learned from these studies, multiple contributors are often
present in combination _in vivo_ , rendering it challenging to identify
relationships among different mechanisms. Nevertheless, it is clear that
characterization of the biomechanics of the aorta is necessary to gain a
deeper knowledge of aortic disease progression. Toward this end, computational
models of aneurysm growth and remodeling can facilitate mechanistic
understanding by isolating the influence of individual biomechanical defects
on subsequent aneurysm progression [16, 17, 18, 19]. A theoretical model with
an ability to generate robust predictions of aneurysm progression from limited
anatomical information can both generate synthetic datasets for analysis and
lay the foundation for enriched diagnosis of aortic function beyond clinically
available measurements.

Current advancements in modeling now provide the opportunity to leverage
machine learning, which has emerged as an effective surrogate model for high-
fidelity solvers, in order to overcome previous computational hurdles that
would otherwise make such modeling intractable for clinically relevant time
frames. Such surrogate models have demonstrated the potential for automated
measurement of aortic geometry, patient risk stratification, and the
prediction of aneurysm growth and rupture [20, 21, 22, 23, 24, 6].
Additionally, physics-informed neural networks (PINNs) [25, 26, 27, 28] have
been a promising advancement in the domain of scientific machine learning.
However, to simulate multiple initial/boundary problems having different
applied loading, one often needs to retrain a PINN. Hence, developing models
that can learn the operator-level mapping between functions (that is,
forecasting the physical system under a variety of initial/boundary
circumstances) is critical [29, 30]. In this case, neural operators can learn
nonlinear mappings between function spaces, providing a novel simulation
paradigm for real-time prediction of complicated dynamics. A deep operator-
based neural network (DeepONet) proposed in [31] is now popular for learning
solutions from labelled input-output datasets consisting of varied
initial/boundary conditions and different forcing functions. The idea of
DeepONet is motivated by the Universal Approximation Theorem for Operators,
which defines a new and relatively underexplored realm for deep neural network
(DNN)-based approaches that map infinite-dimensional functional spaces rather
than finite-dimensional vector spaces (functional regression). The
computational model consists of two classes of DNNs; one encodes the input
function at fixed sensor points (branch net) while the other accounts for the
location of the output function (trunk net).

Here, we present an integrated computational platform for predicting
contributing factors to TAA by melding a constrained mixture model of TAA
enlargement with DeepONet. Specifically, we use a previously established
computationally efficient implementation of our mechanobiologically
equilibrated constrained mixture model to describe the long-term evolution of
TAAs. We propose two surrogate models to approximate aneurysmal initiating
factors: the first framework admits sparse information to be encoded in the
branch network, while the second instead encodes full-field grayscale images
taking into account patient-specific characteristics [32]. The accuracy of
each approach is evaluated for multiple types of simulated aneurysms arising
from multiple contributors investigated in past studies, and several
variations of the neural network design are explored to compare the relative
performance of different architectures. Our findings demonstrate feasibility
of this technique for future analysis of clinical images of the aorta. The key
highlights are:

  * 1.

A novel framework to predict TAA pathology melds a constrained mixture model
for arterial growth and remodeling with a DeepONet-based surrogate model.

  * 2.

3D finite element simulations of TAA progression arise from randomly
distributed losses of elastic fiber integrity and dysfunctional
mechanosensing.

  * 3.

The generalizable DeepONet can predict the solution with sufficient accuracy,
even when provided limited information.

  * 4.

Performance is improved by employing convolutional neural networks rather than
fully connected feed-forward neural networks with sparse information.

  * 5.

The preferred network architecture takes as input grayscale images of
dilatation and distensibility to predict the insult profile.

  * 6.

This approach can be extended to patient-specific medical images to provide
patient-specific solutions.

##  2 Methods

To generate synthetic data for training the surrogate model, we employ a
previously established constrained mixture model within a finite element
platform [33]. Finite element simulations incorporating the initial geometry
of the aorta, bulk mechanical properties of the vessel wall, and _in vivo_
loading conditions (axial pre-stretch and intraluminal blood pressure), along
with a prescribed mechanobiological insult to initiate TAA, allow predictions
of local dilatation and distensibility fields (subsection 2.2). These are
converted into axial–azimuthal maps and randomly categorized into testing and
training sets, which enable training of the DeepONet (subsection 2.4) to
predict the initial insult profile from a given dataset. Finally, prediction
of the spatial distribution and severity of the insult profile is used for the
projection of future TAA growth. The modeling pipeline is summarized in Figure
1.

Figure 1: Pipeline for training the DeepONet to predict mechanobiological
insults for forecasting long-term aneurysm growth. First, 3D finite element
simulations of TAA evolution incorporating vessel geometry, mechanical
properties, _in vivo_ loading conditions, and mechanobiological insult
profiles are used to generate synthetic training and testing data for the
DeepONet surrogate model, which is comprised of a deep neural network with
associated trainable parameters and evaluation points. Subsequently, the
DeepONet is trained using labeled maps of local dilatation and distensibility
corresponding to the prescribed insults, as either axial–azimuthal maps or
converted grayscale images. Finally, the insult profile (distribution and
severity) is predicted from a given set of dilatation and distensibility
information, which can be related to estimation of TAA progression.

###  2.1 Constrained mixture model for arterial growth and remodeling

Over the past several years, we have developed a constrained mixture
formulation to describe growth and remodeling (G&R) of the aortic wall, which
has been used to model the changing composition, structure, and mechanics in
both animal and clinical studies [34, 35, 36, 37]. We have recently
implemented a computationally efficient 3D finite element framework for
modeling the initiation and long-term growth of TAA within a thick-walled
cylindrical vessel segment [33]. In the following, we summarize the pertinent
details of the modeling framework; further details can be found in the
original papers [37, 38].

####  2.1.1 Modeling framework

We model the arterial wall as a mixture of three primary load-bearing
constituents: elastin-dominated matrix $e$, collagen fiber-dominated matrix
$c$, and smooth muscle cells $m$, that can undergo changes in mass (grow) and
microstructure (remodel) in response to biomechanical stimuli. These
constituents $\alpha=e,c,m$ are governed by mass density production rates
$m^{\alpha}(\tau)$ at the time of deposition $\tau$, removal (decay) functions
$q^{\alpha}(s,\tau)$, with $s$ denoting the current time, and stored energy
density functions $\hat{W}^{\alpha}(s,\tau)$ that depend on the constituent-
specific deformations relative to their evolving natural configurations
$n(\tau)$ and thus describe the mechanical behaviors.

Because the extracellular environment is mechanoregulated by cells within the
vessel wall, production and removal rates can be modulated by perturbations
from the homeostatic state (i.e., a quasi-equilibrium where production
balances removal, denoted by the subscript $o$). Production per unit volume at
time $\tau$ is expressed using a potentially evolving nominal rate
$m_{N}^{\alpha}(\tau)>0$ that is modulated by deviations in intramural stress
$\sigma$ (induced by pressure and axial force) and wall shear stress
$\tau_{w}$ (induced by blood flow) from homeostatic values, namely,

$m^{\alpha}(\tau)=m_{N}^{\alpha}(\tau)\bigg{(}1+K_{\sigma}^{\alpha}\Delta\sigma(\tau)-K_{\tau_{w}}^{\alpha}\Delta\tau_{w}(\tau)\bigg{)}=m_{N}^{\alpha}(\tau)\Upsilon^{\alpha}(\tau),$
(1)

where $K_{i}^{\alpha}$ are gain-type parameters controlling the sensitivity to
stress deviations $\Delta\sigma$ and $\Delta\tau_{w}$ from set points
$\sigma_{o}$ and $\tau_{wo}$, summarized in the stimulus function
$\Upsilon^{\alpha}$. Equation 1 is formulated to replicate experimental
observations that increased wall stress and increased flow tend to heighten
and reduce extracellular matrix production, respectively [16, 37]. In this
study, we define $\Delta\sigma=(\sigma-\sigma_{o})/\sigma_{o}$ and
$\Delta\tau_{w}=(\tau_{w}-\tau_{wo})/\tau_{wo}$, where $\sigma$ is one-third
the trace of the in-plane wall stress. Additionally, we assume that the
intraluminal pressure within the aorta remains constant (no pulsatility) when
computing the G&R and that the flow through the aorta remains constant and
laminar.

Constituent removal between times $\tau$ and $s$ is modeled by

$q^{\alpha}(s,\tau)=\exp\bigg{(}-\int_{\tau}^{s}k_{N}^{\alpha}\bigg{(}1+\omega\Big{(}\Delta\sigma(t)\Big{)}^{2}\bigg{)}dt\bigg{)},$
(2)

where $k_{N}^{\alpha}$ is the basal removal rate, and $\omega>0$ is a gain-
type parameter for intramural stress deviations. Note that, while increases in
intramural stress can drive increased production, as stated above, they can
also stimulate the removal of constituents through activation of matrix-
degrading enzymes, as reflected in Equation 2. With the above definitions, we
then express the mass density per unit reference volume of each constituent
$\rho_{R}^{\alpha}$ by

$\rho_{R}^{\alpha}(s)=\int_{-\infty}^{s}m_{R}^{\alpha}(\tau)\,q^{\alpha}(\tau)\,d\tau,$
(3)

where $m_{R}^{\alpha}=Jm^{\alpha}$ is the referential (subscript $R$) mass
density production rate, with $J=\sum\rho_{R}^{\alpha}/\rho$ representing the
volume ratio between the reference and current _in vivo_ configurations of the
tissue constituents.

We employ a hyperelasticity framework to describe the transient biomechanical
properties of the vessel wall, which is modeled as a mixture of nonlinear,
nearly incompressible, anisotropic materials. This allows the intramural
stress-driven changes in production and removal to be defined in terms of
stored energy density relations $\hat{W}^{\alpha}(s,\tau)$, which are
functions of the multiaxial deformations of each constituent. The Cauchy
stress at the tissue (mixture) level $\bm{\sigma}$ includes passive
contributions from each constituent, given by

$\bm{\sigma}(s)=-p(s)\mathbf{I}+\frac{2}{J(s)}\,\mathbf{F}(s)\,\frac{\partial
W_{R}(s)}{\partial\mathbf{C}(s)}\,\mathbf{F}^{\mathrm{T}}(s),$ (4)

where $p$ is a Lagrange multiplier that enforces the transient isochoric
motions, $\mathbf{I}$ is the identity tensor, $\mathbf{F}$ is the deformation
gradient tensor for the tissue from reference to current configurations,
$\mathbf{C}=\mathbf{F}^{\mathrm{T}}\mathbf{F}$ is the right Cauchy-Green
tensor, $J=\det\mathbf{F}$, and $W_{R}=\sum W_{R}^{\alpha}$ is the total
stored energy of the mixture. The constitutive relations are defined
separately for $e$, $c$, and $m$ and are weighted by their respective volume
fractions $\phi^{\alpha}$ ($W_{R}^{\alpha}=\phi^{\alpha}\hat{W}^{\alpha}$).
Further details are described in subsection A.1. Note that, because many TAAs
associate with diminished or absent smooth muscle contractility, we assume
that the stress depends only on passive properties.

####  2.1.2 Mechanobiologically equilibrated constrained mixture model

In cases of TAA, where the characteristic timescale of G&R is frequently
shorter than that of the biomechanical stimulus (e.g., elastic fiber
degradation with gradual hypertension), G&R can be assumed to reach a quasi-
static mechanobiological equilibrium at time $s\gg 0$, thus allowing a time-
independent approach at which production balances removal $\forall s$ and
computation of the heredity integrals is not required. In particular, the
stimulus function in Equation 1 reduces to

$\Upsilon_{h}^{\alpha}(\Delta\sigma_{h},\Delta\tau_{wh})=1,$ (5)

with $h$ denoting the evolved homeostatic state. Accordingly, rule-of-mixtures
expressions can be used to be used for stored energy
($W_{Rh}=\sum\phi_{Rh}^{\alpha}\hat{W}^{\alpha}$, where $\phi_{Rh}^{\alpha}$
are the evolved constituent mass fractions). Similarly, for the Cauchy
stresses (Equation 4),

$\bm{\sigma}_{h}=-p_{h}\mathbf{I}+\sum_{\alpha}^{e,c,m}\phi_{h}^{\alpha}\hat{\bm{\sigma}}_{h}^{\alpha},$
(6)

where $\hat{\bm{\sigma}}_{h}^{\alpha}$ are the constituent-specific Cauchy
stresses, and $p_{h}$ is the equilibrated Lagrange multiplier for the quasi-
static G&R evolution (see subsection A.2 for more details).

###  2.2 Modeling TAA growth from prescribed mechanobiological insults

Although TAAs typically exhibit an irregular diameter, eccentricity, and
thickness, it remains instructive to consider smoothly varying insult profiles
prescribed on an initially straight cylindrical vessel, with _in vivo_
geometry and mechanical properties derived from our previous studies (uniform
wall thickness $h_{o}=40$ µm and luminal radius $r_{o}=647$ µm; see subsection
A.2). Here, aneurysms are initiated by prescribing an insult that emulates
defects such as a localized breakage in elastic fibers or disruption at
cellular integrin binding sites. We adopt a cylindrical coordinate system in
the initial homeostatic state $\\{r_{o},\theta_{o},z_{o}\\}$, and two
approaches for insult profile definitions are considered: analytical and
randomly generated.

####  2.2.1 Analytically defined insult profiles

Insult profiles $\vartheta$ varying in the $z_{o}$–$\theta_{o}$
(axial–azimuthal) plane are defined analytically with the expression [33]

$\vartheta(z_{o},\theta_{o})=\vartheta_{end}\,+\,\bigg{(}\vartheta_{apex}-\vartheta_{end}\bigg{)}\exp\left(-\left|\frac{z_{o}-z_{apex}}{z_{od}}\right|^{\nu_{z}}\right)\exp\left(-\left|\frac{\theta_{o}-\theta_{apex}}{\theta_{od}}\right|^{\nu_{\theta}}\right),$
(7)

where $z_{o}\in[0,l_{o}]$, $l_{o}$ is the initial axial length (15 mm),
$\theta_{o}\in[0,2\pi]$, $z_{od}$ and $\theta_{od}$ are the axial and
circumferential characteristic widths of the insult region, respectively,
$\nu_{z}$ and $\nu_{\theta}$ govern the softness of the boundaries in the
axial and circumferential directions, respectively, and $\vartheta_{end}$ and
$\vartheta_{apex}$ are values of the insult at the ends of the cylinder
($z_{o}=0$, $l_{o}$) and the apex ($z_{o}=z_{apex}$,
$\theta_{o}=\theta_{apex}$) of the profile, respectively. The profile is
normalized to the interval $[0,1]$, with 1 indicating the maximum insult
degree that varies depending on the insult type and severity. Example insult
profiles and their corresponding remodeled vessels are shown in Figure 2, and
parameter ranges are listed in Table 1.

Figure 2: Examples of analytically defined insults and their influence on TAA
geometry. Dilatation (remodeled inner radius normalized by initial inner
radius) maps are shown for variations in the (a) axial location
$z_{apex}\in[6,9]$ mm and (b) circumferential location
$\theta_{apex}\in[0^{\circ},360^{\circ}]$ of the aneurysm apex, as well as the
aneurysm (c) axial characteristic width $z_{od}\in[2,4]$ mm and (d)
circumferential characteristic width $\theta_{od}\in[20^{\circ},360^{\circ}]$.
(e) For each combination of apex location and characteristic width, insults
(e.g., loss of elastic fiber integrity) of increasing severity are modeled.

####  2.2.2 Randomly generated insult profiles

In contrast to the definition above, we also model aneurysms initiated from
insult profiles randomly distributed along the vessel wall to validate the
DeepONet model against more natural (i.e., irregular) insults. These profiles
are initially generated as “latent” (i.e., unobserved) Gaussian random fields
(GRFs), then nonlinearly transformed and censored on $[0,1]$ to match
physically meaningful metrics prescribed by the user. On an unbounded space,
the latent insult profiles are thus sampled according to

$\vartheta^{\ast}(z_{o},\theta_{o})\sim\mathcal{G}\big{(}\mu(z_{o},\theta_{o}),\kappa(z_{o},\theta_{o},z_{o}^{\prime},\theta_{o}^{\prime})\big{)}.$
(8)

The mean $\mu(z_{o},\theta_{o})$ and covariance
$\kappa(z_{o},\theta_{o},z_{o}^{\prime},\theta_{o}^{\prime})$ between two
points $(z_{o},\theta_{o})$ and $(z_{o}^{\prime},\theta_{o}^{\prime})$ can be
parameterized to control the overall propensity of insult ($\varphi$), the
length scale of the insult(s) in the circumferential ($L_{\theta}$) and axial
($L_{z}$) directions, and the softness of the boundaries between normal and
insult regions ($\epsilon$). To enforce periodicity in the circumferential
direction, we define the covariance function as

$\displaystyle\kappa(z_{o},\theta_{o},z_{o}^{\prime},\theta_{o}^{\prime})$
$\displaystyle=\varsigma^{2}\exp\left(-\frac{1}{2}\left[\left(\frac{D_{\theta}\left(\theta_{o},\theta_{o}^{\prime}\right)}{L_{\theta}}\right)^{2}+\left(\frac{D_{z}\left(z_{o},z_{o}^{\prime}\right)}{L_{z}}\right)^{2}\right]\right)$
(9) $\displaystyle\textrm{with}\quad
D_{\theta}\left(\theta_{o},\theta_{o}^{\prime}\right)$
$\displaystyle=2r_{o}\sin{\left(\frac{1}{2}\left|\theta_{o}-\theta_{o}^{\prime}\right|\right)}\quad\text{and}\quad
D_{z}\left(z_{o},z_{o}^{\prime}\right)=\left|z_{o}-z_{o}^{\prime}\right|,$

where $\varsigma^{2}$ is the overall variance of the GRF. To satisfy the
prescribed insult propensity $\varphi$ and boundary softness $\epsilon$, the
mean and variance of the GRF are defined as

$\displaystyle\mu$
$\displaystyle=\frac{1}{2}-\frac{1}{\epsilon\sqrt{\pi}}\,\mathrm{erf}^{-1}\left(1-2\varphi\right)\exp\left(-\left[\mathrm{erf}^{-1}\left(1-2\varphi\right)\right]^{2}\right)$
(10) $\displaystyle\textrm{and}\quad\varsigma^{2}$
$\displaystyle=\frac{1}{2\pi\epsilon^{2}}\exp\left(-2\left[\mathrm{erf}^{-1}\left(1-2\varphi\right)\right]^{2}\right),$

where $\mathrm{erf}^{-1}(\cdot)$ denotes the inverse of the error function.
Note that we choose $\mu$ to be constant herein with respect to $z_{o}$ and
$\theta_{o}$. The insult propensity $\varphi$ corresponds to the fraction of
$\vartheta^{\ast}$ values greater than 0.5, while $\epsilon$ corresponds to
the slope of the cumulative distribution function (CDF) of $\vartheta^{\ast}$
at $\vartheta^{\ast}=0.5$ (see subsection B.1 for the derivation of Equation
10). In practice, to ensure stability of the finite element simulations that
follow, we constrain $\vartheta^{\ast}$ at the vessel boundaries to be low
(e.g., two standard deviations below the mean), such that the censored insult
profile values are zero at the boundaries. Note that $\vartheta^{\ast}$ is
discretized on the finite element mesh of the vessel; thus, the nodal values
$\vartheta_{i}^{\ast}$ jointly follow (and can thus be sampled from) a
multivariate Gaussian distribution with mean vector $\bm{\mu}=\mu\mathbf{1}$
and covariance matrix $\bm{\Sigma}$, where
$\Sigma_{ij}=\Sigma_{ji}=\kappa(z_{o,i},\theta_{o,i},z_{o,j},\theta_{o,j})$.
Partitioning the mesh into the set of interior nodes $a$ and the set of
boundary nodes $b$, it is straightforward to condition the distribution of
$\vartheta_{a}^{\ast}$ on the enforced value of $\vartheta_{b}^{\ast}$ using

$\displaystyle\bm{\mu}_{a}^{\prime}$
$\displaystyle=\bm{\mu}_{a}+\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\left(\vartheta_{b}^{\ast}\bm{1}-\bm{\mu}_{b}\right)=\mu+\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\left(\vartheta_{b}^{\ast}-\mu\right)\bm{1}$
(11) $\displaystyle\bm{\Sigma}_{aa}^{\prime}$
$\displaystyle=\bm{\Sigma}_{aa}-\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}$
$\displaystyle\bm{\mu}_{b}^{\prime}$
$\displaystyle=\vartheta_{b}^{\ast}\bm{1}$
$\displaystyle\bm{\Sigma}_{ab}^{\prime}$
$\displaystyle=\bm{0},\quad\bm{\Sigma}_{ba}^{\prime}=\bm{0},\quad\bm{\Sigma}_{bb}^{\prime}=\bm{0}$

instead for sampling purposes. After $\vartheta_{i}^{\ast}$ is sampled from
$\mathcal{N}\left(\bm{\mu}^{\prime},\bm{\Sigma}^{\prime}\right)$, we perform a
CDF/inverse-CDF transformation, so that the overall distribution of
$\vartheta^{\ast}$ values in each random instance of $\vartheta_{i}^{\ast}$
matches the desired $\mathcal{N}\left(\mu,\varsigma^{2}\right)$ (see Appendix
B for more details). Specifically,

$\left(\vartheta_{i}^{\ast}\right)^{\prime}=\Phi^{-1}\left(F\left(\vartheta_{i}^{\ast}\right);\mu,\varsigma^{2}\right),$
(12)

where $F$ is the CDF of the generated random field values (approximated via
kernel density estimation) and $\Phi^{-1}$ is the inverse CDF (i.e., quantile
function) of the normal distribution with mean $\mu$ and variance
$\varsigma^{2}$. Finally, the insult field values are censored using
$\vartheta_{i}=\min\left(\max\left(\left(\vartheta_{i}^{\ast}\right)^{\prime},0\right),1\right)$.
The randomly generated insult profiles (i.e., loss of elastic fiber integrity
or loss of mechanosensing) correspond in turn to random patterns of dilatation
along the vessel (Figure 3), which are afterward used as input data for the
DeepONet model to perform an inverse prediction of the insult profile. In
subsection B.2, we present a sensitivity analysis to demonstrate how
$\varphi$, $\epsilon$, $L_{\theta}$, and $L_{z}$ jointly control the size,
shape, and appearance of randomly generated insults with high precision.

Figure 3: 3D and axial views of representative dilatation maps of TAAs
initiated with prescribed insult profiles from mapped Gaussian random fields.
Circumferential length scale $L_{\theta}=1.5$ mm and axial length scale
$L_{z}=2.0$ mm with perturbation boundary softness $\epsilon=0.6$. In all
cases, $\varphi=15$% of the total area is perturbed, and the maximum insult
applied is 60% loss in elastic fiber integrity.

####  2.2.3 Insult contributors

Two types of contributors to TAAs are considered: loss of elastic fiber
integrity (analogous to a loss of functional elastin) and dysfunctional
mechanosensing (representing an impaired ability of cells to sense changes in
intramural stress). Insult values ranging from mild to severe are simulated
for each of these functions (Table 1).

##### Elastic fiber integrity

Functional elastin plays a crucial role in the elastic energy storage
capability of the aorta (i.e., compliance and resilience), and genetic
disorders such as Marfan syndrome associate with loss of elastic fiber
integrity from fragmentation or degradation. This is modeled with a user-
defined reduction in the baseline value of the material parameter $c^{e}$,
chosen to achieve a loss in mechanical properties consistent with previous
biaxial testing of aortas from _Fbn1 mgR/mgR_ mouse models for Marfan syndrome
[13, 11, 2].

##### Dysfunctional mechanosensing

Mechanical homeostasis of the aorta is regulated by intramural cells that
sense the local microenvironment, primarily through integrin binding sites.
When these cells cannot accurately detect deviations in stress from their
homeostatic levels, whether through disruptions in extracellular matrix or
impaired actomyosin activity, they tend to drive maladaptations within the
aorta. This is represented with the parameter $\delta\in[0,1]$, with 0
indicating perfect mechanosensing, and the modified expression for deviations
in intramural stress $\Delta\sigma=((1-\delta)\sigma-\sigma_{o})/\sigma_{o}$
[14, 33].

###  2.3 Finite element modeling

####  2.3.1 Normotensive conditions

We employ a 3D finite element model with $1\times 20\times 20$ mesh of 27-node
quadratic hexahedral elements in the radial, circumferential (azimuthal), and
axial directions. Displacement of the ends of the vessel is fixed in the axial
direction, then the vessel is pressurized at normotensive conditions via a
traction boundary condition on the inner surface and axially (pre-)stretched
uniformly according to previous experimental data. The open-source software
FEBio (febio.org) running a custom constrained mixture model plugin enables
both G&R computations and isochoric hyperelasticity computations in which G&R
is arrested [38]. Simulations are performed in the following stages:

  1. 1.

Initialization: uniform pressurization and pre-stretch of the cylindrical
vessel without insult to achieve the initial homeostatic state.

  2. 2.

G&R: computation of growth and remodeling in response to the gradually applied
insult, at fixed pressure.

  3. 3.

Hyperelasticity: post-G&R computation of vessel deformation at normotensive
diastolic and systolic pressures [39].

For analytically defined TAA cases, all combinations of insult profiles with
variations in axial and circumferential location and extent were simulated for
five levels of insult severity and for both insult types, yielding a total of
550 simulations at normotensive conditions. For randomly generated TAAs, 10
unique profiles sharing the same shape parameters are used to yield 100 cases.

####  2.3.2 Superimposed hypertension

Uncontrolled hypertension (elevated blood pressure) is a critical determining
factor in TAA growth [39, 4]. Thus, each insult type in subsubsection 2.2.3 is
also simulated with superimposed hypertension, modeled by a gradual increase
in the intraluminal pressure, concurrent with the prescribed insult. As a
modification to the previously described simulation pipeline, hypertension
cases are simulated as follows:

  1. 1.

Initialization: uniform pressurization and pre-stretch of the cylindrical
vessel without insult to achieve the initial homeostatic state.

  2. 2.

G&R: computation of growth and remodeling in response to the gradually applied
insult, concurrent with gradual increases in pressure by 33%.

  3. 3.

Hyperelasticity: post-G&R computation of vessel deformation at hypertensive
diastolic and systolic pressures [39].

Cases of hypertension highlight how even a mild insult degree with modest
dilatation under normotensive conditions can be exacerbated by increased
pressure, yielding a maximum dilatation comparable to that produced by a
severe insult under normotensive conditions (Table 1). 550 additional
simulations with hypertensive conditions are simulated for analytically
defined profiles. Hypertensive conditions were not considered for randomly
generated insults.

Table 1: G&R and mechanobiological insult profile parameters for normotensive
and hypertensive simulations with analytically defined insults (Equation 7).
To avoid boundary effects, $z_{apex}=7.5$ mm only for $z_{od}=4.0$ mm. For
randomly generated profiles, loss of elastic fiber integrity is varied from 6%
– 60%, and dysfunctional mechanosensing ranges from 2.5% – 25% at normotensive
conditions.

Parameter |  Variable |  Value
---|---|---
Axial characteristic width (mm) |  $z_{od}$ |  2.0, 3.0, 4.0
Circumferential characteristic width (deg) |  $\theta_{od}$ |  $20$, $100$, $180$, $260$, $360$
Axial placement (mm) |  $z_{apex}$ |  6.0, 7.5, 9.0
Circumferential placement (deg) |  $\theta_{apex}$ |  $0$, $90$, $180$, $270$
|  |  Normotensive |  Hypertensive (33%)
Loss of elastic fiber integrity |  |  5.95% – 59.5% |  4.75% – 47.5%
Dysfunctional mechanosensing |  |  1.84% – 18.4% |  1.08% – 10.8%
G&R pressure (mmHg) [33] |  |  105 |  $105\rightarrow 140$
Diastolic pressure (mmHg) [39] |  |  99 |  129
Systolic pressure (mmHg) [39] |  |  121 |  172

####  2.3.3 Post-processing training and testing data

In the Initialization stage, the straight cylindrical vessel with inner radius
$r_{o}$ is uniformly pressurized and pre-stretched according to normotensive
_in vivo_ conditions in one load step while maintaining the same geometry
(Figure 4a). In the G&R stage, an insult (analytically defined or randomly
generated) is gradually applied over 10 subsequent load steps to compute the
evolved post-G&R geometry (Figure 4b) with spatially varying inner radius
$r_{h}$. The G&R evolution is then arrested, and in the Hyperelasticity stage
two additional load steps are simulated in which the luminal pressure is
adjusted to diastolic and systolic values [39] to evaluate the local
distension at diastole $\Lambda_{D}=r_{dia}/r_{o}$ and systole
$\Lambda_{S}=r_{sys}/r_{o}$, where $r_{dia}$ and $r_{sys}$ are the inner
radius at diastolic and systolic pressure, respectively (Figure 4c,d).

We train the DeepONet surrogate model using two quantities of interest from
the finite element-based synthetic data: dilatation (i.e., aneurysmal
enlargement) and distensibility (i.e., cyclic strains). For dilatation, we
utilize a map of $\Lambda_{D}$ for each mesh node along the inner surface of
the geometry. For distensibility, in contrast to the diameter- and pressure-
based “distensibility” defined in previous literature [39], we consider the
actual distensibility $\mathcal{D}=(\Lambda_{S}-\Lambda_{D})/\Lambda_{D}$ for
all nodes on the inner surface. These dilatation and distensibility maps are
provided to the neural network in the form of $z$–$\theta$ maps (Figure 4e) as
well as nondimensionalized grayscale images (Figure 4f). Approximately 10% of
the total number of datasets are used to validate the trained model, and the
remaining are used for training.

Figure 4: Schematic for the generation of synthetic training data for TAAs
arising from localized mechanobiological insults. (a) The initial vessel
geometry with inner radius $r_{o}$ and prescribed insult is evolved over
multiple load steps to compute (b) the final post-G&R geometry, in which the
homeostatic inner radius $r_{h}$ can be calculated with respect to the updated
centerline. The luminal pressure is then adjusted to (c) diastolic and (d)
systolic conditions to evaluate the local distension at diastole $\Lambda_{D}$
and systole $\Lambda_{S}$, respectively. (e) The distensibility is then
computed with the relation
$\mathcal{D}=(\Lambda_{S}-\Lambda_{D})/\Lambda_{D}$. Together, these steps
yield the dilatation $\Lambda_{D}$ and distensibility $\mathcal{D}$ maps
(shown in the flattened $z$–$\theta$ plane). (f) Finally, the maps are
converted into nondimensionalized grayscale intensity maps and contrast
enhanced. For illustrative purposes, this procedure is shown for an
analytically defined insult (Equation 7) under normotensive conditions only.

###  2.4 DeepONet

In this section, we describe the architecture of the surrogate model developed
within the framework of DeepONet to predict mechanobiological insult profiles.
The conventional unstacked DeepONet architecture consists of two deep neural
networks (DNNs): one encodes the input function at fixed sensor points (branch
net) while the other accounts for the locations of the output function (trunk
net). The branch network input is customizable in a generalized setting and
can take the shape of the physical domain, the initial or boundary conditions,
constant or variable coefficients, source terms, and so on, as long as the
input function is discretized at sensor locations $n_{sen}$. For a regularly
spaced discretization of the input function, a convolutional neural network
(CNN) can be used as the branch net, while for a sparse representation, one
may also consider a feed-forward neural network (FNN), or even a recurrent
neural network (RNN) for sequential data. A standard practice is to use an FNN
in the trunk network to take into account the low dimensions of the evaluation
points. The mathematical foundation of DeepONet is provided in Appendix C, and
the interested reader is referred to [31] for additional details.

The cases considered in this work consist of mechanobiologically equilibrated
TAAs developed from an initially cylindrical vessel caused by loss of elastic
fiber integrity or compromised mechanosensing. We generate cases under
normotensive and hypertensive conditions and consider dilatation and
distensibility maps to train the DeepONet. We consider two different network
architectures for the branch network, FNN and CNN, with the aim to identify a
preferred network architecture. With an FNN in the branch network, we show
that even with sparse dilatation and distensibility information, the DeepONet
predicts the insult profile with high accuracy. Subsequently, we use a CNN to
approximate the input function with grayscale images of the dilatation and
distensibility. Both approaches are judged based on two empirical
observations:

  * 1.

the ability of the models to generalize unseen test cases, and

  * 2.

the ability of these models to generalize noisy input functions.

####  2.4.1 FNN-based architectures

The DeepONet in this work considers multiple branch networks to account for
the dilatation map at diastole ($\Lambda_{S}$) and the distensibility map
($\mathcal{D}$), as well as whether the loading condition is normotensive or
hypertensive. The five branch networks $\mathbf{U}^{i}$ for $i=1,2,\ldots,5$
are formulated as FNNs with the following description:

  * 1.

$\mathbf{U}^{1}$: value of $\Lambda_{D}$ at $n_{sen}$.

  * 2.

$\mathbf{U}^{2}$: location of the maximum value in $\Lambda_{D}$.

  * 3.

$\mathbf{U}^{3}$: value of $\mathcal{D}$ at $n_{sen}$.

  * 4.

$\mathbf{U}^{4}$: location of the minimum value in
$\mathcal{D}=(\Lambda_{S}-\Lambda_{D})/\Lambda_{D}$.

  * 5.

$\mathbf{U}^{5}$: binary network to account for hypertension: _normotensive:
0, hypertensive: 1_.

We perform two experiments to decide the optimal number of sensors for
$\mathbf{U}^{1}$ and $\mathbf{U}^{3}$, detailed below.

Experiment 1: $n_{sen}=5\times 5$ single-spaced sensor locations.

In this experiment, a $5\times 5$ lattice of points is employed for
$\mathbf{U}^{1}$ and $\mathbf{U}^{3}$, centered at the location of maximum
dilatation $\mathbf{U}^{2}$ and the location of minimum distensibility
$\mathbf{U}^{4}$. The locations of the sensors are selected such that the
input data is constrained within a fairly localized neighborhood around the
local extremum, dictated by the spacing of the nodes in the finite element
mesh (in other words, the sensor locations correspond to a uniform $5\times 5$
grid of adjacent single-spaced nodes). The arrangement of sensor points is
shown in Figure 5a.

Experiment 2: $n_{sen}=9$ double-spaced sensor locations.

In this experiment, only $9$ sensor locations are employed for for
$\mathbf{U}^{1}$ and $\mathbf{U}^{3}$; however, locations are distributed with
double the spacing in the axial and azimuthal directions compared to the
single-spaced $5\times 5$ locations, encompassing a neighborhood roughly 4
times the size of that in Experiment 1. As above, the sensors for
$\mathbf{U}^{2}$ and $\mathbf{U}^{4}$ are centered at the locations of maximum
dilatation and minimum distensibility, respectively (Figure 5b).

Figure 5: Depictions of the sensor locations for subsubsection 2.4.1 within a
representative TAA case (10.1% loss of mechanosensing, $z_{od}=2$ mm,
$\theta_{od}=20^{\circ}$, $z_{apex}=6$ mm, $\theta_{apex}=270^{\circ}$). (a)
$n_{sen}=5\times 5$ single-spaced sensor locations on $\Lambda_{D}$ and
$\mathcal{D}$ for Experiment 1. (b) Double-spaced sensor locations $n_{sen}=9$
for Experiment 2. The scalar distances $r_{1}$ and $r_{2}$ are computed from
the origin to the maximum dilatation location and the minimum distensibility
location, respectively. Figure 6: Schematic representation of the DeepONet for
subsubsection 2.4.1. The branch network takes as input the dilatation
$\Lambda_{D}$ and distensibility $\mathcal{D}$ at the sensor locations
$n_{sen}$ (branch net: $\mathbf{U}^{1}$ and $\mathbf{U}^{3}$, respectively)
and the locations of maximum dilatation and minimum distensibility (branch
net: $\mathbf{U}^{2}$ and $\mathbf{U}^{4}$, respectively). Additionally, blood
pressure information is fed to the network using a branch net:
$\mathbf{U}^{5}$ as normotensive (input $=0$) or hypertensive (input $=1$).
Each branch net returns features embedding
$[b^{i}_{1},b^{i}_{2},\ldots,b^{i}_{q}]^{\mathrm{T}}\in\mathbb{R}^{q}$,
$i={1,2,\ldots,5}$ as output, and the dot product of these outputs is computed
to obtain $[b_{1},b_{2},\ldots,b_{q}]^{\mathrm{T}}$. The trunk net takes the
continuous coordinates $\bm{y}=(\cos\hat{\theta},\sin\hat{\theta},\hat{z})$ as
inputs and outputs a features embedding
$[t_{1},t_{2},\ldots,t_{q}]^{\mathrm{T}}\in\mathbb{R}^{q}$. The features
embeddings of the branch and trunk networks are merged via an element-wise dot
product to output the solution operator $\mathbb{G}_{\bm{\theta}}$ with
learnable parameters $\bm{\theta}$. Minimization of the loss function
$\mathcal{L}$ determines best-fit parameters $\bm{\theta}^{\ast}$ that enable
estimation of the insult profile.

The trunk network (FNN) considers the locations of evaluation in a cylindrical
coordinate system (i.e.,
$\bm{y}=\\{\bm{y}_{1},\bm{y}_{2},\ldots,\bm{y}_{p}\\}=\\{(\hat{\theta}_{1},\hat{z}_{1}),(\hat{\theta}_{2},\hat{z}_{2}),\ldots,(\hat{\theta}_{p},\hat{z}_{p})\\}$).
Encapsulating previous knowledge into the architecture of DeepONet often
improves the generalization of the output. Here, we replace the $\hat{\theta}$
component of the trunk net with appropriate basis functions in Cartesian
coordinates. Thus, the trunk net input is modified as
$\bm{y}_{i}=(\cos\hat{\theta}_{i},\sin\hat{\theta}_{i},\hat{z}_{i})$. A
schematic representation of the framework is shown in Figure 6.

####  2.4.2 CNN-based architecture

In this section, we discuss the DeepONet architecture considering grayscale
images, converted from the $\Lambda_{D}$ and $\mathcal{D}$ maps and contrast
enhanced, as input functions for the branch network. This approach considers 3
branch networks with the following descriptions:

  * 1.

$\mathbf{U}^{1}$: grayscale image for $\Lambda_{D}$ using a CNN.

  * 2.

$\mathbf{U}^{2}$: grayscale image for $\mathcal{D}$ using a CNN.

  * 3.

$\mathbf{U}^{3}$: binary network to account for hypertension using an FNN.

Note that, in contrast to the sparse sensor point arrays, we use full-field
images for $\mathbf{U}^{1}$ and $\mathbf{U}^{2}$, which have a resolution of
$21\times 20$ pixels. A schematic representation of the DeepONet implemented
with this framework is shown in Figure 7.

Figure 7: Schematic representation of the DeepONet for subsubsection 2.4.2.
Each branch net is a CNN that takes as inputs the grayscale images of
dilatation $\Lambda_{D}$ and distensibility $\mathcal{D}$ (branch nets:
$\mathbf{U}^{1}$ and $\mathbf{U}^{2}$, respectively). Blood pressure
information is fed to the FNN using a branch net: $\mathbf{U}^{3}$ as
normotensive (input $=0$) or hypertensive case (input $=1$). The outputs of
the branch and trunk networks are merged into the solution operator
$\mathbb{G}_{\bm{\theta}}$, and minimization of the loss function
$\mathcal{L}$ enables estimation of the insult profile.

##  3 Results

###  3.1 Dilatation and distensibility in TAA

Locally applied insults in elastic fiber integrity and mechanosensing both
lead to the development of dilatations, as well as aneurysms (defined as a
1.5-fold increase in normalized diameter from baseline) in cases of severe
insults of elastic fiber integrity. In cases with less circumferential
involvement of insult ($\theta_{od}<180^{\circ}$), the unaffected regions of
the vessel help to attenuate dilatation within the insult area. Additionally,
each factor is modeled in combination with superimposed hypertension (elevated
blood pressure), emphasizing how even a mild insult with modest dilatation in
normotensive conditions can be exacerbated by the presence of additional risk
factors. In hypertensive cases, the dilatation also increases significantly,
including in regions where there is no insult applied.

Overall, increased dilatation associated with mechanobiological insults
correlate well with decreases in distensibility. However, there are distinct
differences between distensibility maps resulting from loss of elastic fiber
integrity and those resulting from dysfunctional mechanosensing. Whereas the
location of maximum dilatation and minimum distensibility colocalize in
elastic fiber integrity loss, the maximum dilatation in mechanosensing loss
occurs on the opposite side of the vessel from the minimum distensibility,
although the distensibility also decreases at the location of maximum
dilatation. These findings confirm the need to evaluate both metrics for
robust training of the DeepONet in all cases, especially in experiments
considering multiple types of insults.

###  3.2 Performance of the surrogate models

The effectiveness of the developed surrogate models is demonstrated through
several experiments, put forward in this section. To evaluate performance, we
compute the $L_{2}$ relative error of predictions, and we report its mean and
standard deviation based on five independent training trials. In all cases
presented here, the DeepONet is trained using a combination of Adam [40] and
L-BFGS optimizers [41]. The implementation is carried out using the TensorFlow
framework [42]. Throughout all examples, we initialize the weights and biases
of the DeepONet using Xavier initialization. The experiments carried out in
this work are listed in Table 2.

Table 2: Descriptions of the datasets considered for the experiments carried out in this work. Insult profiles are either analytically defined or randomly generated losses of elastic fiber integrity or mechanosensing, and pressure conditions are either normo- or hypertensive, listed in parentheses. Case # | Description
---|---
Case 1 | Analytically defined elastic fiber integrity (normotensive)
Case 2 | Analytically defined mechanosensing (normotensive)
Case 3 | Analytically defined elastic fiber integrity or mechanosensing (normotensive)
Case 4 | Analytically defined mechanosensing (normotensive & hypertensive)
Case 5 | Analytically defined elastic fiber integrity or mechanosensing (normotensive & hypertensive)
Case 6 | Randomly generated elastic fiber integrity or mechanosensing (normotensive)

The experiments establish the accuracy of the three surrogate models when
trained with data generated for analytically defined insult profiles (Cases
1–5) and for randomly generated insult profiles (Case 6). In Figure 8, we
present the performance of the three surrogate models for representative
analytically defined and randomly generated insults in elastic fiber integrity
(Cases 1 & 6). For each architecture, we compute an error between the true and
predicted insult profiles by normalizing the absolute error by the maximum
insult value throughout the $z$–$\theta$ plane. While each DeepONet
architecture is able to predict the insult profile within 5% error, the design
with $5\times 5$ sensors tends to exhibit the greatest prediction errors,
whereas the design based on grayscale images provides the most accurate
prediction. This discrepancy is most clearly observed in Figure 8b.

Figure 8: Performance of the surrogate models for the prediction of
analytically defined and randomly generated insult profiles. Axial–azimuthal
views of representative true and predicted (a) analytically defined and (b)
randomly generated insult profiles are predicted using $5\times 5$ sensors,
$9$ sensors, and full-field grayscale maps. The insult prediction error,
normalized by the maximum insult value over the axial–azimuthal domain, is
shown for each architecture.

In Table 3, we show the relative $L_{2}$ error of the surrogate models,
averaged over the $z$–$\theta$ domain, for the testing dataset compromised of
mechanosensing insults under normotensive conditions (Case 2). Additionally,
we show the total number of trainable (i.e., learnable) parameters for each
surrogate model, which are the weights and biases of the networks optimized
using back-propagation during the training process. In general, the total
number of parameters learned during the training process for a CNN is less
than that for an FNN. CNNs are very effective in reducing the number of
parameters without losing the quality of models, as the learnable filters of
the network encode the information to reduce the high-dimensional input data
(here, grayscale images) to a lower latent dimension. This is also seen in our
results, which show that the network with full-field grayscale images requires
the lowest number of parameters while simultaneously achieving the highest
predictive accuracy. Finally, we evaluate the performance of the surrogate
model trained with grayscale images on a noisy dataset that consists of the
original testing dataset with $5\%$ added uncorrelated Gaussian noise.
Overall, among the surrogate models considered in this Case, the CNN-based
design again demonstrates improved computational efficiency, prediction
accuracy, and robustness against noise.

Table 3: Relative $L_{2}$ error and the number of trainable parameters of the three surrogate models for $N=545$ training data and $N_{*}=45$ testing data consisting of compromised mechanosensing under normotensive conditions. The noise is added in the inputs of the _testing_ dataset. Method | # parameters |  Relative $L_{2}$ error
---|---|---
Testing data |  $+5\%$ noise
$5\times 5$ sensors | $234800$ | $3.906\pm 0.004\%$ | $44.27\pm 0.008\%$
$9$ sensors | $230784$ | $2.712\pm 0.007\%$ | $32.52\pm 0.009\%$
Grayscale images | $141180$ | $2.340\pm 0.010\%$ | $8.2\pm 0.020\%$

In Table 4, we report the relative $L_{2}$ errors of the three models, along
with the robustness of each model to $5\%$ added uncorrelated Gaussian noise.
Error plots for three representative test cases are shown in Figure 9. Columns
1 and 2 show the 3D geometry and the unfolded axial–azimuthal view of the
aorta, respectively, colored by the insult profile (ground truth). Columns 3–5
depict the prediction errors associated with using $5\times 5$ sensors, $9$
sensors, or grayscale images to train the network, respectively. Error is
normalized with respect to the maximum insult value of each case. These
findings further reinforce the improvements of the CNN-based approach in a
variety of testing scenarios varying by insult type, pressure conditions, and
insult profile generation method.

Table 4: Relative $L_{2}$ error of the three surrogate models for with $N$ training data and $N_{*}$ testing data for the cases listed in Table 2. The noise is added in the grayscale full-field image inputs of the _testing_ dataset. Case # | $N$ | $N_{*}$ |  Relative $L_{2}$ error
---|---|---|---
$5\times 5$ sensors |  $9$ sensors | Grayscale images |  $+5\%$ noise
Case 1 | $545$ | $45$ | $4.79\pm 0.008\%$ | $4.29\pm 0.008\%$ | $3.458\pm 0.011\%$ | $7.80\pm 0.03\%$
Case 2 | $545$ | $45$ | $3.91\pm 0.004\%$ | $2.71\pm 0.007\%$ | $2.340\pm 0.010\%$ | $8.20\pm 0.020\%$
Case 3 | $500$ | $90$ | $7.10\pm 0.001\%$ | $6.75\pm 0.003\%$ | $5.948\pm 0.006\%$ | $10.69\pm 0.015\%$
Case 4 | $500$ | $90$ | $2.88\pm 0.005\%$ | $2.62\pm 0.001\%$ | $2.280\pm 0.005\%$ | $7.30\pm 0.01\%$
Case 5 | $720$ | $180$ | $3.65\pm 0.005\%$ | $2.65\pm 0.003\%$ | $2.540\pm 0.018\%$ | $7.18\pm 0.012\%$
Case 6 | $90$ | $10$ | $7.74\pm 0.011\%$ | $7.20\pm 0.020\%$ | $2.292\pm 0.002\%$ | $15.96\pm 0.029\%$
Figure 9: Error plots for three representative cases of testing samples from
the indicated experiments. Columns 1 and 2 show the 3D geometry and the
unfolded axial–azimuthal view of the aorta, respectively, colored by the
insult profile (ground truth). Columns 3–5 depict the prediction error of
Experiment 1 (using $5\times 5$ sensors), Experiment 2 (using $9$ sensors),
and using grayscale images to train the network, respectively. Error is
normalized with respect to the maximum insult value of each case. The top row
corresponds to a case of analytically defined normotensive elastic fiber
integrity insult (Case 1), the middle row shows an analytically defined
hypertensive mechanosensing insult (Case 4), and the bottom row represents a
randomly generated normotensive elastic fiber integrity insult (Case 6).

##  4 Discussion

###  4.1 Insult profile prediction with sensor point- and image-based
approaches

This study utilizes a constrained mixture model for arterial growth and
remodeling integrated with a DeepONet to predict factors contributing to TAA.
The main challenge in applying deep learning based frameworks to predict TAA
enlargement is the availability of limited information, both in terms of
number of samples as well as information per sample. Additionally, due to the
high complexity of the mechanobiology in patient-specific predictions of TAA
(e.g., age, hypertension, diabetes, prescribed medications), making accurate
predictions with even high-fidelity physical models remains a challenge. To
enable accurate predictions, we have proposed and evaluated three frameworks
of DeepONet that are distinguished based on the information available for
training the network. The results shown in Table 4 indicate that the proposed
frameworks indeed provide effective predictions for both analytically defined
and randomly generated insult profiles. Our observations are summarized as
follows:

  1. 1.

It has been shown that aortic geometry alone is not sufficient to predict TAA
progression [21]. Many have thus sought to incorporate additional information
such as biomechanical properties and patient-level variables to improve
predictive capability. We find that accurate prediction of the insult profile
can be achieved with the inclusion of biomechanical properties along with
dilatation and distensibility fields generated from measurements of diastolic
and systolic phases of the cardiac cycle, similar to the approaches used by
other groups [23, 6].

  2. 2.

Predicting insult profiles with information at $5\times 5$ sensor locations is
sufficiently accurate for cases with $\theta_{od}<260^{\circ}$ and relatively
small $z_{od}$. However, with a wider $\theta_{od}$ and broader $z_{od}$,
limited information within a single-spaced neighborhood around maximum
dilatation and minimum distensibility is insufficient; hence, for such cases,
the relative error increases.

  3. 3.

Overcoming the limitations of $5\times 5$ sensor locations, the double-spaced
arrangement of $9$ sensors more accurately estimates the insult profile within
a wider neighborhood compared to that of the $5\times 5$ sensor array. This
observation is in line with the reduced prediction errors reported in Table 4
(fourth and fifth column), which depicts this improvement in accuracy. This
suggests that the range of variation within the network inputs captured by the
sensor point domain, rather than the absolute number of sensors, is a greater
determining factor in the predictive capability of the network.

  4. 4.

FNNs have proven to work well with limited information and also be robust to
noisy testing inputs. However, in the case of extremely sparse information,
the network fails to generalize well for noisy inputs and reports a relative
$L_{2}$ error of $44.27\pm 0.008\%$ and $32.52\pm 0.009\%$ for $5\times 5$ and
$9$ sensors, respectively, when tested with $5\%$ Gaussian noise added to the
testing inputs of Case 2 (compromised mechanosensing) (last column in Table
3).

  5. 5.

The robustness of the different network architectures is evident in cases of
randomly generated insult profiles. As seen in the third and the fourth
columns of Figure 10, networks trained on sparse sensor point information
often fail to achieve accurate predictions if there is more than one location
of localized dilatation or if the dilatation boundary is irregular. In these
cases, we observe that the model with $9$ double-spaced sensor locations
achieves better prediction accuracy than the $5\times 5$ single-spaced sensor
locations.

  6. 6.

The caveats of the FNN-based frameworks (subsubsection 2.4.1) motivate the
development of a CNN-based framework, which takes as inputs grayscale images
of dilatation and distensibility to more accurately predict the insult
profile. The predictive accuracy of the framework proposed in subsubsection
2.4.2 is clearly depicted in Table 4 (sixth column) and Figure 9.
Additionally, the model is more noise-tolerant than the sparse sensor-based
frameworks (last column of Table 4).

  7. 7.

The computational cost of a network is directly related to the number of
trainable parameters. Neural networks employ back-propagation algorithms to
tune the network parameters, while trying to reflect the best-fit solution to
the training data. Therefore, an FNN trained on real-life images is often
computationally expensive. The choice of a CNN to train on grayscale images
benefits the model not only in terms of predictive accuracy but also
computational efficiency, as this framework requires fewer learnable
parameters compared to the models employing FNNs (second column of Table 3).

  8. 8.

All of the frameworks proposed are capable of handling patient-specific
datasets beyond dilatation and distensibility. In this work, we have
demonstrated that an additional FNN in the branch net is sufficient to handle
information like hypertension. This framework can be easily extended to handle
further information, such as preexisting genetic or pharmacological
conditions, without any additional cost of training the network.

Figure 10: Error plots for two representative cases of loss of elastic fiber
integrity under normotensive conditions for randomly generated insult profiles
(Case 6). Columns 1 and 2 show the 3D geometry and unfolded view of the aorta
with aneurysm (ground truth), respectively, and columns 3–5 depict the
normalized prediction error of Experiment 1 (using $5\times 5$ sensors),
Experiment 2 (using $9$ sensors), and using grayscale images to train the
network, respectively.

###  4.2 Limitations

There are also restrictions to the frameworks proposed and experiments
selected in this study, which will be addressed in future work. The proposed
framework is capable of predicting the insult profile for a given set of
dilatation and distensibility information (that is, for a given time). This
capability represents a scientifically important inverse problem, a class of
problems that are typically ill-posed and yet were tractable herein.
Nevertheless, even though the three surrogate models are capable of handling
patient-specific information, they are not yet capable of determining real-
time TAA evolution, which is a significant unmet need in the assessment of TAA
for diagnosis, prognosis, and follow-up, and also in the determination of
rupture risk. The study has been carried out on synthetic data based on finite
element simulations, which are available in abundance. Future training of the
model on limited clinical information will be challenging in terms of the
number of data points available. In such situations, transfer learning may
need to be integrated with domain adaptation models to leverage the benefits
of a model trained on synthetic data and used for predictions on a different
dataset. Finally, to address the challenge of accounting for additional
genetic and biomechanical factors, we plan to incorporate models with improved
hierarchical relationships (e.g., capsule networks) and quantification of
uncertainties in the local tissue mechanical properties [43]. Multi-fidelity
DeepONet approaches could also be designed, utilizing an additional low-
fidelity dataset to drastically reduce the necessary high-fidelity dataset. To
make use of the available historical data, we plan to use a Bayesian framework
to learn from functional priors (_in silico_) and make posterior estimations
(_in vivo_). Such an approach would reduce the use of high-fidelity data
required for training, potentially enabling improved prognosis of aneurysm
growth and rupture risk in a clinically relevant time frame.

###  4.3 Summary

We have developed a novel framework to predict TAA pathology by integrating a
constrained mixture model for arterial growth and remodeling with a deep
neural network surrogate model. 3D finite element simulations of TAA
development arising from randomly distributed losses of elastic fiber
integrity and dysfunctional mechanosensing provided synthetic training data
for the surrogate model, capable of predicting insult profiles from dilatation
and distensibility information. Finally, we demonstrated improved performance
using convolutional neural networks in our DeepONet construction. This
framework can ultimately be applied to construct patient-specific profiles for
aneurysm growth, which will provide critical information to contextualize the
predicted mechanobiological insult and forecast the short-term evolution of
TAA imaged at one time point in the clinic. Characterization of this
progression could play an integral role in determining future patient risk and
in designing improved therapeutic interventions.

##  Appendix A Constrained mixture model

###  A.1 Constitutive properties

The biomechanical behaviors of the structurally significant constituents of
the aortic wall $\alpha=e,c,m$ (elastin-dominated matrix, collagen-dominated
matrix, and smooth muscle cells, respectively) at current time $s$ are
described by strain energy density functions $W_{R}^{\alpha}$, with subscript
$R$ denoting the reference configuration. The constitutive relations are
weighted by their respective volume fractions $\phi^{\alpha}$. We define
$\mathbf{F}$ as the deformation gradient tensor for the tissue from reference
to current configurations and $J=\det\mathbf{F}$.

Elastin-dominated matrix. A neo-Hookean model is used to describe the
isotropic elastin-dominated mechanical response ($\alpha=e$), defined as

$W_{R}^{e}(s)=\phi_{R}^{e}(s)\hat{W}^{e}(s)=\phi_{R}^{e}(s)\frac{c^{e}}{2}(\mathrm{tr}(\mathbf{F}^{e\mathrm{T}}(s)\mathbf{F}(s))-3),$
(13)

where $\phi_{R}^{e}(s)=J\phi^{e}(s)$ is the referential mass fraction,
$\mathbf{F}^{e}=\mathbf{F}\mathbf{G}^{e}$, with the tensor $\mathbf{G}^{e}$
representing the elastin deposition stretch when produced, and $c^{e}$ is the
shear modulus, which can be estimated by fitting biomechanical testing data.
Because elastin does not turn over continuously in maturity, its natural
configuration is fixed and no convolution integral is required.

Collagen-dominated matrix and passive smooth muscle. For collagen fibers
($\alpha=c$) and smooth muscle fibers ($\alpha=m$), a four-fiber family model
is employed in combination with production and removal, giving the relation

$W_{R}^{c,m}(s)=\frac{1}{\rho}\int_{-\infty}^{s}m_{R}^{c,m}(\tau)q^{c,m}(s,\tau)\hat{W}^{c,m}(\lambda_{n(\tau)}^{c,m}(s))d\tau,$
(14)

where $\rho$ is the mass density of the tissue,
$\phi_{R}^{c,m}(s)=J\phi^{c,m}(s)$, $\lambda_{n(\tau)}^{c,m}(s)$ are the fiber
stretches (relative to their evolving natural configurations $n(\tau)$), and

$\hat{W}^{c,m}(\lambda_{n(\tau)}^{c,m}(s))=\frac{c_{1}^{c,m}}{4c_{2}^{c,m}}\big{(}\exp\big{[}c_{2}^{c,m}(\lambda_{n(\tau)}^{c,m}(s)-1)^{2}\big{]}-1\big{)},$
(15)

where $c_{i}^{c,m}$ are material parameters for each fiber family that can be
estimated similarly from biomechanical testing data. Collagen fibers are
classified into circumferentially, axially, and diagonally oriented (at angle
$\alpha_{0}$) populations with associated fractions $\beta^{\theta}$,
$\beta^{z}$, and $\beta^{d}$, respectively.

###  A.2 Mechanobiologically equilibrated constrained mixture model

In cases of TAA, G&R can be assumed to reach a quasi-static mechanobiological
equilibrium at time $s\gg 0$. Since production balances removal,
$m_{h}^{\alpha}=m_{Nh}^{\alpha}=k_{Nh}^{\alpha}\rho_{h}^{\alpha}$, and
Equation 14 reduces to

$W_{Rh}^{c,m}=\phi_{Rh}^{c,m}\hat{W}^{c,m}(G_{h}^{c,m}),$ (16)

where $G_{h}^{c,m}$ are the homeostatic deposition stretches of collagen and
smooth muscle fibers. Accordingly, rule-of-mixtures expressions can be used to
be used for stored energy ($W_{Rh}=\sum\phi_{Rh}^{\alpha}\hat{W}^{\alpha}$,
where $\phi_{Rh}^{\alpha}$ are the evolved constituent mass fractions).
Similarly, for the Cauchy stresses (Equation 4),

$\bm{\sigma}_{h}=-p_{h}\mathbf{I}+\sum_{\alpha}^{e,c,m}\phi_{h}^{\alpha}\hat{\bm{\sigma}}_{h}^{\alpha},$
(17)

where $\hat{\bm{\sigma}}_{h}^{\alpha}$ are the constituent-specific Cauchy
stresses, and $p_{h}$ is the equilibrated Lagrange multiplier for the quasi-
static G&R evolution.

We assume the turnover rates of smooth muscle and collagen from the initial to
equilibrated states are the equivalent, letting the smooth muscle-to-collagen
turnover ratio $\eta=(k^{m}K_{i}^{m})/(k^{c}K_{i}^{c})=1$
($i=\sigma,\tau_{w}$), which yields
$\rho_{Rh}^{m}/\rho_{o}^{m}=\rho_{Rh}^{c}/\rho_{o}^{c}$. Additionally,
assuming that the cardiac output, and thus blood flow rate $Q$, remains
constant and laminar, we let $\tau_{wh}/\tau_{wo}=Qa_{o}^{3}/(Qa_{h}^{3})$ and
$K_{\tau_{w}}^{\alpha}/K_{\sigma}^{\alpha}=0$. We compute the deviation in
intramural stress using $\sigma_{h}=\frac{1}{3}\mathrm{tr}(\bm{\sigma}_{h})$
and also allow collagen fibers to gradually reorient during G&R. Material
parameters for the model, shown in Table 5, are adapted from [33].

Table 5: Material parameters used in finite element simulations. Superscripts
$e,m,c$ denote elastin-dominated, smooth muscle, and collagen-dominated
matrix, respectively; super/subscripts $r,\theta,z,d$ denote radial,
circumferential, axial, and symmetric diagonal directions, respectively.
Subscript $o$ denotes the initial homeostatic state; subscripts
$i=\sigma,\tau_{w}$ denote intramural and wall shear stress related values,
respectively. For elastic fiber integrity insults, the baseline value of
$c^{e}$ is reduced by a prescribed percentage; for mechanosensing insults,
$\delta$ is increased. The reader is referred to [33] for further details on
these model parameters.

Inner radius, thickness, length | $r_{o},h_{o},l_{o}$ | 0.647 mm, 0.040 mm, 15 mm
---|---|---
Elastin, smooth muscle, collagen mass fractions | $\phi_{o}^{e},\phi_{o}^{m},\phi_{o}^{c}$ | 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_{0o}$ | $29.9^{\circ}$
Elastic material parameters | $c^{e},c_{1}^{m},c_{2}^{m},c_{1}^{c},c_{2}^{c}$ | 89.71 kPa, 261.4 kPa, 0.24, 234.9 kPa, 4.08
Deposition stretches | $G_{\theta}^{e},G_{z}^{e},G_{r}^{e},G^{m},G^{c}$ | 1.90, 1.62, $1/(G_{\theta}^{e}G_{z}^{e})$, 1.20, 1.25
Smooth muscle-to-collagen turnover ratio | $\eta$ | 1.0
Shear-to-intramural gain ratio | $K_{\tau_{w}}/K_{\sigma}$ | 0.00
Dysfunctional mechanosensing | $\delta$ | 0

##  Appendix B Random insult fields

###  B.1 Derivation of mean and variance

Marginally, the values of a GRF with constant mean are themselves normally
distributed with the same mean and with variance equal to the overall variance
of the GRF. Holding $\varphi$ and $\epsilon$ constant, we may thus solve for
the corresponding mean $\mu$ and variance $\varsigma^{2}$ explicitly. To
enforce the overall insult propensity $\varphi$, defined as the fraction of
$\vartheta^{\ast}$ (and thus $\vartheta$) values greater than 0.5, we evaluate
the Gaussian complementary CDF at 0.5 and solve for $\mu$ in terms of
$\varphi$ and $\varsigma$:

$\displaystyle\varphi$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\,\mathrm{erf}\left(\frac{0.5-\mu}{\varsigma\sqrt{2}}\right)$
(18) $\displaystyle\implies\mu$
$\displaystyle=\frac{1}{2}-\varsigma\sqrt{2}\,\mathrm{erf}^{-1}\left(1-2\varphi\right).$

To tune the softness of boundaries between normal and insult regions, we
prescribe the slope $\epsilon$ of the CDF of $\vartheta^{\ast}$ (i.e., its
probability density) at $\vartheta^{\ast}=0.5$. Note that larger $\epsilon$
values ultimately correspond to more $\vartheta^{\ast}$ values lying between 0
and 1 (indicating partial insult), and thus softer normal/insult boundaries
(Figure 11). Evaluating the (Gaussian) probability density at 0.5,
substituting Equation 18 for $\mu$, and solving for $\varsigma$ yields

$\displaystyle\epsilon$
$\displaystyle=\frac{1}{\varsigma\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{0.5-\mu}{\varsigma}\right)^{2}\right)$
(19) $\displaystyle\implies\epsilon$
$\displaystyle=\frac{1}{\varsigma\sqrt{2\pi}}\exp\left(-\left[\mathrm{erf}^{-1}\left(1-2\varphi\right)\right]^{2}\right)$
$\displaystyle\implies\varsigma$
$\displaystyle=\frac{1}{\epsilon\sqrt{2\pi}}\exp\left(-\left[\mathrm{erf}^{-1}\left(1-2\varphi\right)\right]^{2}\right).$

Substituting Equation 19 into Equation 18 and squaring Equation 19, we obtain
final expressions for the mean and variance,

$\displaystyle\mu$
$\displaystyle=\frac{1}{2}-\frac{1}{\epsilon\sqrt{\pi}}\,\mathrm{erf}^{-1}\left(1-2\varphi\right)\exp\left(-\left[\mathrm{erf}^{-1}\left(1-2\varphi\right)\right]^{2}\right)$
(20) $\displaystyle\textrm{and}\quad\varsigma^{2}$
$\displaystyle=\frac{1}{2\pi\epsilon^{2}}\exp\left(-2\left[\mathrm{erf}^{-1}\left(1-2\varphi\right)\right]^{2}\right).$

Figure 11: CDFs of insult field values before ($\vartheta^{\ast}$) and after
($\vartheta$) censoring, with insets illustrating the effect of $\epsilon$ on
the softness of the boundaries between normal and insult regions
($\varphi=50$%, $L_{\theta}=1.5$ mm, and $L_{z}=2$ mm herein). Smaller values
of $\epsilon$ correspond to sharper boundaries, while larger values correspond
to softer boundaries.

###  B.2 Sensitivity analysis

The pipeline detailed in subsubsection 2.2.2 and subsection B.1 is capable of
generating random insult profiles that are greatly variable in appearance
(Figure 12), while precisely controlling desired physically meaningful
metrics, including the overall fraction of the vessel experiencing insult
($\varphi$), the softness of the boundaries between normal and insult regions
($\epsilon$), and the size of the insult region(s) in both the circumferential
and axial directions ($L_{\theta}$ and $L_{z}$, respectively). To demonstrate
and quantify how $\varphi$, $\epsilon$, $L_{\theta}$, and $L_{z}$ jointly
control the size, shape, and appearance of randomly generated insults, we
present here a sensitivity analysis illustrating how each of these affect the
resulting insult profile. Note that while the insult propensity (Figure 13)
and boundary softness (Figure 14) parameters alter the insult profile globally
while mostly preserving the relative location(s) and shape(s) of the insults,
the length scale parameters (Figures 15 and 16) affect the profile more
locally. When varied in tandem (Figure 15), $L_{\theta}$ and $L_{z}$ tend to
uniformly scale the size of insult regions, conjoining them as needed to
preserve the prescribed insult propensity. In contrast, when the length scales
are varied individually (Figure 16), the insult region size is only altered in
the direction of interest, thus tending to change the aspect ratio of the
insult regions that are generated. In the limit, as $L_{\theta}\to\infty$, the
insult profile is constrained to be perfectly symmetric about the axial
direction, thus guaranteeing that only fusiform aneurysms are produced. The
distinction between varying the size and aspect ratio of the insult(s) is
highlighted by examining the corresponding GRF correlation function

$\varrho(z_{o},\theta_{o},z_{o}^{\prime},\theta_{o}^{\prime})=\frac{\kappa(z_{o},\theta_{o},z_{o}^{\prime},\theta_{o}^{\prime})}{\varsigma^{2}}=\exp\left(-\frac{1}{2}\left[\left(\frac{D_{\theta}\left(\theta_{o},\theta_{o}^{\prime}\right)}{L_{\theta}}\right)^{2}+\left(\frac{D_{z}\left(z_{o},z_{o}^{\prime}\right)}{L_{z}}\right)^{2}\right]\right),$
(21)

which governs the pairwise correlation in $\vartheta^{\ast}$ values between
different points on the vessel wall (Figures 15 and 16, bottom panels).

Figure 12: Representative random samples generated using the pipeline
described in subsubsection 2.2.2. Herein, the following were held constant:
$\varphi=35$%, $\epsilon=0.2$, $L_{\theta}=L_{z}=2$ mm. Figure 13: Sensitivity
of representative random insult profiles with respect to the overall
propensity of insult $\varphi$, holding the following constant:
$\epsilon=0.2$, $L_{\theta}=L_{z}=2$ mm. Figure 14: Sensitivity of
representative random insult profiles with respect to the boundary softness
parameter $\epsilon$, holding the following constant: $\varphi=35$%,
$L_{\theta}=L_{z}=2$ mm. Figure 15: (top) Sensitivity of a representative
random insult profile with respect to the insult size $L_{\theta}=L_{z}$,
holding the following constant: $\varphi=35$%, $\epsilon=0.2$. (bottom)
Correlation function $\varrho$ evaluated at a central point, showing how the
bandwidth of spatial correlation between points increases with increasing
insult size. Note that as $L_{\theta}$ and $L_{z}$ (and thus also spatial
correlation) increase, the insult profile transitions from producing several
small insult regions to only a few larger insult regions, while still matching
the prescribed overall perturbed area fraction $\varphi$. Figure 16: (top)
Sensitivity of a representative random insult profile with respect to the
circumferential insult length scale $L_{\theta}$, holding the following
constant: $\varphi=35$%, $\epsilon=0.2$, $L_{z}=2$ mm. (bottom) Correlation
function $\varrho$ evaluated at a central point, showing how the bandwidth of
spatial correlation between points increases in the circumferential direction
with increasing $L_{\theta}$ while remaining unchanged in the axial direction
(since $L_{z}$ is held constant). In the limit, the generated insult profile
is constrained to be perfectly axisymmetric, varying only in the axial
direction.

##  Appendix C DeepONet model

To understand the underlying mathematical foundation of a DNN, we consider a
network with $L$ hidden layers, where the $0$-th layer denotes the input layer
and the $(L+1)$-th layer is the output layer; the weighted input
$\bm{z}^{l}_{i}$ into a $i$th neuron on layer $l$ is a function of weight
$\bm{W}^{l}_{ij}$ and bias $\bm{b}^{l-1}_{j}$ and is represented as

$\bm{z}^{l}_{i}=\mathcal{R}_{l-1}\left(\sum_{j=1}^{m_{l-1}}\left(\bm{W}^{l}_{ij}(\bm{z}^{l-1}_{j})+\bm{b}^{l}_{j}\right)\right),$
(22)

where $\mathcal{R}_{l-1}\left(\cdot\right)$ denotes the activation function of
layer $l$, and $m_{l-1}$ is the number of neurons in layer $l-1$. Based on the
foregoing concepts, the feed-forward algorithm for computing the output
$\bm{Y}^{L}$ is expressed as follows:

$\begin{split}\bm{Y}^{L}&=\mathcal{R}_{L}(\bm{W}^{L+1}\bm{z}^{L}+\bm{b}^{L})\\\
\bm{z}^{L}&=\mathcal{R}_{L-1}\left(\bm{W}^{L}\bm{z}^{L-1}+\bm{b}^{L}\right)\\\
\bm{z}^{L-1}&=\mathcal{R}_{L-2}\left(\bm{W}^{L-1}\bm{z}^{L-2}+\bm{b}^{L-1}\right)\\\
&\;\;\;\;\;\;\;\vdots\\\
\bm{z}^{1}&=\mathcal{R}_{0}\left(\bm{W}^{1}\bm{x}+\bm{b}^{1}\right),\\\
\end{split}$ (23)

where $\bm{x}$ is the input of the neural network. Equation 23 can be encoded
in compressed form as $\bm{Y}=\mathbb{N}(x;\bm{\theta})$, where
$\bm{\theta}=\left(\bm{W},\bm{b}\right)$ includes both the weights and biases
of the neural network $\mathbb{N}$. Taking into account a DeepONet, the branch
network takes as input the function to denote the input realizations
$\bm{U}=\\{\bm{u}_{1},\bm{u}_{2},\ldots,\bm{u}_{N}\\}$ for $N$ samples,
discretized at $n_{sen}$ sensor locations such that
$\bm{u}_{i}=\\{u_{i}(\bm{x}_{1}),u_{i}(\bm{x}_{2}),\ldots,u_{i}(\bm{x}_{n_{sen}})\\}$
and $i\in[1,N]$. The trunk net inputs the location
$\bm{y}=\\{\bm{y}_{1},\bm{y}_{2},\cdots,\bm{y}_{p}\\}=\\{(\hat{x}_{1},\hat{y}_{1}),(\hat{x}_{2},\hat{y}_{2}),\ldots,(\hat{x}_{p},\hat{y}_{p})\\}$
to evaluate the solution operator, where $\hat{x}_{i}$ and $\hat{y}_{i}$
denote the coordinates $x$ and $y$ of the point $\bm{y}_{i}$, respectively.
Let us consider that the branch neural network consists of $l_{br}$ hidden
layers, where the $(l_{br}+1)$th layer is the output layer consisting of $q$
neurons. Considering an input function $\bm{u}_{i}$ in the branch network, the
network returns a feature embedded in
$[b_{1},b_{2},\ldots,b_{q}]^{\mathrm{T}}$ as output. The output
$\bm{z}_{br}^{l_{br}+1}$ of the feed-forward branch neural network is
expressed as

$\begin{split}\bm{z}_{br}^{l_{br}+1}&=\left[b_{1},b_{2},\ldots,b_{q}\right]^{\mathrm{T}}\\\
&=\mathcal{R}_{br}\left(\bm{W}^{l_{br}}\bm{z}^{l_{br}}+\bm{b}^{l_{br}+1}\right),\end{split}$
(24)

where $\mathcal{R}_{br}\left(\cdot\right)$ denotes the nonlinear activation
function for the branch net and
$\bm{z}^{l_{br}}=f_{br}(u_{i}(\bm{x}_{1}),u_{i}(\bm{x}_{2}),\ldots,u_{i}(\bm{x}_{n_{sen}}))$,
where $f_{br}\left(\cdot\right)$ denotes a branch net function. Similarly,
consider a trunk network with $l_{tr}$ hidden layers, where the
$(l_{tr}+1)$-th layer is the output layer consisting of $q$ neurons. The trunk
net outputs a feature embedding $[t_{1},t_{2},\ldots,t_{q}]^{\mathrm{T}}$. The
output of the trunk network can be represented as

$\begin{split}\bm{z}_{tr}^{l_{tr}+1}&=\left[t_{1},t_{2},\ldots,t_{q}\right]^{\mathrm{T}}\\\
&=\mathcal{R}_{tr}\left(\bm{W}^{l_{tr}}\bm{z}^{l_{tr}}+\bm{b}^{l_{tr}+1}\right),\end{split}$
(25)

where $\mathcal{R}_{tr}\left(\cdot\right)$ denotes the non-linear activation
function for the trunk net and
$\bm{z}^{l_{tr}-1}=f_{tr}(\bm{y}_{1},\bm{y}_{2},\ldots,\bm{y}_{p})$. The key
point is that we uncover a new operator $\mathbb{G}_{\bm{\theta}}$ as a neural
network that can infer quantities of interest from unseen and noisy inputs.
The two networks are trained to learn the solution operator such that

$\mathbb{G}_{\bm{\theta}}:\bm{u}_{i}\rightarrow\mathbb{G}_{\bm{\theta}}(\bm{u}_{i}),\;\;\forall\;\;i=\\{1,2,3,\ldots,N\\}.$
(26)

For a single input function $\bm{u}_{i}$, the DeepONet prediction
$\mathbb{G}_{\bm{\theta}}(\bm{u})$ evaluated at any coordinate $\bm{y}$ can be
expressed as

$\begin{split}\mathbb{G}_{\bm{\theta}}(\bm{u}_{i})(\bm{y})&=\sum_{k=1}^{q}\left(\mathcal{R}_{br}(\bm{W}^{l_{br}}_{k}\bm{z}^{l_{br}-1}_{k}+\bm{b}^{l_{br}}_{k})\cdot\mathcal{R}_{tr}(\bm{W}^{l_{tr}}_{k}\bm{z}^{l_{tr}-1}_{k}+\bm{b}^{l_{tr}}_{k})\right)\\\
&=\sum_{k=1}^{q}b_{k}(u_{i}(\bm{x}_{1}),u_{i}(\bm{x}_{2}),\ldots,u_{i}(\bm{x}_{m}))\cdot
t_{k}(\bm{y}).\end{split}$ (27)

DeepONet requires large annotated datasets of paired input-output
observations, but it provides a simple and intuitive model architecture that
is fast to train, allowing a continuous representation of the target output
functions that is resolution-independent. Conventionally, the trainable
parameters of the DeepONet represented by $\bm{\theta}$ in Equation 27 are
obtained by minimizing a loss function. Common loss functions used in the
literature include the $L_{1}$\- and $L_{2}$-loss functions, defined as

$\begin{split}\mathcal{L}_{1}&=\sum_{i=1}^{n}\sum_{j=1}^{p}\big{|}\mathbb{G}(\bm{u}_{i})(\bm{y}_{j})-\mathbb{G}_{\bm{\theta}}(\bm{u}_{i})(\bm{y}_{j})\big{|}\\\
\mathcal{L}_{2}&=\sum_{i=1}^{n}\sum_{j=1}^{p}\big{(}\mathbb{G}(\bm{u}_{i}(\bm{y}_{j})-\mathbb{G}_{\bm{\theta}}(\bm{u}_{i})(\bm{y}_{j})\big{)}^{2},\\\
\end{split}$ (28)

where $\mathbb{G}_{\bm{\theta}}(\bm{u}_{i})(\bm{y}_{j})$ is the predicted
value obtained from the DeepONet, and $\mathbb{G}(\bm{u}_{i})(\bm{y}_{j})$ is
the target value.

## References

  * [1] M. E. Lindsay, H. C. Dietz, The Genetic Basis of Aortic Aneurysm, Cold Spring Harbor Perspectives in Medicine 4 (9) (2014) a015909.
  * [2] F. Ramirez, C. Caescu, E. Wondimu, J. Galatioto, Marfan syndrome; A connective tissue disease at the crossroads of mechanotransduction, TGF signaling and cell stemness, Matrix Biology 71-72 (2018) 82–89.
  * [3] T. F. Vinholo, A. J. Brownstein, B. A. Ziganshin, M. A. Zafar, H. Kuivaniemi, S. C. Body, A. E. Bale, J. A. Elefteriades, Genes Associated with Thoracic Aortic Aneurysm and Dissection: 2019 Update and Clinical Implications, Aorta (Stamford, Conn.) 7 (4) (2019) 99–107.
  * [4] D. M. Milewicz, F. Ramirez, Therapies for Thoracic Aortic Aneurysms and Acute Aortic Dissections, Arteriosclerosis, Thrombosis, and Vascular Biology 39 (2) (2019) 126–136.
  * [5] A. Pinard, G. T. Jones, D. M. Milewicz, Genetics of Thoracic and Abdominal Aortic Diseases, Circulation Research 124 (4) (2019) 588–606.
  * [6] M. Lindquist Liljeqvist, M. Bogdanovic, A. Siika, T. C. Gasser, R. Hultgren, J. Roy, Geometric and biomechanical modeling aided by machine learning improves the prediction of growth and rupture of small abdominal aortic aneurysms, Scientific Reports 11 (1) (2021) 18040.
  * [7] S. Zoli, C. D. Etz, F. Roder, C. S. Mueller, R. M. Brenner, C. A. Bodian, G. Di Luozzo, R. B. Griepp, Long-term survival after open repair of chronic distal aortic dissection, The Annals of Thoracic Surgery 89 (5) (2010) 1458–66.
  * [8] V. Ziza, L. Canaud, N. Molinari, P. Branchereau, C. Marty-Ané, P. Alric, Thoracic endovascular aortic repair: A single center’s 15-year experience, The Journal of Thoracic and Cardiovascular Surgery 151 (6) (2016) 1595–1603.
  * [9] J. P. Pirruccello, M. D. Chaffin, E. L. Chou, S. J. Fleming, H. Lin, M. Nekoui, S. Khurshid, S. F. Friedman, A. G. Bick, A. Arduini, et al., Deep learning enables genetic analysis of the human thoracic aorta, Nature Genetics 54 (1) (2022) 40–51.
  * [10] D. M. Milewicz, S. K. Prakash, F. Ramirez, Therapeutics Targeting Drivers of Thoracic Aortic Aneurysms and Acute Aortic Dissections: Insights from Predisposing Genes and Mouse Models, Annual Review of Medicine 68 (1) (2017) 51–67.
  * [11] C. Bellini, A. Korneva, L. Zilberberg, F. Ramirez, D. B. Rifkin, J. D. Humphrey, Differential Ascending and Descending Aortic Mechanics Parallel Aneurysmal Propensity in a Mouse Model of Marfan Syndrome, Journal of Biomechanics 49 (12) (2016) 2383–2389.
  * [12] C. Cavinato, M. Chen, D. Weiss, M. J. Ruiz-Rodríguez, M. A. Schwartz, J. D. Humphrey, Progressive microstructural deterioration dictates evolving biomechanical dysfunction in the marfan aorta, Frontiers in Cardiovascular Medicine 8 (2021).
  * [13] B. L. Lima, E. J. Santos, G. R. Fernandes, C. Merkel, M. R. Mello, J. P. Gomes, M. Soukoyan, A. Kerkis, S. M. Massironi, J. A. Visintin, et al., A New Mouse Model for Marfan Syndrome Presents Phenotypic Variability Associated with the Genetic Background and Overall Levels of Fbn1 Expression, PloS one 5 (11) (2010) 1–9.
  * [14] J. D. Humphrey, M. A. Schwartz, G. Tellides, D. M. Milewicz, Role of Mechanotransduction in Vascular Biology, Circulation Research 116 (8) (2015) 1448–61.
  * [15] A. C. Estrada, L. Irons, B. V. Rego, G. Li, G. Tellides, J. D. Humphrey, Roles of mTOR in thoracic aortopathy understood by complex intracellular signaling interactions, PloS Computational Biology 17 (12) (2021) e1009683.
  * [16] A. Valentin, J. D. Humphrey, Evaluation of fundamental hypotheses underlying constrained mixture models of arterial growth and remodelling, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367 (1902) (2009) 3585–3606.
  * [17] J. S. Wilson, S. Baek, J. D. Humphrey, Parametric study of effects of collagen turnover on the natural history of abdominal aortic aneurysms, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 469 (2150) (2013) 20120556.
  * [18] L. Liang, M. Liu, C. Martin, J. A. Elefteriades, W. Sun, A machine learning approach to investigate the relationship between shape features and numerically predicted risk of ascending aortic aneurysm, Biomechanics and Modeling in Mechanobiology 16 (5) (2017) 1519–1533.
  * [19] D. Weiss, M. Latorre, B. V. Rego, C. Cavinato, B. J. Tanski, A. G. Berman, C. J. Goergen, J. D. Humphrey, Biomechanical consequences of compromised elastic fiber integrity and matrix cross-linking on abdominal aortic aneurysmal enlargement, Acta Biomaterialia 134 (2021) 422–434.
  * [20] M. Pradella, T. Weikert, J. I. Sperl, R. Kärgel, J. Cyriac, R. Achermann, A. W. Sauter, J. Bremerich, B. Stieltjes, P. Brantner, G. Sommer, Fully automated guideline-compliant diameter measurements of the thoracic aorta on ECG-gated CT angiography using deep learning, Quantitative Imaging in Medicine and Surgery 11 (10) (2021) 4245–4257.
  * [21] N. P. Ostberg, M. A. Zafar, S. K. Mukherjee, B. A. Ziganshin, J. A. Elefteriades, A machine learning approach for predicting complications in descending and thoracoabdominal aortic aneurysms, The Journal of Thoracic and Cardiovascular Surgery (2022).
  * [22] X. He, S. Avril, J. Lu, Estimating aortic thoracic aneurysm rupture risk using tension–strain data in physiological pressure range: an in vitro study, Biomechanics and Modeling in Mechanobiology 20 (2) (2021) 683–699.
  * [23] M. Liu, L. Liang, Y. Ismail, H. Dong, X. Lou, G. Iannucci, E. P. Chen, B. G. Leshnower, J. A. Elefteriades, W. Sun, Computation of a probabilistic and anisotropic failure metric on the aortic wall using a machine learning-based surrogate model, Computers in Biology and Medicine 137 (2021) 104794.
  * [24] M. Zhou, Z. Shi, X. Li, L. Cai, Y. Ding, Y. Si, H. Deng, W. Fu, Prediction of Distal Aortic Enlargement after Proximal Repair of Aortic Dissection Using Machine Learning, Annals of Vascular Surgery 75 (2021) 332–340.
  * [25] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707.
  * [26] S. Goswami, C. Anitescu, T. Rabczuk, Adaptive fourth-order phase field analysis using deep energy minimization, Theoretical and Applied Fracture Mechanics 107 (2020) 102527.
  * [27] S. Goswami, C. Anitescu, S. Chakraborty, T. Rabczuk, Transfer learning enhanced physics informed neural network for phase-field modeling of fracture, Theoretical and Applied Fracture Mechanics 106 (2020) 102447.
  * [28] A. D. Jagtap, E. Kharazmi, G. E. Karniadakis, Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems, Computer Methods in Applied Mechanics and Engineering 365 (2020) 113028.
  * [29] S. Goswami, M. Yin, Y. Yu, G. E. Karniadakis, A physics-informed variational DeepONet for predicting crack path in quasi-brittle materials, Computer Methods in Applied Mechanics and Engineering 391 (2022) 114587.
  * [30] M. Yin, E. Ban, B. V. Rego, E. Zhang, C. Cavinato, J. D. Humphrey, G. Em Karniadakis, Simulating progressive intramural damage leading to aortic dissection using DeepONet: an operator–regression neural network, Journal of the Royal Society Interface 19 (187) (2022) 20210670.
  * [31] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nature Machine Intelligence 3 (3) (2021) 218–229.
  * [32] C. M. Wojnarski, E. E. Roselli, J. J. Idrees, Y. Zhu, T. A. Carnes, A. M. Lowry, P. H. Collier, B. Griffin, J. Ehrlinger, E. H. Blackstone, et al., Machine-learning phenotypic classification of bicuspid aortopathy, The Journal of Thoracic and Cardiovascular Surgery 155 (2) (2018) 461–469.
  * [33] M. Latorre, J. D. Humphrey, Numerical knockouts–In silico assessment of factors predisposing to thoracic aortic aneurysms, PloS Computational Biology 16 (10) (2020) 1–25.
  * [34] 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.
  * [35] S. Baek, K. R. Rajagopal, J. D. Humphrey, A Theoretical Model of Enlarging Intracranial Fusiform Aneurysms, Journal of Biomechanical Engineering 128 (1) (2006) 142–149.
  * [36] A. Valentin, 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) (2009) 293–306.
  * [37] M. Latorre, J. D. Humphrey, A Mechanobiologically Equilibrated Constrained Mixture Model for Growth and Remodeling of Soft Tissues, Zeitschrift für angewandte Mathematik und Mechanik 98 (2018) 2048–2071.
  * [38] 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\.
  * [39] M. R. Bersi, C. Bellini, J. Wu, K. R. Montaniel, D. G. Harrison, J. D. Humphrey, Excessive Adventitial Remodeling Leads to Early Aortic Maladaptation in Angiotensin-Induced Hypertension, Hypertension 67 (5) (2016) 890–896.
  * [40] D. P. Kingma, J. Ba, Adam: A Method for Stochastic Optimization, in: In Proceedings of the 3rd International Conference on Learning Representations, 2015\.
  * [41] D. C. Liu, J. Nocedal, On the limited memory BFGS method for large scale optimization, Mathematical Programming 45 (1) (1989) 503–528.
  * [42] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al., TensorFlow: Large-scale machine learning on heterogeneous systems, Software available from tensorflow. org 1 (2) (2015).
  * [43] B. V. Rego, D. Weiss, M. R. Bersi, J. D. Humphrey, Uncertainty quantification in subject-specific estimation of local vessel mechanical properties, International Journal for Numerical Methods in Biomedical Engineering 37 (12) (2021) e3535.