Aneurysm and biomechanics of vessels

Cardiovascular diseases are the cause of numerous deaths. Among these, abdominal aortic aneurysms are very dangerous for which adequate techniques for diagnosis and prevention have not yet been developed. The implementation of numerical analyzes can help analyze the influence of geometric and material parameters on the risk of aneurysm rupture.

This report considers a simplified model from a geometric, hemodynamic and structural point of view and some numerical fluid dynamics simulations and structural analysis are carried out on a geometry representing an aneurysmatic vessel.

The results of the structural analysis confirm that the maximum stress increases with the diameter and with the reduction in the thickness of the vessel.

The following report is an extract from a project carried out for the course of Modeling and Simulation of Physiological Systems, held for Medical Engineering at the University of Rome Tor Vergata.

Introduction

In recent decades, much of the research on the cardiovascular system has focused on the relationship between local hemodynamics and the onset of pathophysiological phenomena. Cardiovascular diseases are the cause of around 48% of annual deaths in Europe, costing over 100 billion euros [1].

Today it is believed that the onset of some pathologies is linked to particular altered fluid dynamic patterns. For example, it is believed that where there is stagnation of blood linked to low speeds, this has a greater impact on the onset of inflammatory processes.

The hemodynamics of the cardiovascular system is very complex linked to the turbulent nature of the flow where even a small change in blood flow affects the biomechanics of the heart and vessels.

In the case of aneurysms this possibility is strongly supported by the fact that there are no reliable quantitative methods used to concretely validate the risk of rupture and this still remains an open problem. Despite strong developments in surgical and diagnostic imaging techniques, mortality and morbidity rates remain very high.

Currently in the literature there is agreement in looking for a reliable evaluation criterion associated with the potential rupture and not with the maximum size. Furthermore, being a phenomenon that subjects the wall to strong mechanical stress, it is highly probable that a criterion based on the quantification of stresses is optimal [2]. However, there is currently no method capable of obtaining valid in vivo measurements.

This report considers a highly simplified model from a geometric, hemodynamic and structural point of view and some numerical simulations are carried out. Therefore, both a fludio-dynamic simulation is carried out within a geometry representing an aneurysmatic vessel and an analysis of the uncoupled fluid-structure interaction.

Background

An aneurysm is a pathological dilation of a blood vessel, generally of an artery. This results in various health problems and risks, typically linked to the weakening of the dilated wall which risks breaking, leading to massive bleeding and consequent death.

The aorta is the largest and most important arterial vessel in the human body which, originating from the heart, spreads blood to every organ and district. There are several localization sites for an aneurysm but typically we speak of an abdominal, thoracic or intracranial aneurysm.

Abdominal aneurysm

The abdominal aortic aneurysm (AAA) is a localized dilation of the infra-renal aorta. AAA results from changes in the structure of the aortic wall, including the loss of smooth muscle cells and the degradation of the extracellular matrix. The swollen part is fragile and is likely to break relatively easily.

Under normal conditions, the diameter of the abdominal aorta is about 2 cm and we speak of an aneurysm when the swelling is greater than 3 cm. The AAA is considered operable when it reaches a diameter of 5.5 cm [4].

The cause of origin is still unknown but the current literature agrees in recognizing the fundamental role of several factors such as aging, atherosclerosis, hypertension, cigarette smoking, vasculitis and genetic predisposition.

In the last decade, numerous numerical models have been proposed to evaluate and predict AAA rupture showing how different information such as wall thickness, presence of intraluminal thrombi, wall homogeneity and use of an appropriate constitutive model [5].

Analysis

This report addresses some aspects of the biomechanics of the AAA. Several simplifications follow, both from a hemodynamic and structural point of view. It is therefore considered a geometric approximation of a fusiform AAA represented in fig. 1c and described below on which various simulations are carried out aimed at considering the biomechanical behavior of the blood flow and the vascular wall (initially considered rigid).

These numerical analyzes are implemented in the Comsol Multiphysics simulation environment.

Poiseuille flow

Facing the description of blood motion we can talk about the Navier Stokes model for viscous fluid like the incompressible Newtonian where, indicated by 𝐮, the velocity field is:

