
Citation 
 Permanent Link:
 http://digital.auraria.edu/AA00007245/00001
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:
 2019
 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 3element 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 lengthtoradius ratio, lam,
values for each vessel for a tree with asymmetric vessel geometry....... 24
2.4 Upstream radii, Ru, downstream radii, Rd and lengthtoradius 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 tradeoff between computational cost and physiological relevance (i.e. how representative
the model is).......................................................... 4
2.1 HighLevel 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 36 are the terminal vessels. The outflow boundary condition is a 3element 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 2generation RC circuit tree model made in Simulink
showing the inlet flow and outflow boundary condition (BC)............. 19
2.6 An example of a 2generation 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 = le8) 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 = le8) 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 musclelike 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 proinflammatory 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 SwanGanz 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 noninvasively 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 noninvasive ways to ask physiological questions about the body. Generally, they give us the ability to investigate hardtomeasure 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], onedimensional (ID) [12, 13, 14, 15, 16, 17, 18, 19, 20, 21] and threedimensional (3D) [18, 19, 20, 21]. OD models are pointbased 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 tradeoff 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 systemlevel, 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 fluidstructure 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 NavierStokes 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, PHmodeled rats were compared to PHmodeled 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 proinflammatory 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 animalspecific 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 1D0D model. First the 1D0D 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 largescale network.
Specific Aim 2: Execute pulmonary mousespecific 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: HighLevel 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 quasilinear hyperbolic PDEs consists of three equations which are used to solve for flux inside of the vessel lumen, crosssectional area, and pressure. The three equations are solved using a pythonbased, 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/ll. 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 objectoriented 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 NavierStokes equations for conservation of mass and momentum in a onedimensional coordinate system. The details of this solution can be found elsewhere [31] with the resulting system of nondimensionalized 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 36 are the terminal vessels. The outflow boundary condition is a 3element 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 crosssectional 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 crosssectional 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 twostep Lax Wendroff explicit scheme, detailed in [30, 31] which is secondorder accurate in time and space, i.e. 0(At2, Ax2) and must satisfy the CFL condition: At < Ax\^ Â±c1 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: crosssectional 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 zerodimensional, 3element 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 < le6. 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 getdaughters 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 nonlinear 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 1112 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 lengthtoradius 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 2generation 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 2generation 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 2generation 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 postprocessing. 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.2949e4 1.0697e4
AUP Htestgroups 535.2 2.945e5 2.4509e4
Simulink model can be run manually with the variables still being saved in the MATLAB 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 lengthtoradius ratio values for each vessel, taken from data in [23]; however, the radius and lengthtoradius 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 lengthtoradius 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 lengthtoradius 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) le5 le7
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.4e3
PH 1.1683e5 133.8 1.178e3
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 le7 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 2generation 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 2generation 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 le7): 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 1020% [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)
le6 1.366
le7 13.66
le8 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 le6 and for the sake of time, only three dt values were studied: le6, le7 and le8. 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 le8, 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
1e8
Figure 3.8: Maximum Relative Error (compared to dt = le8) 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 = le8) 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 le7 was chosen in order to keep a small dt with a reasonable runtime as the runtime for the le8 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 nonPH), 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 lengthtoradius 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 = le7 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 endcycle pressure was increased compared to the PH despite having a lower peak pressure. The lowered endcycle 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 animalspecific 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 largescale 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 pythonbased 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 circuitbased 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 nonbanded 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:D3441, 2013.
[3] W. Sun and S.Y. Chan. Pulmonary arterial stiffness: An early and pervasive driver of pulmonary arterial hypertension. Frontiers in medicine, 5:204204, 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): 14131435, 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):560580, 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 multiscale, subjectspecific 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:2222, 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 408417, 2008.
61
[12] M.S. Olufsen. Structure tree outflow condition for blood flow in larger systemic arteries. American Journal of Physiology, 276(l):H25768, 1999.
[13] M.S. Olufsen. A onedimensional fluid dynamic model of the systemic arteries. Studies in Health Technologies and Informatics, 71:7998, 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):30473068, 2016.
[16] J.P. Mynard and J.J. Smolich. Onedimensional haemodynamic modeling and wave dynamics in the entire adult circulation. Annals of Biomedical Engineering, 43(6):14431460, 2015.
[17] P. Reymond, F. Merenda, F. Perren, D. Rufenacht, and N. Stergiopulos. Validation of a onedimensional model of the systemic arterial tree. American Journal of PhysiologyHeart and Circulatory Physiology, 297(1):H208H222, 2009.
[18] Jonasova A., Bublik O., and Vimmr J. A comparative study of Id and 3d hemodynamics in patientspecific 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 3d fsi, 3d rigid wall and 1d models. Medical Engineering & Physics, 35(6):784791, 2013.
[21] N. Xiao, J. Alastruey, and C. Alberto Figueroa. A systematic comparison between 1d and 3d 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:193252, 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 modelbased 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 zerod and 1d 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 onedimensional 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 imagebased, computational fluid dynamics study. Pulmonary circulation, 2(4):470476, 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 onedimensional carotid artery bifurcation, part i: Application to stent design. IEEE Transactions on Biomedical Engineering, 54(5):802812, 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 structuredtree 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 375378, 1963.
[35] M.I. Townsley. Structure and composition of pulmonary arteries, capillaries, and veins. Comprehensiue Physiology, 2(l):675709, 2012.
63
[36] N. Rol, E.M. Timmer, T.J.C. Faes, A. Vonk Noordegraaf, K. Grunberg, H. Bogaard, 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, whichappliesaphysiologicalowwaveformattheinletanda2or3elementoutow 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 andlengthtoradiusratio, lam , valuesforeachvesselforatreewithasymmetricvesselgeometry.....24 2.4Upstreamradii, R u ,downstreamradii, R d andlengthtoradiusratio, 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,1Dand0Dmodelsandthetradeobetween computationalcostandphysiologicalrelevancei.e.howrepresentative themodelis..................................4 2.1HighLevelOverviewoftheFullComputationalModel..........8 2.2Anexampleofatreewith3generationswherevessel0representsthe mainpulmonaryarteryMPA,vessel1representstheleftpulmonary arteryLPA,vessel2representstherightpulmonaryarteryRPAand vessels36aretheterminalvessels.Theoutowboundaryconditionisa 3elementwindkesselappliedateachterminalvessel............10 2.3InletowwaveformforcontrolandpulmonaryhypertensionPHtest groups.....................................12 2.4Anexampleofatreewithasymmetricbranchingwherethereareterminal branchesatgeneration2andatgeneration3................17 2.5Anexampleofa2generationRCcircuittreemodelmadeinSimulink showingtheinletowandoutowboundaryconditionBC........19 2.6Anexampleofa2generationRCcircuittreemodelmadeinSimulink 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=1e8forthelargestvessel, Vessel0,forvaryingdtvaluesaswellasthetotaltimeittookforeach simulationtorun...............................39 3.9MaximumRelativeErrorcomparedtodt=1e8forthesmallestvessel, 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 dierentiateintomyobroblastswhicharesmoothmusclelikecells,andmigrateto themediallayer.Fibroblastsintheadventitiaalsochangetheextracellularmatrix ECMcompositionofvascularwallsinPH;morespecically,thereisanincrease 1
PAGE 18
ofcollagentypesIandIII.Anincreaseincollageninthevascularwallleadstoan increaseinstinessofthewallaswell[3,4,5]. EndothelialcellECdysfunctionisknowntooccur,whichcanfurthersupport SMCproliferationandcancauseplexuslesionformationaswellasendothelialchannel formation.ECdysfunctionisalsoassociatedwithlossofpulmonaryarterialvessels [1,6].Aspartofthevascularremodelingprocess,leukocytesarebeingrecruitedby adventitialbroblastsandproinammatorycytokinelevelsarehigherinPHvessels [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 SwanGanzcathetertomeasurepressureandmeanowcardiacoutput[8].A catheterisnotabletotravelfarintothevesselnetworkduetoitssize;thus,thereis 2
PAGE 19
notawaytodirectlymeasureowinthedistalvasculature.Othertechniquesexist tomeasurebloodownoninvasivelysuchasmagneticresonanceimagingMRIand ultrasonictechniques;however,theyalsocomewithlimitationsthatdonowallow forthedirectmeasurementofbloodowindownstreamvessels[9].Duetothese limitations,computationalmodelinghasbecomeanattractivemethodforestimating owindistalvasculature. Computationalmodelsarenoninvasivewaystoaskphysiologicalquestionsabout thebody.Generally,theygiveustheabilitytoinvestigatehardtomeasureproperties, tomakepredictionstocreatemoreeectiveanimalmodels,andtobetterinvestigate treatmentoptions.Therearemanytypesofcomputationalmodels,suchaszerodimensionalD[10,11],onedimensionalD[12,13,14,15,16,17,18,19,20,21] andthreedimensionalD[18,19,20,21].0Dmodelsarepointbasedandonly includethetimedomain,1Dmodelsaddonespatialdomainwhichisthelength alongavessel,and3Dmodelsincludeallthreespatialdimensionsofthevesselas wellasthetimedomain.Modelscanalsobeacombinationofthesetypes. 3
PAGE 20
Figure1.1:Graphicrepresentationof3D,1Dand0Dmodelsandthetradeobetween computationalcostandphysiologicalrelevancei.e.howrepresentativethemodelis. Eachmodelgivesadierentsetofinformationaboutow;3Dmodelsallowfor representationofmorecomplexowbutatahighcomputationalcostwhilea0D modelhasaverylowcomputationalcostbutonlyincludesmoresimpleowfeatures. Thisisrepresentedingure1.1.Typically,0Dmodelsareusedtolookatthesystemlevel,1Dmodelscangiveaccuratebloodpressureandowpropagationinformation foralargenetworkofvessels,and3Dmodelsallowforadeeperinvestigationof uidstructureinteractionsFSI,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,1Dmodelassumptionsallowforagoodrepresentationofrelativelylargerbloodvesselsandusesimplied3DNavierStokesequations.Common numericalschemesexisttosolvethesystemofequationssetupfor1Dmodels,suchas Galerkinvariations,FiniteDierencesmethods,andcharacteristicmethods[24,26]. 1.3Justication InPH,thereisdistalvesselthickeningoccurring.Datafromotherliteraturesuggeststhatproximalstieningcauseshigherowpulsatilitywhichcausesendothelial celldysfunctionandactivatesbroblasts.Activatedbroblastsproliferateanddepositcollagenwhichcausesstieningandthickeningofthedistalvessels.Onestudy providesgreatevidenceforthisproposedmechanismbyshowingthathemodynamic unloadingiseectiveattreatingPHphenotypes[27].Inthisstudy,PHmodeledrats werecomparedtoPHmodeledratswithabandedleftpulmonaryarteryLPAand studiedthehistologyoftherats'arteries.TheyfoundthattheLPAofthebanded ratshadreducedexpressionofproinammatorycytokines,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 studiescouldbedonesimultaneouslywithanimalspecicsimulationstoseeifthere 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 distalowwereinvestigatedusingacoupled1D0Dmodel.Firstthe1D0Dmodel wascreatedinPythonandMATLABandvalidatedagainstpreviouslypublisheddata [23].Then,simulationswereruntosimulatefourdierenttestgroupsbasedon[27] tostudychangesinWSSandpulsepressure.Theaimsaredescribedbelow. SpecicAim1:Developacomputationalmodelofthepulmonaryarterialsystem toallowsolvingofbloodowhemodynamicsinalargescalenetwork. SpecicAim2:Executepulmonarymousespecicsimulationsmatchingproximalpulmonaryarterialstinessandpulmonaryvascularresistancetocomparepulse 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:HighLevelOverviewoftheFullComputationalModel 2.1.1Assumptions 8
PAGE 25
Themodelsassumethevesselsareelasticandaxisymmetrictubes,containing Newtonianuid.Otherassumptionsforthe1Dmodelarethatthebloodpressureis uniformacross ~r ,theradius,andthereispulsatilelaminarow. 2.21DModelVaMpy The1Dcomputationalmodelutilizesuiddynamics,wallmechanics,andboundaryconditions.TheresultingsystemofquasilinearhyperbolicPDEsconsistsofthree equationswhichareusedtosolveforuxinsideofthevessellumen,crosssectional area,andpressure.Thethreeequationsaresolvedusingapythonbased,opensource modelavailableongithub,namedtheVascularModelinginPythontoolkitorVaMpy [30].Thismodelstartswithasingleparentvesselthatbifurcatesintodaughtervessels.Theresultingtwodaughtervesselsalsobifurcate.Thisprocesscontinuesuntil theterminalvesselsarereached.Thetotalnumberofvesselsmadeisdenotedby 2 depth 1.Therefore,ifthearterynetworkhadadepthof3,aka3generations,there wouldbe7vessels.Notethatvesselnumberingbeginswith0.Asanobjectoriented program,eachvesseliscontainedinitsownartery"classwhichareallcontained withinintheArteryNetwork"class.TheVaMpymodelusesacongurationlethat containsthevesselgeometryandvascularpropertyinformation.Anexamplevascular treewith3generationscontainingoutowboundaryconditionsisshowningure2.2 referencingvesselsinthepulmonaryvasculaturestartingwiththemainpulmonary arteryMPA. 2.2.1Equations The1Dmodelsolvesthreeequations.Thersttwoequationsarederivedthrough integrationoftheNavierStokesequationsforconservationofmassandmomentum inaonedimensionalcoordinatesystem.Thedetailsofthissolutioncanbefound elsewhere[31]withtheresultingsystemofnondimensionalizedequationsrepresented bytheequationshownbelowinEq2.1 ~ U t + ~ F x = ~ S .1 9
PAGE 26
Figure2.2:Anexampleofatreewith3generationswherevessel0representsthe mainpulmonaryarteryMPA,vessel1representstheleftpulmonaryarteryLPA, vessel2representstherightpulmonaryarteryRPAandvessels36aretheterminal vessels.Theoutowboundaryconditionisa3elementwindkesselappliedateach 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 ,thecrosssectionalarea,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 withcorrespondingcrosssectionalarea 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'stwostepLaxWendroexplicit scheme,detailedin[30,31]whichissecondorderaccurateintimeandspace,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:crosssectionalarea, 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,azerodimensional,3elementWindkesseloutowboundaryconditionisusedtoclosethecomputationaldomain,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.Theculpritwasthebifurcationboundaryconditionduetoits18nonlinear equationsbeingsolvedbyNewton'smethodwhichincreasesthecomputationalcost. Duetotheslownessofthemodel,the1Dmodelwascombinedwitha0Dmodelin ordertomodelanentirevascularnetworkrelativelyquickly.Oncetheparameters werechangedforvalidation,meshindependence,andactualmodelsimulations,the runtimeforthe1Dmodeldecreasedto1435hours. 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],smalleranimalshaveatotalof1112generations 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 lengthtoradiusratioandstoredintoamatrix, ~ 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.Anexampleofa2generationtreeis showningure2.5.Forthemodeltooutputpressureandowwavesforavessel, Figure2.5:Anexampleofa2generationRCcircuittreemodelmadeinSimulink showingtheinletowandoutowboundaryconditionBC. VoltmetersandAmpMetersareneededinthemodel.Anexampleofthesame2generationtreefromgure2.5isshowningure2.6withthesemetersadded.The VoltmetersandAmpMetersalsoneedtobeconnectedtoascopeblock,whichisjust ablockthatprovidesgraphicvisualization.Inthemodelyoucanseeascopeblock 19
PAGE 36
ontheleftforthevoltageandascopeblockontherightforthecurrent. Figure2.6:Anexampleofa2generationRCcircuittreemodelmadeinSimulink 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 subsequentlyusedinthesameMATLABscriptforpostprocessing.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.2949e41.0697e4 AllPHtestgroups 535.22.945e52.4509e4 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.3showstheradiusandlengthtoradiusratiovaluesforeachvessel,takenfromdatain[23];however,theradiusandlengthtoradius 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 andlengthtoradiusratio, 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 andlengthtoradiusratio, 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 sec1e51e7 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.4e3 PH 1.1683e5133.81.178e3 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. AnNvalueof50andadtvalueof1e7wereusedforthevalidationsimulations.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,asingleRCcircuitwascomparedtothatofa2generationRCcircuit wherebothmodelstotaledtothesameresistanceandcapacitancetoensurethe modelcouldbeexpandedintobifurcationswiththeresultsshowningure3.5.These 32
PAGE 49
resultsingure3.5Ashowgoodagreementbetweenpressureoftheode45solvinga singleRCcircuitcomparedtothetotalpressureofthe2generationRCmodelwhere 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.ThefollowingNvaluesweretestedwithadtof1e7: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 precisionerrortobewithin1020%[37].Therefore,thesimulationsobtainanerror lessthanthatoftheinstrumentthatprovidestheinletboundarycondition.Alarger Nvalueakaasmallerdxismoredesirable,especiallyduetothefactthatalarger 37
PAGE 54
Table3.3:Thetotaltimeeachsimulationranforagivendtvalue. dtsecTimehours 1e61.366 1e713.66 1e8157.2 Nvalueledtoashorterruntime,and,thus,anNvalueof100wasdecidedonfor subsequentsimulations. Next,meshindependencestudywasperformedonvaryingdtvalueswithN= 100aschosenfromthedxmeshindependencestudy.DuetotheCFLcondition,the smallestdtvaluethatcouldbeusedwas1e6andforthesakeoftime,onlythreedt valueswerestudied:1e6,1e7and1e8.Thetimeittookforeachsimulationtorun isshownintable3.3.Theresultsforthechangesinpeakpressurealongthelength ofthevesselforeachdtvalueisshowninAppendixB.Themaximumrelativeerror foreachdtvalue,calculatedagainstthelargestdtvalueof1e8,isshowningure 3.8forVessel0andgure3.9forVessel4.Theseguresshowthatthesolutionis convergingwithincreasingdtvalues.Asshowninbothguresandintable3.3,a decreasingdtvalueleadstoanincreasingruntime,asexpected. 38
PAGE 55
Figure3.8:MaximumRelativeErrorcomparedtodt=1e8forthelargestvessel, Vessel0,forvaryingdtvaluesaswellasthetotaltimeittookforeachsimulationto run. Figure3.9:MaximumRelativeErrorcomparedtodt=1e8forthesmallestvessel, Vessel4,forvaryingdtvaluesaswellasthetotaltimeittookforeachsimulationto run. Again,weseethatthepeakpressureiswithin0.5mmHgseeAppendixBfor 39
PAGE 56
plotsorlessandthattherelativeerrorislessthan0.025%;therefore,anydtvalue couldbeused.Again,asmallerdtvalueisdesirable,soadtvalueof1e7waschosen inordertokeepasmalldtwithareasonableruntimeastheruntimeforthe1e8 simulationtook157.2hoursorroughly6.5dayswhichisanunreasonableamountof time. 3.2Simulations AsstatedintheMethodssection,fourtestgroupswerestudied:controlrepresentingnonPH,PH,PHwithLPAbandingPH+LPAB,andPHwithasynthetic/compliantLPAandRPAPH+replacement. 3.2.11DSimulationTime DuetocomputationalcostofoneofthetestgroupsPH+LPAB,the1D simulationswererunusingatreewithasymmetricgeometrywiththeradiiand lengthtoradiusvaluesshownintable2.4.AnNvalueof46wasusedasthiswas thelargestNvaluethatwouldsatisfytheCFLconditionforadt=1e7tokeepthe 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 isthattheendcyclepressurewasincreasedcomparedtothePHdespitehavinga lowerpeakpressure.TheloweredendcyclepressurealsocontributetothePH+ 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 usinganimalspecicdataforeachsimulation.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 thisprojectweretodevelopacomputationalmodelofalargescalepulmonaryvasculaturenetworkaswellastoinvestigatepulsatility,usingpulsepressurePP,and WSSindistalvasculaturegivenchangestotheproximalpulmonaryvasculature.More specically,simulationofabandedleftpulmonaryarteryLPABinthepresenceof PHandsimulationofareplacementleftandrightpulmonaryarteryinthepresence ofPHwerecomparedtohealthyandPHcontrols. Thecomputationalmodelwasdevelopedusinga1Dand0Dhybridmodelusing anexistingpythonbased1Dmodelasthebase.Themodelwasexpandeduponto allowforinvestigationofthedesiredtestgroupsdescribed,buthadahighcomputationalcost.Thus,a0Dmodelwasaddedtocreatethehybridmodelwherethe 0DmodelusedSimulinktocreateacircuitbasedanalogofthedistalvasculature. ThemodelwasthenusedtoinvestigatePPandWSSinthedistalvessels.PPwas increasedforPHcomparedtothehealthycontrol,asexpected.ThePH+LPAB LshowedaslightdecreasecomparedtoPHwhilethePH+LPABRshowedan increasedPPbeyondthatofthePHcase.TheleftandrightpulmonaryarteryreplacementgroupreturnedPPclosertocontrollevels,whichshowedgoodpromisefor asyntheticleftandrightpulmonaryarteryasatreatmentoptionforPH.TheWSS forPHwasalsoincreasedcomparedtothehealthycontrol,butonlywhenpruning anddecreasedvesselareaweremodeledprovidingevidencethatchangestoWSSare 59
PAGE 76
completelyduetogeometryandnotvesselmechanics.InthePH+LPABcase,the downstreamarteriesofthebandedarteryshowedadecreaseinWSScomparedto PHwhilethearteriesdownstreamfromthenonbandedarteryshowedanincreasein 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 ofmodellingassumptionsinmultiscale,subjectspecicmodelsofaortichaemodynamics. 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.Aonedimensionaluiddynamicmodelofthesystemicarteries. 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.Onedimensionalhaemodynamicmodelingand wavedynamicsintheentireadultcirculation. AnnalsofBiomedicalEngineering , 43:1443{1460,2015. [17]P.Reymond,F.Merenda,F.Perren,D.Rufenacht,andN.Stergiopulos.Validationofaonedimensionalmodelofthesystemicarterialtree. AmericanJournal ofPhysiologyHeartandCirculatoryPhysiology ,297:H208{H222,2009. [18]JonasovaA.,BublikO.,andVimmrJ.Acomparativestudyof1dand3d hemodynamicsinpatientspecichepaticportalveinnetworks. 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 indicesaspredictedby3dfsi,3drigidwalland1dmodels. MedicalEngineering&Physics ,35:784{791,2013. [21]N.Xiao,J.Alastruey,andC.AlbertoFigueroa.Asystematiccomparisonbetween1dand3dhemodynamicsincompliantarterialmodels. 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:amodelbasedanalysisofthediseasemechanism. 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.Reviewofzerodand1dmodelsofbloodow 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 ofnumericalschemesforonedimensionalarterialbloodowmodelling. 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:Animagebased,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 aonedimensionalcarotidarterybifurcation.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 structuredtreeoutowconditions. 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

