Adaptive Remodeling in the Elastase-Induced Rabbit Aneurysms

Abstract

Background

Rupture of brain aneurysms is associated with high fatality and morbidity rates. Through remodeling of the collagen matrix, many aneurysms can remain unruptured for decades, despite an enlarging and evolving geometry.

Objective

Our objective was to explore this adaptive remodeling for the first time in an elastase induced aneurysm model in rabbits.

Methods

Saccular aneurysms were created in 22 New Zealand white rabbits and remodeling was assessed in tissue harvested 2, 4, 8 and 12 weeks after creation.

Results

The intramural principal stress ratio doubled after aneurysm creation due to increased longitudinal loads, triggering a remodeling response. A distinct wall layer with multi-directional collagen fibers developed between the media and adventitia as early as 2 weeks, and in all cases by 4 weeks with an average thickness of 50.6 ± 14.3 μm. Collagen fibers in this layer were multi-directional (AI = 0.56 ± 0.15) with low tortuosity (1.08 ± 0.02) compared with adjacent circumferentially aligned medial fibers (AI = 0.78 ± 0.12) and highly tortuous adventitial fibers (1.22 ± 0.03). A second phase of remodeling replaced circumferentially aligned fibers in the inner media with longitudinal fibers. A structurally motivated constitutive model with both remodeling modes was introduced along with methodology for determining material parameters from mechanical testing and multiphoton imaging.

Conclusions

A new mechanism was identified by which aneurysm walls can rapidly adapt to changes in load, ensuring the structural integrity of the aneurysm until a slower process of medial reorganization occurs. The rabbit model can be used to evaluate therapies to increase aneurysm wall stability.

Appendices

Appendix 1: Fiber dispersion model

The collagen component is treated as a collection of extensible fibers with a distribution of orientations. In general, a single collagen fiber will have a wavy conformation in reference configuration κ_0, Fig. 2(d). Using local Cartesian coordinates with unit base vectors (e₁, e₂) (Fig. 2(c)), a direction is associated with undulated fiber m₀ in reference configuration κ₀ can be decomposed to m₀ = cosθe₁ + sinθe₂. The undulated fiber is assumed to become load bearing when the material element in direction m₀ has undergone a stretch of λ_a, Fig. 2(d). This activation stretch is an additional material property determined from average tortuosity of fibers. Then the fiber commences stretching when λ_f = λ_a . The true stretch of the fiber λ_t is then λ_t = λ_f/λ_a, where λ_f can be obtained from λ_f² = C : m₀ ⨂ m₀.

The distribution of fibers was represented by probability density function (PDF) of fibers at arbitrary orientation m₀ in the reference configuration. Motivated by the MPM data, only collagen fibers in 2D plane were considered. The PDF was represented by a bimodal von Mises distribution [67].

$$ \rho \left(\theta \right)=\frac{1}{2{I}_0\left({b}_1\right)}{e}^{b_1\mathit{\cos}2\left(\theta -{a}_1\right)}+\frac{1}{2{I}_0\left({b}_2\right)}{e}^{b_2\cos 2\left(\theta -{a}_2\right)}\kern1em \theta \in \left[0,\pi \right) $$

(10)

where a₁ and a₂ are symmetrical angles of two peaks, b₁ and b₂ are concentration parameters. I₀(b) is the Bessel function of order zero defined by

$$ {I}_0(b)=\frac{1}{\pi }{\int}_0^{\pi }{e}^{bcos\theta} d\theta $$

(11)

and the distribution is normalized such that [68],

$$ 1=\frac{1}{\pi }{\int}_0^{\pi }{\rho}_i\left(\theta \right) d\theta $$

(12)

The four parameters (a₁, a₂, b₁, b₂) of Eq. (10) of each case were determined using maximum likelihood estimation by function MLE in Matlab (MathWorks Inc. Natick, MA) [51] based on collagen fiber orientations measured from the MPM stacks using ctFIRE. The parameters of ctFIRE were listed in Table 4, with more details on the functions found at https://loci.wisc.edu/software/ctfire.

Table 4 Input parameters of ctFIRE

The strain energy function of collagen in each layer is then constructed by integrating response of fibers over all directions [40, 50]. The strain energy function of collagen in each layer is of the form

$$ {W}_{fi}=\frac{1}{\pi }{\int}_0^{\pi }{\rho}_i\left(\theta \right){w}_{fi} d\theta $$

(13)

where w_fi is the strain energy of fibers in layer i, and defined as

$$ {w}_{fi}=\frac{\beta_i}{2}{\left({\lambda_t}^2-1\right)}^2\cdotp H\left({\lambda}_t-1\right) $$

(14)

where the Heaviside function is used to switch the fiber contribution off when under compression

