Blood Flow Modeling in Cerebral Aneurysm

Mesh	Rupture Site Prediction
Stationary Flow	References
Pulsatile Case

A cerebral aneurysm is an abnormal dilation caused by a weakening of an inner layer of a blood vessel wall. Aneurysms occurs within the brain. The biggest risk is a rupture, which has mortality rate of nearly 40%. The precise prognosis depends on various factors, such as age or the location of the aneurysm.

The goal of this application example is to study a blood flow in a patient-specific aneurysm obtained through a computer tomography (CT) scan and then predict where the aneurysm is most likely to rupture, based on various hemodynamic parameters. The model takes into account the non-Newtonian fluid flow properties of blood. The pulsatile flow data is experimentally obtained by a Doppler ultrasound measurement in the middle cerebral artery. The complexity of the model and geometry make solving the model expensive in both computer time and memory.

This application example is divided into two parts. The first part studies a stationary flow in one patient-specific aneurysm. The second part focuses on the pulsatile flow in the aneurysm and includes the prediction for the rupture site. This application example is based on [Brunátová, 2021].

Many of the animations of the simulation results shown in this application example are generated with a call to Rasterize . This is to reduce the disk space this notebook requires. The downside is that the visual quality of the animations will not be as crisp as without it. To obtain high quality graphics remove or comment out the call to Rasterize .

Units

The standard system of units in the Wolfram Language is the SI system. For blood flow simulations, the centimeter-gram-second (CGS) system is commonly used, and there a mass density is in grams per cubic centimeter. The base unit of pressure is dyne per square centimeter [], which can be converted to standard SI units via 1 =0.1` . The base unit of dynamic viscosity is poise [], which is equal to 1 =0.1` .

More information about units in the context of the finite element method in the Wolfram Language can be found here.

Load the Finite Element package and set the $HistoryLength:

Mesh

The first thing to do is to prepare the mesh. In this case, it is not necessary to create the mesh from scratch, as the mesh already exists and is publically available here. The original mesh is, however, quite large. In the interest of computational time, a downscaled mesh is provided that will speed up the computations. This will, however, negatively impact the accuracy of some results, as discussed later. For a more accurate result, the original mesh should be used.

Import the boundary mesh:

Visualize the boundary mesh:

Show bounds of the mesh:

In order to prescribe boundary conditions on the various parts of the boundary, it is necessary to be able to distinguish between different parts of the boundary. This can be done with element markers. There are three types of element markers, which correspond to different types of elements: mesh element markers, boundary element markers and point element markers.

List all boundary mesh markers:

Show the boundary mesh in different colors corresponding to different element markers:

The red part of the boundary shown in the image on the left represents the inlet. This is where the blood flows into the aneurysm. On the other hand, the green and blue parts of the boundary on the right represent the outlet. This is where the blood flows out of the aneurysm. The white part of the boundary represents the wall.

Set up the Neumann value boundary markers:

The boundary element markers will be used later to prescribe Neumann boundary conditions. The Dirichlet boundary conditions are, however, prescribed at point elements. Thus it is necessary to repeat the procedure for the point elements.

List all point element markers:

Show the point elements of the boundary in different colors corresponding to different point element markers:

Set up the Dirichlet condition boundary markers:

The last thing that is necessary is to pick some reference point for the pressure. If no condition is given for the pressure, then the solution is not unique. To solve this issue, a point is chosen at which a reference value for the pressure is prescribed.

Choose a pressure reference point and show its position in the mesh:

The Dirichlet boundary conditions are prescribed at point elements, but the pressure reference point might not be a point element. Therefore, it is necessary to find the nearest point element.

List all boundary point elements with element marker equal to 1:

Obtain the coordinates of the pointElements:

Find the nearest point element:

Show the position of the new pressure reference point:

As a last step, the full element mesh is created.

Create the full mesh from the boundary mesh:

Stationary Flow

As a first step, to keep things simple, the stationary case is solved. That is easier than to solve for a time-dependent pulsatile inflow case immediately. This has the additional benefit that concepts can be introduced in a simpler setting.

