Citation
Investigating the role of distal flow in progressive pulmonary vascular disease using a computational model

Material Information

Title:
Investigating the role of distal flow in progressive pulmonary vascular disease using a computational model
Creator:
Walter, Rachelle L.
Place of Publication:
Denver, CO
Publisher:
University of Colorado Denver
Publication Date:
Language:
English

Thesis/Dissertation Information

Degree:
Master's ( Master of science)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
Department of Bioengineering, CU Denver
Degree Disciplines:
Bioengineering
Committee Chair:
Kheyfets, Vitaly O.
Committee Members:
Stenmark, Kurt R.
Hunter, Kendall S.

Record Information

Source Institution:
University of Colorado Denver
Rights Management:
All applicable rights reserved by the source institution and holding location.

Downloads

This item has the following downloads:


Full Text
INVESTIGATING THE ROLE OF DISTAL FLOW IN PROGRESSIVE
PULMONARY VASCULAR DISEASE USING A COMPUTATIONAL MODEL.
by
RACHELLE L. WALTER
B.S., University of Colorado Denver, 2017
A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Bioengineering
2019


This thesis for the Master of Science degree by Rachelle L. Walter has been approved for the Bioengineering Program by
Vitaly O. Kheyfets, Chair Kurt R. Stenmark
Kendall S. Hunter


Walter, Rachelle L. (M.S., Bioengineering)
Investigating the role of distal flow in progressive pulmonary vascular disease using a computational model.
Thesis directed by Vitaly O. Kheyfets
ABSTRACT
Pulmonary Arterial Hypertension (PAH), which is associated with pulmonary arterial stiffening and thickening of distal pulmonary vessel walls, leads to right heart failure and sometimes death. Increased thickness of the distal vessel walls primarily occurs in the adventitia and is thought to be due to increased mechanical stress. To provide evidence for this relationship, mechanical stress was simulated for the entire pulmonary vascular network using a combined 1D/0D computational model, which applies a physiological flow waveform at the inlet and a 2 or 3-element outflow boundary condition at the outlets. The computational model was used to evaluate how changes in proximal vascular stiffness impact pulse pressure and wall shear stress in the distal vascular bed in the context of pulmonary hypertension (PH). The four simulated test groups were control, PH, PH with a banded left pulmonary artery (with decreased radius and increased stiffness), and PH with the left and right pulmonary arteries replaced with compliant arteries (decreased stiffness back to control levels). The simulated results showed increased pulse pressure for PH (8.45 mmHg) compared to healthy control (2.39 mmHg). The PH + banded left pulmonary artery showed a slight decrease in pulse pressure compared to PH in the left artery descendants (7.17mmHg) and an increase compared to PH in the right artery descendants(10.12mmHg) artery descendants.The PH + replacement left and right pulmonary arteries had a pulse pressure in between the control and PH levels (4.38mmHg). When only vessel stiffening was modeled, the WSS normalized to control (1.0) showed a decrease for PH (0.79), further decreased for arteries downstream


of the LPA in the banded LPA case (0.67). WSS in the arteries downstream of the RPA in the banded LPA case (0.95) were about equivalent to that of controls, and the WSS in PH with replacement arteries case (0.62) showed a decrease. Once pruning and decreased vessel area were modeled, the WSS normalized to control (1.0) showed an increase for PH (1.5), and a further increase for arteries downstream of the RPA in the banded LPA case (1.80). WSS in the arteries downstream of the LPA in the banded LPA case (1.28) and the WSS in PH with replacement arteries case (1.18) were less than that of the PH case, but still slightly above the control WSS level. These results suggest that a replacement left and right pulmonary artery might be sufficient to return the pulse pressure and WSS to a healthy level given further compliance to the proximal vessels. Further studies should be done to directly investigate the relationship between mechanical stress and vessel remodeling using an animal model. Determining whether or not mechanical stress is causing pulmonary vessel remodeling could lead to the use of mechanical stress as a potential therapeutic target to manage or even reverse PAH disease progression.
The form and content of this abstract are approved. I recommend its publication.
Approved: Vitaly O. Kheyfets
IV


ACKNOWLEDGMENT
I would like to thank:
My Advisors: Dr. Vitaly Kheyfets, Dr. Kendall Hunter, and Dr. Kurt Stenmark My family and friends for providing support
Dr. Michael Yeager, for his mentorship before embarking on this project My fellow Bioengineering Students
Other Bioengineering faculty/staff who have provided invaluable mentorship to me
v


TABLE OF CONTENTS
Tables...................................................................... ix
Figures ..................................................................... x
Chapter
1. Introduction.............................................................. 1
1.1 Pulmonary Hypertension.............................................. 1
1.1.1 Interplay of Mechanics and Biology........................... 2
1.2 Computational Modeling.............................................. 2
1.3 Justification....................................................... 5
1.4 Goals............................................................... 6
2. Methods................................................................... 8
2.1 Model Overview...................................................... 8
2.1.1 Assumptions.................................................. 9
2.2 ID Model (VaMpy).................................................... 9
2.2.1 Equations.................................................... 9
2.2.2 Boundary Conditions......................................... 11
2.2.3 Code Adaptation............................................. 13
2.2.4 Running a ID Simulation..................................... 16
2.2.5 Limitations................................................. 16
2.3 OD Model........................................................... 16
2.3.1 Vessel Network Generation .................................. 17
2.3.2 Simulink Model.............................................. 19
2.3.3 Boundary Conditions......................................... 21
2.3.4 Running a OD Simulation..................................... 21
2.3.5 OD Simulation Data.......................................... 22
2.3.6 Limitations................................................. 23
2.4 Numerical Simulations for Combined Models.......................... 23
vi


2.5 Wall Shear Stress and Pulse Pressure .............................. 26
3. Results..................................................................... 29
3.1 Model Validation and Optimization.................................... 29
3.1.1 ID Model Validation........................................... 29
3.1.2 OD Model Validation........................................... 32
3.1.3 Mesh Independence ............................................ 34
3.2 Simulations.......................................................... 40
3.2.1 ID Simulation Time............................................ 40
3.2.2 Pressure...................................................... 41
3.2.3 Flow.......................................................... 41
3.3 Distal Wall Shear Stress and Pulse Pressure.......................... 41
3.3.1 Effect of proximal stiffening on distal WSS and PP............ 52
3.3.2 Differences in vessel geometry impact distal WSS and PP . . 52
4. Discussion.................................................................. 55
4.1 Pressure Waves ...................................................... 55
4.2 Flow Waves .......................................................... 56
4.3 Distal Wall Shear Stress and Pulse Pressure.......................... 56
4.4 Limitations and Future Work.......................................... 57
5. Conclusion.................................................................. 59
References.................................................................. 61
Appendix
A. Validation Plots for ID Model........................................... 65
A.l Control.............................................................. 65
A. 1.1 Pressure Plots ............................................... 65
A. 1.2 Flow Plots .............................. 67
A.2 PH................................................................... 69
A.2.1 Pressure Plots ............................................... 69
vii


A.2.2 Flow Plots .............................................. 71
B. Mesh Independence Plots.............................................. 73
B.l Changes to Peak Pressure for varying N / dx.................... 73
B. 2 Changes to Peak Pressure for varying dt........................ 74
C. Model Pressure Plots................................................. 76
C. l Control........................................................ 76
C.2 PH............................................................. 80
C.3 PH + LPAB...................................................... 84
C. 4 PH + replacement............................................... 91
D. Model Flow Plots..................................................... 95
D. l Control........................................................ 95
D.2 PH............................................................. 97
D.3 PH + LPAB........................................... 99
D. 4 PH + replacement................................... 104
E. MATLAB Code...................................................... 106
E. l Vessel Network Generation Script ............................. 106
E.2 0D Model Script................................................. 109
E.3 WSS and PP Script............................................. 113
viii


TABLES
Table
2.1 Units comparison between the hydraulic or blood flow model and the
electric or circuit model............................................... 19
2.2 Resistance, Rt, and Capacitance, Ct, at the outlet of the ID model. The scaling factor is applied to the Poiseuille law as shown in equation 2.9 in the OD model when building up the tree such that the total resistance and total capacitance in the OD model are equal to Rt and Ct, respectively. 22
2.3 Upstream radii, RU1 downstream radii, Rd and length-to-radius ratio, lam,
values for each vessel for a tree with asymmetric vessel geometry....... 24
2.4 Upstream radii, Ru, downstream radii, Rd and length-to-radius ratio, lam,
values for each vessel for a tree with symmetric vessel geometry........ 24
2.5 Variables used across all simulations for validation and actual model simulations where rc is the characteristic radius, qc is the characteristic flow, p is blood density, v is the blood viscosity, depth is total number of generations, T is the time period for one cardiac cycle, tc is the number of cardiac cycles to solve for, dt is the time step, and N is the total number
of space steps.......................................................... 25
2.6 Simulation variables used between test groups............................. 26
3.1 The dx values for the largest (dx0) and smallest (daq) vessel sizes for a
given N value........................................................... 34
3.2 The total time each simulation ran for a given N value................... 35
3.3 The total time each simulation ran for a given dt value.................. 38
3.4 The total time each simulation ran for the four test groups.............. 41
IX


FIGURES
Figure
1.1 Graphic representation of 3D, ID and OD models and the trade-off between computational cost and physiological relevance (i.e. how representative
the model is).......................................................... 4
2.1 High-Level Overview of the Full Computational Model..................... 8
2.2 An example of a tree with 3 generations where vessel 0 represents the
main pulmonary artery (MPA), vessel 1 represents the left pulmonary artery (LPA), vessel 2 represents the right pulmonary artery (RPA) and vessels 3-6 are the terminal vessels. The outflow boundary condition is a 3-element windkessel applied at each terminal vessel................... 10
2.3 Inlet flow waveform for control and pulmonary hypertension (PH) test
groups................................................................... 12
2.4 An example of a tree with asymmetric branching where there are terminal
branches at generation 2 and at generation 3........................... 17
2.5 An example of a 2-generation RC circuit tree model made in Simulink
showing the inlet flow and outflow boundary condition (BC)............. 19
2.6 An example of a 2-generation RC circuit tree model made in Simulink
with Voltmeters, Amp Meters, and scope blocks.......................... 20
2.7 The final tree model used for OD simulations. For more information on
what each component or block represents, refer to figure 2.6........... 28
3.1 Example tree of a single bifurcation.................................... 29
3.2 Control test group results for the validation study done on a single bifurcation. The top plots are pressure (mmHg) and the bottom plots are flow
rate (cm3/sec) .......................................................... 30
x


3.3 PH test group results for the validation study done on a single bifurcation.
The top plots are pressure (mmHg) and the bottom plots are flow rate (cm3/sec)........................................................... 31
3.4 Validation between a single RC circuit in Simulink versus solving a single
RC circuit in MATLAB using ODE45....................................... 32
3.5 Validation for a 2 generation RC circuit. (A) shows the total Pressure
of a 2 generation circuit in Simulink versus solving the equivalent single RC circuit in MATLAB using ODE45. (B) shows the pressures for each individual vessel solved in Simulink...................................... 33
3.6 Maximum Relative Error (compared to N = 500) for the largest vessel,
Vessel 0, for varying N / dx values as well as the total time it took for each simulation to run.............................................. 36
3.7 Maximum Relative Error (compared to N = 500) for the smallest vessel,
Vessel 4, for varying N / dx values as well as the total time it took for each simulation to run.............................................. 37
3.8 Maximum Relative Error (compared to dt = le-8) for the largest vessel,
Vessel 0, for varying dt values as well as the total time it took for each simulation to run................................................... 39
3.9 Maximum Relative Error (compared to dt = le-8) for the smallest vessel,
Vessel 4, for varying dt values as well as the total time it took for each simulation to run................................................... 39
3.10 Control test group pressure results for the entire vascular tree. Each vessel at the same generation has the same results due to being completely symmetric, so only one vessel per generation is shown here for simplicity. 42
3.11 PH test group pressure results for the entire vascular tree. Each vessel at
the same generation has the same results due to being completely symmetric, so only one vessel per generation is shown here for simplicity. . . 43
xi


3.12 PH +LPAB test group pressure results for the entire vascular tree. Each vessel at the same generation has the same results due to being completely symmetric, so only one vessel per generation is shown here for simplicity. 44
3.13 PH + replacement test group pressure results for the entire vascular tree.
Each vessel at the same generation has the same results due to being completely symmetric, so only one vessel per generation is shown here for simplicity................................................................... 45
3.14 Pressure results for the initial vessel of all of the test groups of the OD
model simulations............................................................ 46
3.15 Pressure results for the terminal vessels of all of the test groups of the OD
model simulations............................................................ 46
3.16 Control test group pressure results for the entire vascular tree. Each vessel at the same generation has the same results due to being completely symmetric, so only one vessel per generation is shown here for simplicity. 47
3.17 PH test group pressure results for the entire vascular tree. Each vessel at
the same generation has the same results due to being completely symmetric, so only one vessel per generation is shown here for simplicity. . . 48
3.18 PH+LPAB test group pressure results for the entire vascular tree. Each vessel at the same generation has the same results due to being completely symmetric, so only one vessel per generation is shown here for simplicity. 49
3.19 PH + replacement test group pressure results for the entire vascular tree.
Each vessel at the same generation has the same results due to being completely symmetric, so only one vessel per generation is shown here for simplicity................................................................ 50
3.20 Flow results for the initial vessel of all of the test groups of the OD model
simulations............................................................... 51
xii


3.21 Flow results for the initial vessel of all of the test groups of the OD model
simulations................................................................ 51
3.22 WSS results for the OD model simulations for a terminal vessel given
changes to proximal vessel stiffness....................................... 52
3.23 Pulse pressure results for the OD model simulations for a terminal vessel
given changes to vessel stiffness.......................................... 53
3.24 WSS results for the OD model simulations for a terminal vessel given both changes to proximal vessel stiffness, and changes to distal vessel geometry. 53
3.25 Pulse pressure results for the OD model simulations for a terminal vessel
given both changes to proximal vessel stiffness, and changes to distal vessel geometry.............................................................. 54
A.l Validation Pressure Results of the Control test group for the root vessel
0 (MPA).................................................................... 65
A.2 Validation Pressure Results of the Control test group for vessel 1 (LPA) 66
A.3 Validation Pressure Results of the Control test group for vessel 2 (RPA) 66
A.4 Validation Flow Results of the Control test group for the root vessel 0
(MPA)...................................................................... 67
A.5 Validation Flow Results of the Control test group for vessel 1 (LPA) . . 68
A.6 Validation Flow Results of the Control test group for vessel 2 (RPA) . . 68
A.7 Validation Pressure Results of the PH test group for the root vessel 0 (MPA) 69 A.8 Validation Pressure Results of the PH test group for vessel 1 (LPA) ... 70
A.9 Validation Pressure Results of the PH test group for vessel 2 (RPA) ... 70
A. 10 Validation Flow Results of the PH test group for the root vessel 0 (MPA) 71
A. 11 Validation Flow Results of the PH test group for vessel 1 (LPA)........... 72
A. 12 Validation Flow Results of the PH test group for vessel 2 (RPA)........... 72
B. l Changes to Peak Pressure along the length of the largest vessel, Vessel 0,
due to varying N / dx values............................................... 73
xiii


B.2 Changes to Peak Pressure along the length of the smallest vessel, Vessel
4, due to varying N / dx values..................................... 74
B.3 Changes to Peak Pressure along the length of the largest vessel, Vessel 0,
due to varying dt values............................................ 74
B. 4 Changes to Peak Pressure along the length of the smallest vessel, Vessel
4, due to varying dt values......................................... 75
C. l Pressure Results of the Control test group for the root vessel (MPA) . . 76
C.2 Pressure Results of the Control test group for the second generation of
vessels (LPA/RPA)................................................... 77
C.3 Pressure Results of the Control test group for the ID terminal vessels
(Generation 3)...................................................... 78
C.4 Pressure Results of the Control test group for the 0D terminal vessels
(Generation 12)..................................................... 79
C.5 Pressure Results of the PH test group for the root vessel (MPA).. 80
C.6 Pressure Results of the PH test group for the second generation of vessels
(LPA/RPA)........................................................... 81
C.7 Pressure Results of PH test group for the ID terminal vessels (Generation
3).................................................................. 82
C.8 Pressure Results of the PH test group for the 0D terminal vessels (Generation 12)............................................................ 83
C.9 Pressure Results of the PH+LPAB test group for the root vessel (MPA). 84
C.10 Pressure Results of the PH+LPAB (L) test group for the second generation of vessels (LPA) 85
C.ll Pressure Results of the PH+LPAB (R) test group for the second generation of vessels (RPA).................................................. 86
C.12 Pressure Results of the PH+LPAB (L) test group for the ID terminal
vessels that descend from the LPA (Generation 3).................... 87
xiv


C.13 Pressure Results of the PH+LPAB (R) test group for the ID terminal
vessels that descend from the RPA (Generation 3)........................... 88
C.14 Pressure Results of the PH +LPAB (L) test group for the OD terminal
vessels (Generation 12)........................... 89
C.15 Pressure Results of the PH + LPAB (R) test group for the OD terminal
vessels (Generation 12)........................... 90
C.16 Pressure Results of the PH + replacement test group for the root vessel
(MPA)...................................................................... 91
C.17 Pressure Results of the PH + replacement test group for the second generation of vessels (LPA/RPA).......................................... 92
C.18 Pressure Results of the PH+replacement test group for the terminal vessels of the ID model.................................................. 93
C. 19 Pressure Results of the PH +replacement test group for the OD terminal
vessels (Generation 12).................................................... 94
D. l Flow Results of the Control test group for the root vessel (MPA)....... 95
D.2 Flow Results of the Control test group for the second generation of vessels
(LPA/RPA).................................................................. 96
D.3 Flow Results of the Control test group for the ID terminal vessels (Generation 3)............................................................ 96
D.4 Flow Results of the Control test group for the OD terminal vessels (Generation 12)........................................................... 97
D.5 Flow Results of the PH test group for the root vessel (MPA)............. 97
D.6 Flow Results of the PH test group for the second generation of vessels
(LPA/RPA).................................................................. 98
D.7 Flow Results of PH test group for the ID terminal vessels (Generation 3). 98
D.8 Flow Results of the PH test group for the OD terminal vessels (Generation
12)........................................................................ 99
xv


D.9 Flow Results of the PH+LPAB test group for the root vessel (MPA). . . 99
D.10 Flow Results of the PP1+LPAB (L) test group for the second generation
of vessels (LPA)..................................................... 100
D.ll Flow Results of the PP1+LPAB (R) test group for the second generation
of vessels (RPA)..................................................... 100
D.13 Flow Results of the PH+LPAB (R) test group for the ID terminal vessels
that descend from the RPA (Generation 3)........................... 101
D.12 Flow Results of the PH+LPAB (L) test group for the ID terminal vessels
that descend from the LPA (Generation 3)........................... 101
D.14 Flow Results of the PH + LPAB (L) test group for the 0D terminal vessels
(Generation 12)...................................................... 102
D.15 Flow Results of the PH + LPAB (R) test group for the 0D terminal vessels
(Generation 12)...................................................... 103
D.16 Flow Results of the PH + replacement test group for the root vessel (MPA). 104 D.17 Flow Results of the PH + replacement test group for the second generation
of vessels (LPA/RPA)................................................. 104
D.18 Flow Results of the PH+replacement test group for the ID terminal vessels
(Generation 3)....................................................... 105
D.19 Flow Results of the PH + replacement test group for the 0D terminal
vessels (Generation 12).............................................. 105
xvi


1. Introduction
1.1 Pulmonary Hypertension
Pulmonary hypertension (PH) is diagnosed by a mean pulmonary arterial pressure (mPAP) greater than 25 mmHg at rest or 30 mmHg with exercise. There are currently five categories according to the Dana Point classification. Category 1, also called pulmonary arterial hypertension (PAH), is the most common category and includes idiopathic PAH, heritable PAH, drug and toxin induced and other associations such as with connective tissue disease, HIV infection, portal hypertension, congenital heart diseases and schistosomiasis. The other four categories consist of PH due to left heart disease (Cat. 2), PH due to lung diseases and/or hypoxia (Cat. 3), chronic thromboebolic PH (CTEPH, Cat. 4), and others with unclear multifaceted mechanisms (Cat. 5) [1, 2, 3]. Category 1 is focused on here due to its prevalence; however, some information may apply to other categories.
PH is associated with changes to the vascular wall in both the proximal and distal vasculature, which is composed of three layers: the intima, media and adventita, which are primarily composed of endothelial cells, smooth muscle cells, and fibroblasts, respectively [4], During PH, smooth muscle cells (SMCs) proliferate causing the medial layer to grow in towards the lumen, decreasing the vessel radius, resulting in increased pulmonary vascular resistance (PVR), which is seen by many as the only major consequence of PH [5]. Fibroblasts are also proliferating which causes the adventitia to grow out and away from the lumen. The muscularization of the media and stiffening of the adventitia occur in both proximal and distal pulmonary arteries in PH [1, 6].
Adventitial fibroblasts play a role in increasing the medial layer as they can differentiate into myofibroblasts which are smooth muscle-like cells, and migrate to the medial layer. Fibroblasts in the adventitia also change the extracellular matrix (ECM) composition of vascular walls in PH; more specifically, there is an increase
1


of collagen types I and III. An increase in collagen in the vascular wall leads to an increase in stiffness of the wall as well [3, 4, 5].
Endothelial cell (EC) dysfunction is known to occur, which can further support SMC proliferation and can cause plexus lesion formation as well as endothelial channel formation. EC dysfunction is also associated with loss of pulmonary arterial vessels [1, 6]. As part of the vascular remodeling process, leukocytes are being recruited by adventitial fibroblasts and pro-inflammatory cytokine levels are higher in PH vessels [4], which may further act to exacerbate these vascular changes.
1.1.1 Interplay of Mechanics and Biology
The vascular structural changes of the pulmonary vessels in PH are likely to cause changes to the mechanical properties of the vessel; for example, an increase in collagen causes an increase in vascular stiffness. Vascular stiffening is seen by many as a key player of PH disease progression [3, 7, 5]. Stiffening of the large arteries in PH contributes to the increase in mPAP and also can cause an increase in pulse pressure, reflected waves, pulse wave velocity, and pulse wave intensity. All of these factors can lead to increased flow pulsatility. Decreasing lumen radius results in an increase in distal vessel resistance. A decrease in vessel number along with decreasing lumen radius likely causes an increase in distal wall shear stress (WSS). Both an increase in flow pulsatility and WSS in the distal vasculature then contributes to EC dysfunction which further exacerbates the inflammation and remodeling of the vascular walls [7]. Thus WSS and flow pulsatility are important mechanical properties to investigate.
1.2 Computational Modeling
The gold standard for diagnosing pulmonary hypertension (PH) is through right heart catheterization (RHC), which is an invasive procedure that generally uses a Swan-Ganz catheter to measure pressure and mean flow (cardiac output) [8]. A catheter is not able to travel far into the vessel network due to its size; thus, there is
2


not a way to directly measure flow in the distal vasculature. Other techniques exist to measure blood flow non-invasively such as magnetic resonance imaging (MRI) and ultrasonic techniques; however, they also come with limitations that do now allow for the direct measurement of blood flow in downstream vessels [9]. Due to these limitations, computational modeling has become an attractive method for estimating flow in distal vasculature.
Computational models are non-invasive ways to ask physiological questions about the body. Generally, they give us the ability to investigate hard-to-measure properties, to make predictions to create more effective animal models, and to better investigate treatment options. There are many types of computational models, such as zerodimensional (OD) [10, 11], one-dimensional (ID) [12, 13, 14, 15, 16, 17, 18, 19, 20, 21] and three-dimensional (3D) [18, 19, 20, 21]. OD models are point-based and only include the time domain, ID models add one spatial domain which is the length along a vessel, and 3D models include all three spatial dimensions of the vessel as well as the time domain. Models can also be a combination of these types.
3


Most
Computationally
Expensive
3D
ID —
----X
Least
Computation ally Expensive
OD
All include the time domain
Most
Representative
Least
Representative
Figure 1.1: Graphic representation of 3D, ID and OD models and the trade-off between computational cost and physiological relevance (i.e. how representative the model is).
Each model gives a different set of information about flow; 3D models allow for representation of more complex flow but at a high computational cost while a OD model has a very low computational cost but only includes more simple flow features. This is represented in figure 1.1. Typically, OD models are used to look at the system-level, ID models can give accurate blood pressure and flow propagation information for a large network of vessels, and 3D models allow for a deeper investigation of fluid-structure interactions (FSI), which is the relationship between blood flow and vessel wall deformation, and can fully represent blood velocity profiles (ID models can only approximate) [9, 22], Some models exist to study an entire vascular network [12, 13, 16, 21]. Models also exist to study a specific subsystem, such as a pulmonary vasculature model [14, 23], which is what is used here.
For all models, inlet and outlet boundary conditions are required as well as geometric information (e.g. radius, length) and mechanical properties (e.g. stiffness,
4


viscosity) of the blood and vessels. The inlet boundary condition typically comes from clinical measurements such as data from RHC or MRI. Outlet boundary conditions vary based on the modeling approach [9, 24], Common outlet boundary conditions used are Windkessel variations [15, 17, 18, 19], and structured tree [12, 13]. Some 3D models may use ID models for their outflow boundary conditions as well [25, 20, 21]. If a network (i.e. not just a single vessel) is being modeled, a bifurcation boundary condition is also needed at vessel junctions [9, 24],
Due to the ability to have a relatively low computational cost with accurate representation of blood flow and some mechanical properties of the vessels, a ID model is one of the most attractive options for modeling blood vessels and is used here for that reason. Currently, ID model assumptions allow for a good representation of relatively larger blood vessels and use simplified 3D Navier-Stokes equations. Common numerical schemes exist to solve the system of equations set up for ID models, such as Galerkin variations, Finite Differences methods, and characteristic methods [24, 26].
1.3 Justification
In PH, there is distal vessel thickening occurring. Data from other literature suggests that proximal stiffening causes higher flow pulsatility which causes endothelial cell dysfunction and activates fibroblasts. Activated fibroblasts proliferate and deposit collagen which causes stiffening and thickening of the distal vessels. One study provides great evidence for this proposed mechanism by showing that hemodynamic unloading is effective at treating PH phenotypes [27]. In this study, PH-modeled rats were compared to PH-modeled rats with a banded left pulmonary artery (LPA) and studied the histology of the rats’ arteries. They found that the LPA of the banded rats had reduced expression of pro-inflammatory cytokines, and other inflammatory and proliferative markers, which almost matched that of the normal rat. Blood flow to the LPA in the banded rats reduced to less than half of the original blood flow,
while blood flow increased by 40% in the right pulmonary artery [27]. This study
5