$$ H\left({\lambda}_t-1\right)=\left\{\begin{array}{c}0\kern1em if\ {\lambda}_t\le 1\\ {}1\kern1em if\ {\lambda}_t>1\end{array}\right. $$

(15)

Appendix 2: Determination of material parameters from uniaxial tension experiments

Material constants in the structurally motivated constitutive model were determined from multiphoton imaging and uniaxial mechanical testing. Assuming the tissue is incompressible, necessarily the deformation is isochoric with principal stretch values satisfying λ₁λ₂λ₃ = 1. The corresponding deformation gradient tensor and Cauchy Green tensors are

$$ \left[\boldsymbol{F}\right]=\mathit{\operatorname{diag}}\left[{\lambda}_1,{\lambda}_2,\frac{1}{\lambda_1{\lambda}_2}\right],\kern0.75em \left[\boldsymbol{C}\right]=\left[\boldsymbol{b}\right]=\mathit{\operatorname{diag}}\left[{\lambda_1}^2,{\lambda_2}^2,\frac{1}{{\lambda_1}^2{\lambda_2}^2}\right] $$

(16)

where λ_i is the principal stretch in direction e_i. Coordinate axis are chosen such that the loading direction is e₂ and e₃ is orthogonal to the sample. Then the principal components of the Cauchy stress tensor can be written as

$$ {\sigma}_1=-p+{\lambda_1}^2\sum \limits_{i=1}^N{\gamma}_i\left(1-{\mu}_i\right){\alpha}_i+\frac{2}{\pi}\sum \limits_{i=1}^N{\gamma}_i{\mu}_i{\lambda_1}^2\left\{{\int}_0^{\pi }{\rho}_i\frac{\partial {w}_{fi}}{\partial {\lambda_f}^2}{\mathit{\cos}}^2\theta d\theta \right\} $$

(17)

$$ {\sigma}_2=-p+{\lambda_2}^2\sum \limits_{i=1}^N{\gamma}_i\left(1-{\mu}_i\right){\alpha}_i+\frac{2}{\pi}\sum \limits_{i=1}^N{\gamma}_i{\mu}_i{\lambda_2}^2\left\{{\int}_0^{\pi }{\rho}_i\frac{\partial {w}_{fi}}{\partial {\lambda_f}^2}{\mathit{\sin}}^2\theta d\theta \right\} $$

(18)

$$ {\sigma}_3=-p+\frac{1}{{\lambda_1}^2{\lambda_2}^2}\sum \limits_{i=1}^N{\gamma}_i\left(1-{\mu}_i\right){\alpha}_i $$

(19)

By setting σ₃ = 0, an expression for the Lagrange multiplier p can be obtained. An implicit relation between λ₁ and λ₂ can then be obtained by setting Eq. (18) to be zero. Using Eq. (17) along with this implicit relation and the experimental data (relation between σ₁ and λ₁), the parameters α_i and β_i can be obtained from minimizing the least-squares objective function. The interior-point optimization algorithm of FMINCON function in Matlab (R2017a, MathWorks Inc., Natick, MA) was used for this purpose.

Appendix 3: Summary of hemodynamics

Computational fluid dynamic studies were performed in case specific 3D reconstructed geometries of the aneurysm sac and surrounding vasculature to evaluate the consistency in the flow fields across the various models within and across the time points. Specific methods are summarized below and aneurysm metrics (including geometric, flow, and wall shear stress) for each case are given in Table 5.

Table 5 Aneurysm metrics

All aneurysms had a relatively high aspect ratio (AR = height/neck diameter), ranging from 1.79 to 3.56, Table 5, with a bottle neck factor (BF) of 0.80 to 1.53. There was no statistically significant difference of AR or BF at different time points. Three cases, chosen to span a large range of AR and BF, are shown in Fig. 13. As a result of the high AR, the magnitude of wall shear stress was consistently low within the aneurysms, Table 5. For example, the time and space averaged magnitude of wall shear stress (TSAWSSM) had an average of 2.95 ± 1.29 dynes/cm², five times less than the parent artery. The TSAWSSM was less than 7.0 dynes/cm² for all cases.

Fig. 13

Hemodynamic results for three representative cases. Each row represents one of three case chosen to span a large range of AR and BF values found in the rabbit models. (a)-(d) Streamlines at four different time points of the cardiac cycle. (e) Time averaged wall shear stress. (f) Maximum wall shear stress

The qualitative flow featues were similar across all cases, with two transient vortices present in the aneurysm sac, although the duration varied beween cases, Fig. 13. The inflow jet entered at the distal neck, then proceeded to impinge on the proximal wall, splitting into a slower vortex in the dome and a faster vortex near the neck. The flow generally exited through the proximal half of the neck plane. As a result of the concentrated inflow at the distal neck, the time average magnitude of wall shear stress was elevated at this location as exemplified in Fig. 13. The strength and location of jet impingement varied over the aneurysm wall and between cases. The effect of this local feature of the flow field on endothelial cell coverage, for example, is a topic of ongoing investigation.

Methods

Of the 22 aneurysm models created in this study, 3DRA data for 14 cases were available for 3D digital reconstruction. Stereolithography (STL) data of the lumen were produced from the 3DRA data using the commercial software Mimics (Materialise NV, Leuven, Belgium). Finite element models of each lumen were created using a custom-developed software [69]. The finite element domain consisted of linear tetrahedral elements having an average element edge-length of 0.12 mm. The model domain included the ascending and descending aorta, the left subclavian (LS), left common carotid artery (LCCA), and the quadrification distal to the right subclavian arteries.

Blood was approximated as incompressible with constant kinematic viscosity of 3.33e-6 m²/s. The unsteady Navier-Stokes equations were solved over the computational domain assuming a rigid model of the vasculature so the no-slip condition as applied at the vessel wall. Animal specific doppler ultrasound (DUS) waveform data obtained proximal to the rabbit aneurysms in the brachiocephalic artery (BA) were used to set an inlet flow rate for the model. Briefly, the ultrasound velocity was taken as the peak velocity of a zeroeth order Womersley profile, from which the flow rate in the BA could be approximated. Separate DUS data were available in the LCCA and descending aorta and used to similarly approximate those flow rates. The flow split across outlets distal to the aneurysm were assumed to scale by Murray’s Law (i.e., the flow splits scales by the artery radius cubed). At each outlet, a Womersley velocity profile was then applied consistent with this flow division. Finally, the pressure was set via a traction free boundary condition at the ascending aortic inlet. The solution was discretized in time with timesteps of 0.001 s and obtained from FeFlo, the finite element solver described in [69].