Protection of tissue perfusion during micro-infarction by collateral vessels : A modeling study

University of Amsterdam

Academic Medical Center, Biomedical Engineering and

Physics

Report Bachelor Project Physics & Astronomy, size 15 EC

conducted between

04-04-2016 and 07-07-2016

Protection of tissue perfusion during

micro-infarction by collateral vessels

A modeling study

Bobby A. Runderkamp

Student number: 10540083

Supervised by:

Prof. Dr. E.T. van Bavel

Second assessor:

Dr. H.A. Marquering

Abstract

Perfusion of tissue can decline because of occlusion by blood clots in the vascular bed. Vessel collateralization protects perfusion by offering blood ‘sneak routes’ as additional pathways to reach the capillaries where the oxygen exchange with the tissue cells takes place. The degree of collateralization differs between individuals. This paper describes a model built in MATLAB calculating local flows and pressures in two blood flow configurations: The honeycomb and the bifurcating tree. The honeycomb describes a general flow model with high level of heterogeneity containing many collaterals, whereas the bifurcating tree model contains none. The inflow, vessel pressures and capillary flow rates were calculated in different stages of vascular bed occlusion. The honeycomb shows a significantly higher ability than the bifurcating tree to protect inflow and capillary flow rates, as bigger vessel branches are occluded. Also, in the honeycomb model vessel pressures are stabilized increasingly better for bigger branch occlusion, whereas for the bifurcating tree the opposite happens. This indicates the perfusion protecting nature of collaterals in the honeycomb model.

Populair wetenschappelijke samenvatting op 6-VWO niveau

Het weefsel in het lichaam van mens en dier heeft een vaste zuurstoftoevoer nodig om goed te functioneren. Deze zuurstoftoevoer wordt verzorgd door het bloed. De zuurstofuitwisseling vindt plaats in de allerkleinste bloedvaten in het hart- en vaatstelsel: de capillairen, ook wel haarvaten. Deze kunnen tot wel 2.5 µm klein zijn. De capillairen bevinden zich in een vatenbed, waarin arteriolen met zuurstofrijk bloed zich veelvuldig vertakken tot de kleine capillairen, die zich daarna weer op dezelfde manier samenvoegen tot grotere venulen die het zuurstofarme bloed wegvoeren. De precieze structuur van dergelijke vatenbedden verschilt per individu. Dit betekent dat sommige mensen meer wegen hebben voor het bloed om de capillairen te bereiken dan anderen. Mensen met minder vaten als ‘vluchtroutes’ voor het bloed zijn kwetsbaarder voor bloedproppen. Bij deze mensen valt dan een groter gedeelte van het capillairnetwerk uit, waardoor weefsel te weinig zuurstof krijgt en afsterft: een infarct. Tijdens dit onderzoek is er een computermodel gemaakt die de bloedstroom en -druk berekent in elk vat binnen vatenbedden van verschillende vorm, volgens natuurkundige wetten. Het vatenbed kan namelijk worden gezien als een groot netwerk van weerstanden. Er zijn twee vormen vatenbedden gebouwd en bekeken: het honingraatmodel, die veel vluchtroutes voor bloed bevat, en het boomvertakkingsmodel, met geen enkele vluchtroute. Door bloedproppen te simuleren die op verschillende plaatsen in het vatenbed vaten blokkeren, kunnen de gevolgen voor de dieperliggende vaten en capillairen worden bekeken. Als je die gevolgen voor het honingraatmodel vergelijkt met die voor het boomvertakkingsmodel, lijkt het erop dat capillairen in het honingraatmodel over het algemeen minder snel met een zuurstoftekort komen te zitten dan die in het boomvertakkingsmodel. Dit komt waarschijnlijk door het grotere aantal vluchtroutes voor het bloed in het honingraatmodel. Het model is zo ontworpen dat andere soorten berekeningen ook uitgevoerd kunnen worden. Op deze manier kan nauwkeuriger bestudeerd worden hoe vluchtroutes de bloedtoevoer naar de capillairen beschermen.

Contents

I Introduction 4

II Theory 4

III Methods: the model 8

I Honeycomb . . . 9

II Bifurcating tree . . . 14

III Blood clots . . . 16

IV Results 18

V Discussion 24

VI Conclusion 25

I.

Introduction

The human body, and that of almost all ani-mals, has its own mechanism for delivering oxygen to all necessary organs and compart-ments. This is formed by the well known circu-latory system, where oxygen-rich (originating from the lungs) blood is pumped by the heart into the arteries, which divide into smaller arte-rioles and eventually the smallest vessels, cap-illaries. It is here where the oxygen exchange takes place, after which the now oxygen-poor blood travels back to the heart and lungs through the venules and veins. At this smallest level, the vessels divide many times, creating a so called ’capillary bed’, providing better cir-cumstances for oxygen exchange.

Figure 1: Illustration of a branching artery into an

arte-riole and a subsequent capillary bed. Oxygen exchange takes place between these capillaries and tissue cells, followed by anastomosis into a venule and vein [1]

.

This process, of delivering blood to a capillary bed, is called perfusion. However, the configu-ration of these capillary beds differs between individuals. This could make some individu-als more vulnerable than others to blood clots causing decline of organ function, if they have less ’pathways’ from an artery to subsequent

capillary beds. Vessel collateralization would be beneficial for these cases. Collateralization implies the existence of certain ’sneak routes’ (collaterals) through which the tissue can be oxygeniated another way around. The goal of the project presented in this paper is to bet-ter understand how, and to what extent, sneak routes can protect perfusion by giving blood the availability to reach its destination through multiple ways. This is done by building a blood flow (MATLAB) model with a topol-ogy as general as possible for our purpose, namely a honeycomb configuration, contain-ing many collaterals [2]. Three properties are measured, and compared to those of another self-built ’bifurcating tree’ model containing no collaterals. These properties are the volumet-ric rate of inflow, the capillary flow rates and the vessel pressures, and how these develop after the arrival of certain amounts of blood clots. Blood clots reduce the amount of blood flow pathways, so in a non-self-adapting model (constant vessel radii) this means an immedi-ate decline of the inflow and the capillary flow rate. Expected is that the inflow, capillary flow rate and vessel pressures are much less sen-sitive to these blood clots for the honeycomb model than for the tree model, because of the existence of sneak routes [3, 4].

