Simulation of Blood Flow Within the Brain

SIMULATION OF BLOOD FLOW WITHIN THE BRAIN

Sharif Kurdi

This report is produced under the supervision of BIOE310 instructor Prof. Linninger.

Abstract

Cerebrovascular diseases are diseases including strokes that constrict the blood flow within the brain causing major complications for a person [3]. This stoppage of blood flow can happen in any artery of the brain, as well as the Circle of Willis. The goal of this project is to estimate blood flow through a specific network of blood vessels in the brain. This project uses three main equations: the Hagen Poiseuille equation, Conservation of Flow equation, and Species Transport Equation. The first two equations mentioned are used to solve for the flows and pressure of a network while the Species Transport Equation solves for the concentration values. Different boundary conditions (pressure inlet and pressure outlet values) are chosen and the differences between concentrations over time graph curves, as well as the difference in residuals, are compared and analyzed. Overall this paper proposes a model that can help estimate blood flow within the brain in hopes that it will be useful for doctors and potential cerebrovascular diseased patients in the near future.

Keywords: Circle of Willis, Dye Concentration, Blood Flow, Cerebrovascular Disease.

1. Introduction

Cerebral vascular disease is known to be the fourth leading cause of death in the United States; according to the CDC there are approximately 128,932 deaths from cerebrovascular diseases a year in the U.S. [1]. A few types of cerebral vascular disease include strokes, vertebral stenosis, aneurysms, and vascular malformation to name a few [3]. These diseases attack and affect the circulation of blood flow within the brain [3] as well as arteries within the Circle of Willis, which can be seen in Figure 1.

The Circle of Willis’ main function is to help supply blood to the brain and is made up of the anterior cerebral artery (ACA), anterior communicating artery, internal carotid artery (ICA), posterior cerebral artery (PCA), and posterior communicating artery; the basilar artery and middle cerebral artery are also closely related [2]. Without going into too much detail, the posterior cerebral arteries (PCA) streams blood to the back of the brain, while the anterior and middle cerebral artery supply blood to parts of the frontal and parietal lobes of the brain [2].

The goal of this project is attempting to estimate blood flow through a network within the brain. A current technique that is used for this matter is digital subtraction angiography (DSA), though this method is flawed. This method only gives a 2D image of the agent over time. In this proposed method a simulation model to predict blood flow will be created and compared to an estimation model to see if the flow values are correct. If this method can surpass the

current one at hand, doctors could be supplied with the successful models and more effectively and clearly estimate blood flow within the brain whileultimately saving a number of lives in the process.

Figure 1 – The Circle of Willis is a group of arteries in the brain that help supply the brain with blood [1].

2. Methods

2.1 Solving for Flows and Pressures

The first step in creating this model is solving flows and pressures within a network; in this case the flows and pressures were solved first withina bifurcation network, and second in a more complicated Circle of Willis shaped network. In order to solve the flows and pressure the following equations were used:

• (1)

This Conservation of Flow equation is simply saying that the summation of the flow going into the network must be equal to the summation of flow going out of the network.

(2)

This Hagen-Pouiseuille equation is affirming that the change in pressure (∆P) is equaled to the resistance (α) multiplied by the flow (F).

It is necessary to note that while using these two equations to solve for the flows and pressures it is assumed that the resistance is equaled to one. The inlet and outlet pressure values (boundary conditions) are also known values. These equations were subsequently implemented into and solved using MATLAB.

2.2 Solving for Concentration Values

After solving for the flows and pressure, the concentration values through a network over a certain amount of time were solved for. To complete thisa third equation is needed:

(3)

In other words this equation (species-transport equation) is stating that volume (V) multiplied by the derivative of concentration over an amount of time () is equaled to the inflow (FIn) multiplied by concentration going into a network (CIn), subtracted by the outflow (FOut) multiplied by the concentration leaving the network (COut) plus injection rate (Inj). It is important to note that when carrying out this equation that volume and injection rate values are known.

To solve the concentrations using equation (3), Implicit Euler’s Method was applied in MATLAB.

Implicit Euler’s Method will take the differential equation and convert it into an algebraic function; subsequently concentration over time graphs will be plotted with the new algebraic equations.

2.3 Optimization

The next method completed was using different boundary conditions and plugging them into each networkwhich will give uniquepressure, flow, and concentration values. The residual values of each graph compared to the original one were then found and analyzed; this is more clearly explained in the results and the discussion sections.

2.4 Noise

The optimal boundary conditions concentration over time graph for the Circle of Willis Network was given different levels of noise using MATLAB commands, and the residual values of each unique set of boundary conditions were examined.

3. Results

It is critical to mention that in a normal scenario the boundary conditions (PIn and POut values) would initially be altered in an attempt to fit the model that is being dealt with. In other words the values and curves must be as similar as possible, resulting in a residual as close to zero as possible. In this project, however, the ideal boundary conditions are known beforehand as are shown in Figure 2A and 3A.Figure 2B-2F and Figure 3B-3F represent unique boundary conditions, and panels H illustrate the bifurcation network(Figure 2) and Circle of Willis network (Figure 3) worked with.

