Modeling intracranial aneurysm stability and growth: an integrative mechanobiological framework for clinical cases

We present a novel patient-specific fluid-solid-growth framework to model the mechanobiological state of clinically detected intracranial aneurysms (IAs) and their evolution. The artery and IA sac are modeled as thick-walled, non-linear elastic fiber-reinforced composites. We represent the undulation distribution of collagen fibers: the adventitia of the healthy artery is modeled as a protective sheath whereas the aneurysm sac is modeled to bear load within physiological range of pressures. Initially, we assume the detected IA is stable and then consider two flow-related mechanisms to drive enlargement: (1) low wall shear stress; (2) dysfunctional endothelium which is associated with regions of high oscillatory flow. Localized collagen degradation and remodelling gives rise to formation of secondary blebs on the aneurysm dome. Restabilization of blebs is achieved by remodelling of the homeostatic collagen fiber stretch distribution. This integrative mechanobiological modelling workflow provides a step towards a personalized risk-assessment and treatment of clinically detected IAs.


Introduction
Intracranial aneurysm (IA), a localized focal out-pouching of the cerebral vasculature, occurs in 3-5% of the adult population. The risk of rupture is very low; however, if rupture does occur, mortality rates are very high. Physicians select IAs for intervention by weighing up rupture likelihood with risks associated with treatments. Unfortunately, personalized rupture risk assessment remains elusive and, consequently, intervention may be recommended for IAs that are stable.
Despite continuous efforts on identifying the etiology of the disease, the exact causes of aneurysm formation, growth, stabilization and/or rupture are still not fully understood. Correlations with genetic factors [such as tissue disorders, polycistic kidney disease, first-degree family history etc. (Samuel and Radovanovic 2019)], gender (Desai et al. 2019), ethnicity  and external factors (such as high blood pressure, smoking, drugs abuse, head trauma or brain tumor, infection of the arterial wall) have been reported (Toth and Cerejo 2018). Here, we focus on modelling the flow-mediated mechanisms (Cebral et al. 2016) that influence IA pathobiology.
Recent developments in imaging technology have led to asymptomatic IAs being detected more frequently ). Hence there is an increasing need to identify IAs that are stable and not in need of intervention. This will reduce the (increasing) financial burden of the disease and will help to reduce unnecessary intervention. Computational models of IAs that quantify their mechanobiological state (and it's potential evolution) may provide insight into IA stability and thus assist personalized diagnosis.

3
Computational mechano-biological models of IAs need to account for: (1) the biomechanics of the arterial/aneurysm wall; (2) the cellular biology of the tissue and (3) the dynamic interplay between (1) and (2), i.e. the mechanobiology. Over a decade ago, (Humphrey and Taylor 2008) emphasized the need for a new class of Fluid-Solid-Growth models to study aneurysm evolution and proposed the terminology FSG models. These combine fluid and solid mechanics to quantify the flow and structural mechanical stimuli that drive mechanobiological processes. FSG models have been applied by several research groups to model abdominal aortic aneurysm (Zeinali-Davarani and Baek 2012;Watton et al. 2012;Aparício et al. 2014;Wu and Shadden 2015;Grytsan et al. 2015) and IA evolution (Watton et al. 2008(Watton et al. , 2009a(Watton et al. , 2011aSelimovic et al. 2014).
To date, FSG models have explored phenomenological relationships between mechanical stimuli and tissue remodelling that give rise aneurysm evolution. Nevertheless, they provide a framework for integration of cell models and representation of biochemical pathway models of the arterial wall (Harvey et al. 2018;Aparício 2016) which will provide further insight into disease pathophysiology.
Previous FSG models have IA evolution from an initial healthy artery. However, from a clinical perspective characterizing the mechanobiolgical stability (and potential evolution) from the clinical detected geometry is what is important. In this work, we present an integrative FSG modelling framework for the personalized quantification of the mechanobiological state of an IA and it's potential evolution. This has been realized within Sim4Life (Neufeld et al. 2013), a Computational Life Sciences (CLS) simulation platform co-developed by the IT'IS Foundation and ZMT Zurich MedTech AG.
We build on previous studies (e.g. Watton et al. 2011b;Selimovic et al. 2014) to incorporate several new features: integrative pipeline from imaging to simulation; application to patient-specific IA geometries; integration of a novel collagen fiber distribution model into FE code; linking growth and remodelling to pulsatile flow metrics; simulation of secondary blebs; adaption of collagen homeostasis to stabilize IA.
This paper is organized as follows: Section 2 reviews the thick-walled structural model, the computational fluid dynamics solver and the growth and remodeling approach. Simulation results are presented in Sect. 3: first for obtaining the initial homeostatic state in patient-specific geometries; secondly for two different flow-driven remodelling scenarios. Finally, Sect. 5 discusses conclusions and limitations of the model.

Computational framework and mathematical model for the aneurysm evolution
The computational frawework presented is implemented into Sim4Life (Neufeld et al. 2013), a state-of-the-art simulation platform for computational life sciences centered around high-resolution, detailed anatomical phantoms functionalized with dynamic physiology and tissue models; image-based modelling is strongly supported. The platform features a modular graphical user interface (GUI), VTK-and OpenGL-based visualization, a Python scripting interface, and a pipeline framework for advanced analysis. It incorporates: • A robust vessel segmentation software, capable of dealing with diseased vessel geometries and a multitude of imaging modalities, and supporting efficient interactive corrections where automatic segmentation fails. • A comprehensive pipeline supporting all steps from image-based model creation to results analysis, including: tissue property assignment; shape characterization; fluid dynamics modeling (with physiologically realistic boundary conditions); solid mechanics; growth and remodeling simulation; analysis, post-processing and visualisation.
It provides the first complete simulation pipeline for modeling growth of clinically detected IA, see Fig. 1.

Kinematics
The vessel wall experiences large displacement and is made of a nearly incompressible material (Holzapfel 2000). From the structural point of view, it can be decomposed into three layers: intima, media, and adventitia. For a comprehensive analysis of the wall architecture, see Robertson and Watton (2013). This study models the wall with a the three-dimensional, two-layered thick-wall structure (Teixeira 2019). The endothelial layer is assumed to be mechanically irrelevant, while the medial layer is composed of isotropic materials (elastin, passive smooth muscle ground matrix etc.) and two families of anisotropic collagen fibers; the outer adventitial layer is composed of collagen fibers embedded in an isotropic matrix.
The Lagrangian formulation used in this work denotes the deformation gradient as , the Green-Lagrange strain tensor as , and the right Cauchy-Green tensor as = T . The near-incompressible effects are accounted with a multiplicative split of the deformation gradient = J 1∕3̄ , where J = det( ) , see Holzapfel et al. (2000). The modified versions of the right Cauchy-Green tensor and Green-Lagrange (GL) strain tensor are defined by ̄ =̄ T̄ , and ̄ = 1 2 ̄ − , respectively. The modified tensor invariants are simply Ī 1 = tr(̄ ) , The direction of the collagen fibers is denoted by the unit vector m i , i = 1, 2 . The total stretch i 4 expresses the total stretch in the direction m i . The modified version ̄i 4 = J −1∕3 i 4 can be evaluated by i.e., associated with Ī i 4 , a pseudo-invariant of ̄ and m i . We assume the hyperelasticity hypothesis, which implies the existence of an energy function Ψ = Ψ( ) . The strain energy function Ψ (SEF) is split into an additive sum of isochoric (volume-preserving) and dilatational (volumechanging) parts (Simo et al. 1985) In the absence of external body forces, the displacement field must satisfy ∇ ⋅ ( ) = 0 where the Second Piola-Kirchhoff stress tensor is

Strain-energy functions
The isochoric part of each layer L (where L = M denotes medial layer and L = A denotes the adventitial layer) receives contributions from elastin, passive smooth muscle cells and ground matrix (its isotropic components) and collagen fibers (anisotropic components). The volumetric part receives contributions from the dilatational elastic response of the tissue. Hence

Elastinous constitutents
The isotropic components are modeled as a neo-Hookean material (Watton et al. 2009b): with K L,e and K L,sm being stiffness-like material constants, and m L,e and m L,sm the (dimensionless) normalized mass density of elastin and passive smooth muscle cells, respectively.

Collagenous constituents
The constitutive model for the collagen accounts for the experimental observation that collagen fibers have a a distribution of waviness in the unloaded configuration and thus have a distribution of recruitment (Hill et al. 2012). The strain energy function involves integrating the strain energy for a fiber ( Ψ i L,c ) over the distribution of recruitment stretches ( i rec ), (2) (3) Ψ L,e = m L,e ⋅ K L,e + m L,sm ⋅ K L,sm ⋅ Ī 1 − 3 , where m L,c is the (dimensionless) normalized mass density of collagen fibers that can adapt to simulate growth/atrophy (Watton et al. 2004). We follow Chen (2014) and model each individual fiber to have a linear relationship between it's 1st Piola-Kirchoff stress and the collagen fiber stretch (̄i 4c ) where K L,c is a stiffness-like material constant, 4r (Aparício 2016); ̄i ,q 4r relates to the the minimum ( q = min ), modal ( q = mode ) and maximum ( q = max ) recruitment stretches of the distribution (see Fig. 7). More specifically: where Insertion of Eqs. 5 and 6 into 4 and integration yields analytic expressions for the strain energy from which Ψ I i 4 can be derived. More specifically: , see appendix of Aparício et al. (2016) for details. Figure 2 illustrates the 1 st Piola-Kirchhoff stress versus stretch for the fiber-distribution for illustrative recruitment stretch distributions.

Volumetric function
The volumetric/dilatational elastic response adopted here stems from Simo and Armero (1992) and is defined as where the parameter is the bulk modulus.

Structural solvers
A PETSc-based (Balay et al. 2019a(Balay et al. , b, 1997 HPC-enabled FEM framework has been created, that flexibly supports a range of element types (tetrahedral, pyramidal, prismatic, etc.), shape-functions, preconditioning schemes, and multi-physics coupling. Meshes are partitioned with Parmetis (Karypis  A specialized structural solver has been developed based on the FEM solver framework. It supports compressible/incompressible, homogeneous/heterogeneous, linear/non-linear, and isotropic/anisotropic material models, multiple populations of collagen fibers with varying orientations and recruitment stretch distributions.

Computational fluid dynamics
A computational fluid dynamics (CFD) solver was also developed on the basis of the FEM solver framework. It supports P 1 -P 1 elements through Variational Multiscale Stabilization (Bazilevs et al. 2007;Forti and Dedè 2015;Liu and Marsden 2018). The algebraic system is solved with FGMRES (Saad 1993), preconditioned with one iteration of SIMPLE (Patankar and Spalding 1983;Pernice and Tocci 2001). Internally, the velocity and Schur complement blocks are preconditioned with fixed iterations of HYPRE (Falgout and Yang 2002).
The incompressible Navier-Stokes equations are used to model flow. The physical parameters of the simulated fluid resemble those of blood: density f = 1066 kg/m 3 and dynamic viscosity f = 3.5 × 10 −3 Pa s. Three cardiac cycles of 0.8s each were simulated with time-step 0.004s. Figure 3 plots the inlet flow rate and outlet pressure boundary conditions (Villa-Uriol et al. 2011).

Flow metrics
Given the fluid pressure p and velocity , let denote the Newtonian fluid stress tensor, and be an outward unit normal vector at the lumen surface. The wall shear stress vector is given by and represents the shearing component of the traction vector .
The complex geometry of arteries coupled with pulsatile flow gives rise to regions of the vasculature that are subject to oscillatory flow (magnitude and direction). This can be quantified by a novel flow metric, the wall shear stress aspect ratio (WSSAR) which naturally emerges from quantifying the shear stress rosettes, which capture the multidirectionality of flow (Krishna et al. 2020). The WSSAR measures the bidirectionality of the WSS over a given interval [0, T] (Krishna et al. 2020). It is calculated as follows: given a point x on the wall, let = (x, t) denote the WSS at the point x for t ∈ [0, T] , and 1 = 1 (x) and 2 = 2 (x) two orthogonal directions which The directions 1 and 2 are named principal directions at point x. In addition, let be auxiliary values, which by definition, satisfy 1,max ≥ 1,min , 2,max ≥ 2,min , and max ≥ min . Then the WSSAR, denoted AR , is defined as and satisfies 0 ≤ WSSAR ≤ 1 , and

Growth and remodelling of the constituents
Clinical studies indicate that the medial layer degrades during aneurysm development, the intima layer disappears, the internal elastic lamina is lost, and apoptosis of vascular smooth muscle cells is observed (Tong et al. 2014). At the same time, the adventitia layer thickens. For this model, the starting point is an already developed IA. Initial mass-densities of elastin, collagen and smooth muscle cells are prescribed to represent atrophy of the medial layer. Subsequent degradation is driven by linking the mass-densities of constituents to perturations of flow metrics from homeostatic values. Although this approach still does not represent the full complexity of the IA formation with respect to the hemodynamic environment, it is step forward with respect to its application to anatomical IA geometries (Teixeira 2019).

Mass degradation
Additional aneurysm enlargement is modeled as a localized degradation of the layer L through the decay of the elastin ( m L,e ), collagen ( m L,c ), and smooth muscle cells normalized masses ( m L,sm ) (Watton et al. 2011a). The normalized mass degradation is defined as where D max is the maximum rate of degradation, y = e for elastin and y = c for collagen in the layer L; F X represents the influence of hemodynamics on G&R by linking collagen and elastin degradation to WSS-related metrics (Watton et al. 2011b). The hemodynamic environment of an IA is intimately related to it's initiation, growth and rupture (Robertson and Watton 2012;Cebral et al. 2019). It is observed that endothelial cells may still be present in developed aneurysms (Froösen et al. 2004 ) and biochemical pathways associated with functional/dysfunctional endothelium are likely to play a role in enlargement and rupture. In the following examples, we consider simple phenomenological relationships between regions of destructive remodelling and perturbed flow; firstly low WSS and, secondly, high oscillatory flow.
Low WSS: We follow Watton et al. (2009a), Watton et al. (2011b), Watton et al. (2011a), Selimovic et al. (2014) and assume a simple phenomological relationship between the WSS and degenerative processes of an IA, i.e. without explicitly modelling biochemical pathways. More specifically, F WSS in Eq. 10, is taken to be quadratic: where degradation is maximal if || || ≤ L , and there is no degradation if C < || || . For illustration, we follow Watton et al. (2009a) and assume that if the WSS is below 1 Pa, it will drive degenerative processes. We set C = 1 Pa, L = 0.5 Pa and D max = 1.5 . Figure 4 plots the function F WSS with the above-mentioned values.

WSSAR:
We hypothesize that endothelial morphology is related to the oscillatory nature of the flow and suppose that there exists a switch point where endothelial cells have distinct phenotypes, i.e. spindle or irregular morphologies. The morphology of the cells may affect the permeability and functionality of the endothelial layer. Specifically, we assume that regions where endothelial cells have irregular morphologies will drive destructive remodelling processes.
The spatially-variable function F AR ( ) quantifies the WSSAR (defined as AR ), such that, F AR = 0 means no degradation and F AR = 1 yields maximum degradation. The existence of thresholds H,AR is hypothesized, such that the degradation is maximal if H,AR ≤ AR , and absent if AR < C,AR . In this case, F AR is given a simple quadratic form as For illustration, we suppose that a WSSAR magnitude of 0.7 acts as a switch point for endothelial cells (ECs) to be packed irregularly. For the simulations, the critical threshold to activate degradation is set to C,AR = 0.7 Pa, while H,AR = 0.8 Pa. D max is set to 1.5. Figure 5 plots the function F AR with the above-mentioned values.

Collagen remodelling
Collagen fibers are in a continuous state of degradation and deposition and are configured to the artery in the physiological configuration in a state of stretch (Humphrey 1999). A consequence of this is that in the unloaded configuration they may appear wavy and thus the tissue needs to be stretched for a wavy fiber to be recruited to load bearing. Watton et al. 2004 introduced the terminology: attachment stretch, ̄4 c,att , to denote the stretch a fiber is configured to the matrix in the loaded configuration; recruitment stretch, ̄4 r , to denote the amount the unloaded tissue must be stretched (in the fiber direction) for the fiber to straighten out and begin to bear load. The collagen fiber stretch is then related to the tissue stretch ̄i 4 by the simple relationship ̄i 4c =̄i 4 ∕̄i 4r -if the tissue is in homeostasis (at physiological stretch) ̄i 4c =̄4 c,att . In vitro studies show that there is a distribution of collagen fiber recruitment stretches (Hill et al. 2012). Equivalently, this implies that there exists a distribution of attachment stretches: fibers which are recruited first to load bearing will be configured with maximum attachment stretches and fibers which are last to be recruited to load bearing will have minimum attachment stretches. Chen (2014) generalized the remodelling approach of Watton et al. (2004) to account for remodelling of a distribution of fibers. The initial distribution of recruitment ( t = 0 ) is determined by assuming a prescribed range of fiber attachment stretches. Remodelling of the recruitment stretch distribution then acts to maintain the collagen fiber stretch distribution (in the loaded configuration) towards the attachment stretch distribution, simulating the natural consequence of fiber deposition/degradation. Mathematically, this is achieved by evolving the recruitment stretch distribution so that the collagen stretch distribution remodels towards the attachment stretch distribution ): where we assume collagen fibers remodel to achieve a maximum stretch during the cardiac cycle and we assume that this will occur either at systole or diastole (Watton et al. 2011a).
We define three thresholds, namely ̄i ,min 4c,att , ̄i ,mean 4c,att and ̄i ,max 4c,att . The fibers with lowest attachment stretch ̄i ,min 4c,att will be recruited last and therefore have the highest recruitment ̄i ,max 4r , while fibers attached with ̄i ,max 4c,att are recruited first and have the smallest recruitment stretch ̄i ,min 4r . This implies that the collagen fiber recruitment stretches form a distribution, with the thresholds defined as

Collagen growth
The collagen fiber growth is modeled by the evolution of the normalized mass m L,c . We assume that ): • The number of fibroblasts is proportional to the mass of collagen they are maintaining. • Fibroblasts evolve their stress-free configuration as the collagen fabric remodels and they configure to the matrix to achieve a homeostatic stretch ̄i 4f ,att . • Collagen mass growth/atrophy is driven by deviation of fibroblast stretch from the homeostatic levels.
The simplest equation to model this phenomenon is where i L > 0 is a growth rate parameter. We assume that the natural configuration of a fibroblast cell to be identical to the natural configuration of the collagen fibers (with maximum stretch) and so ̄i 4f =̄i 4c and the above equation can be written in terms of the collagen fiber stretch.

Collagen stabilization mechanism
The adventitia acts as a protective sheath against overdistension of an artery due to natural variations in systemic pressure (Schmid et al. 2013). The implication being that it's load-bearing role at physiological pressures may be negligible (age, species, location dependent). Conversely, the parent dome of an intracranial aneurysm can consist of a thin collagenous sac which bears all the pressure load at physiological pressures, i.e. it is not acting as a protective mechanical sheath. These observations imply that the collagen fiber stretch distribution for the sac will differ from that of the protective sheath. Moreover, the fact that IAs generally stabilize in size implies that the collagen fabric is in mechanobiological equilibrium and thus the adventitial collagen fabric homeostasis adapts-suggesting that fibers are configured to bear more load as the load-bearing role of the medial layer is lost.
We simulate this by remodelling the distribution of attachment stretches (Mandaltsi 2016). Here we do not model how this is biologically achieved-rather we simply prescribe its adaption to achieve a new homeostasis where q = max,mode,min and att is a remodelling rate parameter which controls how quickly the attachment stretch adapts. It is set to 0.5 for the analyses in the following sections.

Heterogeneous characterization of the arterial wall
The geometry (obtained from Prof. Robertson, University of Pittsburgh) was extracted from computed tomography (CT). The generation of surface and volumetric meshes were performed in Sim4Life. In order to reduce the computational burden of the structural analysis, the computation was restricted to the neighborhood of the dome region. The planar cuts in the parent artery were selected far enough from the dome to avoid extra stresses caused by the clamped boundaries. The domain of interest was further split into parent artery, neck and dome subregions in order to locally adapt the mechanical characterization of the wall. Figure 6 shows the whole geometry and depicts each subregion.

Physiological parametrization
Tables 1 and 2 specify the list of parameters for the mathematical model (see Sect. 2 and Teixeira (2019) for a detailed description). The values reflect the physiological function of each layer for the wall for distinct healthy (parent artery) and diseased regions (aneurysm dome). For example, in the parent artery, the adventitia is assumed to act as a protective sheath (Schmid et al. 2013), and therefore the collagen fibers in this layer are only recruited at supraphysiological pressure. On the other hand, the adventitia acts as the main load-bearing component in the aneurysm dome (due the degradation of the medial layer), and therefore the collagen fibers must be recruited within the physiological range of pressures so that they are load-bearing (Robertson and Watton 2013). The initial recruitment stretch values are determined so that collagen achieves a prescribed homeostatic fiber stretch distribution at t = 0 ; Fig. 7 illustrates attachment and recruitment stretch distributions using an illustrative stretch of 1.3. The neck region is considered to be a transition region, with parameters being linearly interpolated between the healthy artery and the aneurysm sac.

Collagen fiber orientation
We model the collagen fiber orientations on the aneurysm dome to coincide with the principal curvature directions (Ma et al. 2007), whilst those in the parent artery are modeled as splayed about the circumferential direction (Holzapfel et al. 2000). Figure 8 depicts the resulting fiber orientation control: in the dome region, the directions were based on the principal curvature; in the parent artery, the directions Fig. 7 Illustrative probability density functions for load-bearing (blue) and protective (green) adventitial collagen attachment stretches (left), and the corresponding recruitment stretches (right) assuming a circumferential stretch of 1.3  were modified (rotated) 30 • (in the medial layer) and 60 • (in the adventitial layer) with respect to the circumferential direction. To maintain continuity of the fiber orientations, we interpolate the fiber directions in the neck region.

Simulation setup
Most IAs stabilize and remain asymptomatic. Hence, when an IA is detected, it is more likely to be stable than unstable. The iterative method to define the initial homeostasis in personalized geometries was established in Teixeira (2019). As a next step, the degradation phase implements two degradation approaches to model collagen and elastin evolution: (1) linked to low steady-state WSS, and (2) linked to pulsatile flow metrics (see Sect. 2.4.1). Finally, the stabilization phase proposes a mechanism based on remodelling the attachment stretches so that the collagen fabric can bear more load in homeostasis without increases in mass. Secondary bleb formation is simulated for a period of 3 years after initial homeostasis: 2 years of degradation linked to flow metrics, and 1 year of stabilization (see Sect. 2.5); The applied time-step of 0.02 year requires a total of 150 steps. For each time-step of the disease evolution problem, steady-state structural analysis of the reduced domain (see Fig. 6) is performed for diastolic and systolic pressures. The disease evolution model at each point is driven by the maximum ̄i 4c , see Eq. (5), of the systolic and diastolic configurations. CFD analyses recompute the WSS field over the whole domain every 20 steps after the onset of degradation (at time-step 50). The steady-state CFD analysis uses a constant inlet flow rate of 2.54 × 10 −6 m 3 ∕s and constant outlet pressures of 9671.38 Pa and 9833.19 Pa, respectively, which represents the mean flow rate over the cardiac cycle. For the G&R linked to WSSAR, transient CFD analyses runs for 2.4s (three cardiac cycles) with time-step 0.002s (total of 600 time-steps, 200 per cardiac cycle) and uses the physiological inlet flow rate and outlet pressure waves depicted in Fig. 3.
Each full G&R simulation which is linked to steady-state flow took 7 hours on average to run using 10 MPI processes on a Windows desktop computer equipped with an Intel Core i9-7900X 3.31GHz processor and 64GB of RAM. Those linked to transient flow took 26 hours on average.

Results
In this section, first we overview the numerical approach to model a stabilized IA. From this homeostatic state, development of secondary blebs is subsequently driven by linking collagen degradation to deviation of flow metrics from homeostatic ranges. Lastly, remodelling of the distribution of collagen fiber attachment stretches stabilizes the blebs.

Initial homeostasis
Given an initial prescribed recruitment stretch distribution field and established targets for collagen fiber stretch distribution (i.e. the minimum, mode, maximum attachment stretches), an iterative process searches for a (spatially heterogeneous) recruitment distribution field which yields mechanobiological equilibrium for the loaded IA. For the investigated configuration, a simulation with 50 iterations was sufficient to achieve homeostasis. Figure 9 depicts the displacement evolution from the initial deformed geometry to the (initial) homeostatic configuration. Figure 10 depicts the minimum, mode and maximum values of the collagen fiber stretch distribution. Notice that in the homeostatic state, these values are spatially uniform and equal to prescribe homeostatic values (see Table 2) and the collagen fabric is in mechanobiological equilibrium.
The evolution of the recruitment stretch distribution in Fig. 11, driven by Eq. (13), shows that the recruitment stretch distribution was overestimated and decreases so that collagen fiber stretch distribution remodels to the attachment stretch distribution, i.e. via Eq. (6); Fig. 12 illustrates the evolution of the recruitment stretches over the aneurysm sac from the initial pressurized state to the the initial homeostatic state. Figure 13 illustrates F WSS at t = 0 and the spatial distributions of WSS at t = 0 and t = 2 . In regions where the WSS is below 0.5, maximum degradation occurs (indicated by F WSS =1). It can be seen that the whole aneurysm sac enlarges in size and a prominant secondary bleb develops on the upstream region of the aneurysm sac where a region of low WSS is located. This is facilitated by a feedback mechanism whereby the regions of low WSS enlarge in size as the aneurysm enlarges (compare Fig. 13b, c). Figure 14 illustrates F AR at t = 0 and the spatial distributions of WSSAR at t = 0 and t = 2 . Regions of high WSSAR are associated with oscillatory flow and we hypothesize that in these regions the endothelium will have an irregular morphology. In these locations, i.e. the white regions of Fig. 14a where F AR = 1 , degradation is maximal. For this example, F AR does not have localized maxima and consequently the whole IA sac is seen to enlarge without the development of focal secondary blebs.

Enlargement
The impact on the blood flow velocity field of distinct WSS-related degradation hypotheses is evident in Fig. 15, which compares the velocity streamline from the initial homeostasis (common to all degradation methods) to that of each resulting configuration. Figure 16 illustrates the collagen fiber stretch distributions in the final stabilized state for low-WSS driven enlargement (upper) and high-WSSAR driven enlargement (lower). As material is degraded and the IA sac enlarges, the fiber stretches increase to maintain mechanical equilibrium. As the stabilization mechanism is a point-wise relation, see Eq. 15, the mechanobiological equilibrium at the final homeostatic state contains a spatially-variable attachment stretch field. It can be seen that the magnitudes of the attachment stretch distribution (see Fig. 16) have increased (compare with Fig. 10). More specifically, to stabilize the aneurysm, the attachment stretch distribution evolves from a spatially homogeneous distribution where min/ mode/max attachment stretches of the distribution are 1.01/1.05/1.1 to a spatially heterogeneous distribution where these values increase up to 1.07/1.12/1.17, respectively.

Discussion
We have presented the first Fluid-Solid-Growth framework to model IA growth and stabilization for personalized (image-based) IA geometries. Two illustrative scenarios to drive IA enlargement are presented: low WSS and complex, oscillatory flow. We propose a novel approach to link endothelial morphology (aligned vs. irregular) to a novel pulsatile flow metric (WSSAR) and subsequently localized degradation of the tissue. Moreover, the model integrates a mechanism to account for aneurysm stabilization, i.e. adaption of the adventitial collagen fabric via remodelling of the collagen fiber attachment stretch distribution.
The model is fully implemented into Sim4Life, a state-ofthe-art simulation platform for computational life sciences (Neufeld et al. 2013). To our knowledge, Sim4Life is the first fully integrative framework for modeling IA evolution of its type: it incorporates user-friendly tools which span from image segmentation to simulation of IA enlargement.
We utilized an isochoric split of the deformation gradient for our finite element (FE) model. It is recognized that this can cause problems for fiber-reinforced materials that lead to FE simulations producing unrealistic behaviour (volumetric swelling) of a material (Sansour 2008;Helfenstein et al. 2010;Gültekin et al. 2019). Essentially, from an energy minimization perspective, at a critical fiber stretch, it becomes energetically more favourable to swell the material as opposed to stretch along the fiber direction (Zdunek et al. 2014). Interestingly, the (isochoric) fiber stretch at which this occurs can be calculated by consideration of the relative stiffnesses of the bulk modulus and embedded fiber (Zdunek et al. 2014). If a material is modeled to have a stiffness that increases exponentially, then-for physiological consistent Fig. 9 Evolution of the displacement field from the initial loaded geometry (top) to the (initial) homeostatic state (bottom) material parameters-this issue occurs at relatively small deformation, e.g. fiber stretches < 1.2 . However, if a material has constant stiffness, then a bulk modulus value can be chosen so that the issue does not occur (Zdunek et al. 2014). The constitutive model that we adopt for collagen has a constant stiffness once all fibers are recruited and thus a bulk modulus can be chosen appropriately to maintain incompressibility whilst using an isochoric-split of the deformation gradient; volumetric changes were within 2%, see Fig. 17.
The constitutive model for collagen assumed a linear relationship between the 1st Piola-Kirchoff stress and stretch for individual fibers. A nonlinear relation is applied to model the deformation of the fibers and the total strain-energy of the fiber ensemble requires an integral over the recruitment distribution (Hill et al. 2012). Whilst an exponential strain energy function is often used to capture the gross mechanical response of collagen, the nonlinear behaviour is a consequence of individual fibers being recruited to load over a range of recruitment stresses while the individual fibers have a linear response (see discussion in Zhang et al. (2016). This is found to be suitable for modelling the response of arteries (Bevan et al. 2018). It is of note that the triangular recruitment stretch distribution (and linear fiber response) yields analytic expressions for the total stress of the fiber ensemble,and is therefore straightforward to implement into FE codes.
A typical healthy arterial wall is composed of three layers (intima, media and adventitia-some regions contain internal and external elastic lamina to separate the layers) and several components (ECs, elastin, collagen, proteoglycans, fibrilin, nerves, fibroblast, vasa vasorum, etc.). We assumed that the diseased tissue region features a loss of internal elastin lamina and a degraded extracellular matrix. The constitutive model for reinforced collagenous fiber we assume here is relatively simple to account for the full complexity of the wall constituents and physiological processes which happen in the artery. They could be augmented with models of passive and active SMCs response (Bhogal et al. 2019) (which play an important role in aneurysm formation (Kondo et al. 1998)), angle dispersion of collagen fibers , and macrophage infiltration (proinflammatory signaling-which are linked to flow and drive aneurysm formation (Harvey et al. 2018) and progression (Frösen et al. 2019), and anisotropic volumetric growth models (Grytsan et al. Fig. 10 Evolution of the collagen fiber stretch from the initial loaded geometry (first row), to the initial homeostatic state (second row). The first, second and third columns illustrate the minimum, mean and maximum collagen fiber stretches (of the fiber stretch distribution), respectively Fig. 11 Evolution of the recruitment stretch distribution (from initial prescribed values) to achieve a homeostatic state 2017). In fact, recent experimental research ) observes that the tissue has a much more complex heterogeneous structure, with regions that have an atherosclerotic nature. More sophisticated models are needed which account for this complex pathophysiology.
The clinical geometry is obtained under physiological loading-however here we considered it to be the unloaded Fig. 12 Evolution of the fiber recruitment stretch from the initial deformed geometry (first row), to the initial homeostatic state (second row). The first, second and third columns depict the minimum, mean and maximum collagen recruitment stretches, respectively configuration. A pre-stressing algorithm (Weisbecker et al. 2014) is required to compute the actual unloaded configuration, such that the given input geometry is the loaded configuration. Note also that due to the lack of imaging resolution, the wall thickness has to be artificially defined; here we adopted a uniform thickness.
Experiments show that blood exhibits a Newtonian behaviour in most large arteries (Ku 1997), and a non-Newtonian behavior in vessels with a diameter below 100 m (Popel and Johnson 2005). Several publications debate the impact of the flow model on numerical modeling because of the particulare nature of blood (Blanco et al. 2006;Xiang et al. 2011;Evju and Mardal 2015;Saqr et al. 2019). The general consensus is that Newtonian and non-Newtonian models are useful, and that their applicability depends on the size of the arterial lumen. Even though (Xiang et al. 2012) reports that Newtonian model can overestimate WSS, (Cebral et al. 2011(Cebral et al. , 2016 finds no evidence to prefer the non-Newtonian model over the Newtonian one. Because the arteries of the circle of Willis are still large with respect to the RBCs size, we followed Steinman and Pereira (2019) in this work and modeled blood as a Newtonian flow, with constant viscosity .
ECs morphology is hypothesized to affect wall remodeling (Kaneko et al. 2018). Permeability dysfunction leads to transmigration of leucocytes into the wall and degradation of the extracellular matrix (Babu et al. 2015). Here we proposed a novel approach to predict EC morphology and link it to localized degradation of the wall. For illustration, we implemented this with a novel pulsatile flow metric that quantifies the oscillatory nature of the flow, i.e. WSSAR (Krishna et al. 2020). However other pulsatile flow metrics could be employed for the same purpose, e.g. OSI or high WSS, which were recently correlated to proinflammatory signaling and aneurysm initiation (Frösen et al. 2019). Furthermore, EC morphology is linked to cyclic deformation of the wall (Watton et al. 2011a) and hence the algorithm could be enhanced to include this additional mechanobiological influence. The influence of cyclic deformation on remodelling can be easily included: our framework currently simulate the evolution of systolic and diastolic fields as in Watton et al. (2011a) and already incorporates a monolithic fluid-structure interaction solver.
The mechanisms that lead to IA formation, growth and rupture are incompletely understood. From a CFD perspective, each of these phases appears to be correlated to non-physiological low WSS, high WSS, oscillatory flow and combinations of them (Staarmann et al. 2019). Follow up imaging depicting enlarging IAs obtained from clinical databases ) will assist in generating and testing G&R hypotheses in future applications of the framework.
There is a need to move beyond phenomenological correlations that relate flow metrics to IA stability/growth. More specifically, mechanistic understanding of the governing mechanobiology of IAs is needed to improve management of the disease (Robertson and Watton 2012). We envision that the application of the presented in silico Fig. 16 Final homeostatic state of the collagen fiber stretch under the low-WSS hypothesis (first row), and the WSSAR hypothesis (second row). Minimum, mean and maximum collagen stretches are plotted in the first, second, and third column, respectively. Note: in both scenarios, mechanobiological equilibium is achieved via a spatially heterogeneous attachment stretch distribution Fig. 17 Final pressurized configuration following formation of a secondary bleb. The Jacobian is close to 1 and hence near-incompressibility is achieved framework to in vitro and in vivo (animal/human) models, e.g. (Kaneko et al. 2018;Cebral et al. 2015Cebral et al. , 2016Wang et al. 2018;Cebral et al. 2018Cebral et al. , 2019Frösen et al. 2019), will be a step towards this goal.

Summary
This paper presents the first integrative fluid-solid-growth modelling framework with application to modelling enlargement/stabilization of patient-specific IAs. It will provide new mechanistic insight into the mechanobiology of IA disease.