II.

Theory

The arterial system can be seen as a large net-work of resistances, their values and connectiv-ity determining local flows and pressures. For a given vessel segment, this is described by the Hagen-Poiseuille equation [3]: Q= πr 4_∆P 8ηL = πd4∆P 128ηL (1) where G= Q ∆P = πd4 128ηL (2)

can be seen as the fluid conductance (m3(Pa·s)−1) of the vessel segment, by the analogy with electric circuits (Ohm’s Law) [5]. The assumption is made here that blood flow is laminar. This is reasonable because the Reynolds number is low. Blood also has non-Newtonian properties, but that is ignored here [3]. This equation expresses the volumet-ric flow rate Q (m3s−1) in terms of a number of properties of the particular vessel segment, being (SI-units) [6]:

d, the diameter (m) L, the length (m)

η, the dynamic (fluid) viscosity (Pa·s)

∆P, the pressure difference between the two ends (x=0 to x=L in figure 2) (Pa)

Figure 2: Properties of vessel segments as introduced in

eq. 1,2, where the mother vessel (0) bifurcates into 2 daughter vessels (1,2).

Figure 2 visualizes these properties for a basic vessel bifurcation.

Flow rate calculation in combined vessel sys-tem [7]

As the model serves to resemble large capil-lary beds as found in nature, it will contain thousands of separate vessel segments. This section explains how the flow rate and pres-sures at certain points are calculated, based on the ‘Wheatstone bridge’ configuration in figure 3. This is the simplest non-trivial blood-flow system topology with at least one vessel bifurcation and one collateral.

Figure 3: Blood flow system with the Wheatstone bridge

configuration. Blood flows from the source indicated red, through two dividing and then joining vessels, to the sink indicated blue. In between, there is a collateral through which blood can flow as well, the surrounding flow rates determining its direction. Points where segments meet (nodes) are numbered 1 to 4.

The figure shows the path of blood from a pres-sure source (red), to the end of the route which is called a pressure sink (blue). The points where vessel segments meet are called nodes and are numbered (details in section III). In analogy with electric circuits, the calculation of local flows in such a system follows from Kirchhoff’s first law for a particular node [3]:

n

∑

k=0

Qk=0 (3)

Here, n is the total number of branches with currents flowing towards or away from the node, and k loops over these branches. Hence, it says that the flow entering a node equals the flow going out. Denoting Xabfor the value of

source, this means for node 1 in figure 3 (in combination with eq.2):

Q1S+Q12+Q13+Q14 =

G1S∆P1S+G12∆P12+G13∆P13+G14∆P14 = G1S(PS−P1)+G12(P2−P1)+G13(P3−P1)+G14(P4−P1)=

−(G12+G13+G14)P1+G12P2+G13P3+G14P4+G1S(PS−P1)=0

Then G14vanishes because 1 is not connected

to 4. Extending this to all other nodes, the last line in the derivation can be recognized to be the first row of the matrix equation:

AP+B(PS−P) =0 (4)

This equation holds for any flow system. For a system with n nodes, A is the (n x n)-matrix defining internal connectivity and contains all internal conductances. B is the (n x n) diagonal matrix defining connectivity to sources. Often there are only sources connected to the first and last node, causing B to be a null matrix with only the first and last element nonzero. PS is the n-dimensional column vector

defin-ing the input/output pressure sources/sinks to each node. In the human body, the highest pressure drop occurs at the transition of arte-rioles to capillaries. This is the blood system part we are modeling, and we will use an input pressure source of 100 mmHg (P1S) and

pres-sure sink of 5 mmHg (P4S) (x133Pa in SI-units).

For a general system with 4 nodes, the matrices are shown in figure 4. With all these known, P, the n-dimensional column vector defining the pressure at each node, can be calculated from 4 by some matrix operations:

AP+B(PS−P) =0

AP+BPS−BP=0

(A−B)P= −BPS

P= −(A−B)−1BPS (5)

Combining eq. 5 with 1 gives the desired flow rate Q at every vessel in the system.

Figure 4: The internal connectivity matrix A, source

connectivity matrix B and input/output pres-sure source/sink vector P, for a general flow system containing 4 nodes, from eq. 5. Here, Gabis the conductance between node or source a and b, and PXS is the pressure source/sink value for node X.

Diameter calculation at vessel bifurcation

In our model we are dealing with vessel bifur-cations as illustrated in figure 2. Here, we will call the original vessel the mother vessel, and the branches it bifurcates into are called daugh-ter vessels. It is observed that vessels roughly bifurcate according to the following allometric relation, relating the radii of n daughter vessels to the radius of the mother vessel:

rγ₀ =

n

∑

k=1

rγ_k

Here, rk is the radius of daughter vessel k,

and r0 that of the mother vessel. γ typically

ranges from 2 to 3, in our model it will be taken 3, with which the relation is known as Murray’s Law [3, 8, 9]. As we will only be deal-ing with bifurcations, we get:

rγ 0 =r γ 1+r γ 2, or, in diameters, dγ 0 =d γ 1+d γ 2 (6)

Then, when d0is known, we get one degree of

freedom between d1and d2. To solve this, we