3.1 Bifurcation Network

In Figure 2G, the residualover outlet pressure graph, it displaysthat the low point of the residual seems to be at 5mmHg,and the residual continuously increases the further the outlet pressure is from 5mmHg. The residual values increase more rapidly when choosing outlet pressures larger than 5mmhHg than outlet pressures lower than it.

It can also be observed that a larger pressure outlet value results in more elongated curves while a lower pressure outlet values gives off a slightly more condensed curve. The inlet pressure works in reverse as the higher the inlet pressure value is, the more compressed the curves are; and the lower the inlet pressure is the more drawn-out each curve is.

Figure 2 - Panels A – F represent concentration over time graphs for a bifurcation network given unique boundary conditions; Panel A displays ideal boundary conditions. Panel G represents the relationship between the outlet pressure boundary condition and the residual value for the network; from the graph it’s clear that 5 is the optimal outlet pressure, as the residual equals zero. Panel G illustrates the bifurcation network that is being dealt with including flow direction, and concentration points.

Figure 3 - Panels A – F represent concentration over time graphs for a Circle of Willis network given unique boundary conditions; Panel A displays ideal boundary conditions. Panel G represents the relationship between the inlet pressure boundary condition and the residual value for the network; from the graph it’s clear that 10 is the optimal inlet pressure

The following table, represented as Table1, shows the relationship between the outlet pressure and the residual value within a bifurcation network and depicts the same data as the graph in Figure 2G. The outlet pressure value of 5mmHg is ideal, as the residual is 0 at that point.

POutValue (mmHg) / PIn Value (mmHg) / Residual
1 / 100 / 0.0016
5 / 100 / 0
9 / 100 / 0.0019
13 / 100 / 0.008
17 / 100 / 0.019
21 / 100 / 0.03

3.2 Circle of Willis Network

Figure 3A-3F depictdifferent concentration over time graphs based on their given boundary conditions; Figure 3A contains ideal boundary conditions. From Figure 3B and Figure 3C it seems that a lower inlet pressure value resulted in more elongated and curves while a higher inlet pressure value shown in Figure 3D-3F stemmed a more condensed or compressed curve.

In Figure 3G, the residual vs. inlet pressure graph, it displays that the low point of the residual seems to be at 10, and the residualcontinuously increases from that point whether the inlet pressure increases or decreases; this iswill be explained in the discussion section.

The following table, represented as Table 2, shows the relationship between the inlet pressure and the residual value within a Circle of Willis networkand depicts the same data as Figure 3G.

POutValue (mmHg) / PIn Value (mmHg) / Residual (mmHg)
1 / 2 / 1264.8
1 / 5 / 161.5
1 / 10 / 0
1 / 15 / 29.7
1 / 20 / 68.3
1 / 25 / 97.7

3.3 Circle of Willis Network with Noise

Figure 4 – Optimal Boundary Conditions with a Signal to Noise Ratio of 5dB.

To test the effect that adding noise to a concentration over time graph has on the residual for each graph, noise was added to the Circle of Willis optimal boundary conditions as expressed in Figure 3A. Figure 4 and Figures 6A, 6C, 6E, and 6F shows noise added with signal to noise ratios of 5dB, 10dB, 15dB, 20dB, and 25dB respectively.

Figure 5 – Residual over Inlet Pressure graph comparing each set of boundary conditions’ residuals with the optimal set of boundary conditions with added noise.

Figure 5 and Figures 6B, 6D, 6F, 6Hall show the relationship between the residual and the inlet pressure when a signal to noise ratio of 5dB, 10dB, 15dB, 20dB, and 25dB respectively is added. The graphs are all very identical with one another.

Figure 6 – Panels A,C,E,G represent concentration over time graphs with varied amounts of noise, while panels B,D,F,H represent the relationship between the residual and inlet pressure at each amount of noise.

4. Discussion

4.1 Bifurcation Network

As shown inFigure 2B the POutvalue was reduced to 1mmHg while the PIn value was unchanged at 100mmHg. It’s clear that there isn’t much of a difference between this graph and the optimal graph in Figure 2A. This is furthermore supported from the fact that the residual is the extremely low value of 0.0016. From this graph it can be deduced that merely lowering the POutvalue by 4, to the largest whole number, has a very insignificant effect on the concentration values, though it does very slightly condense the graph’s curves.

