Blood flow modeling in cerebral aneurysm

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].

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 a context of the finite element method in the Wolfram language can be found here.

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 results the original mesh should be used.

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 corresponds to different types of elements - mesh element markers, boundary element markers and point element markers.

The red part of the boundary shown on the picture 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.

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.

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.

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.

As a last step the full element mesh is created.

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. In a second step, in a future version of the Wolfram Language, the pulsatile inflow case will be added. 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 exit two ways to go about this. One is to create a ParametricFunction through, for example, ParametricNDSolveValue and then iteratively increase the amount 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 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 make 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.

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

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.

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

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

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 which has prescribed two different values of pressure at both ends. The pressure difference creates pressure gradient which makes the fluid in the pipe flow.

The formula for the inflow profile is given as:

where r is a radius of the pipe (in our case of the inlet), is the center point of the inlet boundary. Note that this formula for the inlet velocity only works because have the no-slip condition on the wall. If a different boundary condition was 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.

Now the inlet profile could be computed. The vector (x-s_x,y-s_y,z-s_z), that 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 which is tangent to the inlet is obtained and 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 which can be used - the BoundaryUnitNormal function.

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 functions starts from 0, increases until it reaches 1 and is constant thereafter.

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

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

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

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.

Model set up and solution

THe first thing to do is to define variables.

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

The equations can be generated using the FluidFlowPDEComponent.

Now that 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 reduces required memory of storing the returned InterpolatingFunction objects significantly.

The simulation will take some time to finish.

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.

Results

Now that the equations are solved it is possible to take a look at the solution and compute some additional parameters from it.

Use vector plots to visualize how the blood flows.

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 which can not be observed in the stream plot.

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

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

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

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

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

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 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.

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.

Secondly, 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 practise 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.

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

The wall shear stress plays an important role in predicting where the aneurysm rupture is most likely to occur. The prediction of the rupture site will be done in the pulsatile case, which will be added in a future version of the Wolfram Language.

References

1.  Brunátová (formerly Trdlicová), Jana. (2021). Blood Flow Modeling in Cerebral Aneurysm. Master's thesis, Charles University, Faculty of Mathematics and Physics, online available here.