proportionality x. This gives for d1using eq. 6: d2=xd1 dγ 0 =d γ 1+d γ 2  dγ₀ =dγ₁+ (xd1)γ= dγ 0 = (1+xγ)d γ 1 ⇒ d1= dγ 0 1+xγ !_γ1 (7)

The smallest vessels in our model are the capil-laries, so for a capillary diameter p we impose the constraints

d1> p (8)

d2=xd1> p (9)

From here the allowed range of x can be calcu-lated. Writing d1as f(x), we get a first constraint

on x with inequality 8: f(x) >p⇒ dγ₀ 1+xγ !_γ1 >p⇒ 1+xγ_< d γ 0 pγ ⇒ xγ_< d0 p γ −1⇒ x< γ s  d0 p γ −1=b1 (10)

Equation 10 gives the upper limit of x, which we named b1. This means f(b1) = p and we can

get another constraint on x from eq. 9: xd1=x f(x) >p= f(b1) ⇒ x d γ 0 1+xγ !_γ1 > d γ 0 1+bγ 1 !1_γ ⇒ xd0(1+xγ)− 1 γ >d₀(1+bγ 1) −1 γ ⇒ x−γ(1+xγ_{) <1+b}γ 1 ⇒ x−γ<bγ 1 ⇒ x>b−1₁ = 1 b1 =b2 (11) So we get for x: b1= γ s  d0 p γ −1>x>b2=b−1₁ =  d0 p γ −1 −1_γ = 1 dγ pγ −1 !1_γ =  pγ dγ_−pγ _γ1 (12)

For the daughter vessel diameter calculation, x will be chosen random. Because log(b1)

= -log(b−1₁ ) = -log(b2), log(x) will be

cho-sen from the continuous uniform distribution [log(b₁−1),log(b1)]. Then x is achieved by

tak-ing the inverse logarithm. With x known, the daughter vessel diameters are calculated using eq. 7,9.

In the case where d1,d2are known, d0is

calcu-lated by modifying eq. 6 to:

d0= dγ1+d

γ

2

1_γ

(13) The viscosity of blood ranges from 3-4·10−3 Pa·s at T=37°C, depending on its velocity and the vessel radius [3]. Here, it will be taken 3·10−3_{PA·s, constant for blood in all vessels.}

Then the remaining variable to be determined for conductance calculation (eq.2) is L. For this we will use the following relation between L and d [10]:

L=10.2d0.72 (14) This relation was established for porcine ves-sels, and the assumption is that porcine and human vessels behave alike. With conductance G known, P and consequently Q at all vessels can be calculated by the procedure described in the previous paragraph. The pressure P in a vessel is said to be the average of the pressures at its nodes.

III.

Methods: the model

The model was built by programming in MAT-LAB. The purpose was to reconstruct a capil-lary bed as it is found in the body. To get a model resembling the one in figure 1, a capil-lary bed with a honeycomb-shaped configura-tion was chosen. A honeycomb configuraconfigura-tion is constructed by repeating a certain ‘basis ele-ment’. The topology of this element is shown in figure 6.

Figure 6: Blood flow system basis element, multiple of

which the honeycomb configuration is com-posed. Nodes are indicated by dots and are named IN(x), vessels are straight lies named IE(x). Vessels can be connected to the basis el-ement at any node, from/to the node of another basis element or a source/sink. The illustration shows the structure’s topology, vessel lengths and widths are not to scale.

The basis element is a hexagon starting with bifurcation and ending with anastomosis. Con-necting several basis elements amounts to cre-ating bifurcations at nodes. At a larger scale, this gives the availability of a relatively large amount of collaterals, making it a rather gen-eral structure, from which less complicated structures can be extracted, as is the case in our study (subsection II).

Terminology

The terminology in our program code will be explained by means of figure 5. The figure shows a structure where four basis elements are connected symmetrically. All basis ele-ments ‘stacked’ on top of each other vertically, are said to be in the same column. We will call the level of such a structure after the amount of stacked basis elements in the middle column. Hence, the figure shows a honeycomb of level 2. Also drawn are blue dashed lines passing through nodes with equal horizontal distance in the model. These will be called planes and are a measure for how far into the capillary bed a particular node is. IN, IE and SE stand for Internal Node, Internal Element and Source El-ement respectively, meaning node, vessel and source/sink. The desired data lies in these ele-ments, and is stored in structs with the same

Figure 5: Blood flow model for a honeycomb of level 2. Indicated are source elements SE, internal nodes IN and

name. Figure 7 shows the fields of these structs for the level 2 case from figure 5.

Figure 7: MATLAB output: the structs IN, IE and SE

with its fields. All data of the nodes, vessels and sources/sinks are stored here. For an IN, nconnect is the amount of connected nodes, nsources the amount of connected sources and cn, ie and se contain the indices of the con-nected nodes, vessels and sources/sinks respec-tively. nodes in IE gives the nodes IE con-nects. All other parameters are as explained in section II.

I.

Honeycomb

Honeycomb model properties

Before writing code, a number of properties of the theoretical model have been calculated, which have been used to build our model. We have chosen to construct our honeycomb by

symmetrically stacking basis elements, b.e.’s, in the way shown in figure 5. This makes the number of stacked b.e.’s one lower left/right of the middle. Then, for a honeycomb of level x, there are x b.e.’s in the middle column, (x-1)+(x-2)+...+(x-(x-2))+(x-(x-1)) b.e.’s left of the middle and an equal amount on the right, so

Nbe,x=x+2 x−1

∑

i=1 i=x+2(x−1)x 2 = x+x2−x=x2 (15) where Nbe,x is the number of b.e.’s in level x.

For the amount of columns, Nr,x, we can

im-mediately say

Nr,x=2x−1 (16)

Regarding planes, the middle column contains four, where every column to the left or right delivers two. This gives for the number of planes:

middle column ↓

Np,x=4+2(x−1) ·2=4x (17)

% ↑

#columns left/right left and right