almost completely reversed what we think to be associated with PH. Given that banding showed reversal of these disease characteristics, the question of what kind of mechanical changes are occurring to see these biological changes. This paper showed a biological and mechanical relationship between the proximal vasculature, but the relationship between mechanics and biology of the distal vasculature is still missing.
In order to study the distal flow mechanics, a computational model can be used. Given a banded left pulmonary artery as done in [27], flow characteristics in the distal vasculature can be investigated. There is also the potential for replacing the left and right pulmonary arteries with synthetic pulmonary arteries to determine distal flow mechanics. If banding and replacement of the proximal vessels show similar distal flow mechanics, then replacing the left and right pulmonary arteries could be used as a treatment or reversal of PH. Thus four test groups will be simulated here: healthy control, PH, PH with banded LPA, and PH with replacement left and right pulmonary arteries. These simulations will help us determine if proximal vessel stiffness is affecting distal flow. If this relationship is established, then animal model studies could be done simultaneously with animal-specific simulations to see if there is also a reversal of these disease phenotypes as seen in [27] due to changes to proximal vessel stiffness as done in these simulations. Establishing this relationship in these simulations, and confirming them in an animal study, would provide a great starting point to investigate replacing stiffened proximal vessels with synthetic, compliant vessels as a therapeutic treatment for PH.
1.4 Goals
Given a stiffer, vascular network in the pulmonary hypertension (PH) case, we expect to see an increase in flow pulsatility [7], investigated via pulse pressure, and a decrease in wall shear stress (WSS) [28, 29]. We hypothesize that restoring the elasticity of proximal vasculature of a PH mouse model will return pulse pressure and
WSS of distal pulmonary arterioles to control levels.
6


In this study, the relationship between proximal pulmonary artery stiffness and distal flow were investigated using a coupled 1D-0D model. First the 1D-0D model was created in Python and MATLAB and validated against previously published data [23]. Then, simulations were run to simulate four different test groups based on [27] to study changes in WSS and pulse pressure. The aims are described below.
Specific Aim 1: Develop a computational model of the pulmonary arterial system to allow solving of blood flow hemodynamics in a large-scale network.
Specific Aim 2: Execute pulmonary mouse-specific simulations (matching proximal pulmonary arterial stiffness and pulmonary vascular resistance) to compare pulse pressure and wall shear stress of distal vasculature of the following models: control; pulmonary hypertension (PH); PH with banded left pulmonary artery; and PH with a synthetic, compliant left and right pulmonary artery replacement.
7


2. Methods
2.1 Model Overview
To establish evidence that a relationship exists between proximal mechanical stiffness and distal flow, the entire pulmonary vascular network will be modeled using a hybrid 1D/0D model where the ID model is for the larger, proximal vessels and the OD model is used for the smaller, distal vessels. First, the vessel properties and tree are established and put into a configuration file. Then, the ID model sets up the artery network using the configuration file. A physiological waveform is input into the ID model which computes the proximal vessels’ flow, area and pressure waves. The output from the ID model is the input for the OD model including resistance and capacitance (RC) values, and the ID model outlet flow rate. The OD model then sets up a vascular tree using RC values where the RC values of the entire OD model sum to equal the RC values at the outlet of the ID model. The OD model uses the ID model outlet flow as its inlet flow to compute the flow and pressure waves for the distal vessels. The resulting distal pressure and flow are then used to calculate pulse pressure and wall shear stress. A flow chart for the model is shown below:
Figure 2.1: High-Level Overview of the Full Computational Model
2.1.1 Assumptions
8


The models assume the vessels are elastic and axisymmetric tubes, containing Newtonian fluid. Other assumptions for the ID model are that the blood pressure is uniform across r, the radius, and there is pulsatile laminar flow.
2.2 ID Model (VaMpy)
The ID computational model utilizes fluid dynamics, wall mechanics, and boundary conditions. The resulting system of quasi-linear hyperbolic PDEs consists of three equations which are used to solve for flux inside of the vessel lumen, cross-sectional area, and pressure. The three equations are solved using a python-based, open source model available on github, named the Vascular Modeling in Python toolkit or VaMpy [30]. This model starts with a single parent vessel that bifurcates into daughter vessels. The resulting two daughter vessels also bifurcate. This process continues until the terminal vessels are reached. The total number of vessels made is denoted by 2dept/l-l. Therefore, if the artery network had a depth of 3, aka 3 generations, there would be 7 vessels. Note that vessel numbering begins with 0. As an object-oriented program, each vessel is contained in its own “artery” class which are all contained within in the “ArteryNetwork” class. The VaMpy model uses a configuration hie that contains the vessel geometry and vascular property information. An example vascular tree with 3 generations containing outflow boundary conditions is shown in figure 2.2 referencing vessels in the pulmonary vasculature starting with the main pulmonary artery (MPA).
2.2.1 Equations
The ID model solves three equations. The first two equations are derived through integration of the Navier-Stokes equations for conservation of mass and momentum in a one-dimensional coordinate system. The details of this solution can be found elsewhere [31] with the resulting system of non-dimensionalized equations represented by the equation shown below in Eq 2.1
SU SF_ _ -
St 8x
9
(2.1)


Outflow BCfor the ID model
Figure 2.2: An example of a tree with 3 generations where vessel 0 represents the main pulmonary artery (MPA), vessel 1 represents the left pulmonary artery (LPA), vessel 2 represents the right pulmonary artery (RPA) and vessels 3-6 are the terminal vessels. The outflow boundary condition is a 3-element windkessel applied at each terminal vessel.
where
A(:r, t) ( q(x,t) \ M
u = ,F = + f(.r0)y/Ao(x)A(.x,t)J , S = W
Si
2ixR(x,t) q(x,t) 5b Re A(x,t)
pyvwXvV/w + v/TTSA -
U contains the unknowns to be solved: A(x,t), the cross-sectional area, and q(x,t), the flow rate. Elasticity of the vessel is described by f(r0) = and is empirically estimated using Eq 2.2 below:
F,h
---^ = /,' | rA ' + k3 (2.2)
r0
10


R(x,t) represents the radius of the vessel, r0(x) is the initially relaxed vessel radius with corresponding cross-sectional area Aa(x). 6b is the boundary layer thickness estimated by 6b = \JvT/(27r) where v is the kinetmatic viscosity, Re is the Reynold’s number and is calculated using Re = qcJ(vrc) with qc and rc being the characteristic flux (flow) and characteristic radius, respectively, and T is the time period of one cardiac cycle. To model a tapering vessel with upstream and downstream radii, Ru and Rj, the vessel radius can be described as a function down its length L as r0{x) = Ruexp(log( fy)f )■
Eq 2.1 is solved by discretizing using Richtmyer’s two-step Lax Wendroff explicit scheme, detailed in [30, 31] which is second-order accurate in time and space, i.e. 0(At2, Ax2) and must satisfy the CFL condition: At < Ax\^ ±c|-1 where the wave speed is c = \J(A/p)(6p/6A) and Ax and At are the space step and time step, respectively [30].
The third equation needed to solve this system of equations (3 equations with 3 unknowns) is the state equation which is a linear equation relating pressure p(x,t) to area as shown below:
p(x,t) =p0 + f(r0)( 1 - y/Ao(x)/A(x,t)) (2.3)
where pa is the initial (diastolic) pressure [30]. Thus, this model outputs three variables: cross-sectional area, A(x,t), flow rate, q(x,t), and pressure, p(x,t)
2.2.2 Boundary Conditions
To solve these three equations in this ID model, boundary conditions must also be applied at the inlet, at bifurcations, and at the outlets. At the inlet, a physiological flow waveform is given as measured by ultrasound done simultaneously with catheterization. The inflow waveform used for the control, and PH cases are shown in figure 2.3.
A bifurcation boundary condition applies between any parent vessel and its two
daughter vessels. The bifurcation boundary condition consists of pressure continuity,
11


0.6
Figure 2.3: Inlet flow waveform for control and pulmonary hypertension (PH) test groups.
pp(L,t) = pdii0,t), and conservation of flow, qp(L,t) = 5^(0, t). To complete
i= 1,2
the system, these relationships lead to a system of eighteen equations and eighteen unknowns. The system of equations, shown below in Eq 2.4, is solved using Newton’s method.
*fc+i = xk ~ ([J(xk)])~1fj(xk) for k = 1, 2, 3, ... 18 (2.4)
where k is the current iteration from I to 18, [J(;zg.)] is the Jacobian matrix and fj(xk) is the vector of residuals [30].
At the outlet, a zero-dimensional, 3-element Windkessel outflow boundary condition is used to close the computational domain, and is characterized by resistance and capacitance (RCR) values. The same set of RCR values, specified by the configuration file, are applied at all of the outlet vessels. A representation of the RCR circuit as an outlet boundary condition is shown in figure 2.2. These RCR values are used to estimate the outflow of the vessel as shown below in Eq 2.5 where M is the
12