The goal of the stationary flow simulation is to compute the aneurysm's wall shear stress. The main challenge to overcome is the high nonlinearity in the fluid flow model. Not being able to solve highly nonlinear PDEs is not uncommon. In general, there exist two ways to go about this. One is to create a ParametricFunction through, for example, ParametricNDSolveValue and then iteratively increase the amount of nonlinearity. This approach is commonly used to get hyperelastic material models to converge. The second approach, used in this case, is to make the stationary model time dependent and use the time parameter t to ramp up the nonlinearity and integrate until a steady state has been reached. So the equations are set up as time dependent, even though they are used to compute the stationary flow. This allows one to start at the state of rest and slowly ramp up the inflow velocity as a function of time t. The NDSolve family of functions have an adaptive time-step integration method implemented, and NDSolve will automatically increase or decrease the time-step size, which makes it very convenient to use.

Boundary Conditions

There are three parts of the boundary: the wall, the outlets and the inlet. This corresponds to three boundary conditions that have to be specified. Additionally, a condition for the reference pressure is prescribed.

Wall boundary conditions

The so-called no-slip condition is prescribed on the wall. That is, it is required that the velocity on the wall be equal to 0. This is why this condition is called no-slip: the fluid sticks to the boundary and it does not slip.

Define the no-slip condition for the velocity on the boundary:

Outlet conditions

There should be nothing that may interfere with the blood flow at the outlet. Thus 0 traction is prescribed at the outlet. This boundary condition is also called do-nothing boundary condition. Since a Neumann zero boundary condition is the default boundary condition when nothing is specified on a part of the boundary, it is possible to just not specify anything in that case and that would be fine too.

Prescribe the do-nothing condition on the outlet:

Reference pressure

Next, the boundary condition for the pressure is prescribed. This boundary condition is necessary to make the solution unique.

Define the condition for the pressure at the pressure reference point:

Inlet condition

Here the inflow is prescribed, which is composed from several parts.

Inflow profile

The first part is a unit parabolic profile. This is used to prescribe the general profile of the inflow, and it has a shape of a parabola. The parabola has a value of one in the center and is equal to zero at the boundary to match the no-slip condition of the adjacent wall.

It is given as an analytic solution to the Poiseuille flow in a pipe. It is a laminar flow in a cylindrical pipe that has prescribed two different values of pressure at both ends. The pressure difference creates a pressure gradient, which makes the fluid in the pipe flow.

The formula for the inflow profile is given as:

where is a radius of the pipe (in this case of the inlet), and is the center point of the inlet boundary. Note that this formula for the inlet velocity only works because the no-slip condition is on the wall. If a different boundary condition were used, this formula would have to have a different form.

It is necessary to know the radius and the center in order to use the formula. In this case, those values are provided with the mesh.

Define the inlet center and the inlet radius from the provided values:

Now the inlet profile can be computed. The vector , which appears in the definition of the profile, might not lie in the inlet. It is possible that it has a small normal component, which would plague the resulting profile. Thus this potential normal component is subtracted from the vector. This way, a vector that is tangent to the inlet is obtained, which can be used to compute the profile. For this, a unit boundary vector is required. It is possible to once again use a provided vector, but in this case, there already exists a convenient tool that can be used: the NeumannBoundaryUnitNormal function.

Compute the tangential component of the vector by subtracting the normal component from the vector:

Set up the inflow profile:

Inflow damping

The inflow damping function is used to damp the inflow at the beginning, as discussed above. The actual mathematical formulation of the damping function is arbitrary, as long as the function starts from 0, increases until it reaches 1 and is constant thereafter.

Define the damping function:

Plot the damping function on interval [0, 4]:

Inflow scaling

The inflow is scaled with the cycle-averaged cross-section mean velocity. As the name suggests, it is a velocity averaged over the cross section of the inlet and over one cardiac cycle. Its value is 48 .

Define the cycle-averaged cross-section mean velocity:

Full inflow profile

Now it is possible to write down the whole inflow function.

Define the inflow velocity function:

With this function in hand, it is finally possible to prescribe the boundary condition for the inlet.

Prescribe the Dirichlet condition on the inlet:

Wrap all Dirichlet boundary conditions into one variable:

Initial Condition

All boundary conditions are prescribed, but an initial condition has to be given, because the problem is solved as a time-dependent problem. Zero velocity and zero pressure are chosen as the initial condition, as the fluid is desired to be at rest at the beginning of the computation.

Define the initial condition:

Model Setup and Solution

The first thing to do is to define variables.

Define variables:

Next, material parameters are prescribed. The non-Newtonian Carreau model is well suited for describing properties of human blood. It allows the decrease of viscosity when the shear rate increases to be captured, as well as the Newtonian plateau.

Prescribe material parameters:

The equations can be generated using the FluidFlowPDEComponent.

Generate the equations:

Now the equations are solved. For this, the finite element method is used to discretize the equation in the spatial dimensions and then the method of lines is used to deal with the time. The Taylor–Hood element, which is a popular choice for the fluid dynamics, is used. The AccuracyGoal and PrecisionGoal are reduced, and the IDA method with reduced order is used to make the computation faster, but also less accurate.

Furthermore, it is sufficient to store the solution only for the time when the computation ends. For this reason, it is possible to use in the NDSolveValue. This significantly reduces the required memory for storing the returned InterpolatingFunction objects.

Define the time when the computation should end:

The simulation will take some time to finish.

Solve the equations, monitor the progress and measure the time it takes:

The large computational time is expected, as the mesh has many elements and is three dimensional, and the mathematical model is highly nonlinear. All these factors contribute to the long evaluation time. The options given help to reduce the integration time.

Results

Now that the equations are solved, it is possible to take a look at the solution and compute some additional parameters from it. The units of the plotted quantity are written above the legend and are in the CGS system as was discussed at the beginning of this application example.

Velocity

Use vector plots to visualize how the blood flows.

Plot the stationary solution using VectorPlot3D:

It can be seen that the blood flows into the aneurysm through the inlet, then it swirls around the aneurysm until it flows out of the outlets. This behavior can also be observed in the stream plot. Note that it can be seen in the vector plot that the blood flows into the aneurysm through the inlet and flows out through the outlets. This is an expected behavior that cannot be observed in the stream plot.

Plot the streamlines using StreamPlot3D:

A more detailed look can be taken by creating two slices and plotting the norm of the velocity. The slices are defined by a normal vector and a point. This defines unique planes that slice the aneurysm.

Define the point and the normal vector for the first slice:

Define the point and the normal vector for the second slice:

Plot both slices:

Now with both slices defined, it is possible to plot how the norm of the velocity looks on those slices.

Plot the norm of the velocity on the first slice:

Before analysing the plots a list plot of the second slice is created.

Plot the norm of the velocity on the second slice:

It can be seen that the velocity in the center of a large vessel section of the slices is quite small, and at the same time, velocity peaks can be observed near the boundary.

Viscosity

It is now possible to plot the apparent viscosity in the same way the norm of the velocity was plotted. This allows the non-Newtonian properties of blood to be studied. First, however, it is necessary to compute the viscosity, as there is no direct access to it. The first thing to do is to compute the gradient of the velocity, which is then used to compute the strain rate tensor.

Compute the gradient of the velocity:

Compute the strain rate tensor:

Now with the strain rate tensor in hand, the shear rate and then the apparent viscosity can be computed.

Compute the shear rate:

Compute the apparent viscosity:

Plot the apparent viscosity on the first slice:

Plot the apparent viscosity on the second slice:

The plots show the non-Newtonian nature of blood, as the plots show that the viscosity varies from 0.04` to 0.1` . This is a quite significant decrease in the viscosity from the viscosity at zero shear rate of 0.56` . Note that for a Newtonian fluid, both plots would show constant viscosity.

Wall shear stress

The last thing to do with the stationary solution is to compute the wall shear stress from it. The wall shear stress is an important hemodynamic parameter. It represents the tangential part of the traction acting on the wall, and it is considered to be one of the indicators for rupture.

The formula for the wall shear stress is the following:

where is the stress tensor and is the outer unit boundary normal vector. The first part represents the traction acting on the wall, while the second part is the normal component of the traction. If the normal component is subtracted from the traction, the result is the tangential component, which is called the wall shear stress.

Two things should be noted about the wall shear stress. First of all, some sources define the wall shear stress as the norm of the vector , while here it is defined as the vector itself.

Second, the wall shear stress is sensitive to the mesh refinement. This application example focuses more on performance than accuracy, because its goal is to demonstrate concepts and workflow that are associated with solving complex problems arising in computational fluid dynamics. This is why downscaled mesh is used and why the obtained results for the wall shear stress are not as accurate as they could be. For more accurate results that could be used in practice, the problem should be solved on the original, or ideally, even finer mesh.

It is necessary to set up the stress tensor and the unit boundary normal vector to compute the wall shear stress. The stress tensor is given as

where is the pressure, is the identity matrix, is the apparent viscosity and is the strain rate tensor. Most of it was already computed in the previous section, so now it remains to set up the stress tensor and the unit boundary normal.

Define the stress tensor

Obtain the boundary unit normal vector:

Now finally the wall shear stress can be computed and plotted.

Compute the traction acting on the wall:

Compute the wall shear stress:

Plot the wall shear stress:

The wall shear stress plays an important role in predicting where the aneurysm rupture is most likely to occur.

Pulsatile Case

Now it is time to take a look at the pulsatile case. The difference between the pulsatile case and the stationary case is that now the cardiac cycle is taken into account. Most of the things can be reused from the stationary case, because the stationary case was solved as a time dependent problem.

Note that pulsatile case is more expensive to compute both in time and memory than the stationary case.

The stationary solution is no longer needed and therefore it is possible to remove it from memory. This allows to use the same variable names for the solution, namely , and for the velocity and for pressure. It also allows to reuse some parts of the code from the stationary case, which relied on the previously used variable names, such as the material parameters .

Clear the stationary solution from the memory:

Inflow Data

The first thing to do is to import the inflow data. The inflow data used in this application example are real life data obtained by a Doppler ultrasound measurement in the middle cerebral artery. They represent the norm of the velocity of the blood flow during one cardiac cycle. The cardiac cycle is the activity of the human heart from one beat to another.

Import measured inflow data:

Plot the inflow data:

An initial spike can be observed in the plot. This corresponds to the heart beat. The data were obtained by Doppler ultrasound measurement in the middle cerebral artery. Note that this is just a measurement of one specific case. The precise inflow into the aneurysm depends on various factors such as heart rate.

It will be useful to define the length of the measured cardiac cycle.

Define the length of one cardiac cycle:

The inflow data are discrete, hence in order to use them in computations it is necessary to make them continuous. This is possible with either a fit or an interpolation. Both approaches have their advantages and disadvantages. The interpolation is more accurate, but the fit offers a better performance as it is smoother and its evaluation is faster. Since this application example focuses more on performance than accuracy, the fit approach will be used. The interpolation will, however, still be used to measure the accuracy of the fit.

Create the inflow interpolation:

Plot the inflow interpolation together with the inflow data:

On visual inspection, the interpolation nicely represents the data, but it is not as smooth as the fit.

The fit does not suffer from this problem. A fit based on trigonometric functions is used to obtain a periodic function. The accuracy of the fit can be increased by adding additional terms to the fit. This however negatively affects performance as more terms have to be evaluated. This application example uses 10 sine terms and 10 cosine terms which provides reasonable accuracy for this application, as will be shown later.

Create the inflow fit:

Plot the inflow fit together with the the inflow data:

The fit is less accurate than the interpolation, which is expected. It also oscillates quite a lot in the second half of the cardiac cycle. This is caused by the fact that trigonometric functions were used to obtain the fit. The trigonometric nature of the fit also means that the fit is much smoother that the interpolation, which will make the convergence of the model much easier and will positively affect performance.

The fit can be compared to the interpolation to measure the accuracy of the fit.

Maximize the difference between the fit and the interpolation to measure the accuracy of the fit:

Plot the difference between the fit and the interpolation:

The approximation error of the fit varies during the different parts of the cardiac cycle, but it is less than 0.02 most of the time. This is acceptable for this application example.