For the calculation of the amount of nodes, NI N,x, we see that the middle two planes

con-tain x+1 nodes, then the next two planes to the left or right contain x nodes, which repeats itself all the way to 2 planes containing 2 nodes, then followed by one having one node. This gives the sequence:

x+1 x+1 x x x x x−1 x−1 x−1 x−1 x−2 x−2 x−2 x−2 ... ... ... 3 3 3 3 2 2 2 2 1 1

Summing this gives NI N,x: NI N,x =2(x+1) +1+1+4 x

∑

i=2 i= 2x+4+4 _x

∑

i=1 i ! −1 ! = 2x+4+4 x(x+1) 2 −1  = 2x+4+2x(x+1) −4= 2x2+4x= 2x(x+2) (18) Then NIE,x remains to be calculated. Again

looking at figure 5 and recalling our stacking method, we see, working from the middle to the left, there are x+1 IE’s between the middle planes (4,5 in figure), 2x between the next two (3,4), then x, 2x-2, x-1, all the way to 4, 2 and 2 (in general). This gives the following sequence for the left part:

b→ 2x a→ x b→ 2(x−1) a→ x−1 ... b→2(x− (x−2)) a→ x− (x−2) c→2(x− (x−1))

Taking each adjacent couple of b,a together gives: 3x 3(x−1) 3(x−2) ... 3(x− (x−2)) c→2(x− (x−1)) =2

As there are x-1 terms with coefficient 3, we get: NIE le f t,x=3· x−2

∑

i=0 (x−i) +c= 3·1 2(x−1)(x+2) +2

By symmetry NIE le f t,x= NIE right,x, so we

con-clude, adding the middle planes’ IE’s: NIE,x=2  3 2(x−1)(x+2)+2  +x+1= 3(x2+x−2) +4+x+1= 3x2+4x−1 (19) Another useful property to know is the location of the capillaries. The capillaries are the vessels with the smallest radius, so we are looking for the ultimate daughter vessels. For the left part, these are found between the left middle plane and the one left of it, and for the right part like-wise. For the level 2 case in figure 5, these are the IE’s between plane 3 and 4, and 5 and 6 re-spectively. In the figure, one can see that IE(10) is a mother vessel, so the IE’s between the mid-dle planes (4,5) are no capillaries, except the top and bottom ones which do not bifurcate at plane 4,5. We can also illustrate this as a row of circles representing IE 1 to 19 from left to right:

o o o oo o o ooooo o o oo o o o Here, yellow circles are capillaries, red are IE’s in the middle, and blue are both. Black circles are normal IE’s. Now, for a honeycomb of odd level x, NIE,xfrom eq. 19is even, by noting that:

• x2is odd • 3x2_{is odd}

• 4x-1 is odd • 3x2+ 4x-1 is even

Also, the number of middle IE’s, x+1, is even. For x=5 for example, we get the row:

NIE,x 2   y ·· 2x z }| { o o o o o o o o o ooo o o oo | {z } x+1 2x z }| { o o o o o o o o o o·· x   NIE,x 2 - x+12

where yellow circles are followed by black ones. As can be seen, for even NIE,x, NIE,x₂

gives the index of the left middle circle, so the index of the left blue circle is given by NIE,x

2 -x+1

2 + 1. Consequently, the right blue circle has

index NIE,x

2 + x+12 .

For even x, NIE,xis odd, as is x+1:

• x2is even • 3x2is even • 4x-1 is odd • 3x2+ 4x-1 is odd

x=4 gives an example of an odd row:

NIE,x 2 + 12   y ··o o o o o o o ooo o ooo o o o o o o o·· x   NIE,x 2 + 12 -x+1 2 - 12 = NIE,x 2 -x+1 2

We see that for odd NIE,x, NIE,x2 equals the

mid-dle index minus 1₂. However, moving to the index of the left blue circle delivers another 1₂, canceling the previous one. Hence we get the same result for even and odd x. Adding the outer yellow circle indices, we have as capillary indices: Left: NIE,x 2 − x+1 2 +1−2x to NIE,x 2 − x+1 2 +1 Right: NIE,x 2 + x+1 2 to NIE,x 2 + x+1 2 +2x

Substituting NIE,x from eq. 19, we get after

some algebra our final expression for the capil-lary indices in a honeycomb of level x:

1 2x(3x−1) (20a) Lower left 3 2x(x+1) (20b) Upper left 1 2x(3x+5) (20c) Lower right 3 2x(x+3) (20d) Upper right These are the boundary indices, so all left capil-lary indices lie between and including eqs. 20a and 20b, and the right between and including eqs. 20c and 20d.

Building the honeycomb model

With many properties at hand, we were able to start implementing the model in MATLAB. There already was a program available calculat-ing the flows and pressures in the Wheatstone case in figure 3 [7]. This program consisted of three functions:

• wheatstone: This function defines the topology, radii, lengths and conductances of the internal and source elements, for the Wheatstone configuration. For the source elements, the source and sink pres-sures are also defined.

• MakeNodeTable: This function takes IE and SE from wheatstone and defines the topology of the internal nodes.

• solvehemodyn: This function takes IN, IE and SE from the previous functions, de-fines the matrices for local pressure cal-culation and performs it according to eq. 5. Then the local flows at every IN, IE and SE are calculated and assigned, com-pleting the program.

To our purpose, the latter two could be used di-rectly and were not changed. wheatstone had to be changed from the Wheatstone configuration to the honeycomb configuration. Our analog of wheatstone was spread over two functions: build_IE and honeycomb. build_IE defines the topology of the internal elements and honey-comb adds their radii, lengths and conductance. honeycomb also defines all SE properties, and includes MakeNodeTable, after which the model is fully defined for solvehemodyn to calculate Q and P.

build_IE loops over all planes and checks for plane p whether the amount of nodes in plane p+1 differs by 1, -1 or 0. For this, it uses the functions Calculate_nodes_plane and Calcnode-splaneindices. Calculate_nodes_plane calculates the amount of nodes per plane in our honey-comb model as seen in the derivation of eq. 18, and Calcnodesplaneindices gives the appropriate indices. As seen in figure 5, a node difference of 1 means that every node from plane p con-nects to two nodes in plane p+1, and a differ-ence of -1 means the opposite. A differdiffer-ence of 0 gives trivial one-to-one connections. This is equivalent to the base element stacking method as stated earlier. Hence the node-connection of IE is known and defined, and honeycomb can start.

honeycomb first assigns the radii of the cap-illary IE’s, as identified by the indices from eqs. 20a,b,c,d. The capillary radius was set to 2.5·10−6 m [3]. To create an actual capil-lary bed, instead of two columns of capillaries interrupted by the middle column of mother

vessels, these middle IE’s (indices between eq. 20b and 20c) were assigned a capillary radius as well. Then, the function loops over all planes from the middle to the left and right, and gives every IE its radius using eqs. 7,9,13. At bifurca-tions, function x_generator is used to generate the random proportionality factor x from the log of the interval in eq. 12, as explained in section II. One complication arose however in this case, for the first bifurcation. Consider the mother/daughter vessel system in figure 8 from our honeycomb model.

Figure 8: Schematic illustration of anastomosis and

subsequent bifurcation between daughter and mother vessels. 1 and 2 indicate capillaries, 3 is the mother vessel, and 4 and 5 are the next two daughter vessels. A indicates anastomosis, B the bifurcation.

In the figure, 1 and 2 indicate the capillaries which just had their radii assigned. At A they join to form mother vessel 3, which gets a radius of d3= d31+d32

1₃

, eq. 13 (γ = 3). Then at B a bifurcation takes place, hence d3

= d3₄+d3₅
1₃

. We want our capillaries to be the smallest vessels, so d4,d5>d1,d2. But this

gives d3= d34+d35
1₃ > d3 1+d32
1₃ = d3, which

is a contradiction. In other words, Murray’s Law must hold from both sides. Hence the only possibility for d4 and d5 is d4,d5=d1,d2.

But then every program run returns the same model, and we lose all heterogeneity. As vessel systems in nature generally assume a high level of heterogeneity, this is not desired. We solved this issue by allowing d3 to be 1 to 5

times larger than its Murray value. This en-largement factor is chosen random from the

continuous uniform distribution [1,5]. This modification creates a whole column of ran-domized mother vessel radii left and right from the middle, giving a high level of heterogeneity in our model, because now all other IE’s can be calculated with Murray’s Law with random proportionality factors x. The radii of SE also follow from Murray’s Law. Next, honeycomb assigns the lengths and conductances to the IE and SE using eqs. 14 and 2. The topology and source/sink pressures are also defined: SE(1) is connected to node 1 with a pressure of 100 mmHg, SE(2) is connected to node 2x(x+2) from eq. 18 with a pressure of 5 mmHg. In SI-units, these are 13300 and 665 Pa respectively. Now honeycomb is done and solvehemodyn can be called, completing the flow and pressure calculation in the honeycomb.

For visualization, the function draw_graph_currents was built. draw_graph_currents returns a bio-graph in which IN’s are drawn as blocks and IE’s as arrows, where the width of the arrow scales with the value of Q through the IE. Hence, a bigger arrow means a larger flow. IE’s with negative Q are drawn in red instead of backwards, because the latter distorts the honeycomb shape in biograph objects. Also, MATLAB draws it vertically, rather than hor-izontally as we have been doing throughout this paper. Figure 9 shows the biograph object for a honeycomb of level 5.

Figure 9: Biograph object visualizing the flow through

a honeycomb of level 5. IN’s are indicated as blocks, IE’s as arrows. Arrow widths scale with Q, and a negative Q is colored red.

II.

Bifurcating tree

An additional advantage of building such a general (that is, containing extremely many pathways for blood flow) model as the honey-comb is that other simpler models don’t have to be rebuilt separately, but can be ‘extracted’ from the honeycomb model. This can be done by removing IE’s from the honeycomb, until you are left with the IE’s constituting the de-sired model. Through this method our bifur-cating tree model was made, about which the details will now be covered.

Bifurcating tree model construction and properties

The bifurcating tree model is a configuration where all branches bifurcate at certain and the same planes, symmetrically. For example, branch 1 bifurcates in the beginning, plane 1, into branch 2 and 3, which then bifur-cate into branch 4,5 and 6,7 respectively at plane 1+p, etc. That way, at every bifurcat-ing plane the amount of vessels is doubled. While building the bifurcating tree from the honeycomb, we decided that all middle vessels had to be included. Hence, the amount of ul-timate daughter vessels in the middle is equal to 2number of bifurcating planes. This means that a symmetric bifurcating tree can only be made from a honeycomb which number of middle vessels is a power of 2. Because the number of middle vessels for a level x honeycomb equals x+1 we get for the possible bifurcating tree lev-els:

x=2n−1 (21a)

where n is the number of bifurcating planes which can be any positive integer. Conversely, for a tree of level x, we know the number of bifurcating planes:

n=log₂(x+1) (21b)

The first bifurcation of SE(1) into IE(1) and IE(2) is included in n.

We had to find the IE’s from the honeycomb taking part in the bifurcating tree. For that, we needed to know the planes where bifurcations take place such that we ended up including all middle IE’s. For that, a logical bifurcation location is somewhere in the middle of each branch. We consider the example of a honey-comb of level 15, having 4 bifurcating planes by eq. 21b. Writing out the first few cases by hand, we saw the following pattern. The first bifurcation takes place in the beginning, where SE(1) bifurcates into IE(1) and IE(2), cre-ating the two main branches on the edge of the honeycomb. Then the next plane index where bifurcations take place, should be the amount of planes each branch contains, divided by 2, plus the plane index of the previous bifurca-tion, plus 1 (index issue). In our case, the total amount of planes these edge branches contain is 2x=30, as these branches contain all planes from beginning to end. Here, we mean with ‘end’ the middle of the model, where the

ulti-mate daughter vessel reside. The plane index of the previous bifurcation was 1, the begin-ning. So for the second bifurcation plane we get 30₂ + 1 + 1 = 17. Then the same process is repeated for the newly formed branches. The amount of planes each new branch contains is the total amount of planes up to the mid-dle, minus the planes before the bifurcation, so 30-(17-1). The plane index of the previous bifurcation was 17, so for the next bifurcating plane index we get 30−(17−1)₂ + 17 + 1 = 25. We see that the previous bifurcating plane index was used (twice) in the calculation of the next. So the calculation should be done in a loop. We write Np_total,x for the total amount of planes

up to the middle for level x, and Iprev.bif.for the

plane index of the previous bifurcation. Then for level x, the first bifurcating plane index is defined before the loop as Np_total,x₂ + 2, and the

body of the loop is of the form: Iprev.bif.= Np_total,x− (Iprev.bif.−1) 2 +Iprev.bif.+1 =1 2(Np_total,x+Iprev.bif.+3) (22) The loop ends when we have found all bifurcat-ing plane indices accordbifurcat-ing to eq. 21b. Figure 10 shows the process for the case described, level 15. With the bifurcation locations known, all internal elements are as well and we have the bifurcating tree.

Figure 10: Illustration of the bifurcating plane indices

calculation for the level 15 case. Bifurcations are indicated by red circles. Numbers next to vertical arrows indicate the amount of in-cluded planes, where the first and last planes are included. Internal elements taking part in the bifurcating tree are drawn gray, others have zero flow and are drawn yellow.

III.

Blood clots

We wanted to study the reaction of the blood flow systems to the arrival of blood clots. Func-tion blood_clot first defines the radius of a blood clot, and lets it flow through the system. This is done by letting it start at SE(1), and checking whether the clot radius is lower than the vessel radius at every internal element is passes. If it is, the final destination of the clot has been found an the internal element conductance G, hence Q, is set to zero. If it isn’t, the same check is repeated for the next vessel in the clot’s path. This next vessel is determined by the direction in which Q is positive, as the clot flows with the blood. At bifurcations, the blood clot has two path choices. The chosen vessel is determined by chance, where the probability of passing from mother vessel a through daughter vessel b equals Qb

Qa. As the capillaries are the smallest

vessels, choosing a clot radius larger than the capillary radius secures the clot to get stuck somewhere, instead of passing the entire hon-eycomb, giving no information.

Functions execute_blood_clot_honeycomb and ex-ecute_blood_clots_bifurcating_tree use blood_clot, all previously mentioned functions and some minor functions to ultimately gather data of IE, SE and IN after every blood clot for the honeycomb and bifurcating tree respectively. draw_graph_currents draws the last blocked IE black, and all IE’s with zero flow caused by blockage are drawn yellow. The biograph of a level 5 honeycomb after six blood clots is shown in figure 11.

Figure 11: Biograph object visualizing the flow through

a honeycomb of level 5 after the arrival of six blood clots. Here, capillary blockage was chosen by setting the clot radius equal to the capillary radius plus 1·10-10. The most re-cently blocked vessel is drawn black, flowless vessels are yellow. The biograph is a good illustration of how the blood flow chooses its path through the not yet blocked vessels.

Data gathering

To compare the models, three different func-tions were built for both the honeycomb and the bifurcating tree.

In the first, a blood clot with radius equal to the capillary radius plus 1·10−10was chosen. This way, each blood clot gets stuck at a capillary, as these are the smallest vessels. For the second and third functions we decided to abandon the blood clot simulation function blood_clot. Instead, in the second function, for the bifur-cating tree, G of random mother vessels of the last bifurcation is set to zero, until the model is fully occluded. That is, the mother vessels bifurcating at the last bifurcating plane. In the third function, random blockage of mother vessels of one bifurcation earlier was chosen. For the honeycomb, for comparison, the IE’s of the same columns are randomly blocked, of which there are more. Figure 11 shows the pro-cess for capillary blockage for the honeycomb, and figure 12 shows the processes of earlier blockage, for the level 7 case of the tree model. In figure 12 the perfusion declination in the vessels of blocked branches is clearly visible. Some of these still have a negligible flow of the order 10−23due to MATLAB calculation issues, which are not drawn yellow because their flow is nonzero.

For our results, the response to occlusion of three different model properties was measured: the inflow, the pressure and the capillary flow rates. Per property, the same level was cho-sen for each model, for good comparison. For the inflow (SE(1).Q) measurement, level 31 was chosen, and 20 runs were done and averaged over, where one run is from 0% to 100% oc-clusion. Plotting was done as a line graph. For the pressure and capillary flow rate mea-surements, level 63 was chosen, and one run was examined. The data after 0% and 70%, and 10% and 70% respectively were plotted. As there is a discrete amount of blood clots, the data after the amount of clots nearest to 10% and 70% occlusion was used. The

pres-sure meapres-surements were plotted as a scatter plot, and the capillary flow rates in histograms. For the histograms, an amount of bins closest to psamples was used, where samples is the amount of capillaries [3]. For the inflow and pressure measurements, capillary, last bifurca-tion and second last bifurcabifurca-tion occlusion was used, for capillary flow rates only the latter two.

Figure 12: Biograph object visualizing the flow through

a bifurcating tree model of level 7 after ran-dom blockage of bifurcating mother vessels. In the left model, G of a mother vessel of the last bifurcation was set to zero, and in the right model the same was done for a mother vessel one bifurcation earlier. These are indi-cated by a black arrow.

IV.

Results

As explained in subsection III of section III, first the inflow declination was measured caused by occlusion at the capillaries, last bifurcation and second last bifurcation. The resulting graphs are shown in figure 13, 14 and 15 respectively. The inflow is expressed as the fraction of the initial flow, which is where no occlusion has occurred yet and the inflow is maximal. The occlusion is expressed as the percentage of the total amount of blood clots.

Figure 13: Dependence of the volumetric rate of inflow on the occlusion by blood clots at the capillaries.

Figure 14: Dependence of the volumetric rate of inflow on the occlusion by blood clots at the mother vessels of the last

Figure 15: Dependence of the volumetric rate of inflow on the occlusion by blood clots at the mother vessels of the

second last bifurcations.

It can be seen that for capillary and last bifurcation occlusion the honeycomb inflow is lower everywhere than the tree inflow, and for second last bifurcation occlusion the honeycomb inflow is higher. The honeycomb graph moves up, whereas the tree graph moves down.

Secondly, for every vessel in the system, its pressure was plotted against its radius.

Figure 17: Pressure to radius distribution for vessels in the honeycomb, after 0% and 70% occlusion of mother vessels

of the last bifurcations.

Figure 18: Pressure to radius distribution for vessels in the honeycomb, after 0% and 70% occlusion of mother vessels

of the second last bifurcations.

Figure 20: Pressure to radius distribution for vessels in the bifurcating tree, after 0% and 70% occlusion of mother

vessels of the last bifurcations.

Figure 21: Pressure to radius distribution for vessels in the bifurcating tree, after 0% and 70% occlusion of mother

vessels of the second last bifurcations.

Comparing figures 16, 17 and 18, it can be seen that the data start overlapping more and more as earlier bifurcation takes place. Comparing figures 19, 20 and 21, it can be seen that the opposite happens. The data diverges more and more.

Finally, the flow rates of the capillaries were plotted in a histogram, after 10% and 70% occlusion of the mother vessels of the last bifurcation and second last bifurcation respectively.

Figure 22: Dependence of the honeycomb capillary flow rates on the occlusion of 10% and 70% of the mother vessels

of the last bifurcation by blood clots.

Figure 23: Dependence of the honeycomb capillary flow rates on the occlusion of 10% and 70% of the mother vessels

Figure 24: Dependence of the bifurcating tree capillary flow rates on the occlusion of 10% and 70% of the mother

vessels of the last bifurcation by blood clots

Figure 25: Dependence of the bifurcating tree capillary flow rates on the occlusion of 10% and 70% of the mother

Again, comparing figures 22 and 23 it can be seen that with earlier occlusion, the data starts overlapping. The reverse process for the bifurcating tree, as observed at the pressure-radius distributions above, is not visible between figures 24 and 25.

V.

Discussion

For the inflow dependence on vessel occlusion, it can be seen that for earlier occlusion (figure 15), the honeycomb inflow gets higher than the bifurcating tree inflow. This means that the honeycomb configuration has a higher ability than the bifurcating tree to maintain the inflow as earlier occlusion takes place. A higher in-flow means a higher total in-flow, so more in-flow through the capillaries. For the pressure scatter plots, in figure 16 it can be seen that the pres-sure values of all vessels, especially small ones, diverge when capillaries are occluded. Figures 17 and 18 however, show a gradually higher ability of the honeycomb to keep the pressures stable as earlier bifurcation occurs. This shows the power of the sneak routes in the honey-comb. For capillary occlusion, the capillaries itself are getting blocked, so they cannot be perfused. However, when the mother vessels of the last bifurcation are getting blocked, there is still a pathway for neighboring capillaries to perfuse the otherwise dead capillaries. For oc-clusion of one bifurcation earlier, even more of these pathways open up. This is confirmed by these figures. For the bifurcating tree, figures 19, 20 and 21, the opposite happens. There, for each vessel radius the pressures after 70% occlusion diverge more and more. They don’t have collaterals to protect perfusion, and ear-lier occlusion shows more trouble stabilizing the original situation. Between figures 22 and 23 the ability to protect perfusion is observed as well. For occlusion before the second last bifurcation we see a minor change of capillary flow rates after 70% occlusion. The histogram of occlusion before the last bifurcation however shows a much more radical response to 70% occlusion, where many vessels already have Q = 0.

There are several ways to improve and con-tinue this study. It looks like the program pro-duces good, trustworthy calculations for local flows and pressures. However, it doesn’t work as fast as we would want it to. In nature, capil-lary beds are of a size in the order of 104 capil-laries [3]. This is in the order of 100 times larger than our 156 capillaries from the level 32 hon-eycomb we used. That would require waiting several hours until a valuable calculation has been done. The code, containing 15 functions for a tree blood clot calculation alone, should be checked to make it more efficient. However, ≈ 90% of the runtime comes from solvehemo-dyn, especially because of eq. 5. This is an exact matrix calculation which doesn’t show much room for improvement as the matrices are already made sparse. The only possible improvement here seems to be to abandon the exact case and look for approximations. The runtime is especially high for the bifurcating tree model. This is probably because it is built up from the honeycomb model, hence having the same dimensions but a majority of zeros. Therefore, while performing eq. 5 we get the warning: ‘Matrix is singular to working preci-sion’. This increases the runtime drastically for models of high level. For a follow-up study we could consider building a separate model for the bifurcating tree model, solving the prob-lem of big matrices with many zeros. This decreases the runtime, but makes it harder to compare it with the honeycomb.

For the blockage of earlier vessels than cap-illaries, we did not simulate a blood clot, but randomly chose IE’s from a given column. This doesn’t give very reliable results compared to blood clots, as blood clots have higher chance ending up in a vessel with high Q, whereas for the random IE choices the probability dis-tribution is (discrete) uniform. Properties of earlier vessels in the blood clot path aren’t in-cluded in the random case. An improvement would be to change the uniform probability

distribution to one where the chance of an IE getting blocked is proportional to its Q value. For example, prob(IE(n)) = Qn

Q_total where Qtotal

is the sum of the Q’s of all IE’s in the particular column we’re blocking.

This method of choosing random IE’s was cho-sen for the following reason. To cause a blood clot to get stuck at mother vessels of an earlier bifurcation, the clot radius has to be increased. But then there is a high chance that blood clots encounter vessels with lower radius in their path than those of the vessels at our desired destination, causing premature blockage. Then a good comparison of vessel column blockage can not be made. This method still however gives some information regarding the reaction of a flow system to the occlusion of branches. This study set the stage for calculating many model properties. During this study, the whole working program was built and made as user-friendly as possible, but this occupied almost the whole time period for this research project, so a rather low amount of results were gath-ered. This could be done in a follow-up study. For example, occlusion of earlier vessels could be performed, possibly for honeycombs of higher level. Also, a graph showing the per-centage of capillaries having too little flow to function well, caused by occlusion, against the amount of blood clots would be insightful. Ran-dom blood clot radii can be chosen, giving in-formation of vessels having the highest chance of getting occluded. Comparing the models to a third model with an amount of collaterals between that of the honeycomb (many) and the bifurcating tree (none) would be insightful as well. This shows the power of the honeycomb model: Any configuration can be extracted from the general honeycomb, giving good com-parison possibilities, at the cost of high running times.

Another limitation of the model is the one de-scribed below figure 8. The first column of mother vessels got radii above their Murray value. This doesn’t make the model drastically unrealistic, because in and around capillary beds in nature there are many cases where Murray’s Law doesn’t hold [3]. However, this

is a problem arising from the symmetry of the honeycomb, so it might get solved by breaking some symmetry, at the cost of losing a good overview, at least in the biograph drawings. For the pressure and capillary flow rate mea-surements, one run was used. Averaging of the data here would improve the results, as has been done by the inflow calculation. For the capillary flow rates, this could be done by an average of histograms. The improvement is definitely necessary for the capillary flow rate histograms for the bifurcating tree. Quantita-tive data, in numbers, of how well collaterals protect perfusion may be insightful as well. This paper merely presents qualitative data in the form of graphs.

An additional necessary improvement is to make the vessels adaptive. In nature, small arterioles undergo vasoconstriction and vasodi-lation, which changes the blood pressure and radius. Adding this ability to our model would make it more trustworthy.

One last comment should be made about ves-sels in an occluded branch having Q in the or-der 10−23instead of zero, because of some side effects inside MATLAB calculations. This does not cause problems in the data, as it is in the order 109times smaller than the Q in normally functioning capillaries. The biograph however shows them as small gray arrows instead of yellow, because it’s still a nonzero Q. This is especially visible when the honeycomb is fully blocked. This could give a distorted view of the blood flow, so these minor MATLAB errors should be recognized and removed.

VI.

Conclusion

The honeycomb configuration shows a higher ability than the bifurcating tree to protect the inflow and (hence) the capillary flow rates as more branches are occluded. Also, the pressure distribution can be stabilized better. Because the essential difference between the honeycomb

and the bifurcating tree is the amount of collat-erals, this is an indication of sneak routes being able to protect perfusion significantly. Many adjustments can be made however to improve this model, as discussed in section V.

VII.

Acknowledgements and

additional information

I would like to offer special thanks to Prof. Dr. E.T. van Bavel for accepting my Bachelor’s Project proposal, during a very crowded pe-riod at the AMC. Also greatly appreciated was the ability of supervisor Prof. Dr. E.T. van Bavel and 2nd assessor Dr. H.A. Marquering to assess my work at reasonably short notice. Thanks to Sven Boots MSc in collaboration with Prof. Dr. E.T. van Bavel for providing me with the basis of the program code, the Wheatstone case.

In this paper, the built model was broadly discussed using the most important func-tions. The whole code comprises of 20 functions, and can be requested by e-mailing bobby.runderkamp@student.uva.nl or bob-byrunderkamp@hotmail.com if desired. The code is made as user-friendly as possible, with text files explaining every function.

References

[1] National Cancer Institute STM. Classification & structure of blood vessels;. Available from: http://training.seer.cancer.gov/anatomy/ cardiovascular/blood/classification.html. [2] Hacking WJ, Spaan JA, Van Bavel E. Shear stress is

not sufficient to control growth of vascular networks: a model study. American Journal of Physiology. 1996;. [3] Personal communication Prof dr van Bavel ET. Aca-demic Medical Center. Amsterdam. april-june;2016. [4] Wikipedia (high school knowledge) Circulatory

sys-tem;. Available from: https://en.wikipedia.org/ wiki/Circulatory_system.

[5] Dawson CA, Krenz GS. Flow and Pressure Distribu-tions in Vascular Networks Consisting of Distensible Vessels;.

[6] Hagen-Poiseuille equation;. Available from:https: //en.wikipedia.org/wiki/Hagen-Poiseuille_ equation.

[7] Boots S, VanBavel E. wheatstone network solver. 2015;.

[8] Murray CD. The Physiological Principle of Minimum Work: I. The Vascular System and the Cost of Blood Volume. Proceedings of the National Academy of Sciences of the United States of America. 1926;. [9] Murray CD. The Physiological Principle of

Mini-mum Work: II. Oxygen Exchange in Capillaries. Pro-ceedings of the National Academy of Sciences of the United States of America. 1926;.

[10] Spaan JAE, VanBavel E. Branching Patterns in the Porcine Coronary Arterial Tree. Circulation Research. 1992;.