final space point and n is the current time step.
«M+1
„„,p’m'-Pmi AtPit + R2)
9m “T
i?i RiR2Ct RiR2Ct
Area is solved for using a discretized mass conservation equation shown below:
(2.5)
'4"I = '4“- <2'6)
The outlet pressure is then solved using a discretized version of Eq 2.3 and is shown below:
pT = /m+I (l - (2.7)
This process is iterated through after an initial guess for pr^1 until the outlet pressure solution converges, where the error is < le-6. The iteration is limited to a maximum of 1000.
2.2.3 Code Adaptation
The VaMpy model has limitations and therefore had to be expanded upon for this project. The changes made to the model include: solving for multiple bifurcations, using different RCR values for each terminal vessel, allowing dx to vary for each vessel and using different ^ values per vessel than that calculated from Eq 2.2. Other functionality was also included for internal use such as tracking the runtime.
One major limitation of the original VaMpy model, was the fact that it could only solve for a single bifurcation. In order to use this model for a larger network, the model was expanded to solve for multiple bifurcations by changing the indexing of the vessels. Initially, the code indexed the daughter vessels with d\ = p+1 and d2 = p +2 inside of the get-daughters function of the ArteryNetwork class. The get_daughters function was updated to instead return d\ = 2p+ 1 and d2 = 2p + 2 where the parent index starts at 0.
In the initial VaMpy model, RCR values were input into a configuration hie which
were directly used at the terminal vessels. The model has since been changed so that
13


the configuration file instead contains an Rt and Ct value which represent the total vascular resistance and total capacitance, respectively. These values are computed in the MPA and therefore represent the vascular resistance and capacitance for the entire network. These values are then used to iteratively calculate the Rtj and Ctj values for each vessel until reaching the terminal vessels. First, the mean input flow, q, and mean input pressure, p, are calculated. The mean flow is just the average of the inlet waveform and the mean pressure is calculated from RT = where pc is the capillary pressure and drops out of this equation with the assumption that the capillary pressure is 0 and that the mean pressure stays constant for all vessels. The terminal resistance values are then calculated using RTj = A, where qj is the mean flow for vessel j and is calculated by using Poiseuille’s law at each bifurcation until the terminal vessels are reached:
Q l^Gdi \ofiL J d,
i
Gdi is the vessel conductance, qp is the mean flow to the parent vessel and qdi represents the mean flow to the each of the two daughter vessels. Once Rtj has been calculated for the terminal vessels, Ctj can be found using the relationship r = RtCt and r = RtjCtj assuming that r is constant throughout the network. Rtj can then be split up into R\j and i?2j: Rij = ciRtj where a = 0.2 and i?2j = Rtj ~ Ri,j- [23]. This was implemented as a new function inside of the ArteryNetwork class.
Dynamically changing dx based on the vessel was done by applying an N value. This N value specifies the number of space steps to include for all vessels, which means that each vessel will have a different dx based on its length rather than all vessels using the same dx value. In order to include this functionality, a new mesh function was created in the ArteryNetwork class, meshrnx, with the old mesh function being renamed to mesh_dx. The mesh_dx function takes in the dx value and applies the mesh function in the Artery class on each function. The new meshrnx function takes
14


in two variables: nx, the desired number of space steps, and re, the characteristic radius. N is a new variable to include in the configuration hie and is fed in as the nx input for the meshmx function. The nx and re values are then used to calculate dx using the equation: dx = yy3[ which was obtained by rearranging the equation for nx in the mesh function of the Artery class: nx = int( y) + 1.
As described previously, y1 is calculated from Eq 2.2 in the original VaMpy model. In order to fulfill the goals of this project, y1 was adjusted to allow the conhg hie to also input specihc y1 values per vessel. The change was made by adding a simple input argument that checks for the existence of ‘Eh’ in the function arguments just as VaMpy has done with other variables. If ‘Eh’ exists, then the code will use the entered value from the conhg hie, otherwise the program will use Eq 2.2 as originally specihed. This was implemented in the mesh functions of both the artery and ArteryNetwork class.
A Time matrix was added to the Artery class to allow tracking of total simulation runtime as well as the runtime for varying processes. The time tracking was implemented by using the time.clock() function, which required importing the time library to be used. The time spent in the bifurcation boundary condition solve function was tracked as well as the time spent in the solve function as a whole. The solve function in the ArteryNetwork is where the actual solving of the ID model equations occurs and is therefore the majority of the time spent for the model. The bifurcation boundary condition was determined to be the most time consuming part of the solve function at around 80% of the runtime.
An error was also found involving the bifurcation boundary condition in the solve function of the ArteryNetwork class. Originally, the boundary condition solution for the two daughter vessels were accidentally stored in the opposite daughter vessel object, so they were swapped to correct this issue.
2.2.4 Running a ID Simulation
15


The Integrated Development Environment (IDE) used to run the models was Anaconda Navigator with Spyder (Python 3.6). In order to run a simulation with this IDE, the following commands were run: runhle(‘bifurcation_example_Rt_N.py’, args = ‘filename.cfg’) to obtain data and subsequently: runhle(‘plot_example.py’, args = ‘filename.cfg’) to obtain plots where the filename.cfg is the configuration hie for any given simulation.
2.2.5 Limitations
One major limitation of the ID model is its slow run time. Initial runtime was 5 days to solve to a depth of 3 for a total of 7 vessels. The model was investigated to find what part of the model took the longest and then determine if the model could be sped up. The culprit was the bifurcation boundary condition due to its 18 non-linear equations being solved by Newton’s method which increases the computational cost. Due to the slowness of the model, the ID model was combined with a OD model in order to model an entire vascular network relatively quickly. Once the parameters were changed for validation, mesh independence, and actual model simulations, the runtime for the ID model decreased to 14 - 35 hours.
Another limitation of the model is that the changes made to ^ are not compatible with VaMpy’s ability to generate a large tree using a and f3 values that represent the change in radius from a parent vessel to its two daughter vessels. Finally, the ID model code is not formatted to allow for an asymmetric branch distribution in the tree. This means that all vessels must bifurcate except terminal vessels, and all the terminal vessels are located at the same depth. For example, the code cannot accommodate a terminal branch at a depth of 2 (generation 2) with terminal branches at a depth of 3 (generation 3) as shown in figure 2.4.
2.3 OD Model
The OD Model utilizes a circuit model based off of a two element windkessel
model, which is characterized by a single resistance and a single capacitance (RC)
16


Generation 1
Generation 2
Generation 3
Figure 2.4: An example of a tree with asymmetric branching where there are terminal branches at generation 2 and at generation 3.
value. Each vessel in this model is represented by its own RC circuit. The RC circuits are connected in a tree with each parent bifurcating into two daughter vessels, just like with the ID model. The circuit model is graphically created using Simulink with the variables for each R and C value being determined using a program created in MATLAB.
2.3.1 Vessel Network Generation
One of the first steps for the OD model to be able to be created, is to create the
entire vascular network, or tree. There are many considerations when setting up the
tree such as branch symmetry, depth, geometry, and mechanical properties. In this
model, the program starts with a root radius, which will correspond to a terminal
vessel from the ID model, and will use a and /5 values as scaling factors for the
daughter vessels for any given parent vessel. The values chosen for a and /5 were
0.9 and 0.6, respectively [32]. The tree continues to grow until a vessel is less than
the terminal radius of 0.005 cm[33, 34], which represents the capillary bed, or at a
maximum of 9 generations to give a total of 12 generations, including the vessels in
the ID model. According to [35], smaller animals have a total of 11-12 generations
17


for pulmonary artery trees, which is why this number was chosen. If a parent vessel would create a daughter vessel that is less than the terminal radius, then the parent vessel becomes a terminal vessel. Each of these values are represented in a matrix, r, where if a parent is indexed by (i,j) then the two daughter vessels would have an index of (i + 1,2j — 1) and (i + 1, 2j) where i and j represent the row and column positions, respectively. The length for each vessel is then calculated using a constant length-to-radius ratio and stored into a matrix, L.
Once the radius and length matrices are created for the network, the Resistance and Capacitance values are calculated for each vessel. The resistance is calculated by using Poiseulle’s flow with the following equation:
R
8 vL
7rm
(2.9)
where v is the dynamic viscosity. Capacitance is then calculated using the same r value and relationship from the ID model. Once resistance and capacitance matrices are solved, they are summed for the entire tree giving Rsum and Csum. For resistances in parallel, the resistance is summed as: where dl and d2 are
the two daughter vessels. Resistances in series add normally as: Rt = R\ + R2 where R\ and R2 are two resistors in series. For capacitance in parallel, the capacitance adds normally as: Ct = C\ + C2 where C\ and C2 are two capacitors in parallel. Capacitance in series is summed like the resistance in parallel is summed and is given by: ^ where dl and d2 are the two daughter vessels. Resistance and
Capacitance are summed independently of one another, as in the capacitance values do not effect the total resistance values and vice versa.
The total values are then compared to the outlet Rt and Ct values from the ID model and then the tree is adjusted such that Rnew = (scaling factor(*i? and
a

where scaling factor = nRT and r is the same as in the ID model.
o H sum.
2.3.2 Simulink Model
18


Table 2.1: Units comparison between the hydraulic or blood flow model and the electric or circuit model.
Model type Flow or Current Pressure or Volt Resistance Capacitance
Hydraulic mL / s dynj cm2 (mm.Hg) g / scm4 4 2 cm s g
Electric Amps Volts Farads Ohms
Once these values are found for the tree, a physical model of the tree is made in Simulink. A translation between the hydraulic, vessel model and electric, circuit model is shown in table 2.1. Therefore, flow will be treated as current (fluid flow becomes electron flow) and pressure will be treated as voltage (fluid potential becomes electrical potential) inside the Simulink model. An example of a 2-generation tree is shown in figure 2.5. For the model to output pressure and flow waves for a vessel,
Figure 2.5: An example of a 2-generation RC circuit tree model made in Simulink showing the inlet flow and outflow boundary condition (BC).
Voltmeters and Amp Meters are needed in the model. An example of the same 2-
generation tree from figure 2.5 is shown in figure 2.6 with these meters added. The
Voltmeters and Amp Meters also need to be connected to a scope block, which is just
a block that provides graphic visualization. In the model you can see a scope block
19


on the left for the voltage and a scope block on the right for the current.
Figure 2.6: An example of a 2-generation RC circuit tree model made in Simulink with Voltmeters, Amp Meters, and scope blocks.
20


Within the individual component blocks, resistance and capacitance can be hardcoded, which is an extremely manual process especially as the network gets larger. Alternatively, variables named in the MATLAB workspace can be used within the blocks such that the values can be easily changed by updating this variable within the MATLAB workspace. Using this method, a few variables can be changed within the MATLAB script that generates the vessel network and then all of the resistance and capacitance values update automatically in the Simulink model.
2.3.3 Boundary Conditions
In order to solve the entire circuit, boundary conditions must be applied. The inlet flow wave comes from the ID model terminal vessels. CSV hies are automatically saved in folders after running a ID simulation, so the how waveform from the last space step of the terminal vessel is used as the inlet boundary condition and converted into current. At each bifurcation, voltage and current are conserved the same way as pressure and how in the ID model where Vp = and Ip = 1^. The
i= 1,2
outlet boundary consists of another RC circuit that represents all of the downstream vasculature from the terminated vessel. The RC outlet circuits are shown in figures
2.5 and 2.6. It should be noted that if a terminal vessel is representing a capillary bed, then no outlet RC circuit is needed because there is no downstream vasculature.
2.3.4 Running a OD Simulation
To make things easier, MATLAB has a command that can be used inside of a script in order to run a Simulink model. A script was made to create the vessel network parameters and then the command sim(‘modelhle.slx’) can be used to run the Simulink model contained in the modelhle.slx hie. Within the Simulink model, data from the scope blocks are saved as variables in the MATLAB workspace, and can be subsequently used in the same MATLAB script for post-processing. Alternatively, the vessel network parameters can be created with a MATLAB script hie, then the
21


Table 2.2: Resistance, Rt, and Capacitance, Ct, at the outlet of the ID model. The scaling factor is applied to the Poiseuille law as shown in equation 2.9 in the OD model when building up the tree such that the total resistance and total capacitance in the OD model are equal to Rt and Ct, respectively.
Test Group p / mmHg*s\ * t rriL ' CT ( mL„ ) ± v rrimti g ' scaling factor (unitless)
Control 233.6 2.2949e-4 1.0697e-4
AUP Htestgroups 535.2 2.945e-5 2.4509e-4
Simulink model can be run manually with the variables still being saved in the MAT-LAB workspace. The variables are saved in the MATLAB workspace by going to Scope —> Settings —> Logging and checking the box labeled “Log data to workspace” and then a “Variable name” can be specified. The MATLAB scripts used to run the simulations are shown in Appendix E.
The final tree used for OD model is shown in figure 2.7 and only follows one line of vessels down for simplicity.
2.3.5 OD Simulation Data
Outlet flow waveforms from the ID model vessels were used in the OD model as inlet boundary conditions. All of the vessel geometries were symmetric, except for the PH + LPAB condition where vessel 1 (LPA) had a smaller radius than its counterpart, vessel 2. Therefore, only one terminal vessel was used for all of the subsequent OD simulations in Simulink, except the PH + LPAB. In this case, two terminal vessels were tested, one daughter vessel for vessel 1 and one daughter vessel for vessel 2 and are labeled as PH + LPAB (L) and PH + LPAB (R), respectively, to indicate that the vessel is a descendant of the left or right pulmonary artery. The Rt, Ct and scaling factor values obtained for the OD model tree buildup are shown in table 2.2.
22


2.3.6 Limitations
One main limitation of the OD model is that it may not actually ’’end” at the capillaries, since it has a strict cutoff of 9 generations; however, this is a good estimate based on literature showing that small mammals have around 12 generations [35]. The terminal radius value for the control group was 0.0069 cm which is a bit higher than the capillary radius value found in literature of 0.005 cm, so may affect how accurate the results are; however, these simulations were used to look at trends so this error does not have a great effect on what was considered here. Simulations could be done in the future to include more generations, or the a and f3 values may need to be adjusted to achieve more accurate terminal radius values.
Another limitation of the 0D model, is the fact that it does not include wave reflections, and that the pressure and flow waves seem to attenuate more than what is seen in the ID model. While these waves could actually attenuate more in the distal vessels, there is no way to test this due to device limitations described previously. Again, this study was interested in trends of the pressure and flow and since this attenuation would be consistent across all test groups, this attenuation does not greatly affect the results presented here.
2.4 Numerical Simulations for Combined Models
Four test groups were investigated here: control, pulmonary hypertension (PH), PH with a banded left pulmonary artery, and PH with synthetic left and right pulmonary arteries, which will be labeled as: control, PH, PH + LPAB, and PH + replacement, respectively. Table 2.3 shows the radius and length-to-radius ratio values for each vessel, taken from data in [23]; however, the radius and length-to-radius ratio values were averaged for all vessels at the same generation such that each vessel in the same generation has a symmetric vessel geometry. This was done to simplify the model and to decrease the computational cost. The resulting averaged values shown
in table 2.4, which shows Upstream radii, Ru, downstream radii, Rdâ–  These values
23


Table 2.3: Upstream radii, Ru, downstream radii, Rd and length-to-radius ratio, lam, values for each vessel for a tree with asymmetric vessel geometry.
Vessel 0 1 2 3 4 5 6
Ru (cm) 0.047 0.026 0.037 0.024 0.013 0.032 0.017
Rd (cm) 0.047 0.026 0.037 0.024 0.013 0.032 0.017
lam (unitless) 8.723 17.1154 10.0541 10.0417 4.0 6.3125 12.4706
Table 2.4: Upstream radii, Ru, downstream radii, Rd and length-to-radius ratio, lam, values for each vessel for a tree with symmetric vessel geometry.
Vessel 0 1 & 2 3, 4, 5, 6
Ru (cm) 0.047 0.0315 0.0215
Rd (cm) 0.047 0.0315 0.0215
lam (unitless) 8.723 13.5847 8.2062
do not change in these simulations for simplicity, but are shown here to illustrate the original functionality of the code to allow for a changing radius along the length of the vessel. Variables given in the configuration hie to solve a simulation are shown in table 2.5, which shows values used for validating the model as well as running the simulations. Given this configuration hie, the model will solve for 4 cardiac cycles as indicated by tc where T is the time period for one cardiac cycle.
Elasticity (y1), resistance (Rt) and capacitance (Ct) varied between test groups with the values shown in Table 2.6. Generally y1 applies to the proximal vessels, and Rt and Ct apply to the distal vessels. These values came from [23] but y1 values were converted to the units for the ID model before going into the configuration hie. Rt and Ct were entered exactly into the configuration hie and then the ID model code converted them into the desired units. Also, instead of using Eq 2.2 to vary
24


Table 2.5: Variables used across all simulations for validation and actual model simulations where rc is the characteristic radius, qc is the characteristic flow, p is blood density, v is the blood viscosity, depth is total number of generations, T is the time period for one cardiac cycle, tc is the number of cardiac cycles to solve for, dt is the time step, and N is the total number of space steps.
Validation Model
rc (unitless) 1 1
qc (unitless) 10 10
p (gjcm,3) 1.06 1.06
v (cm2/sec) 0.046 0.046
depth (unitless) 2 3
T (sec) 0.109 0.109
tc (unitless) 4 4
dt (sec) le-5 le-7
N (unitless) 50 46
25


Table 2.6: Simulation variables used between test groups.
RUr f S D f mmHg*s \ nm mL 'i
Wrf mL )
Control 4.7306e4 58.4 3.4e-3
PH 1.1683e5 133.8 1.178e-3
PH + banded LPA LPA: 1.1683el0, rQ = 0.5 * rQ, lam = 2 * lam Else same as PH Same as PH Same as PH
PH + replacement LPA/RPA: same as control Else same as PH Same as PH Same as PH
elasticity based on radius value, a constant value was used for the entire network just like in [23]. This was implemented by setting k3 equal to the constant elasticity value with ki,k2 = 0.
In order to test the effect of a decreasing radius and pruning that occurs in PH, the 0D circuit model was adjusted such that the radii values were 20% decreased as given by taking the average of the values shown in [36] and comparing the change between the control and PH cases. There is not much data to quantify the amount of pruning that occurs in PH patients, so a conservative 10% decrease in a and f3 was used to construct the tree, which correlated to about a 13.5% decrease in the number of distal vessels.
2.5 Wall Shear Stress and Pulse Pressure
Wall shear stress (WSS) in the distal vessels is determined by multiplying the
viscosity, q, by the wall shear rate, giving the equation for WSS of tw = The
WSS was normalized to the control group due in order to make comparison easier
due to large resulting WSS values. To investigate pulsatility, pulse pressure (PP) is
investigated. PP is calculated by taking the difference between the maximum and
26


minimum pressure. The MATLAB script used to calculate WSS and PP is shown in Appendix E.
27


Figure 2.7: The final tree model used for OD simulations. For more information what each component or block represents, refer to figure 2.6.


3. Results
3.1 Model Validation and Optimization
A single bifurcation was used for validation of both models.
3.1.1 ID Model Validation
The model was validated using previously published data from [23] for a single bifurcation using the values shown for vessels 0, 1, and 2 shown in table 2.3, where vessel 0 is the parent vessel with two daughter vessels, vessels 1 and 2 as shown in figure 3.1.
Figure 3.1: Example tree of a single bifurcation.
The data from [23] was converted to the VaMpy ID model units when necessary. An N value of 50 and a dt value of le-7 were used for the validation simulations. The pressure and flow results for the control case is shown in figure 3.2. The pressure and flow results for the ph case is shown in figure 3.3. The individual plots shown in the figures can be found in Appendix A.
As shown, the peak pressure for the control is 20 mniHg just as in [23]. The PH model has a peak pressure of 40 mniHg which is a bit higher than [23], but this slight difference is likely due to the fact that this validation used an average of the values reported in [23] for RT, Ct, and Eh/ra as it was unclear what they used as “representative” data.
29


Figure 3.2: Control test group results for the validation study done on a single bifurcation. The top plots are pressure (mmHg) and the bottom plots are flow rate (cm3/sec)
30


Figure 3.3: PH test group results for the validation study done on a single bifurcation. The top plots are pressure (mmHg) and the bottom plots are flow rate (cm3/sec)
31


This data also shows that the two flowrate plots add to equal that of Vessel 0 as expected. It also shows that Vessel 1 has a greater flow than Vessel 2. This is due to the fact that Vessel 1 has a larger diameter than Vessel 2 so has less resistance to flow and therefore has a greater flow rate. The ID model shows expected behavior from the pressure and flowrate, so there is good confidence in the model.
3.1.2 OD Model Validation
The OD Model was validated in a few ways. First, a single RC circuit was made and compared to the same circuit solved with ode45 in MATLAB and compared to ensure the model is solving correctly. The results are shown in figure 3.4 and show agreement between techniques.
Figure 3.4: Validation between a single RC circuit in Simulink versus solving a single RC circuit in MATLAB using ODE45.
Second, a single RC circuit was compared to that of a 2-generation RC circuit where both models totaled to the same resistance and capacitance to ensure the model could be expanded into bifurcations with the results shown in figure 3.5. These
32


results in figure 3.5A show good agreement between pressure of the ode45 solving a single RC circuit compared to the total pressure of the 2-generation RC model where the total resistance and capacitance are equal to the single RC circuit. Another validation shown in figure 3.5B is that the pressure of vessels 1 and 2 are equivalent as they should be since they belong to the same generation just as the voltage drop is equivalent between circuits in parallel.
A. Total Pressure
Figure 3.5: Validation for a 2 generation RC circuit. (A) shows the total Pressure of a 2 generation circuit in Simulink versus solving the equivalent single RC circuit in MATLAB using ODE45. (B) shows the pressures for each individual vessel solved in Simulink.
As shown in both figures 3.4 and 3.5, the model takes a few cycles to reach a
33


Table 3.1: The dx values for the largest (dxo) and smallest (dx4) vessel sizes for a
given N value.
N dx0 (cm) dx4 (cm)
10 0.0456 0.00578
20 0.0216 0.00274
30 0.0141 0.00179
40 0.0105 0.00133
50 0.0084 0.00106
60 0.0069 0.00088
70 0.0059 0.00075
80 0.0052 0.00066
90 0.0046 0.00058
100 0.0041 0.00053
500 0.00082 0.00010
stable solution so only the last cycle or last few cycles are considered in subsequent calculations.
3.1.3 Mesh Independence
Mesh Independence studies were done for the ID model to determine the optimum values for dx and dt to run subsequent simulations, while still satisfying the CFL condition. These simulations were run given a tree with an asymmetric geometry as shown in table 2.3. The following N values were tested (with a dt of le-7): 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 500 with the corresponding dx values for each vessel shown in table 3.1. N represents the number of discretized space steps to use to allow for dynamic dx values based on the vessel length. The time it took for each simulation to run is shown in table 3.2.
34


Table 3.2: The total time each simulation ran for a given N value.
N Time (hours)
10 16.86
20 14.76
30 15.55
40 15.15
50 14.97
60 14.49
70 14.32
80 13.84
90 13.95
too 13.66
500 15.68
35


Changes in peak pressure values along the length was investigated for two vessels: the largest vessel (Vessel 0) and the smallest vessel (Vessel 4). The results for the changes in peak pressure along the length of the vessel for each N value is shown in Appendix B. The maximum relative error for each N value, calculated against the largest N value of 500, is shown here in figure 3.6 for Vessel 0 and figure 3.7 for Vessel
4. These figures show that the solution is converging with increasing N values.
16 14 12
10 If
o
-C
8 a)
E i-
6 4 2 0
10 20 30 40 50 60 70 80 90 100 500
N value
Figure 3.6: Maximum Relative Error (compared to N = 500) for the largest vessel, Vessel 0, for varying N / dx values as well as the total time it took for each simulation to run.
36


0.35
0.3
Co"
^ 0.25
o
L±J
g 0.2
5
CD
^ 0.15
E
X n 1
CO 0.1
2
0.05 0
10 20 30 40 50 60 70 80 90 100 500
N value
Figure 3.7: Maximum Relative Error (compared to N = 500) for the smallest vessel, Vessel 4, for varying N / dx values as well as the total time it took for each simulation to run.
Figures 3.6 and 3.7 also show the time it took for each simulation to run, which is also shown in table 3.2. For some of the cases, increasing the N value actually decreased the time it took for some of the simulations to run. This is likely due to Newton’s method requiring less iterations to converge to a solution due to its starting guess from the previous time step being closer than it would be for a smaller N value.
The results of the mesh independence study done with varying N / dx values showed that any N value could be used due to a couple reasons. First, PH is diagnosed within 1 mniHg, and figures B.l and B.2 in Appendix B show that even with the smallest N value, the peak pressure would be within 0.5 mniHg or less. Second, the maximum relative error is also less than 0.5% as shown in figures 3.6 and 3.7 while a pulmonary artery catheter has been shown to have a random error around 5% and a precision error to be within 10-20% [37]. Therefore, the simulations obtain an error less than that of the instrument that provides the inlet boundary condition. A larger N value (aka a smaller dx) is more desirable, especially due to the fact that a larger
37


Table 3.3: The total time each simulation ran for a given dt value.
dt (sec) Time (hours)
le-6 1.366
le-7 13.66
le-8 157.2
N value led to a shorter runtime, and, thus, an N value of 100 was decided on for subsequent simulations.
Next, mesh independence study was performed on varying dt values with N = 100 as chosen from the dx mesh independence study. Due to the CFL condition, the smallest dt value that could be used was le-6 and for the sake of time, only three dt values were studied: le-6, le-7 and le-8. The time it took for each simulation to run is shown in table 3.3. The results for the changes in peak pressure along the length of the vessel for each dt value is shown in Appendix B. The maximum relative error for each dt value, calculated against the largest dt value of le-8, is shown in figure 3.8 for Vessel 0 and figure 3.9 for Vessel 4. These figures show that the solution is converging with increasing dt values. As shown in both figures and in table 3.3, a decreasing dt value leads to an increasing runtime, as expected.


0.03
160
140
120
100
80
60
40
20
1e-8
Figure 3.8: Maximum Relative Error (compared to dt = le-8) for the largest vessel, Vessel 0, for varying dt values as well as the total time it took for each simulation to run.
0.012
160
140
120
100
Figure 3.9: Maximum Relative Error (compared to dt = le-8) for the smallest vessel, Vessel 4, for varying dt values as well as the total time it took for each simulation to run.
Again, we see that the peak pressure is within 0.5 nnnHg (see Appendix B for
39


plots) or less and that the relative error is less than 0.025%; therefore, any dt value could be used. Again, a smaller dt value is desirable, so a dt value of le-7 was chosen in order to keep a small dt with a reasonable runtime as the runtime for the le-8 simulation took 157.2 hours or roughly 6.5 days which is an unreasonable amount of time.
3.2 Simulations
As stated in the Methods section, four test groups were studied: control (representing non-PH), PH, PH with LPA banding (PH+LPAB), and PH with a synthetic/compliant LPA and RPA (PH + replacement).
3.2.1 ID Simulation Time
Due to computational cost of one of the test groups (PH + LPAB), the ID simulations were run using a tree with a symmetric geometry with the radii and length-to-radius values shown in table 2.4. An N value of 46 was used as this was the largest N value that would satisfy the CFL condition for a dt = le-7 to keep the computational cost down. This N value was kept consistent across all simulations. The time it took each simulation to run is shown in table 3.4. Other variables and data used for these simulations can be found in tables in the methods section. Vessel
0 represents the main pulmonary artery (MPA), vessels 1 and 2 are the two daughter vessels of vessel 0 and represent the left pulmonary artery (LPA) and the right pulmonary artery (RPA), respectively. Vessels 3 and 4 are the daughter vessels for vessel
1 and vessels 5 and 6 are the daughter vessels for vessel 2. This is shown graphically in figure 2.2. The plots for pressure and flow show a time of 0.34 to 0.42 due to the time period for each flow cycle and that the simulations ran for 4 cycles in order to obtain a more stable solution.
3.2.2 Pressure
40


Table 3.4: The total time each simulation ran for the four test groups.
Test Group Time (hours)
Control 16.337
PH 14.523
PH + LPAB 36.373
PH + replacement 14.406
The graphs for pressure along the vessel network for the Control, PH, PH+LPAB, PH+replacement are shown in figures 3.10, 3.11, 3.12, and3.13, respectively. Full graphs of the pressure waves can be found in Appendix C. Consistently among all of the vessels, the control group has a peak pressure around 20 mmHg, the PH and PH+LPAB groups tend to have a peak pressure around 40 mmHg, while the PH+replacement has a slightly lower peak pressure of around 32 mmHg. The plot for pressure of the initial vessel of the 0D model (aka the terminal vessel of the ID model) for all test groups is shown in figure 3.14. A plot showing the distal pressure plots for all of the test groups is shown in figure 3.15.
3.2.3 Flow
The graphs for flow along the vessel network for the Control, PH, PH+LPAB, PH+replacement are shown in figures 3.16, 3.17, 3.18, and3.19, respectively. Full graphs of the pressure waves can be found in Appendix D.
A plot for flow of the initial vessel of the 0D model (aka the terminal vessel of the ID model) for all test groups is shown in figure 3.20. A plot showing the distal flow plots for all of the test groups from the 0D simulations is shown in figure 3.21.
3.3 Distal Wall Shear Stress and Pulse Pressure
Wall shear stress (WSS) and pulse pressure (PP) were calculated for just the
terminal vessel since the goal of this project was to investigate the distal mechanical
41


Figure 3.10: Control test group pressure results for the entire vascular tree. Each
vessel at the same generation has the same results due to being completely symmetric,
so only one vessel per generation is shown here for simplicity.
42


Figure 3.11: PH test group pressure results for the entire vascular tree. Each vessel
at the same generation has the same results due to being completely symmetric, so
only one vessel per generation is shown here for simplicity.
43


A.
0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1
Figure 3.12: PH +LPAB test group pressure results for the entire vascular tree. Each
vessel at the same generation has the same results due to being completely symmetric,
so only one vessel per generation is shown here for simplicity.
44


Figure 3.13: PH + replacement test group pressure results for the entire vascular
tree. Each vessel at the same generation has the same results due to being completely
symmetric, so only one vessel per generation is shown here for simplicity.
45


Figure 3.14: Pressure results for the initial vessel of all of the test groups of the OD model simulations.
Figure 3.15: Pressure results for the terminal vessels of all of the test groups of the OD model simulations.
46


Figure 3.16: Control test group pressure results for the entire vascular tree. Each
vessel at the same generation has the same results due to being completely symmetric,
so only one vessel per generation is shown here for simplicity.
47


2---------â– ---------------â– -------â– -------
0 0.2 0.4 0.6 0.8 1
_________________Time (sec)________________
Figure 3.17: PH test group pressure results for the entire vascular tree. Each vessel
at the same generation has the same results due to being completely symmetric, so
only one vessel per generation is shown here for simplicity.
48


Figure 3.18: PH+LPAB test group pressure results for the entire vascular tree. Each
vessel at the same generation has the same results due to being completely symmetric,
so only one vessel per generation is shown here for simplicity.
49


Figure 3.19: PH + replacement test group pressure results for the entire vascular
tree. Each vessel at the same generation has the same results due to being completely
symmetric, so only one vessel per generation is shown here for simplicity.
50


Figure 3.20: Flow results for the initial vessel of all of the test groups of the 0D model simulations.
Figure 3.21: Flow results for the initial vessel of all of the test groups of the 0D model simulations.
51


Test Group
Figure 3.22: WSS results for the OD model simulations for a terminal vessel given changes to proximal vessel stiffness.
changes.
3.3.1 Effect of proximal stiffening on distal WSS and PP
A bar plot of the WSS, normalized to the control WSS, for each test group is shown in figure 3.22 given changes to proximal vessel stiffness. Only data from the last cycle or two was included in the calculations for WSS and PP to ensure accurate data. A bar plot of the PP is shown in figure 3.23 for all of the test groups given changes to proximal vessel stiffness.
3.3.2 Differences in vessel geometry impact distal WSS and PP
A bar plot of the WSS for each test group is shown in figure 3.24 once changes due to vessel geometry were included. A bar plot of the PP is shown in figure 3.23 for all of the test groups, which includes changes due to vessel geometry.
52


Test Group
Figure 3.23: Pulse pressure results for the OD model simulations for a terminal vessel given changes to vessel stiffness.
1.8
1.6
1.4
E
o
1.2
1
03
03
§
0.8
0.6
0.4
0.2
0
Figure 3.24: WSS results for the OD model simulations for a terminal vessel given both changes to proximal vessel stiffness, and changes to distal vessel geometry.
53


Test Group
Figure 3.25: Pulse pressure results for the OD model simulations for a terminal vessel given both changes to proximal vessel stiffness, and changes to distal vessel geometry.
54


4. Discussion
4.1 Pressure Waves
In the ID model, the pressure changed mostly as expected. The PH cases had higher peak pressure than the control, despite the PH having lower flow rate. There was also a slight increase in peak pressure for the PH + LPAB (45 mmHg) compared to the PH case in Vessel 0, but this makes sense because banding caused a decreased radius and increased stiffness in vessel 1 (LPA) which leads to increased resistance and therefore increased peak pressure. The pressure plot in vessel 1 of the PH + LPAB showed pressure drop off along the length of that vessel, while vessel 2 stayed at the higher peak pressure of 45 mmHg. This is expected as less flow is expected to go through the smaller vessel size. Vessels downstream of the LPA in the PH + LPAB showed a continued decrease in pressure and showed equivalent pressures to that of the PH + replacement case. For the PH + replacement case, there was about a 20% decrease in peak pressure for vessels 1 and 2. The peak pressure did not return to control values; however, the PH case did have a slightly higher peak pressure than expected. If the peak pressure had been closer to what was shown in [23], which was approximately 35 mmHg, then a 20% decrease would have put the peak pressure to 28 mmHg. While this peak pressure is closer to control values, it would still be classified as PH; however, over time this value could decrease even more depending on changes in anatomy and physiology of the distal vasculature due to a decreased pressure.
One interesting difference in the pressure wave of the PH + replacement case is that the end-cycle pressure was increased compared to the PH despite having a lower peak pressure. The lowered end-cycle pressure also contribute to the PH + replacement having a lower pulse intensity. As in the difference between maximum and minimum pressure is only 10 mmHg which is close to the control case (12 mmHg) whereas the PH and PH + LPAB are all around 20 mmHg. This trend continues as
55


you move into the distal vasculature of the OD model and can be seen in the pressure waves in figure 3.15 as well as the pulse pressure measurements shown in figure 3.23. The pulse pressure for the PH + replacement; however, does not decrease to as much as the control case in the distal arterioles which is due to the distal arterioles being much more compliant in the control case allowing for more attenuation.
4.2 Flow Waves
The flow also changed as expected in the ID model. For the PH + LPAB, figure 3.18 shows vessel 1 (LPA) having decreased flow while vessel 2 (RPA) has increased flow, which is what was seen in [27]. These changes give good confidence to the model giving results that are comparable to an animal model. The trends in flow for this case is perpetuated to the vessels downstream from vessels 1 and 2 as shown in subsequent vessel plots. Another trend noticed in vessel 1 compared to vessel 2 in the PH + LPAB case, is that the left daughter vessels show a smoother flow wave than the left daughter vessels as shown in figures 3.20 and 3.21.
In the PH + replacement case, the flow wave seems to become less pulsatile and more sustained. Figure 3.19 shows the pulse of the second generation of vessels diminishing as it travels down the length of the vessel. The ID terminal vessels show a more sustained and less pulsatile flow over time while the wave seems more consistent along the length of the vessels. This trend continues for the distal vessels of the OD model and can be seen in the flow waves shown in figure 3.21 by the elongated flow wave compared to the control and PH cases.
4.3 Distal Wall Shear Stress and Pulse Pressure
Wall shear stress (WSS), shown in figure 3.22, in the PH case decreased in response to increased vessel stiffening; however, this was seen without any changes to geometry as seen with PH. With the decreased radius and pruning being modeled, shown in figure 3.24, WSS in the PH case was increased as expected and gives evidence that the increase in WSS is purely due to geometric changes. With the radius
56


change and pruning being modeled, the WSS for the PH + replacement case returned closer to control levels. The banding of the LPA also decreased WSS in the left lung arteries whereas the WSS for the right lung arteries was slightly increased compared to PH. Due to these changes and the results shown in [27], PP and WSS are likely the cause of the anatomical changes seen in [27].
The pulse pressure (PP), shown in figure 3.23, shows expected trends. Just accounting for changes to vessel stiffness, the PP for the PH case is increased to about four times the control. This increase was seen in the ID model pulse waves, but wasn’t as stark of a difference. Again, due to the downstream vessels being less compliant in the PH case, those pressure waves have less attenuation which is why the PP would be much higher. The PH + LPAB (L) and (R) cases both show increases to PP, but the left side showed a slight decrease compared to PH while the right side had a slight increase compared to PH. This is expected as there was less flow going into the left side, but the difference is very minimal compared to the PP of the PH case. The PH + replacement shows decreased PP back towards the control levels, but is still twice that of the control. Perhaps updating the PH + replacement compliance to be more elastic than the control could further decrease this pulse pressure closer to normal. When pruning and decreased radius values were included, these trends were further exacerbated as shown in figure 3.25.
4.4 Limitations and Future Work
This model used the same geometry for proximal vessels between all groups except for the case of the PH+LPAB. This model also only tested for a single, generalized mouse geometry and a great extension of this model would be to do animal studies using animal-specific data for each simulation. With an animal study, this model could also be used in conjunction with histological data to better connect the fluid dynamics with anatomical changes to the distal blood vessels over time.
Other studies that could be done using this simulation, would be to test varying
57


stiffnesses of replacing the LPA and RPA to determine if there might be a more optimal stiffness. In this model, the same stiffness as control vessels was tested; however, the changes to PP and WSS were not as close to control levels as desired which may mean the stiffness of the replacement vessel would need to be adjusted.
This model also only investigated one subset of the lung and did simulations using a tree with symmetric vessels for simplicity; however, it would be interesting to look at a tree with an asymmetric geometry to investigate what changes occur and if this is a viable simplification for future model use.
58


5. Conclusion
Pulmonary hypertension (PH) is associated with changes to anatomy and physiology of blood vessels [4]. More specifically, PH vessels increase in stiffness, decrease in radius and therefore, increase in vascular resistance. These changes cause an increase in flow pulsatility and wall shear stress (WSS) which further exacerbate these problems [7], and therefore, pulsatility and WSS are important features to study. Unfortunately, these features cannot be studied in distal vessels due to equipment constraints [9], which is where computational models come into play. The goals of this project were to develop a computational model of a large-scale pulmonary vasculature network as well as to investigate pulsatility, using pulse pressure (PP), and WSS in distal vasculature given changes to the proximal pulmonary vasculature. More specifically, simulation of a banded left pulmonary artery (LPAB) in the presence of PH and simulation of a replacement left and right pulmonary artery in the presence of PH were compared to healthy and PH controls.
The computational model was developed using a ID and OD hybrid model using an existing python-based ID model as the base. The model was expanded upon to allow for investigation of the desired test groups described, but had a high computational cost. Thus, a OD model was added to create the hybrid model where the OD model used Simulink to create a circuit-based analog of the distal vasculature. The model was then used to investigate PP and WSS in the distal vessels. PP was increased for PH compared to the healthy control, as expected. The PH + LPAB (L) showed a slight decrease compared to PH while the PH + LPAB (R) showed an increased PP beyond that of the PH case. The left and right pulmonary artery replacement group returned PP closer to control levels, which showed good promise for a synthetic left and right pulmonary artery as a treatment option for PH. The WSS for PH was also increased compared to the healthy control, but only when pruning and decreased vessel area were modeled providing evidence that changes to WSS are
59


completely due to geometry and not vessel mechanics. In the PH + LPAB case, the downstream arteries of the banded artery showed a decrease in WSS compared to PH while the arteries downstream from the non-banded artery showed an increase in WSS beyond that of the PH case. PH + LPAB was shown to decrease inflammatory markers of PH [27] and this model showed that the PH + replacement group showed similar WSS and PP trends as the PH + LPAB case simulated. These results suggest that replacing the left and right pulmonary arteries with synthetic, compliant vessels could be a viable treatment option and should be further investigated. One way this could be investigated, is by testing this treatment option in an animal study and running animal specific simulations simultaneously. Despite many simplifications, this model was able to provide accurate pressure and flow results in pulmonary vessels and is a great tool for future application and investigation. In conclusion, the model was successfully created and used to investigate PP and WSS in distal vessels, and illustrated the importance of mechanical stimuli in the vascular remodeling process of progressive pulmonary vascular disease.
60


REFERENCES
[1] M. Rabinovitch. Molecular pathogenesis of pulmonary arterial hypertension. The Journal of Clinical Investigation, 122(12):4306 4313, 12 2012.
[2] G. Simonneau, M.A. Gatzoulis, I. Adatia, D. Celermajer, C. Denton,
A. Ghofrani, M.A. Gomez Sanchez, K.R. Krishna, M. Landzberg, R.F. Machado, H. Olschewski, I.M. Robbins, and R. Souza. Updated clinical classification of pulmonary hypertension. Journal of the American College of Cardiology, 62:D34-41, 2013.
[3] W. Sun and S.Y. Chan. Pulmonary arterial stiffness: An early and pervasive driver of pulmonary arterial hypertension. Frontiers in medicine, 5:204-204, 2018.
[4] K.R. Stenmark, N. Davie, M. Frid, E. Gerasimovskaya, and M. Das. Role of the adventitia in pulmonary vascular remodeling. Physiology, 21(2):134 145, 2006.
[5] K.S. Hunter, S.R. Lammers, and R. Shandas. Pulmonary vascular stiffness: measurement, modeling, and implications in normal and hypertensive pulmonary circulations. Comprehensive Physiology, 1(3): 1413-1435, 2011.
[6] S. Malenfant, A. Neyron, R. Paulin, F. Potus, J. Meloche, S. Provencher, and
S. Bonnet. Signal transduction in the development of pulmonary arterial hypertension. Pulmonary circulation, 3(2):2T8—293, 2013.
[7] W. Tan, K. Madhavan, K.S. Hunter, D. Park, and K.R. Stenmark. Vascular stiffening in pulmonary hypertension: cause or consequence? (2013 grover conference series). Pulmonary circulation, 4(4):560-580, 2014.
[8] J.C. Grignola. Hemodynamic assessment of pulmonary hypertension. World journal of cardiology, 3(1): 10—IT, 2011.
[9] J. Alastruey, N Xiao, H. Fok, To. Schaeffter, and C.A. Figueroa. On the impact of modelling assumptions in multi-scale, subject-specific models of aortic haemodynamics. Journal of the Royal Society, Interface, 13(119) :20160073, 2016.
[10] K. Hassani, M. Navidbakhsh, and M. Rostami. Modeling of the aorta artery aneurysms and renal artery stenosis using cardiovascular electronic system. Biomedical engineering online, 6:22-22, 2007.
[11] O. Ghasemalizadeh, M.R. Mirzaee, B. Firoozabadi, and K. Hassani. Exact modeling of cardiovascular system using lumped method. International Conference on Bioinformatics and Computational Biology, pages 408-417, 2008.
61


[12] M.S. Olufsen. Structure tree outflow condition for blood flow in larger systemic arteries. American Journal of Physiology, 276(l):H257-68, 1999.
[13] M.S. Olufsen. A one-dimensional fluid dynamic model of the systemic arteries. Studies in Health Technologies and Informatics, 71:79-98, 2000.
[14] G.D.A. Qureshi, M.U.and Vaughan, C. Sainsbury, M. Johnson, C.S. Peskin, M.S. Olufsen, and N.A. Hill. Numerical simulation of blood flow and pressure drop in the pulmonary arterial and venous circulation. Biomechanics and modeling in rnechanobiology, 13(5): 1137—1154, 2014.
[15] J. Flores, J. Alastruey, and E. Corvera Poire. A novel analytical approach to pulsatile blood flow in the arterial network. Annals of Biomedical Engineering, 44(10):3047-3068, 2016.
[16] J.P. Mynard and J.J. Smolich. One-dimensional haemodynamic modeling and wave dynamics in the entire adult circulation. Annals of Biomedical Engineering, 43(6):1443-1460, 2015.
[17] P. Reymond, F. Merenda, F. Perren, D. Rufenacht, and N. Stergiopulos. Validation of a one-dimensional model of the systemic arterial tree. American Journal of Physiology-Heart and Circulatory Physiology, 297(1):H208-H222, 2009.
[18] Jonasova A., Bublik O., and Vimmr J. A comparative study of Id and 3d hemodynamics in patient-specific hepatic portal vein networks. Applied and Computational Mechanics, 8(2): 177—186, 2014.
[19] L. Grinberg, E. Cheever, T. Anor, J.R. Madsen, and G.E. Karniadakis. Modeling blood flow circulation in intracranial arterial networks: A comparative 3d/ld simulation study. Annals of Biomedical Engineering, 39(1):297 309, 2011.
[20] P. Reymond, P. Crosetto, S. Deparis, A. Quarteroni, and N. Stergiopulos. Physiological simulation of blood flow in the aorta: Comparison of hemodynamic indices as predicted by 3-d fsi, 3-d rigid wall and 1-d models. Medical Engineering & Physics, 35(6):784-791, 2013.
[21] N. Xiao, J. Alastruey, and C. Alberto Figueroa. A systematic comparison between 1-d and 3-d hemodynamics in compliant arterial models. International journal for numerical methods in biomedical engineering, 30(2):204—231, 2014.
[22] A. Quarteroni, A. Veneziani, and C. Vergara. Geometric multiscale modeling of the cardiovascular system, between theory and practice. Computer Methods in Applied Mechanics and Engineering, 302:193-252, 2016.
[23] M.U. Qureshi, M.J. Colebank, L.M. Paun, L.E. Fix, N. Chester, M.A. Haider, N.A. Hill, D. Husmeier, and M.S. Olufsen. Hemodynamic assessment of pulmonary hypertension in mice: a model-based analysis of the disease mechanism. Biomechanics and Modeling in Mechanobiology, 18(1):219—243, 2019.
62


[24] N. Bessonov, A. Sequeira, S. Simakov, Y. Vassilevskii, and V. Volpert. Methods of blood flow modeling. Mathematical Modeling of Natural Phenomena, 11(1): 1 25, 2016.
[25] Y. Shi, P. Lawford, and R. Hose. Review of zero-d and 1-d models of blood flow in the cardiovascular system. Biomedical engineering online, 10:33, 2011.
[26] E. Boileau, P. Nithiarasu, P.J. Blanco, L.O. Miiller, F.E. Fossan, L.R. Hellevik, W.P. Donders, W. Huberts, M. Willemet, and J. Alastruey. A benchmark study of numerical schemes for one-dimensional arterial blood flow modelling. International Journal for Numerical Methods in Biomedical Engineering, 31(10):e02732, 2015.
[27] K. Abe, M. Shinoda, M. Tanaka, Y. Kuwabara, K. Yoshida, Y. Hirooka, I.F. McMurtry, M. Oka, and K. Sunagawa. Hemodynamic unloading reverses occlusive vascular lesions in severe pulmonary hypertension. Cardiouascular Research, 1(111): 16 25, 2016.
[28] B.T. Tang, S.S. Pickard, F.P. Chan, P.S. Tsao, C.A. Taylor, and J.A. Feinstein. Wall shear stress is decreased in the pulmonary arteries of patients with pulmonary arterial hypertension: An image-based, computational fluid dynamics study. Pulmonary circulation, 2(4):470-476, 2012.
[29] M. Schafer, V.O. Kheyfets, J.D. Schroeder, J. Dunning, R. Shandas, J.K. Buckner, J. Browning, K.S. Hertzberg, J.and Hunter, and B.E. Fenster. Main pulmonary arterial wall shear stress correlates with invasive hemodynamics and stiffness in pulmonary hypertension. Pulmonary circulation, 6(1):37 45, 2016.
[30] A.K. Diem and N.W. Bresslof. Vampy: A python package to solve Id blood flow problems. Journal Of Open Research Software, 5(17), 2017.
[31] V.B Kolachalama, N.W. Bresslof, P.B. Nair, et al. Predictive haemodynamics in a one-dimensional carotid artery bifurcation, part i: Application to stent design. IEEE Transactions on Biomedical Engineering, 54(5):802-812, 2007.
[32] M.S. Olufsen, C.S. Peskin, W.Y. Kim, E.M. Pedersen, A. Nadim, and J. Larsen. Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions. Annals of Biomedical Enqineerinq, 28:1281 -1299, 2000.
[33] A. Mishra, F.M. O’Farrell, C. Reynell, N.B. Hamilton, C.N. Hall, and D. Attwell. Imaging pericytes and capillary diameter in brain slices and isolated retinae. Nature Protocols, 9:323 EP, 2014.
[34] M.P. Wiedeman. Dimensions of blood vessels from distributing artery to collecting vein. Circulation Research, pages 375-378, 1963.
[35] M.I. Townsley. Structure and composition of pulmonary arteries, capillaries, and veins. Comprehensiue Physiology, 2(l):675-709, 2012.
63


[36] N. Rol, E.M. Timmer, T.J.C. Faes, A. Vonk Noordegraaf, K. Grunberg, H. Bo-gaard, and N. Westerhof. Vascular narrowing in pulmonary arterial hypertension is heterogeneous: rethinking resistance. Physiological reports, 5(6):el3159, 2017.
[37] X.X. Yang, L.A. Critchley, and G.M. Joynt. Determination of the precision error of the pulmonary artery thermodilution catheter using an in vitro continuous flow test rig. Anesthesia and Analgesia, 112(1):70 77, 2011.
64


APPENDIX A. Validation Plots for ID Model
A.l Control
A.1.1 Pressure Plots
Figure A.l: Validation Pressure Results of the Control test group for the root vessel 0 (MPA)
65


Figure A.2: Validation Pressure Results of the Control test group for vessel 1 (LPA)
20
18
16
14
12
10
8
6
â– 
r10
r8 l 6
Figure A.3: Validation Pressure Results of the Control test group for vessel 2 (RPA)
66


A. 1.2 Flow Plots
Figure A.4: Validation Flow Results of the Control test group for the root vessel 0
(MPA)
67


- 0.6
Figure A.5: Validation Flow Results of the Control test group for vessel f (LPA)
Figure A.6: Validation Flow Results of the Control test group for vessel 2 (RPA)
68


A.2 PH
A.2.1 Pressure Plots
Figure A.7: Validation Pressure Results of the PH test group for the root vessel 0 (MPA)
69


Figure A.8: Validation Pressure Results of the PH test group for vessel 1 (LPA)
Ol
x.
£
£
a;
ZI
cn
co
0)
Cl
40
35
30
25
20
15
Figure A.9: Validation Pressure Results of the PH test group for vessel 2 (RPA)
70


A.2.2 Flow Plots
Figure A. 10: Validation Flow Results of the PH test group for the root vessel 0 (MPA)
71


Figure A.ll: Validation Flow Results of the PH test group for vessel 1 (LPA)
Figure A. 12: Validation Flow Results of the PH test group for vessel 2 (RPA)
72


APPENDIX B. Mesh Independence Plots
B.l Changes to Peak Pressure for varying N / dx
Figure B.l: Changes to Peak Pressure along the length of the largest vessel, Vessel 0, due to varying N / dx values.
73


21.02
o N = 10
o N = 20 -
o N = 30
o N = 40
o N = 50
o N = 60
o N = 70
o N = 80
o N = 90 -
o N = 100
o N = 500 -
0.01 0.02 0.03 0.04 0.05
x distance along the length of the vessel (cm)
0.06
Figure B.2: Changes to Peak Pressure along the length of the smallest vessel, Vessel 4, due to varying N / dx values.
B.2 Changes to Peak Pressure for varying dt
Figure B.3: Changes to Peak Pressure along the length of the largest vessel, Vessel 0, due to varying dt values.
74


Figure B.4: Changes to Peak Pressure along the length of the smallest vessel, Vessel 4, due to varying dt values.
75


APPENDIX C. Model Pressure Plots
C.l Control
Figure C.l: Pressure Results of the Control test group for the root vessel (MPA)
76


Figure C.2: Pressure Results of the Control test group for the second generation of vessels (LPA/RPA)
77


Figure C.3: Pressure Results of the Control test group for the ID terminal vessels (Generation 3).
78


Figure C.4: Pressure Results of the Control test group for the OD terminal vessels (Generation 12).
C.2 PH
79


Figure C.5: Pressure Results of the PH test group for the root vessel (MPA)
40
35
30
25
20
15
80


Figure C.6: Pressure Results of the PH test group for the second generation of vessels (LPA/RPA).
81


82


15
Figure C.8: Pressure Results of the PH test group for the OD terminal vessels (Generation 12).
C.3 PH + LPAB
83


Figure C.9: Pressure Results of the PH+LPAB test group for the root vessel (MPA).
84


Full Text

PAGE 1

INVESTIGATINGTHEROLEOFDISTALFLOWINPROGRESSIVE PULMONARYVASCULARDISEASEUSINGACOMPUTATIONALMODEL. by RACHELLEL.WALTER B.S.,UniversityofColoradoDenver,2017 Athesissubmittedtothe FacultyoftheGraduateSchoolofthe UniversityofColoradoinpartialfulllment oftherequirementsforthedegreeof MasterofScience Bioengineering 2019

PAGE 2

ThisthesisfortheMasterofSciencedegreeby RachelleL.Walter hasbeenapprovedforthe BioengineeringProgram by VitalyO.Kheyfets,Chair KurtR.Stenmark KendallS.Hunter June21,2019 ii

PAGE 3

Walter,RachelleL.M.S.,Bioengineering Investigatingtheroleofdistalowinprogressivepulmonaryvasculardiseaseusinga computationalmodel. ThesisdirectedbyVitalyO.Kheyfets ABSTRACT PulmonaryArterialHypertensionPAH,whichisassociatedwithpulmonary arterialstieningandthickeningofdistalpulmonaryvesselwalls,leadstorightheart failureandsometimesdeath.Increasedthicknessofthedistalvesselwallsprimarilyoccursintheadventitiaandisthoughttobeduetoincreasedmechanicalstress. Toprovideevidenceforthisrelationship,mechanicalstresswassimulatedforthe entirepulmonaryvascularnetworkusingacombined1D/0Dcomputationalmodel, whichappliesaphysiologicalowwaveformattheinletanda2or3-elementoutow boundaryconditionattheoutlets.Thecomputationalmodelwasusedtoevaluatehowchangesinproximalvascularstinessimpactpulsepressureandwallshear stressinthedistalvascularbedinthecontextofpulmonaryhypertensionPH. Thefoursimulatedtestgroupswerecontrol,PH,PHwithabandedleftpulmonary arterywithdecreasedradiusandincreasedstiness,andPHwiththeleftand rightpulmonaryarteriesreplacedwithcompliantarteriesdecreasedstinessback tocontrollevels.ThesimulatedresultsshowedincreasedpulsepressureforPH .45mmHgcomparedtohealthycontrol.39mmHg.ThePH+bandedleft pulmonaryarteryshowedaslightdecreaseinpulsepressurecomparedtoPHinthe leftarterydescendants.17mmHgandanincreasecomparedtoPHintheright arterydescendants.12mmHgarterydescendants.ThePH+replacementleftand rightpulmonaryarterieshadapulsepressureinbetweenthecontrolandPHlevels .38mmHg.Whenonlyvesselstieningwasmodeled,theWSSnormalizedtocontrol.0showedadecreaseforPH.79,furtherdecreasedforarteriesdownstream iii

PAGE 4

oftheLPAinthebandedLPAcase.67.WSSinthearteriesdownstreamofthe RPAinthebandedLPAcase.95wereaboutequivalenttothatofcontrols,and theWSSinPHwithreplacementarteriescase.62showedadecrease.Oncepruninganddecreasedvesselareaweremodeled,theWSSnormalizedtocontrol.0 showedanincreaseforPH.5,andafurtherincreaseforarteriesdownstreamof theRPAinthebandedLPAcase.80.WSSinthearteriesdownstreamoftheLPA inthebandedLPAcase.28andtheWSSinPHwithreplacementarteriescase .18werelessthanthatofthePHcase,butstillslightlyabovethecontrolWSS level.Theseresultssuggestthatareplacementleftandrightpulmonaryarterymight besucienttoreturnthepulsepressureandWSStoahealthylevelgivenfurther compliancetotheproximalvessels.Furtherstudiesshouldbedonetodirectlyinvestigatetherelationshipbetweenmechanicalstressandvesselremodelingusingan animalmodel.Determiningwhetherornotmechanicalstressiscausingpulmonary vesselremodelingcouldleadtotheuseofmechanicalstressasapotentialtherapeutic targettomanageorevenreversePAHdiseaseprogression. Theformandcontentofthisabstractareapproved.Irecommenditspublication. Approved:VitalyO.Kheyfets iv

PAGE 5

ACKNOWLEDGMENT Iwouldliketothank: MyAdvisors:Dr.VitalyKheyfets,Dr.KendallHunter,andDr.KurtStenmark Myfamilyandfriendsforprovidingsupport Dr.MichaelYeager,forhismentorshipbeforeembarkingonthisproject MyfellowBioengineeringStudents OtherBioengineeringfaculty/stawhohaveprovidedinvaluablementorshiptome v

PAGE 6

TABLEOFCONTENTS Tables........................................ix Figures.......................................x Chapter 1.Introduction...................................1 1.1PulmonaryHypertension........................1 1.1.1InterplayofMechanicsandBiology...............2 1.2ComputationalModeling........................2 1.3Justication...............................5 1.4Goals...................................6 2.Methods.....................................8 2.1ModelOverview.............................8 2.1.1Assumptions...........................9 2.21DModelVaMpy...........................9 2.2.1Equations.............................9 2.2.2BoundaryConditions.......................11 2.2.3CodeAdaptation.........................13 2.2.4Runninga1DSimulation....................16 2.2.5Limitations............................16 2.30DModel.................................16 2.3.1VesselNetworkGeneration...................17 2.3.2SimulinkModel..........................19 2.3.3BoundaryConditions.......................21 2.3.4Runninga0DSimulation....................21 2.3.50DSimulationData.......................22 2.3.6Limitations............................23 2.4NumericalSimulationsforCombinedModels.............23 vi

PAGE 7

2.5WallShearStressandPulsePressure.................26 3.Results......................................29 3.1ModelValidationandOptimization..................29 3.1.11DModelValidation.......................29 3.1.20DModelValidation.......................32 3.1.3MeshIndependence.......................34 3.2Simulations................................40 3.2.11DSimulationTime.......................40 3.2.2Pressure..............................41 3.2.3Flow................................41 3.3DistalWallShearStressandPulsePressure..............41 3.3.1EectofproximalstieningondistalWSSandPP......52 3.3.2DierencesinvesselgeometryimpactdistalWSSandPP..52 4.Discussion....................................55 4.1PressureWaves.............................55 4.2FlowWaves...............................56 4.3DistalWallShearStressandPulsePressure..............56 4.4LimitationsandFutureWork......................57 5.Conclusion....................................59 References ......................................61 Appendix A.ValidationPlotsfor1DModel.........................65 A.1Control..................................65 A.1.1PressurePlots..........................65 A.1.2FlowPlots............................67 A.2PH....................................69 A.2.1PressurePlots..........................69 vii

PAGE 8

A.2.2FlowPlots............................71 B.MeshIndependencePlots............................73 B.1ChangestoPeakPressureforvaryingN/dx.............73 B.2ChangestoPeakPressureforvaryingdt................74 C.ModelPressurePlots..............................76 C.1Control..................................76 C.2PH....................................80 C.3PH+LPAB...............................84 C.4PH+replacement............................91 D.ModelFlowPlots................................95 D.1Control..................................95 D.2PH....................................97 D.3PH+LPAB...............................99 D.4PH+replacement............................104 E.MATLABCode.................................106 E.1VesselNetworkGenerationScript...................106 E.20DModelScript.............................109 E.3WSSandPPScript...........................113 viii

PAGE 9

TABLES Table 2.1Unitscomparisonbetweenthehydraulicorbloodowmodelandthe electricorcircuitmodel............................19 2.2Resistance, R T ,andCapacitance, C T ,attheoutletofthe1Dmodel.The scalingfactorisappliedtothePoiseuillelawasshowninequation2.9in the0Dmodelwhenbuildingupthetreesuchthatthetotalresistanceand totalcapacitanceinthe0Dmodelareequalto R T and C T ,respectively.22 2.3Upstreamradii, R u ,downstreamradii, R d andlength-to-radiusratio, lam , valuesforeachvesselforatreewithasymmetricvesselgeometry.....24 2.4Upstreamradii, R u ,downstreamradii, R d andlength-to-radiusratio, lam , valuesforeachvesselforatreewithsymmetricvesselgeometry......24 2.5Variablesusedacrossallsimulationsforvalidationandactualmodelsimulationswhere r c isthecharacteristicradius, q c isthecharacteristicow, isblooddensity, isthebloodviscosity, depth istotalnumberofgenerations, T isthetimeperiodforonecardiaccycle, t c isthenumberof cardiaccyclestosolvefor, dt isthetimestep,and N isthetotalnumber ofspacesteps..................................25 2.6Simulationvariablesusedbetweentestgroups................26 3.1The dx valuesforthelargest dx 0 andsmallest dx 4 vesselsizesfora givenNvalue..................................34 3.2ThetotaltimeeachsimulationranforagivenNvalue...........35 3.3Thetotaltimeeachsimulationranforagivendtvalue...........38 3.4Thetotaltimeeachsimulationranforthefourtestgroups........41 ix

PAGE 10

FIGURES Figure 1.1Graphicrepresentationof3D,1Dand0Dmodelsandthetrade-obetween computationalcostandphysiologicalrelevancei.e.howrepresentative themodelis..................................4 2.1High-LevelOverviewoftheFullComputationalModel..........8 2.2Anexampleofatreewith3generationswherevessel0representsthe mainpulmonaryarteryMPA,vessel1representstheleftpulmonary arteryLPA,vessel2representstherightpulmonaryarteryRPAand vessels3-6aretheterminalvessels.Theoutowboundaryconditionisa 3-elementwindkesselappliedateachterminalvessel............10 2.3InletowwaveformforcontrolandpulmonaryhypertensionPHtest groups.....................................12 2.4Anexampleofatreewithasymmetricbranchingwherethereareterminal branchesatgeneration2andatgeneration3................17 2.5Anexampleofa2-generationRCcircuittreemodelmadeinSimulink showingtheinletowandoutowboundaryconditionBC........19 2.6Anexampleofa2-generationRCcircuittreemodelmadeinSimulink withVoltmeters,AmpMeters,andscopeblocks..............20 2.7Thenaltreemodelusedfor0Dsimulations.Formoreinformationon whateachcomponentorblockrepresents,refertogure2.6........28 3.1Exampletreeofasinglebifurcation.....................29 3.2Controltestgroupresultsforthevalidationstudydoneonasinglebifurcation.ThetopplotsarepressuremmHgandthebottomplotsareow rate cm 3 /sec................................30 x

PAGE 11

3.3PHtestgroupresultsforthevalidationstudydoneonasinglebifurcation. ThetopplotsarepressuremmHgandthebottomplotsareowrate cm 3 /sec...................................31 3.4ValidationbetweenasingleRCcircuitinSimulinkversussolvingasingle RCcircuitinMATLABusingODE45....................32 3.5Validationfora2generationRCcircuit.AshowsthetotalPressure ofa2generationcircuitinSimulinkversussolvingtheequivalentsingle RCcircuitinMATLABusingODE45.Bshowsthepressuresforeach individualvesselsolvedinSimulink.....................33 3.6MaximumRelativeErrorcomparedtoN=500forthelargestvessel, Vessel0,forvaryingN/dxvaluesaswellasthetotaltimeittookfor eachsimulationtorun.............................36 3.7MaximumRelativeErrorcomparedtoN=500forthesmallestvessel, Vessel4,forvaryingN/dxvaluesaswellasthetotaltimeittookfor eachsimulationtorun.............................37 3.8MaximumRelativeErrorcomparedtodt=1e-8forthelargestvessel, Vessel0,forvaryingdtvaluesaswellasthetotaltimeittookforeach simulationtorun...............................39 3.9MaximumRelativeErrorcomparedtodt=1e-8forthesmallestvessel, Vessel4,forvaryingdtvaluesaswellasthetotaltimeittookforeach simulationtorun...............................39 3.10Controltestgrouppressureresultsfortheentirevasculartree.Each vesselatthesamegenerationhasthesameresultsduetobeingcompletely symmetric,soonlyonevesselpergenerationisshownhereforsimplicity.42 3.11PHtestgrouppressureresultsfortheentirevasculartree.Eachvesselat thesamegenerationhasthesameresultsduetobeingcompletelysymmetric,soonlyonevesselpergenerationisshownhereforsimplicity...43 xi

PAGE 12

3.12PH+LPABtestgrouppressureresultsfortheentirevasculartree.Each vesselatthesamegenerationhasthesameresultsduetobeingcompletely symmetric,soonlyonevesselpergenerationisshownhereforsimplicity.44 3.13PH+replacementtestgrouppressureresultsfortheentirevasculartree. Eachvesselatthesamegenerationhasthesameresultsduetobeing completelysymmetric,soonlyonevesselpergenerationisshownherefor simplicity....................................45 3.14Pressureresultsfortheinitialvesselofallofthetestgroupsofthe0D modelsimulations...............................46 3.15Pressureresultsfortheterminalvesselsofallofthetestgroupsofthe0D modelsimulations...............................46 3.16Controltestgrouppressureresultsfortheentirevasculartree.Each vesselatthesamegenerationhasthesameresultsduetobeingcompletely symmetric,soonlyonevesselpergenerationisshownhereforsimplicity.47 3.17PHtestgrouppressureresultsfortheentirevasculartree.Eachvesselat thesamegenerationhasthesameresultsduetobeingcompletelysymmetric,soonlyonevesselpergenerationisshownhereforsimplicity...48 3.18PH+LPABtestgrouppressureresultsfortheentirevasculartree.Each vesselatthesamegenerationhasthesameresultsduetobeingcompletely symmetric,soonlyonevesselpergenerationisshownhereforsimplicity.49 3.19PH+replacementtestgrouppressureresultsfortheentirevasculartree. Eachvesselatthesamegenerationhasthesameresultsduetobeing completelysymmetric,soonlyonevesselpergenerationisshownherefor simplicity....................................50 3.20Flowresultsfortheinitialvesselofallofthetestgroupsofthe0Dmodel simulations...................................51 xii

PAGE 13

3.21Flowresultsfortheinitialvesselofallofthetestgroupsofthe0Dmodel simulations...................................51 3.22WSSresultsforthe0Dmodelsimulationsforaterminalvesselgiven changestoproximalvesselstiness......................52 3.23Pulsepressureresultsforthe0Dmodelsimulationsforaterminalvessel givenchangestovesselstiness........................53 3.24WSSresultsforthe0Dmodelsimulationsforaterminalvesselgivenboth changestoproximalvesselstiness,andchangestodistalvesselgeometry.53 3.25Pulsepressureresultsforthe0Dmodelsimulationsforaterminalvessel givenbothchangestoproximalvesselstiness,andchangestodistalvessel geometry....................................54 A.1ValidationPressureResultsoftheControltestgroupfortherootvessel 0MPA....................................65 A.2ValidationPressureResultsoftheControltestgroupforvessel1LPA66 A.3ValidationPressureResultsoftheControltestgroupforvessel2RPA66 A.4ValidationFlowResultsoftheControltestgroupfortherootvessel0 MPA.....................................67 A.5ValidationFlowResultsoftheControltestgroupforvessel1LPA..68 A.6ValidationFlowResultsoftheControltestgroupforvessel2RPA..68 A.7ValidationPressureResultsofthePHtestgroupfortherootvessel0MPA69 A.8ValidationPressureResultsofthePHtestgroupforvessel1LPA...70 A.9ValidationPressureResultsofthePHtestgroupforvessel2RPA...70 A.10ValidationFlowResultsofthePHtestgroupfortherootvessel0MPA71 A.11ValidationFlowResultsofthePHtestgroupforvessel1LPA.....72 A.12ValidationFlowResultsofthePHtestgroupforvessel2RPA.....72 B.1ChangestoPeakPressurealongthelengthofthelargestvessel,Vessel0, duetovaryingN/dxvalues.........................73 xiii

PAGE 14

B.2ChangestoPeakPressurealongthelengthofthesmallestvessel,Vessel 4,duetovaryingN/dxvalues........................74 B.3ChangestoPeakPressurealongthelengthofthelargestvessel,Vessel0, duetovaryingdtvalues............................74 B.4ChangestoPeakPressurealongthelengthofthesmallestvessel,Vessel 4,duetovaryingdtvalues..........................75 C.1PressureResultsoftheControltestgroupfortherootvesselMPA..76 C.2PressureResultsoftheControltestgroupforthesecondgenerationof vesselsLPA/RPA..............................77 C.3PressureResultsoftheControltestgroupforthe1Dterminalvessels Generation3.................................78 C.4PressureResultsoftheControltestgroupforthe0Dterminalvessels Generation12................................79 C.5PressureResultsofthePHtestgroupfortherootvesselMPA.....80 C.6PressureResultsofthePHtestgroupforthesecondgenerationofvessels LPA/RPA..................................81 C.7PressureResultsofPHtestgroupforthe1DterminalvesselsGeneration 3........................................82 C.8PressureResultsofthePHtestgroupforthe0DterminalvesselsGeneration12...................................83 C.9PressureResultsofthePH+LPABtestgroupfortherootvesselMPA.84 C.10PressureResultsofthePH+LPABLtestgroupforthesecondgenerationofvesselsLPA..............................85 C.11PressureResultsofthePH+LPABRtestgroupforthesecondgenerationofvesselsRPA.............................86 C.12PressureResultsofthePH+LPABLtestgroupforthe1Dterminal vesselsthatdescendfromtheLPAGeneration3.............87 xiv

PAGE 15

C.13PressureResultsofthePH+LPABRtestgroupforthe1Dterminal vesselsthatdescendfromtheRPAGeneration3.............88 C.14PressureResultsofthePH+LPABLtestgroupforthe0Dterminal vesselsGeneration12............................89 C.15PressureResultsofthePH+LPABRtestgroupforthe0Dterminal vesselsGeneration12............................90 C.16PressureResultsofthePH+replacementtestgroupfortherootvessel MPA.....................................91 C.17PressureResultsofthePH+replacementtestgroupforthesecondgenerationofvesselsLPA/RPA.........................92 C.18PressureResultsofthePH+replacementtestgroupfortheterminalvesselsofthe1Dmodel..............................93 C.19PressureResultsofthePH+replacementtestgroupforthe0Dterminal vesselsGeneration12............................94 D.1FlowResultsoftheControltestgroupfortherootvesselMPA.....95 D.2FlowResultsoftheControltestgroupforthesecondgenerationofvessels LPA/RPA..................................96 D.3FlowResultsoftheControltestgroupforthe1DterminalvesselsGeneration3....................................96 D.4FlowResultsoftheControltestgroupforthe0DterminalvesselsGeneration12...................................97 D.5FlowResultsofthePHtestgroupfortherootvesselMPA.......97 D.6FlowResultsofthePHtestgroupforthesecondgenerationofvessels LPA/RPA..................................98 D.7FlowResultsofPHtestgroupforthe1DterminalvesselsGeneration3.98 D.8FlowResultsofthePHtestgroupforthe0DterminalvesselsGeneration 12.......................................99 xv

PAGE 16

D.9FlowResultsofthePH+LPABtestgroupfortherootvesselMPA...99 D.10FlowResultsofthePH+LPABLtestgroupforthesecondgeneration ofvesselsLPA................................100 D.11FlowResultsofthePH+LPABRtestgroupforthesecondgeneration ofvesselsRPA................................100 D.13FlowResultsofthePH+LPABRtestgroupforthe1Dterminalvessels thatdescendfromtheRPAGeneration3.................101 D.12FlowResultsofthePH+LPABLtestgroupforthe1Dterminalvessels thatdescendfromtheLPAGeneration3.................101 D.14FlowResultsofthePH+LPABLtestgroupforthe0Dterminalvessels Generation12................................102 D.15FlowResultsofthePH+LPABRtestgroupforthe0Dterminalvessels Generation12................................103 D.16FlowResultsofthePH+replacementtestgroupfortherootvesselMPA.104 D.17FlowResultsofthePH+replacementtestgroupforthesecondgeneration ofvesselsLPA/RPA.............................104 D.18FlowResultsofthePH+replacementtestgroupforthe1Dterminalvessels Generation3.................................105 D.19FlowResultsofthePH+replacementtestgroupforthe0Dterminal vesselsGeneration12............................105 xvi

PAGE 17

1.Introduction 1.1PulmonaryHypertension PulmonaryhypertensionPHisdiagnosedbyameanpulmonaryarterialpressuremPAPgreaterthan25mmHgatrestor30mmHgwithexercise.Thereare currentlyvecategoriesaccordingtotheDanaPointclassication.Category1,also calledpulmonaryarterialhypertensionPAH,isthemostcommoncategoryandincludesidiopathicPAH,heritablePAH,drugandtoxininducedandotherassociations suchaswithconnectivetissuedisease,HIVinfection,portalhypertension,congenital heartdiseasesandschistosomiasis.TheotherfourcategoriesconsistofPHdueto leftheartdiseaseCat.2,PHduetolungdiseasesand/orhypoxiaCat.3,chronic thromboebolicPHCTEPH,Cat.4,andotherswithunclearmultifacetedmechanismsCat.5[1,2,3].Category1isfocusedonhereduetoitsprevalence;however, someinformationmayapplytoothercategories. PHisassociatedwithchangestothevascularwallinboththeproximalanddistalvasculature,whichiscomposedofthreelayers:theintima,mediaandadventita, whichareprimarilycomposedofendothelialcells,smoothmusclecells,andbroblasts,respectively[4].DuringPH,smoothmusclecellsSMCsproliferatecausing themediallayertogrowintowardsthelumen,decreasingthevesselradius,resultinginincreasedpulmonaryvascularresistancePVR,whichisseenbymanyasthe onlymajorconsequenceofPH[5].Fibroblastsarealsoproliferatingwhichcausesthe adventitiatogrowoutandawayfromthelumen.Themuscularizationofthemedia andstieningoftheadventitiaoccurinbothproximalanddistalpulmonaryarteries inPH[1,6]. Adventitialbroblastsplayaroleinincreasingthemediallayerastheycan dierentiateintomyobroblastswhicharesmoothmuscle-likecells,andmigrateto themediallayer.Fibroblastsintheadventitiaalsochangetheextracellularmatrix ECMcompositionofvascularwallsinPH;morespecically,thereisanincrease 1

PAGE 18

ofcollagentypesIandIII.Anincreaseincollageninthevascularwallleadstoan increaseinstinessofthewallaswell[3,4,5]. EndothelialcellECdysfunctionisknowntooccur,whichcanfurthersupport SMCproliferationandcancauseplexuslesionformationaswellasendothelialchannel formation.ECdysfunctionisalsoassociatedwithlossofpulmonaryarterialvessels [1,6].Aspartofthevascularremodelingprocess,leukocytesarebeingrecruitedby adventitialbroblastsandpro-inammatorycytokinelevelsarehigherinPHvessels [4],whichmayfurtheracttoexacerbatethesevascularchanges. 1.1.1InterplayofMechanicsandBiology ThevascularstructuralchangesofthepulmonaryvesselsinPHarelikelyto causechangestothemechanicalpropertiesofthevessel;forexample,anincreasein collagencausesanincreaseinvascularstiness.Vascularstieningisseenbymanyas akeyplayerofPHdiseaseprogression[3,7,5].StieningofthelargearteriesinPH contributestotheincreaseinmPAPandalsocancauseanincreaseinpulsepressure, reectedwaves,pulsewavevelocity,andpulsewaveintensity.Allofthesefactors canleadtoincreasedowpulsatility.Decreasinglumenradiusresultsinanincrease indistalvesselresistance.Adecreaseinvesselnumberalongwithdecreasinglumen radiuslikelycausesanincreaseindistalwallshearstressWSS.Bothanincreasein owpulsatilityandWSSinthedistalvasculaturethencontributestoECdysfunction whichfurtherexacerbatestheinammationandremodelingofthevascularwalls[7]. ThusWSSandowpulsatilityareimportantmechanicalpropertiestoinvestigate. 1.2ComputationalModeling ThegoldstandardfordiagnosingpulmonaryhypertensionPHisthroughright heartcatheterizationRHC,whichisaninvasiveprocedurethatgenerallyusesa Swan-Ganzcathetertomeasurepressureandmeanowcardiacoutput[8].A catheterisnotabletotravelfarintothevesselnetworkduetoitssize;thus,thereis 2

PAGE 19

notawaytodirectlymeasureowinthedistalvasculature.Othertechniquesexist tomeasurebloodownon-invasivelysuchasmagneticresonanceimagingMRIand ultrasonictechniques;however,theyalsocomewithlimitationsthatdonowallow forthedirectmeasurementofbloodowindownstreamvessels[9].Duetothese limitations,computationalmodelinghasbecomeanattractivemethodforestimating owindistalvasculature. Computationalmodelsarenon-invasivewaystoaskphysiologicalquestionsabout thebody.Generally,theygiveustheabilitytoinvestigatehard-to-measureproperties, tomakepredictionstocreatemoreeectiveanimalmodels,andtobetterinvestigate treatmentoptions.Therearemanytypesofcomputationalmodels,suchaszerodimensionalD[10,11],one-dimensionalD[12,13,14,15,16,17,18,19,20,21] andthree-dimensionalD[18,19,20,21].0Dmodelsarepoint-basedandonly includethetimedomain,1Dmodelsaddonespatialdomainwhichisthelength alongavessel,and3Dmodelsincludeallthreespatialdimensionsofthevesselas wellasthetimedomain.Modelscanalsobeacombinationofthesetypes. 3

PAGE 20

Figure1.1:Graphicrepresentationof3D,1Dand0Dmodelsandthetrade-obetween computationalcostandphysiologicalrelevancei.e.howrepresentativethemodelis. Eachmodelgivesadierentsetofinformationaboutow;3Dmodelsallowfor representationofmorecomplexowbutatahighcomputationalcostwhilea0D modelhasaverylowcomputationalcostbutonlyincludesmoresimpleowfeatures. Thisisrepresentedingure1.1.Typically,0Dmodelsareusedtolookatthesystemlevel,1Dmodelscangiveaccuratebloodpressureandowpropagationinformation foralargenetworkofvessels,and3Dmodelsallowforadeeperinvestigationof uid-structureinteractionsFSI,whichistherelationshipbetweenbloodowand vesselwalldeformation,andcanfullyrepresentbloodvelocityprolesDmodels canonlyapproximate[9,22].Somemodelsexisttostudyanentirevascularnetwork [12,13,16,21].Modelsalsoexisttostudyaspecicsubsystem,suchasapulmonary vasculaturemodel[14,23],whichiswhatisusedhere. Forallmodels,inletandoutletboundaryconditionsarerequiredaswellasgeometricinformatione.g.radius,lengthandmechanicalpropertiese.g.stiness, 4

PAGE 21

viscosityofthebloodandvessels.Theinletboundaryconditiontypicallycomesfrom clinicalmeasurementssuchasdatafromRHCorMRI.Outletboundaryconditions varybasedonthemodelingapproach[9,24].Commonoutletboundaryconditions usedareWindkesselvariations[15,17,18,19],andstructuredtree[12,13].Some3D modelsmayuse1Dmodelsfortheiroutowboundaryconditionsaswell[25,20,21]. Ifanetworki.e.notjustasinglevesselisbeingmodeled,abifurcationboundary conditionisalsoneededatvesseljunctions[9,24]. Duetotheabilitytohavearelativelylowcomputationalcostwithaccuraterepresentationofbloodowandsomemechanicalpropertiesofthevessels,a1Dmodel isoneofthemostattractiveoptionsformodelingbloodvesselsandisusedherefor thatreason.Currently,1Dmodelassumptionsallowforagoodrepresentationofrelativelylargerbloodvesselsandusesimplied3DNavier-Stokesequations.Common numericalschemesexisttosolvethesystemofequationssetupfor1Dmodels,suchas Galerkinvariations,FiniteDierencesmethods,andcharacteristicmethods[24,26]. 1.3Justication InPH,thereisdistalvesselthickeningoccurring.Datafromotherliteraturesuggeststhatproximalstieningcauseshigherowpulsatilitywhichcausesendothelial celldysfunctionandactivatesbroblasts.Activatedbroblastsproliferateanddepositcollagenwhichcausesstieningandthickeningofthedistalvessels.Onestudy providesgreatevidenceforthisproposedmechanismbyshowingthathemodynamic unloadingiseectiveattreatingPHphenotypes[27].Inthisstudy,PH-modeledrats werecomparedtoPH-modeledratswithabandedleftpulmonaryarteryLPAand studiedthehistologyoftherats'arteries.TheyfoundthattheLPAofthebanded ratshadreducedexpressionofpro-inammatorycytokines,andotherinammatory andproliferativemarkers,whichalmostmatchedthatofthenormalrat.Bloodow totheLPAinthebandedratsreducedtolessthanhalfoftheoriginalbloodow, whilebloodowincreasedby40%intherightpulmonaryartery[27].Thisstudy 5

PAGE 22

almostcompletelyreversedwhatwethinktobeassociatedwithPH.Giventhat bandingshowedreversalofthesediseasecharacteristics,thequestionofwhatkindof mechanicalchangesareoccurringtoseethesebiologicalchanges.Thispapershowed abiologicalandmechanicalrelationshipbetweentheproximalvasculature,butthe relationshipbetweenmechanicsandbiologyofthedistalvasculatureisstillmissing. Inordertostudythedistalowmechanics,acomputationalmodelcanbeused. Givenabandedleftpulmonaryarteryasdonein[27],owcharacteristicsinthedistal vasculaturecanbeinvestigated.Thereisalsothepotentialforreplacingtheleft andrightpulmonaryarterieswithsyntheticpulmonaryarteriestodeterminedistal owmechanics.Ifbandingandreplacementoftheproximalvesselsshowsimilar distalowmechanics,thenreplacingtheleftandrightpulmonaryarteriescould beusedasatreatmentorreversalofPH.Thusfourtestgroupswillbesimulated here:healthycontrol,PH,PHwithbandedLPA,andPHwithreplacementleftand rightpulmonaryarteries.Thesesimulationswillhelpusdetermineifproximalvessel stinessisaectingdistalow.Ifthisrelationshipisestablished,thenanimalmodel studiescouldbedonesimultaneouslywithanimal-specicsimulationstoseeifthere isalsoareversalofthesediseasephenotypesasseenin[27]duetochangestoproximal vesselstinessasdoneinthesesimulations.Establishingthisrelationshipinthese simulations,andconrmingtheminananimalstudy,wouldprovideagreatstarting pointtoinvestigatereplacingstienedproximalvesselswithsynthetic,compliant vesselsasatherapeutictreatmentforPH. 1.4Goals Givenastier,vascularnetworkinthepulmonaryhypertensionPHcase,we expecttoseeanincreaseinowpulsatility[7],investigatedviapulsepressure,and adecreaseinwallshearstressWSS[28,29].Wehypothesizethatrestoringthe elasticityofproximalvasculatureofaPHmousemodelwillreturnpulsepressureand WSSofdistalpulmonaryarteriolestocontrollevels. 6

PAGE 23

Inthisstudy,therelationshipbetweenproximalpulmonaryarterystinessand distalowwereinvestigatedusingacoupled1D-0Dmodel.Firstthe1D-0Dmodel wascreatedinPythonandMATLABandvalidatedagainstpreviouslypublisheddata [23].Then,simulationswereruntosimulatefourdierenttestgroupsbasedon[27] tostudychangesinWSSandpulsepressure.Theaimsaredescribedbelow. SpecicAim1:Developacomputationalmodelofthepulmonaryarterialsystem toallowsolvingofbloodowhemodynamicsinalarge-scalenetwork. SpecicAim2:Executepulmonarymouse-specicsimulationsmatchingproximalpulmonaryarterialstinessandpulmonaryvascularresistancetocomparepulse pressureandwallshearstressofdistalvasculatureofthefollowingmodels:control; pulmonaryhypertensionPH;PHwithbandedleftpulmonaryartery;andPHwith asynthetic,compliantleftandrightpulmonaryarteryreplacement. 7

PAGE 24

2.Methods 2.1ModelOverview Toestablishevidencethatarelationshipexistsbetweenproximalmechanical stinessanddistalow,theentirepulmonaryvascularnetworkwillbemodeledusing ahybrid1D/0Dmodelwherethe1Dmodelisforthelarger,proximalvesselsand the0Dmodelisusedforthesmaller,distalvessels.First,thevesselpropertiesand treeareestablishedandputintoacongurationle.Then,the1Dmodelsetsupthe arterynetworkusingthecongurationle.Aphysiologicalwaveformisinputinto the1Dmodelwhichcomputestheproximalvessels'ow,areaandpressurewaves. Theoutputfromthe1Dmodelistheinputforthe0Dmodelincludingresistance andcapacitanceRCvalues,andthe1Dmodeloutletowrate.The0Dmodelthen setsupavasculartreeusingRCvalueswheretheRCvaluesoftheentire0Dmodel sumtoequaltheRCvaluesattheoutletofthe1Dmodel.The0Dmodelusesthe 1Dmodeloutletowasitsinletowtocomputetheowandpressurewavesforthe distalvessels.Theresultingdistalpressureandowarethenusedtocalculatepulse pressureandwallshearstress.Aowchartforthemodelisshownbelow: Figure2.1:High-LevelOverviewoftheFullComputationalModel 2.1.1Assumptions 8

PAGE 25

Themodelsassumethevesselsareelasticandaxisymmetrictubes,containing Newtonianuid.Otherassumptionsforthe1Dmodelarethatthebloodpressureis uniformacross ~r ,theradius,andthereispulsatilelaminarow. 2.21DModelVaMpy The1Dcomputationalmodelutilizesuiddynamics,wallmechanics,andboundaryconditions.Theresultingsystemofquasi-linearhyperbolicPDEsconsistsofthree equationswhichareusedtosolveforuxinsideofthevessellumen,cross-sectional area,andpressure.Thethreeequationsaresolvedusingapython-based,opensource modelavailableongithub,namedtheVascularModelinginPythontoolkitorVaMpy [30].Thismodelstartswithasingleparentvesselthatbifurcatesintodaughtervessels.Theresultingtwodaughtervesselsalsobifurcate.Thisprocesscontinuesuntil theterminalvesselsarereached.Thetotalnumberofvesselsmadeisdenotedby 2 depth -1.Therefore,ifthearterynetworkhadadepthof3,aka3generations,there wouldbe7vessels.Notethatvesselnumberingbeginswith0.Asanobject-oriented program,eachvesseliscontainedinitsownartery"classwhichareallcontained withinintheArteryNetwork"class.TheVaMpymodelusesacongurationlethat containsthevesselgeometryandvascularpropertyinformation.Anexamplevascular treewith3generationscontainingoutowboundaryconditionsisshowningure2.2 referencingvesselsinthepulmonaryvasculaturestartingwiththemainpulmonary arteryMPA. 2.2.1Equations The1Dmodelsolvesthreeequations.Thersttwoequationsarederivedthrough integrationoftheNavier-Stokesequationsforconservationofmassandmomentum inaone-dimensionalcoordinatesystem.Thedetailsofthissolutioncanbefound elsewhere[31]withtheresultingsystemofnon-dimensionalizedequationsrepresented bytheequationshownbelowinEq2.1 ~ U t + ~ F x = ~ S .1 9

PAGE 26

Figure2.2:Anexampleofatreewith3generationswherevessel0representsthe mainpulmonaryarteryMPA,vessel1representstheleftpulmonaryarteryLPA, vessel2representstherightpulmonaryarteryRPAandvessels3-6aretheterminal vessels.Theoutowboundaryconditionisa3-elementwindkesselappliedateach terminalvessel. where ~ U = 0 B @ A x;t q x;t 1 C A , ~ F = 0 B @ q x;t q x;t 2 A x;t + f r o p A o x A x;t 1 C A , ~ S = 0 B @ 0 S 1 1 C A S 1 = )]TJ/F17 7.9701 Tf 10.494 5.698 Td [(2 R x;t b Re q x;t A x;t +[2 p A x;t p f r o + p A o x df r o dr o )]TJ/F19 11.9552 Tf 11.955 0 Td [(A x;t df r o dr o ] dr o x dx ~ U containstheunknownstobesolved: A x;t ,thecross-sectionalarea,and q x;t , theowrate.Elasticityofthevesselisdescribedby f r o = 4 3 Eh r o ,andisempirically estimatedusingEq2.2below: Eh r o = k 1 e k 2 r o + k 3 .2 10

PAGE 27

R x;t representstheradiusofthevessel, r o x istheinitiallyrelaxedvesselradius withcorrespondingcross-sectionalarea A o x . b istheboundarylayerthickness estimatedby b = p T= where isthekinetmaticviscosity, Re istheReynold's numberandiscalculatedusing Re = q c = r c with q c and r c beingthecharacteristic uxowandcharacteristicradius,respectively,and T isthetimeperiodofone cardiaccycle.Tomodelataperingvesselwithupstreamanddownstreamradii, R u and R d ,thevesselradiuscanbedescribedasafunctiondownitslength L as r o x = R u exp log R d R u x L . Eq2.1issolvedbydiscretizingusingRichtmyer'stwo-stepLaxWendroexplicit scheme,detailedin[30,31]whichissecond-orderaccurateintimeandspace,i.e. O t 2 ; x 2 andmustsatisfytheCFLcondition: t x j q A c j )]TJ/F17 7.9701 Tf 6.587 0 Td [(1 wherethewave speedis c = p A= p=A and x and t arethespacestepandtimestep, respectively[30]. Thethirdequationneededtosolvethissystemofequationsequationswith3 unknownsisthestateequationwhichisalinearequationrelatingpressure p x;t to areaasshownbelow: p x;t = p o + f r o )]TJ/F25 11.9552 Tf 11.955 10.806 Td [(p A o x =A x;t .3 where p o istheinitialdiastolicpressure[30].Thus,thismodeloutputsthreevariables:cross-sectionalarea, A x;t ,owrate, q x;t ,andpressure, p x;t 2.2.2BoundaryConditions Tosolvethesethreeequationsinthis1Dmodel,boundaryconditionsmustalso beappliedattheinlet,atbifurcations,andattheoutlets.Attheinlet,aphysiologicalowwaveformisgivenasmeasuredbyultrasounddonesimultaneouslywith catheterization.Theinowwaveformusedforthecontrol,andPHcasesareshown ingure2.3. Abifurcationboundaryconditionappliesbetweenanyparentvesselanditstwo daughtervessels.Thebifurcationboundaryconditionconsistsofpressurecontinuity, 11

PAGE 28

Figure2.3:InletowwaveformforcontrolandpulmonaryhypertensionPHtest groups. p p L;t = p d i ;t ,andconservationofow, q p L;t = P i =1 ; 2 q d i ;t .Tocomplete thesystem,theserelationshipsleadtoasystemofeighteenequationsandeighteen unknowns.Thesystemofequations,shownbelowinEq2.4,issolvedusingNewton's method. x k +1 = x k )]TJ/F15 11.9552 Tf 11.955 0 Td [([ J x k ] )]TJ/F17 7.9701 Tf 6.586 0 Td [(1 ~ f J x k fork =1 ; 2 ; 3 ;::: 18.4 where k isthecurrentiterationfrom1to18,[ J x k ]istheJacobianmatrixand ~ f J x k isthevectorofresiduals[30]. Attheoutlet,azero-dimensional,3-elementWindkesseloutowboundaryconditionisusedtoclosethecomputationaldomain,andischaracterizedbyresistance andcapacitanceRCRvalues.ThesamesetofRCRvalues,speciedbythecongurationle,areappliedatalloftheoutletvessels.ArepresentationoftheRCR circuitasanoutletboundaryconditionisshowningure2.2.TheseRCRvaluesare usedtoestimatetheoutowofthevesselasshownbelowinEq2.5whereMisthe 12

PAGE 29

nalspacepointandnisthecurrenttimestep. q n +1 M = q n M + p n +1 M )]TJ/F19 11.9552 Tf 11.955 0 Td [(p n M R 1 + tp n M R 1 R 2 C T )]TJ/F15 11.9552 Tf 13.151 8.088 Td [( tq n M R 1 + R 2 R 1 R 2 C T .5 Areaissolvedforusingadiscretizedmassconservationequationshownbelow: A n +1 M = A n o )]TJ/F15 11.9552 Tf 14.363 8.088 Td [( t x q n +1 M )]TJ/F19 11.9552 Tf 11.955 0 Td [(q n +1 M )]TJ/F17 7.9701 Tf 6.586 0 Td [(1 .6 TheoutletpressureisthensolvedusingadiscretizedversionofEq2.3andisshown below: p n +1 M = f n +1 M 1 )]TJ/F25 11.9552 Tf 11.955 21.099 Td [(s A o M A n +1 M ! .7 Thisprocessisiteratedthroughafteraninitialguessfor p n +1 M untiltheoutletpressure solutionconverges,wheretheerroris < 1 e )]TJ/F17 7.9701 Tf 6.587 0 Td [(6 .Theiterationislimitedtoamaximum of1000. 2.2.3CodeAdaptation TheVaMpymodelhaslimitationsandthereforehadtobeexpandeduponforthis project.Thechangesmadetothemodelinclude:solvingformultiplebifurcations, usingdierentRCRvaluesforeachterminalvessel,allowingdxtovaryforeach vesselandusingdierent Eh r o valuespervesselthanthatcalculatedfromEq2.2. Otherfunctionalitywasalsoincludedforinternalusesuchastrackingtheruntime. OnemajorlimitationoftheoriginalVaMpymodel,wasthefactthatitcould onlysolveforasinglebifurcation.Inordertousethismodelforalargernetwork,the modelwasexpandedtosolveformultiplebifurcationsbychangingtheindexingofthe vessels.Initially,thecodeindexedthedaughtervesselswith d 1 = p +1and d 2 = p +2 insideoftheget daughtersfunctionoftheArteryNetworkclass.Theget daughters functionwasupdatedtoinsteadreturn d 1 =2 p +1and d 2 =2 p +2wheretheparent indexstartsat0. IntheinitialVaMpymodel,RCRvalueswereinputintoacongurationlewhich weredirectlyusedattheterminalvessels.Themodelhassincebeenchangedsothat 13

PAGE 30

thecongurationleinsteadcontainsan R T and C T valuewhichrepresentthetotal vascularresistanceandtotalcapacitance,respectively.Thesevaluesarecomputed intheMPAandthereforerepresentthevascularresistanceandcapacitanceforthe entirenetwork.Thesevaluesarethenusedtoiterativelycalculatethe R T;j and C T;j valuesforeachvesseluntilreachingtheterminalvessels.First,themeaninputow, q ,andmeaninputpressure, p ,arecalculated.Themeanowisjusttheaverageof theinletwaveformandthemeanpressureiscalculatedfrom R T = p )]TJ/F20 7.9701 Tf 6.586 0 Td [(p c q where p c is thecapillarypressureanddropsoutofthisequationwiththeassumptionthatthe capillarypressureis0andthatthemeanpressurestaysconstantforallvessels.The terminalresistancevaluesarethencalculatedusing R T;j = p q j ,where q j isthemean owforvesseljandiscalculatedbyusingPoiseuille'slawateachbifurcationuntil theterminalvesselsarereached: q d i = G d i P i G d i q p ;whereG d i = r 4 o 8 L d i fori =1 ; 2.8 G d i isthevesselconductance, q p isthemeanowtotheparentvesseland q d i represents themeanowtotheeachofthetwodaughtervessels.Once R T;j hasbeencalculated fortheterminalvessels, C T;j canbefoundusingtherelationship = R T C T and = R T;j C T;j assumingthat isconstantthroughoutthenetwork. R T;j canthenbe splitupinto R 1 ;j and R 2 ;j : R 1 ;j = aR T;j where a =0 : 2and R 2 ;j = R T;j )]TJ/F19 11.9552 Tf 12.083 0 Td [(R 1 ;j .[23]. ThiswasimplementedasanewfunctioninsideoftheArteryNetworkclass. DynamicallychangingdxbasedonthevesselwasdonebyapplyinganNvalue. ThisNvaluespeciesthenumberofspacestepstoincludeforallvessels,whichmeans thateachvesselwillhaveadierentdxbasedonitslengthratherthanallvessels usingthesamedxvalue.Inordertoincludethisfunctionality,anewmeshfunction wascreatedintheArteryNetworkclass,mesh nx,withtheoldmeshfunctionbeing renamedtomesh dx.Themesh dxfunctiontakesinthedxvalueandappliesthe meshfunctionintheArteryclassoneachfunction.Thenewmesh nxfunctiontakes 14

PAGE 31

intwovariables:nx,thedesirednumberofspacesteps,andrc,thecharacteristic radius.Nisanewvariabletoincludeinthecongurationleandisfedinasthenx inputforthemesh nxfunction.Thenxandrcvaluesarethenusedtocalculatedx usingtheequation: dx = Lr c nx )]TJ/F17 7.9701 Tf 6.586 0 Td [(1 whichwasobtainedbyrearrangingtheequationfor nxinthemeshfunctionoftheArteryclass: nx = int L dx +1. Asdescribedpreviously, Eh r o iscalculatedfromEq2.2intheoriginalVaMpymodel. Inordertofulllthegoalsofthisproject, Eh r o wasadjustedtoallowthecongle toalsoinputspecic Eh r o valuespervessel.Thechangewasmadebyaddingasimple inputargumentthatchecksfortheexistenceof`Eh'inthefunctionargumentsjustas VaMpyhasdonewithothervariables.If`Eh'exists,thenthecodewillusetheentered valuefromthecongle,otherwisetheprogramwilluseEq2.2asoriginallyspecied. ThiswasimplementedinthemeshfunctionsofboththearteryandArteryNetwork class. ATimematrixwasaddedtotheArteryclasstoallowtrackingoftotalsimulationruntimeaswellastheruntimeforvaryingprocesses.Thetimetrackingwas implementedbyusingthetime.clockfunction,whichrequiredimportingthetime librarytobeused.Thetimespentinthebifurcationboundaryconditionsolvefunctionwastrackedaswellasthetimespentinthesolvefunctionasawhole.Thesolve functionintheArteryNetworkiswheretheactualsolvingofthe1Dmodelequations occursandisthereforethemajorityofthetimespentforthemodel.Thebifurcation boundaryconditionwasdeterminedtobethemosttimeconsumingpartofthesolve functionataround80%oftheruntime. Anerrorwasalsofoundinvolvingthebifurcationboundaryconditioninthesolve functionoftheArteryNetworkclass.Originally,theboundaryconditionsolution forthetwodaughtervesselswereaccidentallystoredintheoppositedaughtervessel object,sotheywereswappedtocorrectthisissue. 2.2.4Runninga1DSimulation 15

PAGE 32

TheIntegratedDevelopmentEnvironmentIDEusedtorunthemodelswas AnacondaNavigatorwithSpyderPython3.6.Inordertorunasimulationwith thisIDE,thefollowingcommandswererun:runle`bifurcation example Rt N.py', args=`lename.cfg'toobtaindataandsubsequently:runle`plot example.py',args =`lename.cfg'toobtainplotswherethelename.cfgisthecongurationlefor anygivensimulation. 2.2.5Limitations Onemajorlimitationofthe1Dmodelisitsslowruntime.Initialruntimewas5 daystosolvetoadepthof3foratotalof7vessels.Themodelwasinvestigatedto ndwhatpartofthemodeltookthelongestandthendetermineifthemodelcouldbe spedup.Theculpritwasthebifurcationboundaryconditionduetoits18non-linear equationsbeingsolvedbyNewton'smethodwhichincreasesthecomputationalcost. Duetotheslownessofthemodel,the1Dmodelwascombinedwitha0Dmodelin ordertomodelanentirevascularnetworkrelativelyquickly.Oncetheparameters werechangedforvalidation,meshindependence,andactualmodelsimulations,the runtimeforthe1Dmodeldecreasedto14-35hours. Anotherlimitationofthemodelisthatthechangesmadeto Eh r o arenotcompatible withVaMpy'sabilitytogeneratealargetreeusing and valuesthatrepresentthe changeinradiusfromaparentvesseltoitstwodaughtervessels.Finally,the1D modelcodeisnotformattedtoallowforanasymmetricbranchdistributioninthe tree.Thismeansthatallvesselsmustbifurcateexceptterminalvessels,andall theterminalvesselsarelocatedatthesamedepth.Forexample,thecodecannot accommodateaterminalbranchatadepthof2generation2withterminalbranches atadepthof3generation3asshowningure2.4. 2.30DModel The0DModelutilizesacircuitmodelbasedoofatwoelementwindkessel model,whichischaracterizedbyasingleresistanceandasinglecapacitanceRC 16

PAGE 33

Figure2.4:Anexampleofatreewithasymmetricbranchingwherethereareterminal branchesatgeneration2andatgeneration3. value.EachvesselinthismodelisrepresentedbyitsownRCcircuit.TheRCcircuits areconnectedinatreewitheachparentbifurcatingintotwodaughtervessels,just likewiththe1Dmodel.ThecircuitmodelisgraphicallycreatedusingSimulinkwith thevariablesforeachRandCvaluebeingdeterminedusingaprogramcreatedin MATLAB. 2.3.1VesselNetworkGeneration Oneoftherststepsforthe0Dmodeltobeabletobecreated,istocreatethe entirevascularnetwork,ortree.Therearemanyconsiderationswhensettingupthe treesuchasbranchsymmetry,depth,geometry,andmechanicalproperties.Inthis model,theprogramstartswitharootradius,whichwillcorrespondtoaterminal vesselfromthe1Dmodel,andwilluse and valuesasscalingfactorsforthe daughtervesselsforanygivenparentvessel.Thevalueschosenfor and were 0.9and0.6,respectively[32].Thetreecontinuestogrowuntilavesselislessthan theterminalradiusof0.005cm[33,34],whichrepresentsthecapillarybed,orata maximumof9generationstogiveatotalof12generations,includingthevesselsin the1Dmodel.Accordingto[35],smalleranimalshaveatotalof11-12generations 17

PAGE 34

forpulmonaryarterytrees,whichiswhythisnumberwaschosen.Ifaparentvessel wouldcreateadaughtervesselthatislessthantheterminalradius,thentheparent vesselbecomesaterminalvessel.Eachofthesevaluesarerepresentedinamatrix, ~r ,whereifaparentisindexedby i;j thenthetwodaughtervesselswouldhavean indexof i +1 ; 2 j )]TJ/F15 11.9552 Tf 12.328 0 Td [(1and i +1 ; 2 j where i and j representtherowandcolumn positions,respectively.Thelengthforeachvesselisthencalculatedusingaconstant length-to-radiusratioandstoredintoamatrix, ~ L . Oncetheradiusandlengthmatricesarecreatedforthenetwork,theResistance andCapacitancevaluesarecalculatedforeachvessel.Theresistanceiscalculatedby usingPoiseulle'sowwiththefollowingequation: ~ R = 8 ~ L ~r 4 .9 where isthedynamicviscosity.Capacitanceisthencalculatedusingthesame valueandrelationshipfromthe1Dmodel.Onceresistanceandcapacitancematrices aresolved,theyaresummedfortheentiretreegiving R sum and C sum .For resistancesinparallel,theresistanceissummedas: 1 R T = 1 R d 1 + 1 R d 2 where d 1and d 2are thetwodaughtervessels.Resistancesinseriesaddnormallyas: R T = R 1 + R 2 where R 1 and R 2 aretworesistorsinseries.Forcapacitanceinparallel,thecapacitance addsnormallyas: C T = C 1 + C 2 where C 1 and C 2 aretwocapacitorsinparallel. Capacitanceinseriesissummedliketheresistanceinparallelissummedandisgiven by: 1 C T = 1 C d 1 + 1 C d 2 where d 1and d 2arethetwodaughtervessels.Resistanceand Capacitancearesummedindependentlyofoneanother,asinthecapacitancevalues donoteectthetotalresistancevaluesandviceversa. Thetotalvaluesarethencomparedtotheoutlet R T and C T valuesfromthe 1Dmodelandthenthetreeisadjustedsuchthat ~ R new =scalingfactor* ~ R and ~ C new = ~ R new ,wherescalingfactor= R T R sum and isthesameasinthe1Dmodel. 2.3.2SimulinkModel 18

PAGE 35

Table2.1:Unitscomparisonbetweenthehydraulicorbloodowmodelandthe electricorcircuitmodel. ModeltypeFloworCurrentPressureorVoltResistanceCapacitance Hydraulic mL=sdyn=cm 2 mmHg g=scm 4 cm 4 s 2 g Electric AmpsVoltsFaradsOhms Oncethesevaluesarefoundforthetree,aphysicalmodelofthetreeismade inSimulink.Atranslationbetweenthehydraulic,vesselmodelandelectric,circuit modelisshownintable2.1.Therefore,owwillbetreatedascurrentuidow becomeselectronowandpressurewillbetreatedasvoltageuidpotentialbecomes electricalpotentialinsidetheSimulinkmodel.Anexampleofa2-generationtreeis showningure2.5.Forthemodeltooutputpressureandowwavesforavessel, Figure2.5:Anexampleofa2-generationRCcircuittreemodelmadeinSimulink showingtheinletowandoutowboundaryconditionBC. VoltmetersandAmpMetersareneededinthemodel.Anexampleofthesame2generationtreefromgure2.5isshowningure2.6withthesemetersadded.The VoltmetersandAmpMetersalsoneedtobeconnectedtoascopeblock,whichisjust ablockthatprovidesgraphicvisualization.Inthemodelyoucanseeascopeblock 19

PAGE 36

ontheleftforthevoltageandascopeblockontherightforthecurrent. Figure2.6:Anexampleofa2-generationRCcircuittreemodelmadeinSimulink withVoltmeters,AmpMeters,andscopeblocks. 20

PAGE 37

Withintheindividualcomponentblocks,resistanceandcapacitancecanbehardcoded,whichisanextremelymanualprocessespeciallyasthenetworkgetslarger. Alternatively,variablesnamedintheMATLABworkspacecanbeusedwithinthe blockssuchthatthevaluescanbeeasilychangedbyupdatingthisvariablewithin theMATLABworkspace.Usingthismethod,afewvariablescanbechangedwithin theMATLABscriptthatgeneratesthevesselnetworkandthenalloftheresistance andcapacitancevaluesupdateautomaticallyintheSimulinkmodel. 2.3.3BoundaryConditions Inordertosolvetheentirecircuit,boundaryconditionsmustbeapplied.The inletowwavecomesfromthe1Dmodelterminalvessels.CSVlesareautomaticallysavedinfoldersafterrunninga1Dsimulation,sotheowwaveformfromthe lastspacestepoftheterminalvesselisusedastheinletboundaryconditionand convertedintocurrent.Ateachbifurcation,voltageandcurrentareconservedthe samewayaspressureandowinthe1Dmodelwhere V p = V d i and I p = P i =1 ; 2 I d i .The outletboundaryconsistsofanotherRCcircuitthatrepresentsallofthedownstream vasculaturefromtheterminatedvessel.TheRCoutletcircuitsareshowningures 2.5and2.6.Itshouldbenotedthatifaterminalvesselisrepresentingacapillary bed,thennooutletRCcircuitisneededbecausethereisnodownstreamvasculature. 2.3.4Runninga0DSimulation Tomakethingseasier,MATLABhasacommandthatcanbeusedinsideofa scriptinordertorunaSimulinkmodel.Ascriptwasmadetocreatethevesselnetworkparametersandthenthecommandsim`modelle.slx'canbeusedtorunthe Simulinkmodelcontainedinthemodelle.slxle.WithintheSimulinkmodel,data fromthescopeblocksaresavedasvariablesintheMATLABworkspace,andcanbe subsequentlyusedinthesameMATLABscriptforpost-processing.Alternatively, thevesselnetworkparameterscanbecreatedwithaMATLABscriptle,thenthe 21

PAGE 38

Table2.2:Resistance, R T ,andCapacitance, C T ,attheoutletofthe1Dmodel.The scalingfactorisappliedtothePoiseuillelawasshowninequation2.9inthe0Dmodel whenbuildingupthetreesuchthatthetotalresistanceandtotalcapacitanceinthe 0Dmodelareequalto R T and C T ,respectively. TestGroup R T mmHg s mL C T mL mmHg scalingfactor unitless Control 233.62.2949e-41.0697e-4 AllPHtestgroups 535.22.945e-52.4509e-4 SimulinkmodelcanberunmanuallywiththevariablesstillbeingsavedintheMATLABworkspace.ThevariablesaresavedintheMATLABworkspacebygoingto Scope ! Settings ! LoggingandcheckingtheboxlabeledLogdatatoworkspace" andthenaVariablename"canbespecied.TheMATLABscriptsusedtorunthe simulationsareshowninAppendixE. Thenaltreeusedfor0Dmodelisshowningure2.7andonlyfollowsoneline ofvesselsdownforsimplicity. 2.3.50DSimulationData Outletowwaveformsfromthe1Dmodelvesselswereusedinthe0Dmodel asinletboundaryconditions.Allofthevesselgeometriesweresymmetric,except forthePH+LPABconditionwherevessel1LPAhadasmallerradiusthanits counterpart,vessel2.Therefore,onlyoneterminalvesselwasusedforallofthe subsequent0DsimulationsinSimulink,exceptthePH+LPAB.Inthiscase,two terminalvesselsweretested,onedaughtervesselforvessel1andonedaughtervessel forvessel2andarelabeledasPH+LPABLandPH+LPABR,respectively, toindicatethatthevesselisadescendantoftheleftorrightpulmonaryartery.The R T , C T andscalingfactorvaluesobtainedforthe0Dmodeltreebuildupareshown intable2.2. 22

PAGE 39

2.3.6Limitations Onemainlimitationofthe0Dmodelisthatitmaynotactually"end"atthe capillaries,sinceithasastrictcutoof9generations;however,thisisagoodestimate basedonliteratureshowingthatsmallmammalshavearound12generations[35].The terminalradiusvalueforthecontrolgroupwas0.0069cmwhichisabithigherthan thecapillaryradiusvaluefoundinliteratureof0.005cm,somayaecthowaccurate theresultsare;however,thesesimulationswereusedtolookattrendssothiserror doesnothaveagreateectonwhatwasconsideredhere.Simulationscouldbedone inthefuturetoincludemoregenerations,orthe and valuesmayneedtobe adjustedtoachievemoreaccurateterminalradiusvalues. Anotherlimitationofthe0Dmodel,isthefactthatitdoesnotincludewave reections,andthatthepressureandowwavesseemtoattenuatemorethanwhatis seeninthe1Dmodel.Whilethesewavescouldactuallyattenuatemoreinthedistal vessels,thereisnowaytotestthisduetodevicelimitationsdescribedpreviously. Again,thisstudywasinterestedintrendsofthepressureandowandsincethis attenuationwouldbeconsistentacrossalltestgroups,thisattenuationdoesnot greatlyaecttheresultspresentedhere. 2.4NumericalSimulationsforCombinedModels Fourtestgroupswereinvestigatedhere:control,pulmonaryhypertensionPH, PHwithabandedleftpulmonaryartery,andPHwithsyntheticleftandrightpulmonaryarteries,whichwillbelabeledas:control,PH,PH+LPAB,andPH+ replacement,respectively.Table2.3showstheradiusandlength-to-radiusratiovaluesforeachvessel,takenfromdatain[23];however,theradiusandlength-to-radius ratiovalueswereaveragedforallvesselsatthesamegenerationsuchthateachvessel inthesamegenerationhasasymmetricvesselgeometry.Thiswasdonetosimplifythe modelandtodecreasethecomputationalcost.Theresultingaveragedvaluesshown intable2.4,whichshowsUpstreamradii, R u ,downstreamradii, R d .Thesevalues 23

PAGE 40

Table2.3:Upstreamradii, R u ,downstreamradii, R d andlength-to-radiusratio, lam , valuesforeachvesselforatreewithasymmetricvesselgeometry. Vessel0123456 R u cm0.0470.0260.0370.0240.0130.0320.017 R d cm0.0470.0260.0370.0240.0130.0320.017 lam unitless8.72317.115410.054110.04174.06.312512.4706 Table2.4:Upstreamradii, R u ,downstreamradii, R d andlength-to-radiusratio, lam , valuesforeachvesselforatreewithsymmetricvesselgeometry. Vessel01&23,4,5,6 R u cm0.0470.03150.0215 R d cm0.0470.03150.0215 lam unitless8.72313.58478.2062 donotchangeinthesesimulationsforsimplicity,butareshownheretoillustratethe originalfunctionalityofthecodetoallowforachangingradiusalongthelengthof thevessel.Variablesgiveninthecongurationletosolveasimulationareshown intable2.5,whichshowsvaluesusedforvalidatingthemodelaswellasrunningthe simulations.Giventhiscongurationle,themodelwillsolvefor4cardiaccyclesas indicatedby t c where T isthetimeperiodforonecardiaccycle. Elasticity Eh r o ,resistance R T andcapacitance C T variedbetweentestgroups withthevaluesshowninTable2.6.Generally Eh r o appliestotheproximalvessels,and R T and C T applytothedistalvessels.Thesevaluescamefrom[23]but Eh r o values wereconvertedtotheunitsforthe1Dmodelbeforegoingintothecongurationle. R T and C T wereenteredexactlyintothecongurationleandthenthe1Dmodel codeconvertedthemintothedesiredunits.Also,insteadofusingEq2.2tovary 24

PAGE 41

Table2.5:Variablesusedacrossallsimulationsforvalidationandactualmodelsimulationswhere r c isthecharacteristicradius, q c isthecharacteristicow, isblood density, isthebloodviscosity, depth istotalnumberofgenerations, T isthetime periodforonecardiaccycle, t c isthenumberofcardiaccyclestosolvefor, dt isthe timestep,and N isthetotalnumberofspacesteps. ValidationModel r c unitless11 q c unitless1010 g=cm 3 1.061.06 cm 2 =sec 0.0460.046 depth unitless23 T sec0.1090.109 t c unitless44 dt sec1e-51e-7 N unitless5046 25

PAGE 42

Table2.6:Simulationvariablesusedbetweentestgroups. Eh = r o g cm s 2 R T mmHg s mL C T mL mmHg Control 4.7306e458.43.4e-3 PH 1.1683e5133.81.178e-3 PH+ LPA:1.1683e10,SameasPHSameasPH bandedLPA r o =0 : 5 r o , lam =2 lam ElsesameasPH PH+ LPA/RPA:sameascontrolSameasPHSameasPH replacement ElsesameasPH elasticitybasedonradiusvalue,aconstantvaluewasusedfortheentirenetworkjust likein[23].Thiswasimplementedbysetting k 3 equaltotheconstantelasticityvalue with k 1 ;k 2 =0. InordertotesttheeectofadecreasingradiusandpruningthatoccursinPH, the0Dcircuitmodelwasadjustedsuchthattheradiivalueswere20%decreasedas givenbytakingtheaverageofthevaluesshownin[36]andcomparingthechange betweenthecontrolandPHcases.Thereisnotmuchdatatoquantifytheamountof pruningthatoccursinPHpatients,soaconservative10%decreasein and was usedtoconstructthetree,whichcorrelatedtoabouta13.5%decreaseinthenumber ofdistalvessels. 2.5WallShearStressandPulsePressure WallshearstressWSSinthedistalvesselsisdeterminedbymultiplyingthe viscosity, ,bythewallshearrate,givingtheequationforWSSof w = 4 Q rA .The WSSwasnormalizedtothecontrolgroupdueinordertomakecomparisoneasier duetolargeresultingWSSvalues.Toinvestigatepulsatility,pulsepressurePPis investigated.PPiscalculatedbytakingthedierencebetweenthemaximumand 26

PAGE 43

minimumpressure.TheMATLABscriptusedtocalculateWSSandPPisshownin AppendixE. 27

PAGE 44

Figure2.7:Thenaltreemodelusedfor0Dsimulations.Formoreinformationon whateachcomponentorblockrepresents,refertogure2.6. 28

PAGE 45

3.Results 3.1ModelValidationandOptimization Asinglebifurcationwasusedforvalidationofbothmodels. 3.1.11DModelValidation Themodelwasvalidatedusingpreviouslypublisheddatafrom[23]forasingle bifurcationusingthevaluesshownforvessels0,1,and2shownintable2.3,where vessel0istheparentvesselwithtwodaughtervessels,vessels1and2asshownin gure3.1. Figure3.1:Exampletreeofasinglebifurcation. Thedatafrom[23]wasconvertedtotheVaMpy1Dmodelunitswhennecessary. AnNvalueof50andadtvalueof1e-7wereusedforthevalidationsimulations.The pressureandowresultsforthecontrolcaseisshowningure3.2.Thepressureand owresultsforthephcaseisshowningure3.3.Theindividualplotsshowninthe gurescanbefoundinAppendixA. Asshown,thepeakpressureforthecontrolis20mmHgjustasin[23].The PHmodelhasapeakpressureof40mmHgwhichisabithigherthan[23],butthis slightdierenceislikelyduetothefactthatthisvalidationusedanaverageofthe valuesreportedin[23]for R T , C T ,and Eh=r o asitwasunclearwhattheyusedas representative"data. 29

PAGE 46

Figure3.2:Controltestgroupresultsforthevalidationstudydoneonasinglebifurcation.ThetopplotsarepressuremmHgandthebottomplotsareowrate cm 3 /sec 30

PAGE 47

Figure3.3:PHtestgroupresultsforthevalidationstudydoneonasinglebifurcation. ThetopplotsarepressuremmHgandthebottomplotsareowrate cm 3 /sec 31

PAGE 48

ThisdataalsoshowsthatthetwoowrateplotsaddtoequalthatofVessel0as expected.ItalsoshowsthatVessel1hasagreaterowthanVessel2.Thisisdue tothefactthatVessel1hasalargerdiameterthanVessel2sohaslessresistanceto owandthereforehasagreaterowrate.The1Dmodelshowsexpectedbehavior fromthepressureandowrate,sothereisgoodcondenceinthemodel. 3.1.20DModelValidation The0DModelwasvalidatedinafewways.First,asingleRCcircuitwasmade andcomparedtothesamecircuitsolvedwithode45inMATLABandcomparedto ensurethemodelissolvingcorrectly.Theresultsareshowningure3.4andshow agreementbetweentechniques. Figure3.4:ValidationbetweenasingleRCcircuitinSimulinkversussolvingasingle RCcircuitinMATLABusingODE45. Second,asingleRCcircuitwascomparedtothatofa2-generationRCcircuit wherebothmodelstotaledtothesameresistanceandcapacitancetoensurethe modelcouldbeexpandedintobifurcationswiththeresultsshowningure3.5.These 32

PAGE 49

resultsingure3.5Ashowgoodagreementbetweenpressureoftheode45solvinga singleRCcircuitcomparedtothetotalpressureofthe2-generationRCmodelwhere thetotalresistanceandcapacitanceareequaltothesingleRCcircuit.Another validationshowningure3.5Bisthatthepressureofvessels1and2areequivalent astheyshouldbesincetheybelongtothesamegenerationjustasthevoltagedrop isequivalentbetweencircuitsinparallel. Figure3.5:Validationfora2generationRCcircuit.AshowsthetotalPressureof a2generationcircuitinSimulinkversussolvingtheequivalentsingleRCcircuitin MATLABusingODE45.Bshowsthepressuresforeachindividualvesselsolvedin Simulink. Asshowninbothgures3.4and3.5,themodeltakesafewcyclestoreacha 33

PAGE 50

Table3.1:The dx valuesforthelargest dx 0 andsmallest dx 4 vesselsizesfora givenNvalue. N dx 0 cm dx 4 cm 100.04560.00578 200.02160.00274 300.01410.00179 400.01050.00133 500.00840.00106 600.00690.00088 700.00590.00075 800.00520.00066 900.00460.00058 1000.00410.00053 5000.000820.00010 stablesolutionsoonlythelastcycleorlastfewcyclesareconsideredinsubsequent calculations. 3.1.3MeshIndependence MeshIndependencestudiesweredoneforthe1Dmodeltodeterminetheoptimum valuesfordxanddttorunsubsequentsimulations,whilestillsatisfyingtheCFL condition.Thesesimulationswererungivenatreewithanasymmetricgeometryas shownintable2.3.ThefollowingNvaluesweretestedwithadtof1e-7:10,20,30, 40,50,60,70,80,90,100,500withthecorrespondingdxvaluesforeachvesselshown intable3.1.Nrepresentsthenumberofdiscretizedspacestepstousetoallowfor dynamicdxvaluesbasedonthevessellength.Thetimeittookforeachsimulation torunisshownintable3.2. 34

PAGE 51

Table3.2:ThetotaltimeeachsimulationranforagivenNvalue. NTimehours 1016.86 2014.76 3015.55 4015.15 5014.97 6014.49 7014.32 8013.84 9013.95 10013.66 50015.68 35

PAGE 52

Changesinpeakpressurevaluesalongthelengthwasinvestigatedfortwovessels: thelargestvesselVessel0andthesmallestvesselVessel4.Theresultsforthe changesinpeakpressurealongthelengthofthevesselforeachNvalueisshownin AppendixB.ThemaximumrelativeerrorforeachNvalue,calculatedagainstthe largestNvalueof500,isshownhereingure3.6forVessel0andgure3.7forVessel 4.TheseguresshowthatthesolutionisconvergingwithincreasingNvalues. Figure3.6:MaximumRelativeErrorcomparedtoN=500forthelargestvessel, Vessel0,forvaryingN/dxvaluesaswellasthetotaltimeittookforeachsimulation torun. 36

PAGE 53

Figure3.7:MaximumRelativeErrorcomparedtoN=500forthesmallestvessel, Vessel4,forvaryingN/dxvaluesaswellasthetotaltimeittookforeachsimulation torun. Figures3.6and3.7alsoshowthetimeittookforeachsimulationtorun,which isalsoshownintable3.2.Forsomeofthecases,increasingtheNvalueactually decreasedthetimeittookforsomeofthesimulationstorun.Thisislikelydueto Newton'smethodrequiringlessiterationstoconvergetoasolutionduetoitsstarting guessfromtheprevioustimestepbeingcloserthanitwouldbeforasmallerNvalue. TheresultsofthemeshindependencestudydonewithvaryingN/dxvalues showedthatanyNvaluecouldbeusedduetoacouplereasons.First,PHisdiagnosed within1mmHg,andguresB.1andB.2inAppendixBshowthatevenwiththe smallestNvalue,thepeakpressurewouldbewithin0.5mmHgorless.Second,the maximumrelativeerrorisalsolessthan0.5%asshowningures3.6and3.7whilea pulmonaryarterycatheterhasbeenshowntohavearandomerroraround5%anda precisionerrortobewithin10-20%[37].Therefore,thesimulationsobtainanerror lessthanthatoftheinstrumentthatprovidestheinletboundarycondition.Alarger Nvalueakaasmallerdxismoredesirable,especiallyduetothefactthatalarger 37

PAGE 54

Table3.3:Thetotaltimeeachsimulationranforagivendtvalue. dtsecTimehours 1e-61.366 1e-713.66 1e-8157.2 Nvalueledtoashorterruntime,and,thus,anNvalueof100wasdecidedonfor subsequentsimulations. Next,meshindependencestudywasperformedonvaryingdtvalueswithN= 100aschosenfromthedxmeshindependencestudy.DuetotheCFLcondition,the smallestdtvaluethatcouldbeusedwas1e-6andforthesakeoftime,onlythreedt valueswerestudied:1e-6,1e-7and1e-8.Thetimeittookforeachsimulationtorun isshownintable3.3.Theresultsforthechangesinpeakpressurealongthelength ofthevesselforeachdtvalueisshowninAppendixB.Themaximumrelativeerror foreachdtvalue,calculatedagainstthelargestdtvalueof1e-8,isshowningure 3.8forVessel0andgure3.9forVessel4.Theseguresshowthatthesolutionis convergingwithincreasingdtvalues.Asshowninbothguresandintable3.3,a decreasingdtvalueleadstoanincreasingruntime,asexpected. 38

PAGE 55

Figure3.8:MaximumRelativeErrorcomparedtodt=1e-8forthelargestvessel, Vessel0,forvaryingdtvaluesaswellasthetotaltimeittookforeachsimulationto run. Figure3.9:MaximumRelativeErrorcomparedtodt=1e-8forthesmallestvessel, Vessel4,forvaryingdtvaluesaswellasthetotaltimeittookforeachsimulationto run. Again,weseethatthepeakpressureiswithin0.5mmHgseeAppendixBfor 39

PAGE 56

plotsorlessandthattherelativeerrorislessthan0.025%;therefore,anydtvalue couldbeused.Again,asmallerdtvalueisdesirable,soadtvalueof1e-7waschosen inordertokeepasmalldtwithareasonableruntimeastheruntimeforthe1e-8 simulationtook157.2hoursorroughly6.5dayswhichisanunreasonableamountof time. 3.2Simulations AsstatedintheMethodssection,fourtestgroupswerestudied:controlrepresentingnon-PH,PH,PHwithLPAbandingPH+LPAB,andPHwithasynthetic/compliantLPAandRPAPH+replacement. 3.2.11DSimulationTime DuetocomputationalcostofoneofthetestgroupsPH+LPAB,the1D simulationswererunusingatreewithasymmetricgeometrywiththeradiiand length-to-radiusvaluesshownintable2.4.AnNvalueof46wasusedasthiswas thelargestNvaluethatwouldsatisfytheCFLconditionforadt=1e-7tokeepthe computationalcostdown.ThisNvaluewaskeptconsistentacrossallsimulations. Thetimeittookeachsimulationtorunisshownintable3.4.Othervariablesand datausedforthesesimulationscanbefoundintablesinthemethodssection.Vessel 0representsthemainpulmonaryarteryMPA,vessels1and2arethetwodaughter vesselsofvessel0andrepresenttheleftpulmonaryarteryLPAandtherightpulmonaryarteryRPA,respectively.Vessels3and4arethedaughtervesselsforvessel 1andvessels5and6arethedaughtervesselsforvessel2.Thisisshowngraphically ingure2.2.Theplotsforpressureandowshowatimeof0.34to0.42duetothe timeperiodforeachowcycleandthatthesimulationsranfor4cyclesinorderto obtainamorestablesolution. 3.2.2Pressure 40

PAGE 57

Table3.4:Thetotaltimeeachsimulationranforthefourtestgroups. TestGroupTimehours Control16.337 PH14.523 PH+LPAB36.373 PH+replacement14.406 ThegraphsforpressurealongthevesselnetworkfortheControl,PH,PH+LPAB, PH+replacementareshowningures3.10,3.11,3.12,and3.13,respectively.Full graphsofthepressurewavescanbefoundinAppendixC.Consistentlyamong allofthevessels,thecontrolgrouphasapeakpressurearound20mmHg,thePH andPH+LPABgroupstendtohaveapeakpressurearound40mmHg,whilethe PH+replacementhasaslightlylowerpeakpressureofaround32mmHg.Theplot forpressureoftheinitialvesselofthe0Dmodelakatheterminalvesselofthe1D modelforalltestgroupsisshowningure3.14.Aplotshowingthedistalpressure plotsforallofthetestgroupsisshowningure3.15. 3.2.3Flow ThegraphsforowalongthevesselnetworkfortheControl,PH,PH+LPAB, PH+replacementareshowningures3.16,3.17,3.18,and3.19,respectively.Full graphsofthepressurewavescanbefoundinAppendixD. Aplotforowoftheinitialvesselofthe0Dmodelakatheterminalvesselof the1Dmodelforalltestgroupsisshowningure3.20.Aplotshowingthedistal owplotsforallofthetestgroupsfromthe0Dsimulationsisshowningure3.21. 3.3DistalWallShearStressandPulsePressure WallshearstressWSSandpulsepressurePPwerecalculatedforjustthe terminalvesselsincethegoalofthisprojectwastoinvestigatethedistalmechanical 41

PAGE 58

Figure3.10:Controltestgrouppressureresultsfortheentirevasculartree.Each vesselatthesamegenerationhasthesameresultsduetobeingcompletelysymmetric, soonlyonevesselpergenerationisshownhereforsimplicity. 42

PAGE 59

Figure3.11:PHtestgrouppressureresultsfortheentirevasculartree.Eachvessel atthesamegenerationhasthesameresultsduetobeingcompletelysymmetric,so onlyonevesselpergenerationisshownhereforsimplicity. 43

PAGE 60

Figure3.12:PH+LPABtestgrouppressureresultsfortheentirevasculartree.Each vesselatthesamegenerationhasthesameresultsduetobeingcompletelysymmetric, soonlyonevesselpergenerationisshownhereforsimplicity. 44

PAGE 61

Figure3.13:PH+replacementtestgrouppressureresultsfortheentirevascular tree.Eachvesselatthesamegenerationhasthesameresultsduetobeingcompletely symmetric,soonlyonevesselpergenerationisshownhereforsimplicity. 45

PAGE 62

Figure3.14:Pressureresultsfortheinitialvesselofallofthetestgroupsofthe0D modelsimulations. Figure3.15:Pressureresultsfortheterminalvesselsofallofthetestgroupsofthe 0Dmodelsimulations. 46

PAGE 63

Figure3.16:Controltestgrouppressureresultsfortheentirevasculartree.Each vesselatthesamegenerationhasthesameresultsduetobeingcompletelysymmetric, soonlyonevesselpergenerationisshownhereforsimplicity. 47

PAGE 64

Figure3.17:PHtestgrouppressureresultsfortheentirevasculartree.Eachvessel atthesamegenerationhasthesameresultsduetobeingcompletelysymmetric,so onlyonevesselpergenerationisshownhereforsimplicity. 48

PAGE 65

Figure3.18:PH+LPABtestgrouppressureresultsfortheentirevasculartree.Each vesselatthesamegenerationhasthesameresultsduetobeingcompletelysymmetric, soonlyonevesselpergenerationisshownhereforsimplicity. 49

PAGE 66

Figure3.19:PH+replacementtestgrouppressureresultsfortheentirevascular tree.Eachvesselatthesamegenerationhasthesameresultsduetobeingcompletely symmetric,soonlyonevesselpergenerationisshownhereforsimplicity. 50

PAGE 67

Figure3.20:Flowresultsfortheinitialvesselofallofthetestgroupsofthe0Dmodel simulations. Figure3.21:Flowresultsfortheinitialvesselofallofthetestgroupsofthe0Dmodel simulations. 51

PAGE 68

Figure3.22:WSSresultsforthe0Dmodelsimulationsforaterminalvesselgiven changestoproximalvesselstiness. changes. 3.3.1EectofproximalstieningondistalWSSandPP AbarplotoftheWSS,normalizedtothecontrolWSS,foreachtestgroupis showningure3.22givenchangestoproximalvesselstiness.Onlydatafromthe lastcycleortwowasincludedinthecalculationsforWSSandPPtoensureaccurate data.AbarplotofthePPisshowningure3.23forallofthetestgroupsgiven changestoproximalvesselstiness. 3.3.2DierencesinvesselgeometryimpactdistalWSSandPP AbarplotoftheWSSforeachtestgroupisshowningure3.24oncechanges duetovesselgeometrywereincluded.AbarplotofthePPisshowningure3.23 forallofthetestgroups,whichincludeschangesduetovesselgeometry. 52

PAGE 69

Figure3.23:Pulsepressureresultsforthe0Dmodelsimulationsforaterminalvessel givenchangestovesselstiness. Figure3.24:WSSresultsforthe0Dmodelsimulationsforaterminalvesselgiven bothchangestoproximalvesselstiness,andchangestodistalvesselgeometry. 53

PAGE 70

Figure3.25:Pulsepressureresultsforthe0Dmodelsimulationsforaterminalvessel givenbothchangestoproximalvesselstiness,andchangestodistalvesselgeometry. 54

PAGE 71

4.Discussion 4.1PressureWaves Inthe1Dmodel,thepressurechangedmostlyasexpected.ThePHcaseshad higherpeakpressurethanthecontrol,despitethePHhavinglowerowrate.There wasalsoaslightincreaseinpeakpressureforthePH+LPABmmHgcompared tothePHcaseinVessel0,butthismakessensebecausebandingcausedadecreased radiusandincreasedstinessinvessel1LPAwhichleadstoincreasedresistance andthereforeincreasedpeakpressure.Thepressureplotinvessel1ofthePH+ LPABshowedpressuredropoalongthelengthofthatvessel,whilevessel2stayed atthehigherpeakpressureof45mmHg.Thisisexpectedaslessowisexpected togothroughthesmallervesselsize.VesselsdownstreamoftheLPAinthePH+ LPABshowedacontinueddecreaseinpressureandshowedequivalentpressuresto thatofthePH+replacementcase.ForthePH+replacementcase,therewasabout a20%decreaseinpeakpressureforvessels1and2.Thepeakpressuredidnotreturn tocontrolvalues;however,thePHcasedidhaveaslightlyhigherpeakpressurethan expected.Ifthepeakpressurehadbeenclosertowhatwasshownin[23],which wasapproximately35mmHg,thena20%decreasewouldhaveputthepeakpressure to28mmHg.Whilethispeakpressureisclosertocontrolvalues,itwouldstillbe classiedasPH;however,overtimethisvaluecoulddecreaseevenmoredepending onchangesinanatomyandphysiologyofthedistalvasculatureduetoadecreased pressure. OneinterestingdierenceinthepressurewaveofthePH+replacementcase isthattheend-cyclepressurewasincreasedcomparedtothePHdespitehavinga lowerpeakpressure.Theloweredend-cyclepressurealsocontributetothePH+ replacementhavingalowerpulseintensity.Asinthedierencebetweenmaximum andminimumpressureisonly10mmHgwhichisclosetothecontrolcasemmHg whereasthePHandPH+LPABareallaround20mmHg.Thistrendcontinuesas 55

PAGE 72

youmoveintothedistalvasculatureofthe0Dmodelandcanbeseeninthepressure wavesingure3.15aswellasthepulsepressuremeasurementsshowningure3.23. ThepulsepressureforthePH+replacement;however,doesnotdecreasetoasmuch asthecontrolcaseinthedistalarterioleswhichisduetothedistalarteriolesbeing muchmorecompliantinthecontrolcaseallowingformoreattenuation. 4.2FlowWaves Theowalsochangedasexpectedinthe1Dmodel.ForthePH+LPAB,gure 3.18showsvessel1LPAhavingdecreasedowwhilevessel2RPAhasincreased ow,whichiswhatwasseenin[27].Thesechangesgivegoodcondencetothe modelgivingresultsthatarecomparabletoananimalmodel.Thetrendsinowfor thiscaseisperpetuatedtothevesselsdownstreamfromvessels1and2asshownin subsequentvesselplots.Anothertrendnoticedinvessel1comparedtovessel2inthe PH+LPABcase,isthattheleftdaughtervesselsshowasmootherowwavethan theleftdaughtervesselsasshowningures3.20and3.21. InthePH+replacementcase,theowwaveseemstobecomelesspulsatile andmoresustained.Figure3.19showsthepulseofthesecondgenerationofvessels diminishingasittravelsdownthelengthofthevessel.The1Dterminalvesselsshowa moresustainedandlesspulsatileowovertimewhilethewaveseemsmoreconsistent alongthelengthofthevessels.Thistrendcontinuesforthedistalvesselsofthe0D modelandcanbeseenintheowwavesshowningure3.21bytheelongatedow wavecomparedtothecontrolandPHcases. 4.3DistalWallShearStressandPulsePressure WallshearstressWSS,showningure3.22,inthePHcasedecreasedinresponsetoincreasedvesselstiening;however,thiswasseenwithoutanychangesto geometryasseenwithPH.Withthedecreasedradiusandpruningbeingmodeled, showningure3.24,WSSinthePHcasewasincreasedasexpectedandgivesevidencethattheincreaseinWSSispurelyduetogeometricchanges.Withtheradius 56

PAGE 73

changeandpruningbeingmodeled,theWSSforthePH+replacementcasereturned closertocontrollevels.ThebandingoftheLPAalsodecreasedWSSintheleftlung arterieswhereastheWSSfortherightlungarterieswasslightlyincreasedcompared toPH.Duetothesechangesandtheresultsshownin[27],PPandWSSarelikely thecauseoftheanatomicalchangesseenin[27]. ThepulsepressurePP,showningure3.23,showsexpectedtrends.Just accountingforchangestovesselstiness,thePPforthePHcaseisincreasedto aboutfourtimesthecontrol.Thisincreasewasseeninthe1Dmodelpulsewaves, butwasn'tasstarkofadierence.Again,duetothedownstreamvesselsbeingless compliantinthePHcase,thosepressurewaveshavelessattenuationwhichiswhythe PPwouldbemuchhigher.ThePH+LPABLandRcasesbothshowincreases toPP,buttheleftsideshowedaslightdecreasecomparedtoPHwhiletherightside hadaslightincreasecomparedtoPH.Thisisexpectedastherewaslessowgoing intotheleftside,butthedierenceisveryminimalcomparedtothePPofthePH case.ThePH+replacementshowsdecreasedPPbacktowardsthecontrollevels,but isstilltwicethatofthecontrol.PerhapsupdatingthePH+replacementcompliance tobemoreelasticthanthecontrolcouldfurtherdecreasethispulsepressurecloser tonormal.Whenpruninganddecreasedradiusvalueswereincluded,thesetrends werefurtherexacerbatedasshowningure3.25. 4.4LimitationsandFutureWork Thismodelusedthesamegeometryforproximalvesselsbetweenallgroupsexcept forthecaseofthePH+LPAB.Thismodelalsoonlytestedforasingle,generalized mousegeometryandagreatextensionofthismodelwouldbetodoanimalstudies usinganimal-specicdataforeachsimulation.Withananimalstudy,thismodel couldalsobeusedinconjunctionwithhistologicaldatatobetterconnecttheuid dynamicswithanatomicalchangestothedistalbloodvesselsovertime. Otherstudiesthatcouldbedoneusingthissimulation,wouldbetotestvarying 57

PAGE 74

stinessesofreplacingtheLPAandRPAtodetermineiftheremightbeamore optimalstiness.Inthismodel,thesamestinessascontrolvesselswastested; however,thechangestoPPandWSSwerenotasclosetocontrollevelsasdesired whichmaymeanthestinessofthereplacementvesselwouldneedtobeadjusted. Thismodelalsoonlyinvestigatedonesubsetofthelunganddidsimulationsusing atreewithsymmetricvesselsforsimplicity;however,itwouldbeinterestingtolook atatreewithanasymmetricgeometrytoinvestigatewhatchangesoccurandifthis isaviablesimplicationforfuturemodeluse. 58

PAGE 75

5.Conclusion PulmonaryhypertensionPHisassociatedwithchangestoanatomyandphysiologyofbloodvessels[4].Morespecically,PHvesselsincreaseinstiness,decrease inradiusandtherefore,increaseinvascularresistance.ThesechangescauseanincreaseinowpulsatilityandwallshearstressWSSwhichfurtherexacerbatethese problems[7],andtherefore,pulsatilityandWSSareimportantfeaturestostudy. Unfortunately,thesefeaturescannotbestudiedindistalvesselsduetoequipment constraints[9],whichiswherecomputationalmodelscomeintoplay.Thegoalsof thisprojectweretodevelopacomputationalmodelofalarge-scalepulmonaryvasculaturenetworkaswellastoinvestigatepulsatility,usingpulsepressurePP,and WSSindistalvasculaturegivenchangestotheproximalpulmonaryvasculature.More specically,simulationofabandedleftpulmonaryarteryLPABinthepresenceof PHandsimulationofareplacementleftandrightpulmonaryarteryinthepresence ofPHwerecomparedtohealthyandPHcontrols. Thecomputationalmodelwasdevelopedusinga1Dand0Dhybridmodelusing anexistingpython-based1Dmodelasthebase.Themodelwasexpandeduponto allowforinvestigationofthedesiredtestgroupsdescribed,buthadahighcomputationalcost.Thus,a0Dmodelwasaddedtocreatethehybridmodelwherethe 0DmodelusedSimulinktocreateacircuit-basedanalogofthedistalvasculature. ThemodelwasthenusedtoinvestigatePPandWSSinthedistalvessels.PPwas increasedforPHcomparedtothehealthycontrol,asexpected.ThePH+LPAB LshowedaslightdecreasecomparedtoPHwhilethePH+LPABRshowedan increasedPPbeyondthatofthePHcase.TheleftandrightpulmonaryarteryreplacementgroupreturnedPPclosertocontrollevels,whichshowedgoodpromisefor asyntheticleftandrightpulmonaryarteryasatreatmentoptionforPH.TheWSS forPHwasalsoincreasedcomparedtothehealthycontrol,butonlywhenpruning anddecreasedvesselareaweremodeledprovidingevidencethatchangestoWSSare 59

PAGE 76

completelyduetogeometryandnotvesselmechanics.InthePH+LPABcase,the downstreamarteriesofthebandedarteryshowedadecreaseinWSScomparedto PHwhilethearteriesdownstreamfromthenon-bandedarteryshowedanincreasein WSSbeyondthatofthePHcase.PH+LPABwasshowntodecreaseinammatory markersofPH[27]andthismodelshowedthatthePH+replacementgroupshowed similarWSSandPPtrendsasthePH+LPABcasesimulated.Theseresultssuggest thatreplacingtheleftandrightpulmonaryarterieswithsynthetic,compliantvessels couldbeaviabletreatmentoptionandshouldbefurtherinvestigated.Onewaythis couldbeinvestigated,isbytestingthistreatmentoptioninananimalstudyandrunninganimalspecicsimulationssimultaneously.Despitemanysimplications,this modelwasabletoprovideaccuratepressureandowresultsinpulmonaryvessels andisagreattoolforfutureapplicationandinvestigation.Inconclusion,themodel wassuccessfullycreatedandusedtoinvestigatePPandWSSindistalvessels,and illustratedtheimportanceofmechanicalstimuliinthevascularremodelingprocess ofprogressivepulmonaryvasculardisease. 60

PAGE 77

REFERENCES [1]M.Rabinovitch.Molecularpathogenesisofpulmonaryarterialhypertension. TheJournalofClinicalInvestigation ,122:4306{4313,122012. [2]G.Simonneau,M.A.Gatzoulis,I.Adatia,D.Celermajer,C.Denton, A.Ghofrani,M.A.GomezSanchez,K.R.Krishna,M.Landzberg,R.F.Machado, H.Olschewski,I.M.Robbins,andR.Souza.Updatedclinicalclassicationofpulmonaryhypertension. JournaloftheAmericanCollegeofCardiology ,62:D34{41, 2013. [3]W.SunandS.Y.Chan.Pulmonaryarterialstiness:Anearlyandpervasive driverofpulmonaryarterialhypertension. Frontiersinmedicine ,5:204{204, 2018. [4]K.R.Stenmark,N.Davie,M.Frid,E.Gerasimovskaya,andM.Das.Roleofthe adventitiainpulmonaryvascularremodeling. Physiology ,21:134{145,2006. [5]K.S.Hunter,S.R.Lammers,andR.Shandas.Pulmonaryvascularstiness: measurement,modeling,andimplicationsinnormalandhypertensivepulmonary circulations. ComprehensivePhysiology ,1:1413{1435,2011. [6]S.Malenfant,A.Neyron,R.Paulin,F.Potus,J.Meloche,S.Provencher,and S.Bonnet.Signaltransductioninthedevelopmentofpulmonaryarterialhypertension. Pulmonarycirculation ,3:278{293,2013. [7]W.Tan,K.Madhavan,K.S.Hunter,D.Park,andK.R.Stenmark.Vascular stieninginpulmonaryhypertension:causeorconsequence?groverconferenceseries. Pulmonarycirculation ,4:560{580,2014. [8]J.C.Grignola.Hemodynamicassessmentofpulmonaryhypertension. World journalofcardiology ,3:10{17,2011. [9]J.Alastruey,NXiao,H.Fok,To.Schaeter,andC.A.Figueroa.Ontheimpact ofmodellingassumptionsinmulti-scale,subject-specicmodelsofaortichaemodynamics. JournaloftheRoyalSociety,Interface ,13:20160073,2016. [10]K.Hassani,M.Navidbakhsh,andM.Rostami.Modelingoftheaortaartery aneurysmsandrenalarterystenosisusingcardiovascularelectronicsystem. Biomedicalengineeringonline ,6:22{22,2007. [11]O.Ghasemalizadeh,M.R.Mirzaee,B.Firoozabadi,andK.Hassani.Exactmodelingofcardiovascularsystemusinglumpedmethod. InternationalConference onBioinformaticsandComputationalBiology ,pages408{417,2008. 61

PAGE 78

[12]M.S.Olufsen.Structuretreeoutowconditionforbloodowinlargersystemic arteries. AmericanJournalofPhysiology ,276:H257{68,1999. [13]M.S.Olufsen.Aone-dimensionaluiddynamicmodelofthesystemicarteries. StudiesinHealthTechnologiesandInformatics ,71:79{98,2000. [14]G.D.A.Qureshi,M.U.andVaughan,C.Sainsbury,M.Johnson,C.S.Peskin,M.S. Olufsen,andN.A.Hill.Numericalsimulationofbloodowandpressuredropin thepulmonaryarterialandvenouscirculation. Biomechanicsandmodelingin mechanobiology ,13:1137{1154,2014. [15]J.Flores,J.Alastruey,andE.CorveraPoire.Anovelanalyticalapproachto pulsatilebloodowinthearterialnetwork. AnnalsofBiomedicalEngineering , 44:3047{3068,2016. [16]J.P.MynardandJ.J.Smolich.One-dimensionalhaemodynamicmodelingand wavedynamicsintheentireadultcirculation. AnnalsofBiomedicalEngineering , 43:1443{1460,2015. [17]P.Reymond,F.Merenda,F.Perren,D.Rufenacht,andN.Stergiopulos.Validationofaone-dimensionalmodelofthesystemicarterialtree. AmericanJournal ofPhysiology-HeartandCirculatoryPhysiology ,297:H208{H222,2009. [18]JonasovaA.,BublikO.,andVimmrJ.Acomparativestudyof1dand3d hemodynamicsinpatient-specichepaticportalveinnetworks. Appliedand ComputationalMechanics ,8:177{186,2014. [19]L.Grinberg,E.Cheever,T.Anor,J.R.Madsen,andG.E.Karniadakis.Modeling bloodowcirculationinintracranialarterialnetworks:Acomparative3d/1d simulationstudy. AnnalsofBiomedicalEngineering ,39:297{309,2011. [20]P.Reymond,P.Crosetto,S.Deparis,A.Quarteroni,andN.Stergiopulos.Physiologicalsimulationofbloodowintheaorta:Comparisonofhemodynamic indicesaspredictedby3-dfsi,3-drigidwalland1-dmodels. MedicalEngineering&Physics ,35:784{791,2013. [21]N.Xiao,J.Alastruey,andC.AlbertoFigueroa.Asystematiccomparisonbetween1-dand3-dhemodynamicsincompliantarterialmodels. International journalfornumericalmethodsinbiomedicalengineering ,30:204{231,2014. [22]A.Quarteroni,A.Veneziani,andC.Vergara.Geometricmultiscalemodelingof thecardiovascularsystem,betweentheoryandpractice. ComputerMethodsin AppliedMechanicsandEngineering ,302:193{252,2016. [23]M.U.Qureshi,M.J.Colebank,L.M.Paun,L.E.Fix,N.Chesler,M.A.Haider, N.A.Hill,D.Husmeier,andM.S.Olufsen.Hemodynamicassessmentofpulmonaryhypertensioninmice:amodel-basedanalysisofthediseasemechanism. BiomechanicsandModelinginMechanobiology ,18:219{243,2019. 62

PAGE 79

[24]N.Bessonov,A.Sequeira,S.Simakov,Y.Vassilevskii,andV.Volpert.Methods ofbloodowmodeling. MathematicalModelingofNaturalPhenomena ,11:1{ 25,2016. [25]Y.Shi,P.Lawford,andR.Hose.Reviewofzero-dand1-dmodelsofbloodow inthecardiovascularsystem. Biomedicalengineeringonline ,10:33,2011. [26]E.Boileau,P.Nithiarasu,P.J.Blanco,L.O.Muller,F.E.Fossan,L.R.Hellevik, W.P.Donders,W.Huberts,M.Willemet,andJ.Alastruey.Abenchmarkstudy ofnumericalschemesforone-dimensionalarterialbloodowmodelling. InternationalJournalforNumericalMethodsinBiomedicalEngineering ,31:e02732, 2015. [27]K.Abe,M.Shinoda,M.Tanaka,Y.Kuwabara,K.Yoshida,Y.Hirooka,I.F. McMurtry,M.Oka,andK.Sunagawa.Hemodynamicunloadingreversesocclusivevascularlesionsinseverepulmonaryhypertension. CardiovascularResearch , 1:16{25,2016. [28]B.T.Tang,S.S.Pickard,F.P.Chan,P.S.Tsao,C.A.Taylor,andJ.A.Feinstein. Wallshearstressisdecreasedinthepulmonaryarteriesofpatientswithpulmonaryarterialhypertension:Animage-based,computationaluiddynamics study. Pulmonarycirculation ,2:470{476,2012. [29]M.Schafer,V.O.Kheyfets,J.D.Schroeder,J.Dunning,R.Shandas,J.K.Buckner,J.Browning,K.S.Hertzberg,J.andHunter,andB.E.Fenster.Mainpulmonaryarterialwallshearstresscorrelateswithinvasivehemodynamicsand stinessinpulmonaryhypertension. Pulmonarycirculation ,6:37{45,2016. [30]A.K.DiemandN.W.Bresslof.Vampy:Apythonpackagetosolve1dbloodow problems. JournalOfOpenResearchSoftware ,5,2017. [31]V.BKolachalama,N.W.Bresslof,P.B.Nair,etal.Predictivehaemodynamicsin aone-dimensionalcarotidarterybifurcation.parti:Applicationtostentdesign. IEEETransactionsonBiomedicalEngineering ,54:802{812,2007. [32]M.S.Olufsen,C.S.Peskin,W.Y.Kim,E.M.Pedersen,A.Nadim,andJ.Larsen. Numericalsimulationandexperimentalvalidationofbloodowinarterieswith structured-treeoutowconditions. AnnalsofBiomedicalEngineering ,28:1281{ 1299,2000. [33]A.Mishra,F.M.O'Farrell,C.Reynell,N.B.Hamilton,C.N.Hall,andD.Attwell. Imagingpericytesandcapillarydiameterinbrainslicesandisolatedretinae. NatureProtocols ,9:323EP,2014. [34]M.P.Wiedeman.Dimensionsofbloodvesselsfromdistributingarterytocollectingvein. CirculationResearch ,pages375{378,1963. [35]M.I.Townsley.Structureandcompositionofpulmonaryarteries,capillaries,and veins. ComprehensivePhysiology ,2:675{709,2012. 63

PAGE 80

[36]N.Rol,E.M.Timmer,T.J.C.Faes,A.VonkNoordegraaf,K.Grunberg,H.Bogaard,andN.Westerhof.Vascularnarrowinginpulmonaryarterialhypertension isheterogeneous:rethinkingresistance. Physiologicalreports ,5:e13159,2017. [37]X.X.Yang,L.A.Critchley,andG.M.Joynt.Determinationoftheprecisionerror ofthepulmonaryarterythermodilutioncatheterusinganinvitrocontinuousow testrig. AnesthesiaandAnalgesia ,112:70{77,2011. 64

PAGE 81

APPENDIXA.ValidationPlotsfor1DModel A.1Control A.1.1PressurePlots FigureA.1:ValidationPressureResultsoftheControltestgroupfortherootvessel 0MPA 65

PAGE 82

FigureA.2:ValidationPressureResultsoftheControltestgroupforvessel1LPA FigureA.3:ValidationPressureResultsoftheControltestgroupforvessel2RPA 66

PAGE 83

A.1.2FlowPlots FigureA.4:ValidationFlowResultsoftheControltestgroupfortherootvessel0 MPA 67

PAGE 84

FigureA.5:ValidationFlowResultsoftheControltestgroupforvessel1LPA FigureA.6:ValidationFlowResultsoftheControltestgroupforvessel2RPA 68

PAGE 85

A.2PH A.2.1PressurePlots FigureA.7:ValidationPressureResultsofthePHtestgroupfortherootvessel0 MPA 69

PAGE 86

FigureA.8:ValidationPressureResultsofthePHtestgroupforvessel1LPA FigureA.9:ValidationPressureResultsofthePHtestgroupforvessel2RPA 70

PAGE 87

A.2.2FlowPlots FigureA.10:ValidationFlowResultsofthePHtestgroupfortherootvessel0MPA 71

PAGE 88

FigureA.11:ValidationFlowResultsofthePHtestgroupforvessel1LPA FigureA.12:ValidationFlowResultsofthePHtestgroupforvessel2RPA 72

PAGE 89

APPENDIXB.MeshIndependencePlots B.1ChangestoPeakPressureforvaryingN/dx FigureB.1:ChangestoPeakPressurealongthelengthofthelargestvessel,Vessel 0,duetovaryingN/dxvalues. 73

PAGE 90

FigureB.2:ChangestoPeakPressurealongthelengthofthesmallestvessel,Vessel 4,duetovaryingN/dxvalues. B.2ChangestoPeakPressureforvaryingdt FigureB.3:ChangestoPeakPressurealongthelengthofthelargestvessel,Vessel 0,duetovaryingdtvalues. 74

PAGE 91

FigureB.4:ChangestoPeakPressurealongthelengthofthesmallestvessel,Vessel 4,duetovaryingdtvalues. 75

PAGE 92

APPENDIXC.ModelPressurePlots C.1Control FigureC.1:PressureResultsoftheControltestgroupfortherootvesselMPA 76

PAGE 93

FigureC.2:PressureResultsoftheControltestgroupforthesecondgenerationof vesselsLPA/RPA 77

PAGE 94

FigureC.3:PressureResultsoftheControltestgroupforthe1Dterminalvessels Generation3. 78

PAGE 95

FigureC.4:PressureResultsoftheControltestgroupforthe0Dterminalvessels Generation12. C.2PH 79

PAGE 96

FigureC.5:PressureResultsofthePHtestgroupfortherootvesselMPA. 80

PAGE 97

FigureC.6:PressureResultsofthePHtestgroupforthesecondgenerationofvessels LPA/RPA. 81

PAGE 98

FigureC.7:PressureResultsofPHtestgroupforthe1DterminalvesselsGeneration 3. 82

PAGE 99

FigureC.8:PressureResultsofthePHtestgroupforthe0DterminalvesselsGeneration12. C.3PH+LPAB 83

PAGE 100

FigureC.9:PressureResultsofthePH+LPABtestgroupfortherootvesselMPA. 84

PAGE 101

FigureC.10:PressureResultsofthePH+LPABLtestgroupforthesecondgenerationofvesselsLPA. 85

PAGE 102

FigureC.11:PressureResultsofthePH+LPABRtestgroupforthesecondgenerationofvesselsRPA. 86

PAGE 103

FigureC.12:PressureResultsofthePH+LPABLtestgroupforthe1Dterminal vesselsthatdescendfromtheLPAGeneration3. 87

PAGE 104

FigureC.13:PressureResultsofthePH+LPABRtestgroupforthe1Dterminal vesselsthatdescendfromtheRPAGeneration3. 88

PAGE 105

FigureC.14:PressureResultsofthePH+LPABLtestgroupforthe0Dterminal vesselsGeneration12. 89

PAGE 106

FigureC.15:PressureResultsofthePH+LPABRtestgroupforthe0Dterminal vesselsGeneration12. C.4PH+replacement 90

PAGE 107

FigureC.16:PressureResultsofthePH+replacementtestgroupfortherootvessel MPA. 91

PAGE 108

FigureC.17:PressureResultsofthePH+replacementtestgroupforthesecond generationofvesselsLPA/RPA. 92

PAGE 109

FigureC.18:PressureResultsofthePH+replacementtestgroupfortheterminal vesselsofthe1Dmodel. 93

PAGE 110

FigureC.19:PressureResultsofthePH+replacementtestgroupforthe0Dterminal vesselsGeneration12. 94

PAGE 111

APPENDIXD.ModelFlowPlots D.1Control FigureD.1:FlowResultsoftheControltestgroupfortherootvesselMPA. 95

PAGE 112

FigureD.2:FlowResultsoftheControltestgroupforthesecondgenerationofvessels LPA/RPA. FigureD.3:FlowResultsoftheControltestgroupforthe1DterminalvesselsGeneration3. 96

PAGE 113

FigureD.4:FlowResultsoftheControltestgroupforthe0DterminalvesselsGeneration12. D.2PH FigureD.5:FlowResultsofthePHtestgroupfortherootvesselMPA. 97

PAGE 114

FigureD.6:FlowResultsofthePHtestgroupforthesecondgenerationofvessels LPA/RPA. FigureD.7:FlowResultsofPHtestgroupforthe1DterminalvesselsGeneration 3. 98

PAGE 115

FigureD.8:FlowResultsofthePHtestgroupforthe0DterminalvesselsGeneration 12. D.3PH+LPAB FigureD.9:FlowResultsofthePH+LPABtestgroupfortherootvesselMPA. 99

PAGE 116

FigureD.10:FlowResultsofthePH+LPABLtestgroupforthesecondgeneration ofvesselsLPA. FigureD.11:FlowResultsofthePH+LPABRtestgroupforthesecondgeneration ofvesselsRPA. 100

PAGE 117

FigureD.13:FlowResultsofthePH+LPABRtestgroupforthe1Dterminal vesselsthatdescendfromtheRPAGeneration3. FigureD.12:FlowResultsofthePH+LPABLtestgroupforthe1Dterminal vesselsthatdescendfromtheLPAGeneration3. 101

PAGE 118

FigureD.14:FlowResultsofthePH+LPABLtestgroupforthe0Dterminal vesselsGeneration12. 102

PAGE 119

FigureD.15:FlowResultsofthePH+LPABRtestgroupforthe0Dterminal vesselsGeneration12. D.4PH+replacement 103

PAGE 120

FigureD.16:FlowResultsofthePH+replacementtestgroupfortherootvessel MPA. FigureD.17:FlowResultsofthePH+replacementtestgroupforthesecondgenerationofvesselsLPA/RPA. 104

PAGE 121

FigureD.18:FlowResultsofthePH+replacementtestgroupforthe1Dterminal vesselsGeneration3. FigureD.19:FlowResultsofthePH+replacementtestgroupforthe0Dterminal vesselsGeneration12. 105

PAGE 122

APPENDIXE.MATLABCode E.1VesselNetworkGenerationScript %closeall; %clearall; %clc; %%BuildupTree %DiameterofMPA %Diameterofterminalvessels %ifoneofthetwodaughtervesselsislessthantheterminalvessel %diameter,theparentbecomesaterminalvesselsindicatedby0fordau %Matrixpositionsetupofradiusofvesselswhereif %ifparent=row,col:daughters=row+1,2 * col )]TJ0 g 0 G.133 .545 .133 rg .133 .545 .133 RG0 g 0 G.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(1&2 * col %ifdau=row,col:parent=row )]TJ/F51 9.9626 Tf 7.395 0 Td [(1,floorcol+1/2 %function[d1,d2]=get daucol %d1=2 * col )]TJ0 g 0 G.133 .545 .133 rg .133 .545 .133 RG0 g 0 G.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(1; %d2=2 * col; %end %%Inputs r root=0.0215;%rootradiusakatheendofthe1Dmodel )]TJ0 g 0 G.133 .545 .133 rg .133 .545 .133 RG0 g 0 G.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(cm r term=0.0005;%terminalradiusakacapillaryradius )]TJ0 g 0 G.133 .545 .133 rg .133 .545 .133 RG0 g 0 G.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(cm ESR=0;%forsimulinkmodel alpha=0.85;%rightdaughterscalingfactor beta=0.62;%leftdaughterscalingfactor lam=10;%length )]TJ0 g 0 G.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 6.864 0 Td [(to )]TJ0 g 0 G.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 6.864 0 Td [(radiusratio N=9;%Maxnumberofgenerations 106

PAGE 123

%%EndInputs r matrix=zerosN,2 * N; r matrix,1=r root; fori=1:N )]TJ/F51 9.9626 Tf 7.158 0 Td [(1%rows forj=1:sizer matrix,2%columns [d1,d2]=get dauj; r d1=alpha * r matrixi,j; r d2=beta * r matrixi,j; ifr d1 > =r term&&r d2 > =r term r matrixi+1,d1=r d1; r matrixi+1,d2=r d2; end ifr matrixi,j+1:end==0 %disp[Fori=,num2stri,Columnsstoppedatj=,num2strj]; break; end end ifr matrixi+1,:==0 disp[Depthstoppedati=,num2stri]; break; end end l matrix=lam. * r matrix;%lengthcalculatedusinganaveragelength )]TJ0 g 0 G.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 6.864 0 Td [(to )]TJ0 g 0 G.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 6.863 0 Td [(radiusratio. %%RCCircuit %FindResistanceusingpoiseuillelawandcomparetotalresistanceto %literaturevaluestoscaleaccordinglyassumeenoughgenerationswhere %P=0fromcapillarybed %C=tau/Rtau=RCandisconstantforentiretree qc=10; 107

PAGE 124

rho=1.06; %%Inputs nu=0.046;%dynamicviscosity ifstrcmpistate,valid R tot=58.4;%TotalResistanceoftree C tot=3.4e )]TJ/F51 9.9626 Tf 7.306 0 Td [(3;%TotalCapacitanceoftree scale=2.6744e )]TJ/F51 9.9626 Tf 7.394 0 Td [(05;%factortoensuretotalresistanceofthetreeequalsR tot %R tot=133.84;%PHvalue %C tot=3.392e )]TJ/F51 9.9626 Tf 7.306 0 Td [(3;%PHvalue %fudge=4.131013234515888e )]TJ/F51 9.9626 Tf 7.394 0 Td [(05;%PHvalue elseifstrcmpistate,control R tot=29381.230626415083 * qc * rho/1333.22365;%Control C tot=.760895451037044e )]TJ/F51 9.9626 Tf 7.495 0 Td [(06/qc/rho * 1333.22365;%Control scale=1.069748781321109e )]TJ/F51 9.9626 Tf 7.395 0 Td [(04;%Control elseifstrcmpistate,ph R tot=67315.21674339622 * qc * rho/1333.22365;%PH C tot=.3414676149796783e )]TJ/F51 9.9626 Tf 7.495 0 Td [(06/qc/rho * 1333.22365;%PH scale=2.450897036656925e )]TJ/F51 9.9626 Tf 7.395 0 Td [(04;%PH else R tot=58.4;%TotalResistanceoftree C tot=3.4e )]TJ/F51 9.9626 Tf 7.306 0 Td [(3;%TotalCapacitanceoftree scale=1.802534166883801e )]TJ/F51 9.9626 Tf 7.395 0 Td [(05;%factortoensuretotalresistanceofthetreeequalsR tot end %%EndInputs tau=R tot * C tot; R=scale * 8 * nu * l matrix./pi * r matrix.; RisnanR=0; C=tau./R; CisinfC=0;%convertInfto0 108

PAGE 125

%Thenneedtosumalloftheresistorsandcapacitorstogetoutflow %BoundaryConditionsfor1Dmodel.AlsotobeabletoconvertRandC %valuesbasedonTotalResistanceandcapacitance. R sum=R; C sum=C; [num rows,num cols]=sizer matrix; fori=num rows )]TJ/F51 9.9626 Tf 7.355 0 Td [(1: )]TJ/F51 9.9626 Tf 7.355 0 Td [(1:1 ifRi,:==0 else forj=1:1:num cols/2 [d1,d2]=get dauj; ifRi+1,d1=0&&Ri+1,d2=0 R sumi,j=R sumi,j+1//R sumi+1,d1+1/R sumi+1,d2; C sumi,j=1//C sumi,j+1/C sumi+1,d1+C sumi+1,d2; end end end end function[d1,d2]=get daucol %d1andd2arethecolumnpositionsfordaughtervessels. d1=2 * col )]TJ/F51 9.9626 Tf 13.136 0 Td [(1; d2=2 * col; end E.20DModelScript clearall; closeall; clc; %%Gettingthesimulationdata state=control; 109

PAGE 126

runTreeBuildup.m simFinal model control.slx; state=ph; runTreeBuildup.m simFinal model ph.slx; simFinal model ph LPAB3.slx; simFinal model ph LPAB5.slx; simFinal model ph repl.slx; %%Storingthedata P c 1:,1=Pressure control3 f 1 g .Values.Time; P c 1:,2=Pressure control3 f 1 g .Values.Data; P c 9:,1=Pressure control3 f 9 g .Values.Time; P c 9:,2=Pressure control3 f 9 g .Values.Data; Q c 1:,1=Flow control3 f 1 g .Values.Time; Q c 1:,2=Flow control3 f 1 g .Values.Data; Q c 9:,1=Flow control3 f 9 g .Values.Time; Q c 9:,2=Flow control3 f 9 g .Values.Data; P ph 1:,1=Pressure ph f 1 g .Values.Time; P ph 1:,2=Pressure ph f 1 g .Values.Data; P ph 9:,1=Pressure ph f 9 g .Values.Time; P ph 9:,2=Pressure ph f 9 g .Values.Data; Q ph 1:,1=Flow ph f 1 g .Values.Time; Q ph 1:,2=Flow ph f 1 g .Values.Data; Q ph 9:,1=Flow ph f 9 g .Values.Time; Q ph 9:,2=Flow ph f 9 g .Values.Data; P ph lpab 1 3:,1=Pressure ph lpab3 f 1 g .Values.Time; 110

PAGE 127

P ph lpab 1 3:,2=Pressure ph lpab3 f 1 g .Values.Data; P ph lpab 9 3:,1=Pressure ph lpab3 f 9 g .Values.Time; P ph lpab 9 3:,2=Pressure ph lpab3 f 9 g .Values.Data; Q ph lpab 1 3:,1=Flow ph lpab3 f 1 g .Values.Time; Q ph lpab 1 3:,2=Flow ph lpab3 f 1 g .Values.Data; Q ph lpab 9 3:,1=Flow ph lpab3 f 9 g .Values.Time; Q ph lpab 9 3:,2=Flow ph lpab3 f 9 g .Values.Data; P ph lpab 1 5:,1=Pressure ph lpab5 f 1 g .Values.Time; P ph lpab 1 5:,2=Pressure ph lpab5 f 1 g .Values.Data; P ph lpab 9 5:,1=Pressure ph lpab5 f 9 g .Values.Time; P ph lpab 9 5:,2=Pressure ph lpab5 f 9 g .Values.Data; Q ph lpab 1 5:,1=Flow ph lpab5 f 1 g .Values.Time; Q ph lpab 1 5:,2=Flow ph lpab5 f 1 g .Values.Data; Q ph lpab 9 5:,1=Flow ph lpab5 f 9 g .Values.Time; Q ph lpab 9 5:,2=Flow ph lpab5 f 9 g .Values.Data; P ph repl 1:,1=Pressure ph repl f 1 g .Values.Time; P ph repl 1:,2=Pressure ph repl f 1 g .Values.Data; P ph repl 9:,1=Pressure ph repl f 9 g .Values.Time; P ph repl 9:,2=Pressure ph repl f 9 g .Values.Data; Q ph repl 1:,1=Flow ph repl f 1 g .Values.Time; Q ph repl 1:,2=Flow ph repl f 1 g .Values.Data; Q ph repl 9:,1=Flow ph repl f 9 g .Values.Time; Q ph repl 9:,2=Flow ph repl f 9 g .Values.Data; %%PlotsofPressureandFlowwaves leg= f Control,PH,PH+LPABL,PH+LPABR,PH+Replacement g ; figure; 111

PAGE 128

%titleInitialVessel holdon; plotP c 1:,1,P c 1:,2; plotP ph 1:,1,P ph 1:,2; plotP ph lpab 1 3:,1,P ph lpab 1 3:,2; plotP ph lpab 1 5:,1,P ph lpab 1 5:,2; plotP ph repl 1:,1,P ph repl 1:,2; legendleg,location,southeast; xlabelTimesec ylabelPressuremmHg figure; %titleInitialVessel holdon; plotQ c 1:,1,Q c 1:,2; plotQ ph 1:,1,Q ph 1:,2; plotQ ph lpab 1 3:,1,Q ph lpab 1 3:,2; plotQ ph lpab 1 5:,1,Q ph lpab 1 5:,2; plotQ ph repl 1:,1,Q ph repl 1:,2; legendleg,location,northwest; xlabelTimesec ylabelFlowRatemL/sec figure; %titleTerminalVesselBranch holdon; plotP c 9:,1,P c 9:,2; plotP ph 9:,1,P ph 9:,2; plotP ph lpab 9 3:,1,P ph lpab 9 3:,2; plotP ph lpab 9 5:,1,P ph lpab 9 5:,2; plotP ph repl 9:,1,P ph repl 9:,2; legendleg,location,southeast; 112

PAGE 129

xlabelTimesec ylabelPressuremmHg figure; %titleTerminalVesselBranch holdon; plotQ c 9:,1,Q c 9:,2; plotQ ph 9:,1,Q ph 9:,2; plotQ ph lpab 9 3:,1,Q ph lpab 9 3:,2; plotQ ph lpab 9 5:,1,Q ph lpab 9 5:,2; plotQ ph repl 9:,1,Q ph repl 9:,2; legendleg,location,northwest; xlabelTimesec ylabelFlowmL/sec E.3WSSandPPScript %%WallShearStress mu=0.035;%g/cm * s r=r matrix,1;%cm A=pi * r;%cm q ctl=maxFlow control3 f 9 g .Values.Data:end/2;%cm/sec WSS ctl=0.1 * 4 * mu * q ctl/A * r;%Pascals q ph=maxFlow ph f 9 g .Values.Data:end/2;%cm/sec WSS ph=0.1 * 4 * mu * q ph/A * r;%Pascals q ph lpab3=maxFlow ph lpab3 f 9 g .Values.Data:end/2;%cm/sec WSS ph lpab3=0.1 * 4 * mu * q ph lpab3/A * r;%Pascals q ph lpab5=maxFlow ph lpab5 f 9 g .Values.Data:end/2;%cm/sec WSS ph lpab5=0.1 * 4 * mu * q ph lpab5/A * r;%Pascals 113

PAGE 130

q ph repl=maxFlow ph repl f 9 g .Values.Data:end/2;%cm/sec WSS ph repl=0.1 * 4 * mu * q ph repl/A * r;%Pascals WSS ctl n=WSS ctl/WSS ctl; WSS ph n=WSS ph/WSS ctl; WSS ph lpab3 n=WSS ph lpab3/WSS ctl; WSS ph lpab5 n=WSS ph lpab5/WSS ctl; WSS ph repl n=WSS ph repl/WSS ctl; leg c=categorical f Control,PH,PH+LPABL,PH+LPABR,... PH+Replacement g ; figure; barleg c,[WSS ctl n,WSS ph n,WSS ph lpab3 n,WSS ph lpab5 n,WSS ph repl n]; xlabelTestGroup ylabelWSSnorm %%PulsePressure PP ctl=maxPressure control3 f 9 g .Values.Data:end )]TJ.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(... minPressure control3 f 9 g .Values.Data:end; PP ph=maxPressure ph f 9 g .Values.Data:end )]TJ.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(... minPressure ph f 9 g .Values.Data:end; PP ph lpab3=maxPressure ph lpab3 f 9 g .Values.Data:end )]TJ.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(... minPressure ph lpab3 f 9 g .Values.Data:end; PP ph lpab5=maxPressure ph lpab5 f 9 g .Values.Data:end )]TJ.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(... minPressure ph lpab5 f 9 g .Values.Data:end; PP ph repl=maxPressure ph repl f 9 g .Values.Data:end )]TJ.133 .545 .133 rg .133 .545 .133 RG/F51 9.9626 Tf 13.136 0 Td [(... minPressure ph repl f 9 g .Values.Data:end; figure; barleg c,[PP ctl,PP ph,PP ph lpab3,PP ph lpab5,PP ph repl]; 114

PAGE 131

xlabelTestGroup ylabelPulsePressuremmHg 115