\begin{equation}
\left\{\begin{array}{l}
\rho \left(\frac{\partial \mathbf{u}}{\partial t}+ \nabla \mathbf{u} \cdot \mathbf{u}\right)=-\nabla p+\mu \Delta_{2} \mathbf{u}+\mathbf{f} \\
\nabla \cdot \mathbf{u}=0
\end{array}\right.
\end{equation}

These equations represent precisely the second principle of dynamics and the conservation of mass.

Hollow cylinder

In Poiseuille motion we consider a unidirectional motion along the axis of the cylinder of the type:

\begin{equation}
\mathbf{u}=(u, 0,0)
\end{equation}

This motion, together with the continuity equation, express how the derivative of speed in the longitudinal direction is zero:

\begin{equation}
\frac{\partial u}{\partial x}=0
\end{equation}

Furthermore, the volume forces 𝐟 are neglected and a stationary case is considered. These assumptions lead to considerably simplifying the Navier-Stokes equations. In particular, the convective contribution, a source of non-linearity, is no longer present. The spatial dependencies are greatly simplified and, considering cylindrical coordinates with geometry and symmetrical loads with respect to the anomaly, the variables remain:

\begin{equation}
\begin{aligned}
p(x, y, z,) & \rightarrow p(x) \\
\mathbf{u}(x, y, z) & \rightarrow u(r)
\end{aligned}
\end{equation}

Or rather:

\begin{equation}
\mu \Delta_{2} u=\frac{\partial p}{\partial x}
\end{equation}

Velocity profile

This equation, integrated twice and considering the condition of zero velocity at the boundary, translates into the analytical solution of the Poiseuille motion:

\begin{equation}
u(r)=\frac{1}{4 \mu} \frac{P_{L}-P_{0}}{L}\left(r^{2}-R^{2}\right)
\end{equation}

That is, there is a parabolic trend where the speed is maximum on the axis of the channel:

\begin{equation}
u(r)=\frac{1}{4 \mu} \frac{\Delta p}{L} R^{2}
\end{equation}

This problem is tackled both analytically, graphing its trend, and numerically. The comparison is present in fig. 2a.

In CFD modeling, a 2 cm wide channel with a length of 20 cm [2] and a pressure of approximately 90 mmHg [6], or 12000 Pa, with a pressure drop of 2 Pa [7] is considered. These parameters are representative of the intrarenal aorta.

Aneurysmal vessel

The fluid domain is varied by approximating an AAA by means of an axisymmetric domain. This domain is defined starting from the edge of the previous fluid domain by adding, within the Comsol environment, a sphere of such a radius as to return the aneurysmatic radius overall. That is, let 𝚍 be the diameter of the vessel and 𝑎𝑛𝑒𝑢𝑟𝑦𝑠𝑚 _ 𝐷 the diameter of the aneurysm, the sphere is drawn through a semicircle of radius (𝚊𝚗𝚎𝚞𝚛𝚢𝚜𝚖_𝙳 – 𝚍) ∕ 𝟸 centered on the outside of the vase, halfway up.

Internal pressure

Subsequently, the two geometries are joined through a Boolean operation and the connection points are rounded off by means of a joint with a radius equal to 40% of the radius of the aneurysm. This allows both to have a more realistic geometry and to avoid numerical singularities linked to discontinuity.

The same conditions apply to the previous case for the inlet, the outlet and the wall with no slip condition. However, the inlet is considered in the upper section and the outlet in the lower section, referring to the anatomical situation.

The flow lines are present in fig. 3 and a three-dimensional render is present in the appendix.

We also extract further data such as the pressure along the radius of the vessel at different heights, in fig. 4 and the pressure on the wall, fig. 5. This information can be extracted both via the graphical interface and within Matlab (see appendix).

Structural analysis

Once the pressure load is known, a structural analysis is carried out, reporting this load as a condition at the edge. Therefore, the two simulations remain decoupled but the pressure loads are reported.

To implement the analysis it is necessary to define a structural domain. Therefore, a procedure is implemented, within Comsol, such as to maintain the edge of the internal geometry (previously described for the fluid domain) but to return a shell of constant thickness. The domain is then opened leaving only the wall and then the latter is thickened towards the outside with a thickness of 1.5 mm [8].

The load is defined using a point function starting from the numerical data extracted from the previous simulation. Therefore, the tabulated data is loaded and then interpolated via Comsol providing a function 𝑝𝑤𝑎𝑙𝑙(z) used as an edge condition on the inner wall of the vessel.

An almost incompressible Neo-Hookean hyperelastic model is then used with parameters 𝜇=1.7⋅10⁶ Pa, 𝑘=34,7⋅10¹⁰ Pa [9] and 𝜌 = 104 kg / m3. The Neo-Hookean model is a hyperelastic material model that can be used to describe the stress-strain relationship for a material with non-linear behavior in finite elasticity. In particular, this model can describe the behavior from blood vessels with sufficient accuracy [10].

Neo Hooke

In finite elasticity the material is defined by the deformation energy density, a function of the elastic elastic state typically such as to be described by the right Cauchy-Green tensor.

This allows us to calculate the stress in local coordinates using the second Piola-kirchoff tensor:

\begin{equation}
\mathbb{S}=2 \frac{\partial W_{s}}{\partial \mathbb{C}}
\end{equation}

Which is related to stresses such as:

\begin{equation}
\boldsymbol{\sigma}=\boldsymbol{J}^{-1} \mathbb{F} S \mathbb{F}^{T}
\end{equation}

Where:

\begin{equation}
J=\operatorname{det} \mathbb{F}
\end{equation}

Or as a function of the Green-Lagrange tensor:

\begin{equation}
\mathbb{E}=\frac{\mathbb{C}-\mathbb{I}}{2}
\end{equation}

Where:

\begin{equation}
\mathbb{E}=\frac{1}{2}\left(\mathbb{F}^{T} \mathbb{F}-\mathbb{I}\right)
\end{equation}

And F is the strain gradient:

\begin{equation}
\mathbb{F}=\mathbb{I}+\nabla \mathbf{u}
\end{equation}

Strain energy

In the Neo-Hookean model [11], therefore (10) is defined the deformation energy density as a function of a volumetric contribution and an isochoric contribution:

\begin{equation}
W_{\mathrm{s}}=W_{\text {iso }}+W_{\mathrm{vol}}
\end{equation}

Where the volumetric contribution is defined as a function of the Bulk modulus and the volumetric elastic deformation 𝐽𝑒𝑙 as:

\begin{equation}
W_{\text {vol }}\left(J_{{el}}\right)=\frac{1}{2} \kappa\left(J_{\mathrm{el}}-\mathbb{I}\right)^{2}
\end{equation}

While the isobaric contribution is defined as a function of the first isochore invariant and the second Lamè parameter:

\begin{equation}
W_{\text {iso }}=\frac{1}{2} \mu\left(\overline{\mathbb{I}_{1}}-3\right)
\end{equation}

Vascular wall

This constitutive model allows to describe the shell representing the vessel wall. Known the pressure field acting on the internal wall of the vessel, obtained from the fluid dynamics analysis (Fig. 5) is applied as a condition to the edge in the structural analysis.

Therefore we analyze the displacement present in fig. 7b and the wall stress shown in fig. 6. Stress is calculated along the three directions along a curvilinear abscissa passing through the inner edge of the wall. The result is therefore a maximum value in the tangential direction of 0.2 MPa while in the radial direction it is less than 0.5 MPa.

The maximum displacement is less than a millimeter and is concentrated in the connection area between the cylindrical geometry and the sphere.

From fig. 7a it is possible to observe the resultant of the stresses according to Von Mises and it is observed that there is a greater concentration of stresses on the connection area between the circumference and the vessel wall. This area is in fact the one most subjected to elongation as can be seen from the false color diagram in fig. 7a.

These results, although the result of a decoupled and very simplified analysis, show how the swollen geometry has a strong impact on the amplification of stresses.

Geometric variability

From the Comsol model you then pass to the corresponding model in Matlab through which it is possible to extract values and parameters of interest. In particular, the fluid dynamic analysis is repeated starting from a variation in the diameter of the vessel and the aneurysm according to pathophysiological values [12, 13].

The pressure field on the inner wall is then extracted and passed to a structural simulation based on a similarly varied geometry. The results are present in fig. 9 where the maximum stress on the wall is considered as a function of the geometric parameters.

In particular, these results are in agreement with the literature showing how the value of the maximum stress increases as the diameter of the aneurysm increases. For the wall thickness the trend is inverse, that is, the thinner the wall is, the higher the stress will be.

Change in blood pressure

A non-stationary fluid dynamics simulation is then implemented, the result of an attempt at a non-convergent full-coupled FSI. In this simulation we consider a time-varying inlet and outlet. For the inlet a velocity curve is chosen and for the outlet a pressure curve [14, 2] typical for the arterial compartment of interest (indicated in the appendix).

Time dependent analysis

These curves are defined in Comsol by points and then the interpolation and transformation into a function of time is entrusted to Comsol. The result, in terms of pressure on the wall, is shown in fig. 10 where it shows itself at different moments in time. Although the curves appear flat, they show, observing them in the correct scale (in Pascal), a trend similar to the stationary case.

The one at the instant of time where the maximum pressure is reached (0.5 s) is selected and is used again for structural analysis. The results (fig. 11) show how, having increased the average pressure, the stress on the wall increases, even though the trend in space is very similar to the previous case.

Conclusions

The results of the analysis confirm that the maximum stress increases with the diameter of the aneurysm and as the thickness of the vessel decreases. Moreover, this stress is strongly linked to the pressure curve. The geometric irregularities also provide an amplification of the tensions on the aneurysmal wall.

Furthermore, the implementation of an FSI model would lead to more precise results as it would allow to take into account the high deformability of the vessels which confers compliance with the cardiovascular system.

References

• [1]  Michal Tendera. “How much does Europe invest in the treatment of cardiovascular diseases?The opinions ex- pressed in this article are not necessarily those of the Editors of the European Heart Journal or of the Eu- ropean Society of Cardiology.”
• [2]  Christine M. Scotti et al. “Wall stress and flow dynam- ics in abdominal aortic aneurysms: finite element anal- ysis vs. fluid–structure interaction”. en.
• [3]  Nathan M. Wilson, Ana K. Ortiz, and Allison B. Johnson. “The Vascular Model Repository: A Pub- lic Resource of Medical Imaging Data and Blood Flow Simulation Results”.
• [4]  Natzi Sakalihasan et al. “Abdominal aortic aneurysms”. en. In: Nature Reviews Disease Primers 
• [5]  Mamadou Toungara, Gregory Chagnon, and Chris- tian Geindreau. “NUMERICAL ANALYSIS OF THE WALL STRESS IN ABDOMINAL AORTIC ANEURYSM: INFLUENCE OF THE MATERIAL MODEL NEAR-INCOMPRESSIBILITY”.
• [6] C. J. Mills et al. “Pressure-flow relationships and vascular impedance in man”. en. In: Cardiovascular Re- search 
• [7]  Peter C. Luchsinger et al. “Instantaneous Pressure Dis- tribution Along the Human Aorta”.
• [8]  José F. Rodríguez et al. “Mechanical Stresses in Abdominal Aortic Aneurysms: Influence of Diame- ter, Asymmetry, and Material Anisotropy”.
• [9]  Mohammed Yahya. “Three dimensional finite-element modeling of blood flow in elastic vessels: effects of arterial geometry and elasticity on aneurysm growth and rupture”.
• [10]  Benjamin Owen et al. “Structural modelling of the car- diovascular system”.
• [11]  Comsol Multiphysics. “Nearly incompressible hyper- elastic materials”.
• [12]  A Thapar et al. “Internal or external wall di- ameter for abdominal aortic aneurysm screening?” [13]  Todd E. Rasmussen and John W. Hallett. “Inflam- matory Aortic Aneurysms: A Clinical Review with New Perspectives in Pathogenesis”.
• [14] M.L. Raghavan and David A. Vorp. “Toward a biome- chanical tool to evaluate rupture potential of abdom- inal aortic aneurysm: identification of a finite strain constitutive model and evaluation of its applicabil- ity”.

Appendix

Poiseuille

% load parameters
parameters.length=0.2; 
parameters.diam=0.02; 
parameters.R=parameters.diam/2; 
parameters.viscosity=3.5*1e-3; 
parameters.P0=12000; 
parameters.PL=11998; %[Pa]
import com.comsol.model.* ; 
import com.comsol.model.util.* % load comsol utility 
model = ModelUtil.create(’Model’); 
path=pwd; % load current directory path
addpath(strcat(path,’\subroutine\’)); % add subrotuine
model.modelPath(path); model.label(’poiseuille.mph’); model=setup_parameters(model,parameters);
model=create_geometry(model); 
model=create_physics(model); model=solve_study(model);
[model,speed, radius]=post_process(model); 
mphsave(model,’poiseuille_mat.mph’) 
pouiseuille_analytic=@(r) (1/(4*parameters.viscosity))*((parameters.PL-parameters.P0)/parameters.length)*(r.^2-(parameters.R*ones(size(r))).^2); % analytic solution: r_span=-parameters.R:0.001:parameters.R;
figure; hold on; 
scale=1e3; 
plot([0,0],[0,max(max(speed))],’--k’,’LineWidth’,0.5) %plot central axis
ylabel(strcat(’Speed’,’ [m/s]’)) xlabel(strcat(’r’,’ [mm]’)) 
plot(scale*r_span,pouiseuille_analytic(r_span),’-’,’Color’, ’#A2142F’,’LineWidth’,2)
plot(scale*radius,speed,’o’,’Color’,’#0072BD’,’LineWidth’,1) 
legend({’Cyl axis’,’Analytc’,’Numerical’})

Post processing

Postprocessing subroutine snippet to extract pressure data on the inner edge of the vessel (𝚋𝚘𝚞𝚗𝚍𝚊𝚛𝚢 4,5,6,7,8 and 9) and pressure on the cut lines. The code also rearranges the curvilinear abscissas updated with the respective indices also the ordinates.

wall_pressure = mpheval(model,’p’,’selection’,[4 5 6 7 8 9],’edim’,’boundary’);
% extract wall pressure on selected boundary % sort longitudinal value
[line,index]=sort(wall_pressure.p(2,:),’ascend’); pressure_wall=wall_pressure.d1(index);
figure; 
plot(line,pressure_wall,’Color’,’#77AC30’,’LineWidth’,1.5) ylabel(strcat(’Pressure’,’ [’,wall_pressure.unit,’]’)); xlabel(’z [m]’); title(’Wall pressure’)
% extract pressure over cutline 1
[U, radius,unit]=mphinterp(model,{’spf.U’,’cln1x’},’dataset’,’cln1’); 
clear index
[space,index]=sort(radius’,’ascend’); x_to_plot=space; pressure_cln1=U’; 
for j=1:length(pressure_cln1(1,:))
pressure_cln1(:,j)=pressure_cln1(index(:,j),j); 
end
% extract pressure over cutline 2
[U1, radius1,unit]=mphinterp(model,{’spf.U’,’cln2x’},’dataset’,’cln2’);
clear index; clear space
[space,index]=sort(radius1’,’ascend’); x_to_plot1=space; pressure_cln2=U1’;
for j=1:length(pressure_cln2(1,:)) pressure_cln2(:,j)=pressure_cln2(index(:,j),j);
end
% compute symmetrization
x_to_plot=[-fliplr(x_to_plot),x_to_plot]; pressure_cln1=[fliplr(pressure_cln1), pressure_cln1]; x_to_plot1=[-fliplr(x_to_plot1),x_to_plot1]; pressure_cln2=[fliplr(pressure_cln2),
pressure_cln2];
figure; scale=1e3; % scale to [mm]
plot(scale*x_to_plot ,pressure_cln1 ,’Color’,’#0072BD’,’LineWidth’,2,’DisplayName’,’cos(x)’); 
hold on; 
plot(scale*x_to_plot1 ,pressure_cln2 ,’Color’,’#A2142F’,’LineWidth’,2)
plot([0,0],[0,max(max(pressure_cln2))],’--k’)%plot central axis ylabel(strcat(’Speed’,’[’,unit(1),’]’)); xlabel(strcat(’r’,’ [mm]’)); legend({’z=0.5*L’,’’,’z=0.98*L’})

Streamline 3D

Three-dimensional flow lines for fluid dynamics analysis with stationary inlets and outlets.

Time dependent boundary conditions

Pressure (a) and speed (b) curves used respectively for the outlet and the pressure inlet referring to the physiological parameters measured in Scotti et al. [2].