As represented in Figure 2C the next set of boundary conditions uses the POutvalues of 50mmHg, and PIn value of 100mmHg. From raising the POutvalue tenfold and leaving the PIn value the same, it is clear that the curves reach a higher concentration value than the original’s, and subside at a later time. This proves that a larger POutvalue will give a larger concentration value at each node, and will take the concentration a longer amount of time to leave the network. This makes sense, because a larger POutvalue will result in smaller flow values in bifurcation networks. Smaller flow values will subsequently keep the dye in each node for a longer amount of time. Theresidual value is at 0.5764for this graph which makes sense, as there visually is a higher amount of error from the ideal graph in panel A than there was in panel B.

In the case of Figure 2D the POutis kept at 5mmHg (from the original model), and instead the PIn value is halved to 50mmHg. These curves look very similar to the previous graph’s curves, because the concentration values are larger at each node, along with the time elapsed for the dye to leave the network being longer. This makes complete sense: since the PIn value is lesser than the ideal, it takes a longer time for the dye to enter the network resulting in a longer amount of time for the dye to leave the network. The larger elapsed time can also be explained by the smaller flow values that are obtained by the lower pressure inlet values. After analyzing these last few graphs it’s clear that a rise in the POutvalue OR a drop in the PIn value give off the same result: higher concentration values and more time being elapsed before the dye completely leaves the network. The residual of the graph is 0.8193 which signifies that we are getting farther from the solution compared to the last few.

The graph inFigure 2E was created using the POutvalue of 5mmHg (unchanged from original) and a PIn of 200mmHg (doubled from the original). Compared to the original this graph has more condensed curves; the curves not only have overall smaller concentration values than our original, but a smaller amount of time passed for all of the dye to leave the network. This is the effect drastically raising the PIn value gives which is logical: lowering the PInvalue lowers flow, but raising PIn raises the flow value; this results in the dye being pushed through the network more quickly. The residual value of this graph is 0.3068 which is closer to 0 than our last. This signifies that doubling our PIn gives is a more correct solution than halving our PIn value.

The last concentration over dye graph for this section, labeled as Figure 2F, had PIn and POutvalues 10x larger than our optimal boundary condition values. This graph has the lowest concentration and least amount of time elapsed thus far by a large amount as well as the greatest residual (1.0209). This graph helps to understand the relationship between boundary conditions well in the bifurcation network: the farther apart the boundary conditions are from one another, the smaller the concentration values and time elapsed will be for each graph (assuming POutis lower than PIn). The residual values will also be greater as the distance between the two initial conditions differ assuming PIn > 100mmHg (the optimal output value).

The first table describes the relationship between the outlet pressure value and the residual value; the ideal outlet pressure value appears to be 5mmHg as already known; this is signified by the residual equaling zero at that boundary condition. It is also clear that the farther away the inlet pressure is from the ideal inlet value, the larger the residual will be; in other words the farther away the boundary conditions are from each other the higher the percent error will be. The same data can be looked at visually within Figure 2G.

4.2 Circle of Willis Network

Concerning the Circle of Willis network the inlet pressures were the focus and subsequently were the only boundary condition changed. The main point to be taken out of these concentration over time graphs for the Circle of Willis is simple: the larger the inlet pressure boundary condition becomes, the higher the flow will be. On the other side of things the lower the inlet pressure value becomes, the lower the flow will become. Since flow has a major effect on the time it take for all of the dye to completely leave the network, it is able to be said that the more flow a network has, the less time it takes to flow out of the network; and the smaller the flow values are the longer it takes to leave the network. This can be proven by Figures 3A-3F. It can also be deduced from these graphs that merely changing inlet pressure (at the values used) does not affect the value of the peak concentration. The peak point looks to be the same for each graph at each individual node.

Figure 3G and table 2 show the relationship between the inlet pressure boundary condition and the residual. From looking at either demonstration it is confirmed that an inlet pressure of 10mmHg is an ideal boundary condition as that is when the residual reaches 0. It is also shown by the graph and table that the farther the inlet pressure is from the perfect pressure, the larger the residual will become. Therefore it can be deduced that the farther away the inlet pressure is from 10mmHg, the higher the error.

4.3 Circle of Willis Network with Noise

After adding noise to the optimal boundary conditions of Circle of Willis, one thing is clear: the more noise you add to a signal, the worse the residual value is all around. In other words the lower the signal to noise ratio is, the higher the residual values will be. Noise is generally unwanted data which worsens and interferes with the results [4] so it is logical that the more noise present would result in a larger residual for every single set of boundary conditions tested.

5. Conclusion

To conclude the goal was to simulate blood flow within the brain, specifically the Circle of Willis to help supply doctors with a more advanced and efficient way to predict blood flow. The Conservation of Flow Equation, Hagen Poiseuille Equation, and the Species transport equation were used to find the Flows, Pressures, and Concentration values respectively within a bifurcation network as well as a Circle of Willis network. Pressure boundary conditions then must be continuously altered in an attempt to find the most accurate set of boundary conditions; in other words to find the boundary conditions that give the closest residual value to zero. The residual values were found and analyzed in conditions with and without noise. It was deduced that the more noise present, the higher the residual value would be.