A multiscale model of vascular function in chronic thromboembolic pulmonary hypertension
Mitchel J. Colebank, M. Umar Qureshi, Sudarshan Rajagopal, Richard A. Krasuski, Mette S. Olufsen
AA multiscale model of vascular function in chronic thromboembolicpulmonary hypertension
Mitchel J. Colebank a , M. Umar Qureshi a , Sudarshan Rajagopal b , Richard A. Krasuski c , Mette S.Olufsen a, ∗ a Department of Mathematics, North Carolina State University, Raleigh, NC, NC 27695, USA b Division of Cardiology, Department of Medicine, Duke University Medical Center, Durham, NC 27710, USA c Department of Cardiovascular Medicine, Duke University Health System, Durham, NC 27710, USA
Abstract
Pulmonary hypertension (PH) is a debilitating disease that includes five main subgroups, but only one iscurable: chronic thromboembolic pulmonary hypertension (CTEPH). CTEPH is caused by recurrent orunresolved pulmonary thromboemboli, leading to perfusion defects and increased arterial wave reflection.Treatment for CTEPH depends on lesion location and severity, and patients with distal lesions are inoperableby standard surgical intervention. Instead, these patients are treated with balloon pulmonary angioplasty(BPA), a multi-session, minimally invasive surgery that disrupts the thromboembolic material blockingthe vessel lumen using a catheter balloon. The goal of BPA is to reduce pulmonary arterial pressure andreestablish adequate blood flow to all regions of the lung. However, there still lacks an integrative, holistictool for identifying target lesions for an optimal treatment. To address this insufficiency, we provide a proof ofconcept study simulating CTEPH hemodynamics and BPA therapy using a multiscale, one-dimensional (1D)fluid dynamics model. The large pulmonary arterial geometry is derived from a computed tomography (CT)image, whereas a fractal, structured tree represents the small vessels. We model ring- and web-like lesions,common in CTEPH, using two pressure-loss terms. We simulate normotensive conditions and four CTEPHdisease scenarios; the latter includes 20 large artery lesions and introduces both local and global vascularremodeling. BPA therapy is simulated by simultaneously reducing lesion severity in three locations. Ourpredictions mimic severe CTEPH, manifested by an increase in mean proximal pulmonary arterial pressureabove 20 mmHg and prominent wave reflections. Both flow and pressure decrease in vessels distal to thelesions and increase in unobstructed vascular regions. We use the main pulmonary artery (MPA) pressure,a wave reflection index, and a measure of flow heterogeneity to select optimal target lesions for BPA. Insummary, this study provides a preliminary multiscale, image-to-hemodynamics pipeline for BPA therapyplanning for inoperable CTEPH patients.
Keywords: elsarticle.cls , L A TEX, Elsevier, template
1. Introduction
Pulmonary hypertension (PH), defined as a mean pulmonary arterial pressure (mPAP) >
20 mmHg atrest [1], degrades pulmonary vascular function. The disease is associated with vascular wall thickening anddecreased vessel compliance, elevating pulmonary arterial pressure. If untreated, PH progressively leadsto vascular-ventricular decoupling and right ventricular dysfunction [2]. PH encompasses five etiologies: ∗ Corresponding author
Email addresses: [email protected] (Mitchel J. Colebank), [email protected] (M. Umar Qureshi), [email protected] (Sudarshan Rajagopal), [email protected] (Richard A. Krasuski), [email protected] (Mette S. Olufsen )
Preprint submitted to Computer Methods in Applied Mechanics and Engineering February 15, 2021 a r X i v : . [ q - b i o . T O ] F e b
1) pulmonary arterial hypertension (PAH); (2) PH due to left heart disease; (3) PH due to lung disease,including hypoxia and chronic obstructive pulmonary disease; (4) chronic thromboembolic pulmonary hyper-tension (CTEPH); and (5), PH due to unclear or multifactorial mechanisms [1]. Though vasodilator therapycan be used to treat several PH subclasses, only CTEPH (group 4) is curable via surgical intervention [3, 4].Recently, a less invasive approach has emerged, known as balloon pulmonary angioplasty (BPA), that opensthe blood vessels using cardiac catheterization [5]. However, it is unclear if this more novel treatment hasthe same hemodynamic improvements as the gold-standard surgical methods. To this extent, this studydevelops a multi-scale, computational fluid dynamics (CFD) model to predict hemodynamics in CTEPHand simulate the impact of balloon pulmonary angioplasty (BPA) treatment.CTEPH affects 3–30 people per million worldwide [6]. However, its prevalence is likely underestimateddue to the diagnostic complexity and generic nature of its two vital symptoms: dyspnea and fatigue [7].CTEPH etiology involves recurrent obstruction of the pulmonary arteries, leading to perfusion deficitsand pulmonary vascular remodeling. In addition to elevated mPAP, clinical manifestations of CTEPHinclude increased flow heterogeneity and elevated right ventricular afterload [8, 9]. Diagnosis requires acombination of medical imaging (e.g., chest computed tomography (CT) or magnetic resonance imaging(MRI)), ventilation-perfusion (V/Q) scanning, and invasive right heart catheterization. CT images provideinformation regarding pulmonary vascular structure, while V/Q scans display blood flow and oxygen deliveryto the alveoli. Right heart catheterization measures dynamic pulmonary arterial blood pressure in the main,left, and right pulmonary arteries (MPA, LPA, and RPA, respectively). Catheterization also estimates leftatrial pressure by calculating pulmonary arterial wedge pressure (PAWP), which is ≤
15 mmHg for CTEPHand other pre-capillary PH [1]. Cardiac output is measured during catheterization via thermodilution orFick’s method. CTEPH diagnosis requires a pulmonary vascular resistance (PVR) ≥ > Z ) can be adjusted and matched topressure data during respiration. While these studies account for the distal vascular structure to determineimpedance boundary conditions, they do not predict hemodynamics within the distal vasculature itself.Several studies predict multiscale hemodynamics, imperative for understanding the interactions betweenthe large and small vessels during PH [26, 33, 34, 35]. Yang et al. [26] use the impedance boundary conditionpreviously implemented by Spilker et al. [29] in a 3D pulmonary circulation model. They predict meanWSS in a Strahler-based distal vascular structure, and illustate that severe PAH increases mean WSS inthe smallest arteriolar branches. Clark and Tawhai [33] use a wave-transmission model of the pulmonaryarterial and venous trees and couple it to a capillary sheet model to predict hemodynamics throughout thepulmonary circulation. They use imaging to obtain arteries and veins to the subsegmental level, whereasthe arterioles and venules are generated using a space-filling algorithm. Their multiscale study also providesqualitative predictions of lung perfusion, showing increased flow heterogeneity with reduced pulmonaryvascular compliance. The study by Olufsen et al. [35] solves nonlinear, inertia-driven 1D flow equations in3he large pulmonary arteries and a linearized wave-equation in the arterioles of the structured tree model.The authors use the model framework to investigate pulmonary microvascular rarefaction and show thatreducing the number of arterioles increases mPAP. A similar modeling approach is used in Chambers etal.[34], which compares control and hypertensive (PH induced by hypoxia - PH group 3) mice. Their studyuses large artery geometry from micro-CT imaging, whereas the small arterial and arteriolar networks areconstructed from fractal trees with scaling parameters extracted from high-resolution imaging data.To our knowledge, no previous studies incorporate pressure loss models to study the effect of lesions inCTEPH, nor do they attempt to predict large and small vessel fluid dynamics in CTEPH. To address thiscritical gap in the literature, we propose a novel multiscale approach that predicts large and small vesselhemodynamics. We conduct the former simulations in a patient-specific geometry extracted from CT images,whereas the latter simulations are done in asymmetric, structured trees. We combine this computationalmodel with two energy loss models for ring- and web-like thromboembolic lesions to predict pressure, flow,and area in normotensive and CTEPH scenarios. To understand how lesions and pulmonary vascularremodeling affect pulmonary hemodynamics, we compute mPAP, FFR, and a wave reflection index, as wellas lung perfusion. We compute these quantities under normotensive and four physiologically-driven CTEPHconditions. The flow distribution within the lung is converted to a quantitative flow heterogeneity indexutilizing the Kullback-Leibler (KL) divergence [36]. Finally, we simulate “virtual BPA” by modulating lesionseverity in multiple vessels and identify the most significant set of lesions for percutaneous intervention. Thisstudy develops a proof-of-concept modeling tool, which can assist the interdisciplinary team of cardiologists,pulmonologists, and radiologists in determining the optimal intervention strategy for BPA therapy.
2. Materials and Methods
In this study, the pulmonary arterial geometry consists of an image-based large artery network generatedfrom a CT image and a self-similar, fractal-tree representing the small arteries and arterioles.
We generate a network of large pulmonary arteries from a publicly available chest CT scan (from Sim-vascular ). The CT image is from a healthy, 67-year-old female volunteer, captured using Omnipaque 350contrast agent. Similar to our previous study [34], we segment and reconstruct the pulmonary arterialnetwork and lung cavity using the open-source segmentation software 3DSlicer [37]. When segmenting thegeometry, we first identify the MPA, LPA, and RPA, and then manually segment all lobar, segmental, andsubsegmental vessels. We generate the 1D network by reducing the 3D vessel geometry to centerlines usingthe Vascular Modeling ToolKit (VMTK)[38] and use our network extraction algorithm to generate a directedgraph from the centerlines [39]. This algorithm converts the network to a series of arcs (vessels) and nodes(junctions), encoded by a connectivity matrix. For non-terminal arteries, each vessel’s outlet is connected toeither two or three vessel inlets, deemed the parent and daughter vessels, respectively. Terminal vessels areprescribed a structured tree model at their outlet, as defined below. Each vessel has a constant radius and alength. Radii are computed as the mean within the center 80% of the vessel, and the lengths are computedas the Euclidean distance between each point in the vessel. Figure 1 illustrates this process showing threeplanes in the CT image, the segmented network and lung cavity, and the 1D network. It is difficult to obtain high-fidelity measurements of the small arteries and arterioles from CT images.To circumvent this, we represent the small vessels with fractal, self-similar structured trees [32, 35]. Weassume that each daughter vessel is related to its parent via the radius scaling factors α and β , and a length http://simvascular.github.io/ essel and lung segmentation Skeletonization and model constructionImage analysis LesionsModel of arterioles Figure 1: Schematic of model development, starting from image analysis and image segmentation and ending in constructionof the skeletonized network and computational domain. to radius ratio L rr . The dimensions of any vessel in the structured tree can be expressed in terms of theterminal large artery radius r term proximal to the beginning of the structured tree r n,m = r term α n β m , L n,m = L rr r n,m . (1)Following [40], who use a combination of pulmonary artery literature data, the length to radius ratio usedin this study is L rr = (cid:40) . r . r ≥ . . r . r ≤ . . (2)The structured tree is terminated after a specified minimal radius, r min , has been reached. We use param-eter values from a previous investigation from high-resolution micro-CT scans of control and hypertensivemice [34], and set α = 0 .
88 and β = 0 .
68. The terminal radius of the structured tree is based on anatomicalliterature and set to 7 . µ m [41]. The 1D model consists of the large pulmonary arteries segmented from the CT image, and the smallpulmonary arteries and arterioles represented by structured trees. In the large arteries, we predict inertia-driven hemodynamics using a system of nonlinear partial differential equations (PDEs), while fluid dynamicsin the network of small vessels are predicted by solving a linearized wave equation.
Similar to previous studies [40, 42, 43], we compute flow q ( x, t ) (mL/s), pressure p ( x, t ) (mmHg), andcross-sectional area A ( x, t ) (cm ) in each large vessel by enforcing conservation of volume ∂A∂t + ∂q∂x = 0 , (3)and the balance of momentum ∂q∂t + ∂∂x (cid:18) q A (cid:19) + Aρ ∂p∂x = − πνRδ qA , (4)5here ρ = 1 .
055 (g/mL) is the blood density and ν (cm /s) is the kinematic viscosity. The right hand sideof equation (4) is derived under the assumption that the fluid velocity profile, u x ( r, x, t ) = q/A (cm/s), hasa linearly decreasing boundary layer thickness δ (cm), i.e., u x = (cid:40) ¯ u x , ≤ r ≤ R − δ ¯ u x / ( R − δ ) , R − δ ≤ r ≤ R. (5)The boundary layer δ is approximated by (cid:112) νT / π , where T (s) is the cardiac cycle length [44]. Under thisassumption, the shear stress in the large vessels is τ ( x, t ) = − µ ∂u x ∂r (cid:12)(cid:12)(cid:12) r = R = µ ¯ u x δ . (6)To close the system of equations, we impose a pressure-area relation governing the interaction between thevascular wall and the blood. Under the assumption that vessels are thin-walled, linearly elastic, homogenous,and isotropic [32, 27], we get p ( x, t ) = p + 43 Ehr (cid:32)(cid:114) AA − (cid:33) , (7)where p (mmHg) is the reference pressure at which A = A = πr , E (mmHg) is Young’s modulus in thecircumferential direction, and h (cm) is the wall thickness. For clarity, simulations are conducted in cgsunits but the results are provided in mmHg using the conversion 1mmHg = 1333 .
22 g/cm/s .Previous investigations [45] suggest that pulmonary vascular wall properties, e.g., Young’s modulus ( E )and wall thickness ( h ), change with vessel caliber. To account for these varying wall properties, we assumethat the wall-stiffness Eh/r can be represented by Ehr = k la e r k la + k la (8)as described by Olufsen [32]. Here k la (mmHg), k la (cm − ), and k la (mmHg) are stiffness parameters forthe large arteries. Normotensive values of k la = 2 . × , k la = −
15, and k la = 6 . × are based on ourprevious investigation [40, 42].The nonlinear system of PDEs is hyperbolic with characteristics pointing in the opposite direction; thus,we require boundary conditions at each vessel inlet and outlet ( x = 0 and x = L , respectively). We digitize anaverage control flow waveform, available from Simvascular, using GraphClick to provide an inflow boundarycondition at the inlet of the MPA. The cardiac output is scaled to match the flow values used in [35]. Atnon-stenotic vessel junctions, we prescribe continuity of static pressure and conservation of flow p p ( L p , t ) = p d (0 , t ) = p d (0 , t ) , (9) q p ( L p , t ) = q d (0 , t ) + q d (0 , t ) (10)which holds ∀ t ∈ [0 , T ]. Lastly, we couple the terminal large arteries to the structured tree model, asdescribed later. As discussed earlier, CTEPH progression is attributed to multiple emboli and vascular proliferation,which are grouped into five types [3]: ring-like, web-like, subtotal occlusion, total occlusive, and tortuouslesions. In this work, we model ring-like lesions (accounting for approximately 30% of all lesions) andweb-like lesions (accounting for approximately 64% of all lesions) [3].
6e model the former lesion type using the pressure loss term derived from Young and Tsai’s work [46]∆ p ring = µK v π ( r p ) q + ρK t A p ) (cid:18) A p A s − (cid:19) | q | q + ρK u L s A p ∂q∂t , (11)where K v , K t , and K u (dimensionless) are the loss coefficients due to viscous, turbulent, and inertial forces,respectively. Moreover, r p (cm) and A p (cm ) denote the unobstructed radius and area proximal to thelesion, and L s is the length of the ring-like lesion. This stenosis model has been implemented in 1D modelsand validated against 3D models in the systemic circulation [27, 28, 47]. The viscous loss K v , related to thestenosis geometry, is given by [27] K v = 16 L s r p (cid:18) A p A s (cid:19) . (12)Ring-like lesions are typically located proximal to vascular junctions [3], hence we prescribe these lesionsand their energy loss proximal to vessel junctions.We model web-like occlusions as semi-porous structures, which only partially occlude the vascular lumen[11]. The pressure drop in web-like lesions is modeled as∆ p web = L w (cid:18) µK perm qA + ρK w (cid:16) qA (cid:17) (cid:19) , (13)where L w (cm) is the length of the obstructed area, K perm (dimensionless) is the permeability of theweb-like lesion, and K w (dimensionless) is the pressure drop due to inertial effects. Equation (13),theDarcy-Forchheimer equation, is used for fluid flow through porous media at moderate Reynold’s numbers(i.e. >
10) [48]. Previous work [48, 49] assumes that K w ≈ G/K / perm , where G (dimensionless) depends onthe material in question. Since material properties are unknown, we assume G = 0 .
10 for all simulationsinvolving web-like lesion and place them at the vessel midpoint as these lesions can extend throughout asingle vessel [11].We compute the FFR for each lesion to quantify the pressure-losses due to lesion type and severity. Asdone in prior coronary studies [50], the FFR is calculated asFFR = ¯ p dist ¯ p prox (14)where ¯ p prox and ¯ p dist are the mean pressure proximal to and distal to the lesion, respectively. While inertial effects dominate large vessel hemodynamics, blood flow in the small arteries and arteriolesis dominated by viscous forces. For this reason, we solve linearized forms of equations 3 and 4. The linearizedequations are converted to a discrete, frequency domain Fourier series, resulting in a wave-equation. This isderived in Appendix A, and leads to analytical pressure and flow solutions, P ( x, ω ) and Q ( x, ω ), for eachfrequency ω , and are related via the impedance, Z ( x, ω ) = P ( x, ω ) /Q ( x, ω ). The input impedance at x = 0can then be written as Z (0 , ω ) = ig − sin( ωL/c ) + Z ( L, ω ) cos( ωL/c )cos( ωL/c ) + igZ ( L, ω ) sin( ωL/c ) (15)for any frequency ω (cid:54) = 0, whereas the zero frequency value, Z (0 , ω → Z (0 , ω ) = 8 µ(cid:96) rr πr + Z ( L, . (16)The above equations depend on the frequency domain quantities g ω and c , which are defined and relegatedto Appendix A. These terms depend on Eh/r , which has the same form as equation 8, but with parameters k sa , k sa , and k sa . 7sing the impedance relation, the pressure and flow solutions at x = L are [35] P ( L, ω ) = P (0 , ω ) cos ( ωL/c ) − Q (0 , ω ) ig ω sin ( ωL/c ) , (17) Z ( L, ω ) = ig − ω sin( ωL/c ) − Z (0 , ω ) cos( ωL/c ) g ω cos( ωL/c ) Z (0 , ω ) − cos( ωL/c ) , (18) Q ( L, ω ) = P ( L, ω ) /Z ( L, ω ) , (19)for ω (cid:54) = 0. For the zeroth frequency, the solution is P ( L,
0) = P (0 , − Q (0 ,
0) 8 µLπr , (20) Z ( L,
0) = Z (0 , − µLπr (21) Q ( L,
0) = P ( L, /Z ( L, . (22)This enables us to predict pressure and flow in any daughter vessel provided the pressure and flow of theparent are known.The non-Newtonian effects of blood viscosity µ (and hence kinematic viscosity ν ) become more influentialas the vessel radii decrease towards the arterioles. Following [51], and our previous study [34], we model theviscosity in the small arteries and arterioles by a nonlinear radius dependent function of the form µ ∗ ( r ) = (cid:104) µ . − (cid:18) − Hct ) C − − . C − (cid:19) D (cid:105) D , (23) µ . ( r ) = 6 e − . r + 3 . − . e − . r . , (24)where Hct is the blood hematocrit, D = (cid:16) r r − . (cid:17) , µ . ( r ) is the relative viscosity at an average hemat-ocrit level of 0 .
45, and C is a shape parameter derived from a nonlinear, empirical formula [51]. We assumea typical healthy blood hematocrit here, Hct = 0 . Z p ( L, ω k ) = (cid:18) Z d (0 , ω k ) + 1 Z d (0 , ω k ) (cid:19) − . (25)At the terminating end of the structured tree (i.e., Z ( L,
0) at the last generation branches), the terminalimpedance Z trm =constant is prescribed. Since the structured tree results in thousands of vessels, the flowat the end of the structured tree’s are on the order of 10 − − − . For this reason, the capillary pressure, P trm = Z trm · Q trm , requires the magnitude of Z trm to be on the order of 10 − to achieve non-negligiblepressure values. Once we reach the terminal structured tree branch, all proximal impedance values in thetree are computed. Knowing the impedance at the root of the tree provides an interface between the largeand small vessels. The large artery predictions are then used to solve for pressure and flow throughout thestructured tree using equations 17 and 19. Large vessel hemodynamics (equations (3, 4, 7)) are solved using the Ritchmeyer two-step Lax-Wendroffscheme [32], which is second order accurate in space and time. At non-stenotic junctions, pressure continuityand flow conservation are enforced via a Newton-Raphson root finding routine. Pressure losses due toring-like lesions, located proximal to bifurcations, are also solved using a Newton-Raphson routine. Forweb-like lesions, located at vessels’ midpoint, we minimize residual equations using gradient-descent with abacktracking line search algorithm [52]. To increase numerical accuracy at the junctions and lesion locations,we derive and utilize the analytical Jacobian of the residuals.8o couple the large vessels to the structured tree, the frequency impedance is convolved with the pressureat the end of the large arteries. This is numerically approximated as q nL = ∆ t n (cid:88) k =1 p n − kL (cid:0) Z k (cid:1) − , (26)where q nL and p n − kL are the numerical solutions of flow and pressure at the end of the large terminal arteries,and Z k is the vascular impedance of the whole structured tree calculated using equation (15). Wave intensity analysis (WIA) is a time-domain decomposition of forward and backward waves and isutilized in understanding the progression of cardiovascular disease [40, 42, 53]. The so-called pressure andvelocity ‘wavefronts’ [54] δp ± and δu ± are calculated as δp ± = δp ± ρc ( p ) δu i , δu ± = 12 (cid:18) δu ± δpρc ( p ) (cid:19) , (27)where the notation ± indicates the forward (‘+’) or backward (‘-’) running wave components, and δp and δu are the incremental changes in pressure and velocity, respectively. The wave speed c ( p ) is calculated fromthe pressure-area relation c ( p ) = (cid:115) Aρ ∂p∂A = (cid:115) Ehr ρ (cid:18) AA (cid:19) / . (28)Both forward and backward waves are classified as compression ( δp ± >
0) or decompression ( δp ± <
0) waves.Forward compression waves (FCW) originate from the right ventricle, increasing pressure and acceleratingflow as blood travels towards the periphery. Forward expansion waves (FEW) also originate in the proximalarteries but decrease the pressure and flow velocity. In contrast, backward compression waves (BCW) andexpansion waves (BEW) propagate from the periphery, the former increasing pressure while deceleratingflow, while the latter does the opposite. The wave intensity W I ± is calculated as a product of the velocityand pressure wavefronts. Since wavefronts are not continuous but discrete approximations, we use thetime-normalized wave intensity W I ± = δp ± δt δu ± δt . (29)After determining the type and nature of the waves, a wave reflection coefficient, R , can be calculated. Inthis study, we use R = ∆ p − ∆ p + , (30)where ∆ p + and ∆ p − are the sums of the backward and forward compression waves, respectively[55]. We map flow predictions from the large arteries to the 3D geometry [34] to predict lung perfusion. Todo so, we first determine which vessels are closest to the lung tissue using a nearest-neighbor algorithm inMATLAB (Mathworks, Natick, MA). Then, we project the 1D flow predictions onto the lung tissue andgenerate a lung perfusion map. These predictions are converted to a probability distribution function (PDF)using kernel density estimation; we use the ksdensity function in MATLAB. To quantify flow heterogeneityin control, CTEPH, and post-surgical intervention simulations, we compute the KL divergence, a measureof relative entropy between two probability distributions [36], for the flow PDFs.The KL divergence is a probabilistic measure between two distributions [36]. Let A and B be discreteprobability distributions defined on some probability space W , where A is our reference or target distribution.Then, the KL-divergence is defined as D KL (cid:16) A (cid:12)(cid:12)(cid:12) B (cid:17) = (cid:88) w ∈W A ( w ) log (cid:18) A ( w ) B ( w ) (cid:19) , (31)9here D KL → D KL indicates a bigger mismatchbetween the two distributions.In summary, we compute the perfusion map and quantify flow heterogeneity as follows:(1) project the 1D flow predictions to the centerline network in 3D;(2) determine the terminal vessels closest to the segmented lung tissue via a nearest-neighbor algorithm;(3) prescribe the flow to the entire portion of tissue;(4) construct the flow PDF, B , via kernel density estimation; and(5) compute D KL (cid:16) A (cid:12)(cid:12)(cid:12) B (cid:17) . To simulate CTEPH, we add ring-like (concentric stenosis) and web-like (longitudinal narrowing withina vessel) [3] lesions to the proximal pulmonary arteries. We also introduce arterial vasoconstriction andprogressive remodeling of the distal vasculature [12]. The former is done by decreasing the radii of vesselssegmented from the CT image, while the latter is done by increasing terminal PVR ( Z trm ) and vessel stiffnessin the structured trees. We set up a normotensive simulation and four disease cases and examine changes inhemodynamic predictions. These disease cases are summarized in Table 1 and parameter values are given inTable 2. To simulate the effects of BPA, we reduce both stenosis and downstream vasoconstriction severityin three lesions and quantify improvements in hemodynamics after BPA. This scenario represents a CTEPH patient without pulmonary vascular remodeling,illustrating the sole effects of arterial obstructions (ring- and web-like lesions). We add a total of 20 lesionsin the network: 9 ring-like and 11 web-like. Ring-like lesions have an area reduction of 90% and spreadover 25% of the original vessel. These lesions are placed proximal to vessel junctions [3]. Web-like lesionshave a permeability parameter K perm = 0 .
009 and a lesion length parameter, which is 50% of the vessel’slength [3]. Vessels that contain either lesion type are ten fold stiffer to account for wall stiffening due tothromboembolic material.
Case b: Lesions with local remodeling.
In this scenario, the disease has progressed downstream to thelesions. Remodeling narrows the vascular lumen and moderately thickens the vascular walls, decreasingvessel compliance. Parameters k , k , and k are increased in all vessels (both large and small) distal to thelesion while reducing the reference vessel radius r by 40%. We also increase pulmonary vascular resistanceby increasing Z trm in terminal vessels that are narrowed. Case c: Lesions with small vessel arteriopathy.
Prolonged exposure to CTEPH is associated with secondaryarteriopathy throughout the pulmonary vasculature, affecting vessels with cross-sectional areas ≤ r in all vessels with a cross-sectional area < by 40% and increasing vessel stiffness.We also increase Z trm in the structured trees. Case d: Lesions with proximal dilation and global small vessel arteriopathy.
In the most severe cases ofCTEPH, small vessel arteriopathy causes large artery dilation and stiffening [3, 5]. This scenario includesthe global small vessel remodeling described in case c, with dilation and thickening of arteries with cross-sectional areas ≥ by 10%, as documented previously [12]. We further increase stiffness and raisePVR in all terminal vessels. 10isease scenario Description Modeling approachControl Typical pulmonary vasculature - Low stiffness- Minimal PVR(a) Lesions only Mechanical obstructions due to - 90% area reductionring-like and web-like lesions in ring lesions-0 .
05 permeability inweb-like occlusions(b) Lesions with local Increased vascular stiffening -Increased stiffness inremodeling and narrowing in all vessels distal vessels distal to lesionto lesion locations, and increased -Reduce radius by 40% invascular resistance vessels distal to lesion-Increase PVR in terminalvessels distal to lesion(c) Lesions with global small Narrowing, stiffening, and - Reduce radius by 40% invessel arteriopathy increased resistance of all vessels with A < small vessels - Increased stiffness innarrowed vessels- Increased PVR inaffected terminal vessels.(d) Lesions with proximal Narrowing and stiffening of - Small vessel remodeling from (c)remodeling and small small vessels, dilation and - Dilation of proximal arteriesvessel arteriopathy stiffening of large proximal with A ≥ vessels, and increased PVR. - Increased network stiffness- Increased PVR in all terminal vessels Table 1: Summary of disease scenarios simulated in this study. isease case Parameter Normotensive a b c d k la (g/(cm s )) 2 . × . × . × . × . × k la (1/cm) − − − − − k la (g/(cm s )) 6 . × . × . × . × . × − − (cid:0) . × (cid:1) (cid:0) . × (cid:1) (cid:0) . × (cid:1) k sa (g/(cm s )) 2 . × . × . × . × . × k sa (1/cm) − − − − − k sa (g/(cm s )) 8 . × . × . × . × . × − − (cid:0) . × (cid:1) (cid:0) . × (cid:1) (cid:0) . × (cid:1) k lesion (g/(cm s )) − . × . × . × . × Z trm (g/(cm s) 1 . × . × . × . × . × − − (cid:0) . × (cid:1) (cid:0) . × (cid:1) (cid:0) . × (cid:1) Table 2: Parameter values for simulating the five different disease scenarios described in section 2.7. k and Z trm values arelisted for vessels both in unobstructed and obstructed paths (in parentheses). Stenosis parameters Description Value K t (ND) Turbulent pressure loss coefficient 1 . K u (ND) Inertial pressure loss coefficient 1 . L s (cm) Ring-like stenosis length L i / K perm (cm ) Darcy loss coefficient for web-like lesion permeability 0 . G (ND) Inertial loss coefficient for web-like lesions 0 . L w (cm) Length of web-like lesion L i / Table 3: Parameter describing ring- and web-like lesions, and their values. .7.2. BPA therapy Inoperable CTEPH patients with inaccessible distal lesions for PTE are ideal candidates for BPA. Wesimulate BPA intervention by decreasing lesion severity simultaneously in three locations at a time. BPAdoes not completely remove the lesions; hence, we still include them in the model but reduce their severity.Specifically, we reduce the stenosed area A s and lesion length L s from 90% to 50% for ring-like lesions andincrease K perm (equation (13)) from 0 .
009 to 0 .
05 for web-like lesions. We also return vessel reference radii, r , to 90% of their original value in the obstructed pathways, mimicking decreased vasoconstriction due tolung reperfusion.Since multiple interventions must be performed to improve hemodynamic function, the most severe lesionsare prioritized for treatment to optimize increases in lung perfusion and decreases in mPAP, PVR, and wavereflections. This study uses simulations to select the “optimal” interventional option by quantifying flowheterogeneity, MPA pressure, wave reflection coefficients, along with post-BPA perfusion maps. We performall virtual procedures in disease case (d), and compare hemodynamic indices pre- and post-BPA.
3. Results
The large pulmonary arterial network includes 226 vessels including the MPA, LPA, and RPA, andall visible segmental and subsegmental branches. Centerlines constructed by VMTK are processed andimplemented in the 1D model. With the exception of one trifurcation, all network junctions are bifurcations.Vessels with a length ≤ .
25 mm (3 of the 226 vessels) are elongated to match the spatial resolution of thefluid dynamics solver.We simulate normotensive hemodynamics and four disease cases with increasing CTEPH severity (casesa-d). CTEPH is induced by including 20 lesions in the segmental and subsegmental vessels and by increasingthe stiffness of the large and small arteries. We compare large and small vessel hemodynamics, networkperfusion, and wave intensities for all simulations. Using case (d) as our baseline for BPA, we performvirtual treatments by reducing lesion severity and stenosis induced vasoconstriction in three lesions at atime. We select an optimal procedure based on improvements in flow heterogeneity, MPA pressure, and thewave reflection coefficient.
To understand how lesions affect pulmonary hemodynamics, we increase the stenosis area of ring-likelesions and reduce the permeability of the web-like lesions. We increase the relative stenosis area, A s /A p ,from 0% to 90%, and decrease K perm from 1 . . ∼
18 mmHg. Decreasing the permeability of web-like lesions is less influentail on MPApressure, raising it by approximately 0.5 mmHg. The minimum FFR for ring-like lesions is 0 .
16 when thearea reduction is 90%, while the minimum FFR in web-like lesions is 0 .
66 for K perm = 0 . K perm = 0 . K perm in web-like lesions alone does not alter perfusionpredictions substantially (results not shown). The flow PDF’s for each simulation are shown in Figure 3.These PDFs are constructed using kernel density estimation with a bandwidth of 0 .
009 and a reflectionboundary correction.
Normotensive and CTEPH predictions of pressure, flow, and shear stress at the midpoint of the MPA,LPA, and RPA are displayed in Figure 4. Simulations from case (a) (normotensive parameters with lesionsonly) result in a mean MPA pressure of 22.5 mmHg. The downstream stiffening and narrowing in case (b)13 FF R Area reduction K perm % % % % % (a) (b) P r e ss u r e ( mm H g ) .
003 0 .
005 0 . . . . . . . Figure 2: Predictions of mean MPA pressure (top), and fractional flow reserve (FFR, bottom) when increasing lesion severityin ring-like (left) or web-like (right) lesions. FFR predictions correspond to 9 ring-like lesions and 11 web-like lesions. Overall,ring-like lesions have a greater effect on MPA pressure and FFR than web-like lesions. % % % % K L D i v e r gen c e %
0% 10% 20% 30%40% 50% 60% 70%80% 90% R e l a t i v e C a r d i a c O u t pu t P D F Relative Cardiac Output
Figure 3: Perfusion maps for increasing ring-like stenosis severity. The top left panel (0%) shows the normotensive perfusionmap, while the bottom right perfusion map shows 90% severity. The Kullback–Leibler (KL) divergence, constructed from theflow probability density function (PDF), quantifies the flow heterogeneity. have a negligible effect on pressure from case (a). Including global small vessel stiffening in case (c) elevatessystolic and diastolic pressures slightly. In contrast, severe CTEPH, case (d), raises the pulse pressure (thedifference between systolic and diastolic values), with systolic pressures ≈
45 mmHg and diastolic pressures ≈ P r e ss u r e ( mm H g ) F l o w ( m l / s ) S hea r s t r e ss ( g / c m s ) MPA LPA RPA controlcase (a)case (b)case (c)case (d)
Figure 4: Pressure, flow, and shear stress in the main, left, and right pulmonary arteries in each disease scenario. stiffness, cases (b)-(d), dampens the oscillations. All simulations have a larger flow distribution to the RPAthan the LPA; normotensive predictions have a 43/57 flow split to the LPA/RPA, while in severe CTEPH,case (d), the ratio is 22/78. MPA and LPA Shear stress magnitude is greater in normotension than in anyof the CTEPH simulations, while RPA shear stress is largest in case (d).Figure 5 shows pressure predictions throughout the large pulmonary arteries in normotension, case (a),and case (d) at different times within the cardiac cycle. Yellow vessel segments denote lesion locations.Results show a dramatic rise in pressure proximal to lesions for cases (a) and (d), while pressure dropsdownstream to the lesions. For case (a), the average FFR is 0.33 (range [0 .
16 0 . .
88 0 . .
19 0 . .
81 0 . ∼
45 mmHg in vessels in the unobstructed path, whereas the terminalsystolic pressure is 10 mmHg in the obstructed pathway.15 ontrol Case (a) Case (d)
Pressure (mmHg) 400 10 20 30 50 MPA Inflow
Figure 5: Pressure predictions in normotension and disease cases (a) and (d) throughout the network at different time points.Yellow segments indicate regions where lesions are present. P r e ss u r e ( mm H g ) F l o w ( m l / s ) Obstructed pathUnobstructed path P r e ss u r e ( mm H g ) F l o w ( m l / s ) F l o w ( m l / s ) F l o w ( m l / s ) controlcase (a)case (b)case (c)case (d) MPA MPA TerminalTerminal
Figure 6: Representative obstructed (red) and unobstructed (blue) pathways in the pulmonary vasculature, along with predic-tions of pressure and flow for the normotensive and four CTEPH disease scenarios. A zoom in of the flow predictions in themost distal vessels of the tree are shown to the right. Results are consistent with other pathways in the network.
The structured tree model predicts pressure and flow in the smaller arteries and arterioles. Figure 7shows normotensive predictions as well as CTEPH simulations in six large terminal arteries: three thatare downstream to a lesion and three that are not. Mean pressure and flow are shown along the α and β pathways, the longest and shortest pathways in the structured tree.In the obstructed pathway, we see a significant drop in mean pressure and mean flow in case (a). Thisis consistent in both the α and β pathways. In contrast, flow redistribution to the unobstructed vesselscause an increase in pressure and flow. Vasoconstriction of the terminal arteries and stiffening of themicrovasculature in case (d) increases pressure, yet flow in this pathway does not increase. Moreover, themean pressure magnitudes in the unobstructed pathway are similar in both case (a) and case (d), thoughthe systolic pressures (not shown) increase from 30 mmHg to 40 mmHg. Figure 8 presents wave intensity profiles in the MPA, LPA, and RPA for the normotensive and fourdisease cases. The wave reflection coefficient for each case is reported in Table 4. Normotensive simulationshave a predominant FCW in all three branches, followed by a minor BCW. Overall, normotensive expansionwaves are negligible. The wave reflection coefficient in the MPA is 0.47, smaller than the LPA and RPA(0.76 and 0.89). Adding lesions in case (a) leads to a prominent second FCW in the MPA and LPA, whileRPA forward and backward compression waves have a bimodal profile that occur consecutively. Case (a)also produces a distinct bimodal BCW in both the MPA and LPA. In all three branches, both the forwardand backward expansion waves are more oscillatory. The wave reflection coefficient increases in case (a) inthe MPA and LPA, but decreases in the RPA. 17 b s t r u c t e d p a t h w ay Alpha pathway Beta pathway Alpha pathway Beta pathway Alpha pathway Beta pathway U nob s t r u c t e d p a t h w ay T1T2 T3T4 T5T6
Control Case (a) Case (d)
Terminal β -pathway α -pathway αα α β β β P r e ss u r e ( mm H g ) F l o w ( m l / s ) P r e ss u r e ( mm H g ) F l o w ( m l / s ) Radius (cm) T4 T5 T6 T1 T2 T3 (a) (b)(c)(d) Figure 7: Hemodynamic predictions in the α and β pathways of the structured tree in the normotensive case and CTEPHcases (a) and (d). (a) Large vessel network with terminal vessels in obstructed (T1-T3, red) and unobstructed (T4-T6, blue)pathways. (b) Schematic of the structured tree (grey branches) and the α and β pathways (red and blue paths, respectively.(c) Mean pressure and flow predictions in the α and β pathways in terminal vessels T1-T3 as a function of radius. (d) Meanpressure and flow predictions in the α and β pathways in terminal vessels T4-T6 as a function of radius. Small vessel stiffening in case (b) dampens expansion wave oscillations. The wave intensity profile issimilar to case (a) but shows an elevated FCW across the three branches. A diminished BCW magnitudedecreases the MPA and RPA reflection coefficients; the LPA reflection coefficient is the same in case (a) and(b). Downstream narrowing in case (c) has a negligible effect on the wave intensities, but slightly intensifiesBCW magnitudes, increasing the reflection coefficient in all three vessels. For case (d), the most severeCTEPH, the peak wave intensity increases for all four wave types and the bimodal FCW and BEW aregone. Though the compression wave peak increases, the wave reflection coefficient decreases in all threevessels, suggesting a greater increase in magnitude for the FCWs than the BCWs.18 -50510 10 -50510 10 controlcase (a)case (b)scenario (c)case (d)MPA LPA RPAFCW BEWFEW BCW W ave i n t e n s i t y ( W m s ) - - Time (s) -50510 10 case (c) Figure 8: Wave intensity analysis (WIA) at the midpoint of the main, left, and right pulmonary arteries (MPA, LPA, and RPA,respectively). Foward and backward running compression waves (FCW, BCW) are distinguished from forward and backwardrunning expansion waves (FEW, BEW). ave reflection coefficient Disease case MPA LPA RPAcontrol 0.47 0.76 0.89case (a) 0.74 0.88 0.86case (b) 0.68 0.88 0.83case (c) 0.74 0.90 0.86case (d) 0.49 0.83 0.67
Table 4: Wave reflection coefficient in the main, left, and right pulmonary arteries (MPA, LPA, and RPA) in each disease case.The coefficient is calculated using equation 30.
BPA therapy for inoperable CTEPH typically addresses 3-5 lesions per session; therefore, we target 3lesions for the virtual BPA. This provides C = 1140 possible lesion combinations, which are all analyzed.We choose the optimal BPA using a combination of MPA pressure, WIA results and flow heterogeneity(KL-divergence), and want to minimize all three quantities. These three metrics are shown in Figure 9for both the pre- (case (d)) and post-BPA simulations. Results show that the optimal treatment targetsthree ring-like lesions: two in the right lower lobe and one in the left upper lobe (here, anatomical left andright are opposite to the left and right sides of the figure). Lung tissue perfusion maps show reperfusion toboth regions, and the flow PDF shifts leftward towards a more uniform network flow. The KL-divergencedecreases from 3.05 to 2.95, signifying a smaller distance between the control and post-BPA flow PDF.Systolic MPA pressure decreases from 45 to 42 mmHg after BPA, and mean pressure decreases from 22mmHg to 19 mmHg, below the diagnostic PH cutoff. Lastly, wave intensity results illustrate a decreasein BCWs and an increase in FEWs post-BPA. The wave reflection index is reduced in all three proximalarteries after BPA, changing from 0.49, 0.83, and 0.67 to 0.48, 0.76, and 0.61 in the MPA, LPA, and RPA,respectively. The pre-BPA FFR in these lesions are 0.44, 0.32, and 0.30, which improve to 0.88, 0.76, and0.59, respectively, post-BPA.
4. Discussion
Computational hemodynamics is emerging as a useful clinical tool for cardiovascular disease. Whileseveral previous studies employ modeling to quantify treatment strategies for systemic artery disease [27, 50],less attention has been given to pulmonary vascular disease. In this study, we solve a patient-specific, 1Dhemodynamics model in the large pulmonary arteries constructed from a CT image. The large vessels arecoupled to a fractal-based small vessel model, enabling hemodynamic predictions at multiple scales. Ourprevious studies utilize a similar approach [43, 39, 34], yet this study is the first to utilize these methodswith pressure loss models to generate perfusion maps and study CTEPH and BPA treatment strategies.
We consider five scenarios, ranging from normotensive to severe CTEPH. Figure 4 shows that lesionsalone, case (a), elevate pressures above 20 mmHg and into the PH range. This agrees with results reported byBurrowes et al. [57], who use a steady flow, Poiseuille flow model to predict hemodynamics with and withoutpulmonary emboli. They argue that pulmonary vasoconstriction and stiffening are necessary to elevate meanarterial pressures above 25 mmHg. We observe similar trends in cases (b) and (c), which introduce narrowingand stiffening of the small vessels. Though cases (a)-(c) increase the mean MPA pressure, typical CTEPHpatients have a systolic MPA pressure >
55 mmHg [53]. We approach this in case (d), which also dilates andstiffens the large pulmonary arteries. In addition to elevated pressure, our results show that the LPA/RPAflow imbalance increases with disease severity. This parallels results by Spazzapan et al. [25], which comparesimulations between stenosed and non-stenosed geometries.WSS in the MPA and LPA decreases in case (a), while the RPA shear stress increases. An elevatedpressure without a change in wall stiffness leads to a larger area deformation, decreasing the flow velocity20 e r c en t o f c a r d i a c ou t pu t MPA pre-BPA post-BPARPA (a)(b) (c)(d) controlpre-BPApost-BPA
FCWBEWFEWBCW W ave i n t e n s i t y ( W m s ) - - F l o w P D F Percent of cardiac output P r ess u r e ( mm H g ) Time (s)Time (s)
LPA
Figure 9: Pre and post-BPA metrics used to select the most effective treatment strategies. A total of 1140 interventions areconsidered. (a) Network of vessels with untreated lesions (cyan) and lesions treated in the best BPA simulation (red) alongwith pre and post-BPA perfusion maps. The region of the tissue most affected by BPA is identified within the black, dash-dotboxes. (b) Flow probability distribution functions (PDFs) for the normotensive (black), pre-BPA (red), and post-BPA (blue).(c) Pressure in the MPA pre-BPA (red) and post-BPA (black). (d) Wave intensities in the MPA, LPA, and RPA, pre-BPA(solid line) and post-BPA (dotted line). and WSS in the MPA. Case (a) leads to a larger flow and WSS in the RPA but reduces the WSS in the LPA.Shear stress does not increase significantly in case (d), a result of the proximal artery dilation which increasesarea (in the denominator of equation (6). Results by Yang et al. [26] show a reduced time-averaged wallshear stress in pediatric patients with severe PH, congruent with our findings. Another study by Tang etal. [58] shows that proximal and distal WSS magnitude is smaller in PAH when compared to normotensivemeasurements, a consequence of MPA, LPA, and RPA dilation, also shown in this study. Though themechanisms of PAH and CTEPH are different, both lead to proximal artery dilation [12, 26] and distalvascular narrowing [4, 12], hence a comparison of our results with PAH and CTEPH is appropriate.A benefit of our model is the ability to project 1D network predictions onto surrounding lung tissue.Figure 3 shows a conversion of 1D model simulations to a 3D flow map, providing a meaningful and qualitativetool for physicians studying the effects of CTEPH. A recent study by Clark et al. [33] use a transmissionline model to simulate large vessel hemodynamics and subsequently construct a lung perfusion map, similarto our work. However, their study did not account for perfusion deficits due to pulmonary lesions nor did21t quantify flow heterogeneity, both of which are done here. Our model couples nonlinear, large arterialhemodynamics to both a linearized small vessel fluid dynamics model and a pressure loss term, capturingthe complex relationship between pulmonary lesions and pulmonary hemodynamics. Figure 3 demonstratesthe capability of our methodology in predicting large perfusion deficits in obstructed regions, as seen in vivo[11].Our ring-like lesion model is based on the work by Young and Tsai [46] utilized in prior studies. Thestudies by Bezerra et al. [28] and Yin et al. [27] show that the combination of 1D CFD modeling withthe aforementioned energy loss model is nearly identical to 3D stenosis model predictions in the coronarycirculation. Only one prior pulmonary study (Spilker et al. [29]) utilized this energy loss model with a 1Dframework in pigs; their results show that RPA stenosis leads to a drastic reduction of flow to the right lungand 25 mmHg increase in systolic pressure. However, none of these studies consider the effects of multipledistal lesions in an expansive pulmonary tree. In addition, we introduce a novel web-like lesion pressure lossterm based on the Darcy-Forchheimer equation [59], which is a significant contribution to modeling CTEPH.These lesion types have not been modeled before, but are characteristic of a majority of CTEPH lesions [3].Hemodynamics simulations throughout the pulmonary vasculature are feasible using our 1D CFD modelframework. Figure (5) shows the mean pressure in each large pulmonary artery in both normotension andCTEPH, cases (a) and (d), during the cardiac cycle. While control simulations shows a relatively consistentpressure throughout the vasculature at each time point, cases (a) and (d) illustrate a significant pressure dropdownstream to lesions. This pressure drop, quantified by the FFR, is measured in CTEPH patients in twoprevious studies [11, 60] using optical coherence tomography. These studies report FFR values ranging from0.22 to 0.90, similar to the FFR magnitudes provided here. Our results in Figure 2 show that ring like lesionshave a larger pressure drop than web-like lesions, but that both lesion types have varied effects on MPApressure. Pressure drops across different types of pulmonary lesions have only recently been documented[3, 11], and these studies suggest that web-like lesions are more common and easier to treat with BPA thanother lesion types. These studies also show that the changes in pressure and flow downstream of web-likelesions are more volatile than other lesion types. To date, the properties of web-like lesions have not beeninvestigated, but more information regarding the structure and composition of web-like lesions may be usedto inform our loss model. Overall, our model shows that both lesion types contribute to a decrease inpressure downstream, with ring-like lesions playing a larger role.Current hypotheses suggest that secondary arteriopathy is a determining factor in CTEPH recovery[4, 5, 13]. As shown in Figure 6, model predictions in both obstructed and unobstructed pathways areaffected in the different disease cases. Clearly, blood pressure drops downstream from the lesion whileresistance to flow increases, leading to underperfusion in these vessels. Conversely, flow and pressure bothincrease in the unobstructed pathways. The study by Lang et al. [61] claims that decreased flow, andhence decreased shear stress, promotes pulmonary arterial remodeling downstream to a lesion. The authorsalso argue that increased flow and shear stress in unobstructed pathways leads to flow vasculopathy and,consequentially, pulmonary vascular remodeling. This remodeling of both obstructed and unobstructedarterial pathways explains why all small arterial vessels, i.e. with an area ≤ , narrow in CTEPH[5, 12]. Our results agree with these hypotheses and illustrate underperfusion in obstructed pathways andhyperperfusion in unobstructed pathways. Predictions in the obstructed pathdo not return to their originalnormotensive values, but they gradually approach these values as we introduce stiffening and narrowingin cases (b) through (d). Simulations from case (d) are consistent with the physiological conclusions thatstiffening and narrowing of the pulmonary tree partially corrects pressure and flow imbalances due to CTEPHlesions [61]. Our multiscale approach utilizes the structured tree model to simulate arteriolar hemodynamics distalto the large subsegmental arteries. In CTEPH, lesions in segmental and subsegmental arteries decrease flowto the microvasculature, leading to small vessel remodeling [5, 23]. We simulate this in cases (b), (c), and(d) by decreasing the area of the terminal vessels and, consequently, the radii and number of branches in thestructured tree. Results in Figure 6 show hemodynamic predictions in three structured trees in an obstructedpath (T1, T2, and T3) and three structured trees in an unobstructed path (T4, T5, and T6). As expected,22ow and pressure decrease substantially in the obstructed pathway with the addition of vascular lesionsin case (a), whereas flow redistribution leads to an increase in flow and pressure in unobstructed terminalvessels and their arteriolar beds. Flow within the obstructed pathways does not increase when stiffness andPVR ( Z term ) are increased, though flow in the unobstructed pathway decreases as a result of stiffening.These results and methods agree with previous physiological studies. For example, the study by Stam et al.[22] induced CTEPH in swine, showing that animals with pulmonary lesions have increased wall-thicknessand decreased microvascular luminal area compared to control. We model both these phenotypes of CTEPHin cases (b), (c), and (d), increasing MPA pressure. The review by Lang et. al [5] also supports the ideathat an increase in flow in unobstructed pathways leads to secondary arteriopathy and wall thickening inthe pulmonary arterioles.Decreasing the radii of terminal arteries leads to a smaller arteriolar tree, shown in Figure 7 for case (d),as r min is held constant. A reduced microvascular density, also called microvascular rarefaction, is a knownconsequence of both CTEPH and PAH [62]. An imaging study by Come et al. [12] concludes that CTEPHpatients have a reduced volume distribution of small arteries when compared to control, and suggests pruningof the pulmonary arterioles. This phenomenon is modeled by Olufsen et al. [35] who simulate rarefactionof the pulmonary microvasculature and consequently predict increased mean and diastolic pressure in theMPA. The study by Yang et al. [26] uses a 3D-0D coupled model of the proximal and distal vasculatureand predicts increased right ventricular pressure and distal WSS as a consequence of the pruned pulmonarymicrovasculature. However, these distal vascular predictions, as well as our own here, are hypothetical, ascurrent imaging technologies cannot capture dynamics in vessels of this magnitude. This model type cantest large and small vessel disease hypotheses, and can be validated if measurements are available in boththe large and small pulmonary arterial vessels. WIA is increasingly recognized as a measure of pulmonary vascular function [15, 53]. Results in Figure8 show how the addition of lesions and vascular stiffening affects forward and backward traveling waves.Results from the normotensive simulations show a large FCW in the MPA during systole followed by asmaller magnitude BCW and FEW, which rapidly decay in magnitude during diastole, and that the BEWis relatively negligible. In contrast, wave intensity magnitudes in scenario (d) nearly double in all four wavetypes. This agrees with the study by Su et al. [53], who measure area-velocity signals in both normotensiveand PH patients, and show that wave intensity magnitudes nearly double in PH. Similarly, they show thatFCWs in PH are approximately twice the magnitude of FEWs and BCWs, and that BEWs do not dissipateas quickly as in the normotensive patients. A similar conclusion is drawn in the study by Lau et al. [63];PAH patients in the study have a much larger BCW than control subjects, ultimately leading to an increasedwave reflection coefficient. In this study, we define the wave reflection coefficient as the sum of forward andbackward compression waves in the time domain. However, the wave reflection coefficient can be calculatedusing other methods [53, 55], including impedance analysis, either in the frequency or time domain [64].Our results suggest that wave reflection coefficients increase with disease severity. However, dilatingproximal vessels leads to a decreased reflection coefficient in the MPA, LPA, and RPA, as shown in Table 4.Though the reflection index is decreased in case (d), BCW magnitude increases significantly. This is shownto correlate with right ventricular dysfunction in the study by Schafer et al. [65]. The authors’ of this studyalso argue that BCWs are more indicative of increased proximal stiffening, agreeing with our results in case(d) after stiffening the large proximal arteries. The alignment between our disease scenarios and previouslypublished in-vivo WIA suggest that this model framework is suitable for addressing wave-propagation inPH and, more specifically, CTEPH. It is reported that proximal arteries both stiffen and dilate in longterm PH [12], which is accounted for in case (d) in this study. Our simulations show that this physiologicalphenomenon decreases the wave reflection coefficient. This suggests that dilation of the MPA attempts toreduce the load on the heart by minimizing the mismatch between right ventricular ejection and reflectedwaves from the arterial periphery. 23 .4. Computational treatment planning
For inoperable CTEPH patients, a combination of drug therapy and BPA is the best alternative forimproving vascular function. BPA is a physician dependent strategy and varies with both lesion locationand severity of the disease [10]. One common treatment strategy is to target the lobe with the largestperfusion deficit [66], yet there may be multiple lesions within each lobe. For this reason, an integratedmathematical model of CTEPH hemodynamics with patient-specific imaging can rank lesion importance inprocedural planning. Image-based computational hemodynamics modeling is already recognized as a usefultool in surgical planning for coronary artery disease [27, 50], and this study is a first step in utilizing thesame framework for understanding CTEPH and hemodynamic improvement after BPA intervention.We utilize three indices to determine the best treatment strategy: the mean MPA pressure, the wavereflection coefficient in the MPA, and the KL-divergence of the flow field. The former two quantify how BPAaffects vascular-ventricular coupling and ventricular afterload, while the latter quantifies improvements inperfusion. Our results show that a combination of all three indices leads to a best treatment for improvingproximal hemodynamics and lung perfusion. Within our treatment framework, we reduce stenosis severitythrough A s for ring-like lesions or K perm for web-like lesions, but do not remove lesions. This is important toconsider, as BPA does not remove the lesions completely [5, 11]. In addition to targeting the lesions directly,we also decrease the degree of vasoconstriction in vessels downstream from a lesion, returning their radiito 90% of their original value. This is supported by the study by Boulate et al. [67] who show that distalartery thickness in pigs can regress to typical values seen in shams after removing lesions. Their study alsoanalyzes arterial thickness 6 weeks after reperfusing the obstructed lung tissue, demonstrating that surgicalremoval of clots reestablishes a normal lumen area after surgery.Results shown in Figure 9 illustrate that the combination of pressure predictions, WIA, and flow het-erogeneity lead to an optimal BPA treatment. The optimal treatment here reperfuses the right lower lobeand restores some flow to the periphery in the right middle lobe and left upper lobe. Though physicianstypically target the lobes of the lung with the largest perfusion deficit seen on V/Q scans [5], advances innew imaging technologies (such as optical coherence tomography) also assist in determining which lesionswill best improve hemodynamics [3, 14]. Our virtual BPA targets a region that is underperfused, but doesnot identify the most flow deficient region (the left lower lobe) as the optimal location. There are severallesions located in this lobe, hence addressing only 3 lesions has little effect on redistributing the flow. Ouroptimal BPA treatment also identifies proximal lesions rather than distal lesions, consistent with clinicalpractice [5]. Panel (b) in Figure 9 displays the relative cardiac output PDF for the normotensive controland pre- and post-BPA simulations, which are used to compute the KL-divergence. The control PDF has abimodal structure, signifying two levels of perfused tissue as seen in Figure 3. Both pre- and post-BPA havea similar trend, but show a small area of lung that is hyperperfused relative to normotensive simulations.Our simulation results show that both MPA pressure and wave intensities change only minutely inresponse to the single BPA treatment. However, improvements in hemodynamics post-BPA are typicallyreported over the course of months or years and after multiple sessions [5, 68]. Improvements immediatelyafter BPA have only been reported in terms of FFR improvement [11, 66]. Our optimal BPA treatmentincreased FFR from 0.44, 0.32, and 0.30 to 0.88, 0.76, and 0.59, respectively. These improvements aresimilar to recordings from the study by Ishiguro et al. [11], which discover an improved FFR after BPA.Their study shows that FFR increases from 0.22 to 0.59 and 0.34 to 0.86 in two different lesions, similar inmagnitude to our results. The long-term effects of BPA likely include decreased vasoconstriction in smallarteries and normalized wall thickness [67], but only in the case of reversible PH. Recent investigations intosuccess rates in CTEPH treatment find that the severity of small vessel disease correlates negatively withsuccessful surgical outcomes, suggesting that the degree of small vessel arteriopathy dictates whether PH ispersistent after intervention [24]. Histological data on arterioles after BPA would provide insight into theremodeling process after surgery, and could be reflected in our simulations of post-BPA hemodynamics.
5. Limitations and future work
We acknowledge several limitation in this study. We compute the wave reflection coefficient as theratio of the cumulative forward and backward compression pressure waves. This is only one possible way of24omputing the wave-reflection coefficient in the time-domain [55], and can also be calculated in the frequencydomain [64]. It is unclear which index best captures features of the vasculature, yet all these metrics wouldsupport our conclusion that large vessel stiffening and the presence of pulmonary lesions increase BCWsand wave reflections.This study implements a novel web-like lesion energy loss model; however, it is unclear if this modelcaptures proximal and distal hemodynamic features in human web-like lesions. Future studies will be devotedto validating this model type by obtaining catheter measurements prior to and after web-like lesions [3].Moreover, we plan to investigate the sensitivity of model outputs, such as the wave reflection coefficient,mean MPA pressure, and flow heterogeneity index, to lesion and boundary condition parameters. In thismanner, we will determine which parameters influence outputs after BPA and use this to guide parameterinference in patient-specific modeling studies for patients with CTEPH.Though we provide hemodynamic predictions immediately after BPA and hypothesize some degree ofvasodilation after an intervention, it is unclear what the immediate effects of a single BPA session are. Datacollection immediately after single BPA procedures may provide more insight into the immediate remodelingprocess. To this point, we do not attempt to predict the long term remodeling of the pulmonary vasculatureafter intervention, which will be investigated in future studies via a growth and remodeling framework [69].Lastly, this proof-of-concept study models CTEPH in a control geometry by placing lesions in the segmentaland subsegmental branches. Future studies will integrate CTEPH images, including lesion location, type,and severity.
6. Conclusions
This study provides a framework that integrates patient-specific CT imaging with both large and smallvessel fluid dynamics to predict multiscale hemodynamics in CTEPH. We model normotensive hemody-namics and four CTEPH disease cases, the latter driven by physiological hypotheses including large andsmall vessel remodeling. We utilize two lesion models representing ring-like and web-like lesions, commonin CTEPH and imperative for simulation studies. Our results show that combining thromboembolic lesionswith pulmonary vascular remodeling increases pulmonary pressures to the CTEPH range and mimics clinicalobservations. We predict perfusion in the lung tissue and provide a novel, quantitative metric for measuringperfusion heterogenity, which is essential in understanding the link between flow deficits in the lung and dis-ease severity. WIA results from the model framework are akin to prior measurements of wave intensities inCTEPH patients. Small vessel predictions are in agreement with clinical knowledge of CTEPH progression,and illustrate that lesions lead to downstream underperfusion as well as hyperperfusion in unobstructedterritories. We propose a combination of indices predicted by the model, and utilize these in priortizinglesions for BPA therapy. Our modeling framework shows improvements in hemodynamics post-BPA, layingthe foundation for future patient-specific investigations and validation studies using CTEPH data.
7. Acknowledgements
This study was funded in part by American Heart Association Grant
Appendix A. Structured tree model
Writing pressure and flow as periodic, discrete, frequency domain functions, P ( x, ω k ) and Q ( x, ω k ), thelinearized mass conservation and momentum balance equations for each frequency ω k are given by iω k CP + ∂Q∂x = 0 , (A.1)25nd iω k Q + A (1 − F J ) ρ ∂P∂x = 0 . (A.2)The quantity F J is the quotient of first and zeroth order Bessel functions F J = 2 J (cid:16) rw r (cid:17) w J ( w ) , w = √ i w , (A.3)where w = (cid:112) r ω k /ν is the non-dimensional Womersley number.Equation (A.1) depends on the vessel compliance C , approximated by C = ∂p∂A ≈ r Eh , (A.4)where
Eh/r has the same form as equation 8, but with parameters k sa , k sa , and k sa .Differentiating equation (A.1) with respect to x and solving for P in equation (A.2) gives ω k c Q + ∂ Q∂x = 0 , c = (cid:115) A (1 − F J ) ρC , (A.5)where c (cm/s) is the pulse wave propagation velocity. Solving the above wave equation for Q and plugginginto equation (A.2) gives Q ( x, ω k ) = a cos( ω k x/c ) + b sin( ω k x/c ) ,P ( x, ω k ) = ig ( − a sin( ω k x/c ) + b cos( ω k x/c )) , (A.6)where a and b are unknown integration constants and g = cC = (cid:112) CA (1 − F J ) /ρ .Using these definitions, vascular impedance is computed as Z ( x, ω ) = P ( x, ω ) Q ( x, ω ) = i ( b cos( ωx/c ) − a sin( ωx/c )) g ( a cos( ωx/c ) + b sin( ωx/c )) . (A.7)The above formulation is used to compute impedance at the inlet ( x = 0) and outlet ( x = L ) of each smallartery in the structured tree Z (0 , ω ) = ig ba ,Z ( L, ω ) = i ( b cos( ωL/c ) − a sin( ωL/c )) g ( a cos( ωL/c ) + b sin( ωL/c )) . (A.8)Combining these equations gives the impedance at x = 0 Z (0 , ω ) = ig − sin( ωL/c ) + Z ( L, ω ) cos( ωL/c )cos( ωL/c ) + igZ ( L, ω ) sin( ωL/c ) . (A.9)At the end of the structured tree, the distal impedance is set equal to the terminal impedance, Z ( L,
0) = Z trm . This impedance is then used to back solve the impedance values all the way through the root-impedance at the terminal end of the large arteries. 26 eferences [1] G. Simonneau, D. Montani, D. S. Celermajer, C. P. Denton, M. A. Gatzoulis, M. Krowka, P. G. Williams, R. Souza,Haemodynamic definitions and updated clinical classification of pulmonary hypertension, The European respiratory journal53 (1) (2019). doi:10.1183/13993003.01913-2018 .[2] Z. Wang, N. C. Chesler, Pulmonary vascular wall stiffness: An important contributor to the increased right ventricularafterload with pulmonary hypertension, Pulmonary Circulation 1 (2011) 212–223. doi:10.4103/2045-8932.83453 .[3] T. Kawakami, A. Ogawa, K. Miyaji, H. Mizoguchi, H. Shimokawahara, T. Naito, T. Oka, K. Yunoki, M. Munemasa,H. Matsubara, Novel Angiographic Classification of Each Vascular Lesion in Chronic Thromboembolic Pulmonary Hy-pertension Based on Selective Angiogram and Results of Balloon Pulmonary Angioplasty, Circulation: CardiovascularInterventions 9 (10) (2016) 1–8. doi:10.1161/CIRCINTERVENTIONS.115.003318 .[4] D. T. Matthews, A. R. Hemnes, Current concepts in the pathogenesis of chronic thromboembolic pulmonary hypertension,Pulmonary Circulation 6 (2) (2016) 145–154. doi:10.1086/686011 .[5] I. Lang, B. C. Meyer, T. Ogo, H. Matsubara, M. Kurzyna, H. A. Ghofrani, E. Mayer, P. Brenot, Balloon pulmonaryangioplasty in chronic thromboembolic pulmonary hypertension, European Respiratory Review 26 (143) (2017). doi:10.1183/16000617.0119-2016 .[6] I. M. Lang, R. Pesavento, D. Bonderman, J. X.-J. Yuan, Risk factors and basic mechanisms of chronic thromboembolicpulmonary hypertension: a current understanding, European Respiratory Journal 41 (2) (2013) 462–468. arXiv:https://erj.ersjournals.com/content/41/2/462.full.pdf , doi:10.1183/09031936.00049312 .[7] H. A. Memon, H. C. Lin, A. Guha, Chronic thromboembolic pulmonary hypertension: Pearls and pitfalls of diagnosis,Methodist Debakey Cardiovasc J 12 (4) (2016) 199–204. doi:10.14797/mdcj-12-4-199 .[8] M. Sarah, S. Zeenat, Epidemiology and pathophysiology of chronic thromboembolic pulmonary hypertension: Risk factorsand mechanisms, Methodist Debakey Cardiovasc J 12 (4) (2016) 195–198. doi:10.14797/mdcj-12-4-195 .[9] P. F. Fedullo, W. R. Auger, K. M. Kerr, L. J. Rubin, Chronic thromboembolic pulmonary hypertension, New EnglandJournal of Medicine 345 (20) (2001) 1465–1472, pMID: 11794196. arXiv:https://doi.org/10.1056/NEJMra010902 , doi:10.1056/NEJMra010902 .[10] H. Wilkens, S. Konstantinides, I. M. Lang, A. C. Bunck, M. Gerges, F. Gerhardt, A. Grgic, C. Groh´e, S. Guth, M. Held,J. B. Hinrichs, M. M. Hoeper, W. Klepetko, T. Kramm, U. Kr¨uger, M. Lankeit, B. C. Meyer, K. M. Olsson, H. J. Sch¨afers,M. Schmidt, H. J. Seyfarth, S. Ulrich, C. B. Wiedenroth, E. Mayer, Chronic thromboembolic pulmonary hypertension(CTEPH): Updated Recommendations from the Cologne Consensus Conference 2018, International Journal of Cardiology272 (2018) 69–78. doi:10.1016/j.ijcard.2018.08.079 .[11] H. Ishiguro, M. Kataoka, T. Inami, N. Shimura, R. Yanagisawa, T. Kawakami, K. Fukuda, H. Yoshino, T. Satoh, Diversityof Lesion Morphology in CTEPH Analyzed by OCT, Pressure Wire, and Angiography, JACC: Cardiovascular Imaging9 (3) (2016) 324–325. doi:10.1016/j.jcmg.2015.02.015 .[12] C. E. Come, A. A. Diaz, F. N. Rahaghi, J. C. Ross, M. Agarwal, G. Gonz´alez, C. E. Come, A. A. Diaz, G. Vegas-S´anchez-Ferrero, A. Hunsaker, R. San Jos´e Est´epar, A. B. Waxman, G. R. Washko, R. S. J. Est´epar, A. B. Waxman, G. R.Washko, Pulmonary Vascular Morphology as an Imaging Biomarker in Chronic Thromboembolic Pulmonary Hypertension,Pulmonary circulation 6 (1) (2016) 0–11. doi:10.1086/685081.Chronic .[13] M. Delcroix, I. Lang, J. Pepke-Zaba, P. Jansa, A. M. D’Armini, R. Snijder, P. Bresser, A. Torbicki, S. Mellemkjaer,J. Lewczuk, et al., Long-term outcome of patients with chronic thromboembolic pulmonary hypertension: results from aninternational prospective registry, Circulation 133 (9) (2016) 859–871.[14] M. S. Khan, E. Amin, M. M. Memon, N. Yamani, T. J. Siddiqi, S. U. Khan, M. H. Murad, F. Mookadam, V. M.Figueredo, R. Doukky, et al., Meta-analysis of use of balloon pulmonary angioplasty in patients with inoperable chronicthromboembolic pulmonary hypertension, International journal of cardiology 291 (2019) 134–139.[15] P. Segers, S. Moledina, D. S. Knight, M. A. Quail, A. M. Taylor, G. J. Coghlan, J. A. Steeden, L. Taelman, V. Muthurangu,Noninvasive pulmonary artery wave intensity analysis in pulmonary hypertension, American Journal of Physiology-Heartand Circulatory Physiology 308 (12) (2015) H1603–H1611. doi:10.1152/ajpheart.00480.2014 .[16] F. Vanden Eynden, B. El-Oumeiri, T. Bov´e, G. Van Nooten, P. Segers, Proximal pressure reducing effect of wave reflectionin the pulmonary circulation disappear in obstructive disease: insight from a rabbit model, American Journal of Physiology-Heart and Circulatory Physiology 316 (5) (2019) H992–H1004. doi:10.1152/ajpheart.00635.2018 .[17] J. M. Wells, A. S. Iyer, F. N. Rahaghi, S. P. Bhatt, H. Gupta, T. S. Denney, S. G. Lloyd, L. J. Dell’Italia, H. Nath, R. S. J.Estepar, G. R. Washko, M. T. Dransfield, Pulmonary artery enlargement is associated with right ventricular dysfunctionand loss of blood volume in small pulmonary vessels in chronic obstructive pulmonary disease, Circulation: CardiovascularImaging 8 (4) (2015) 1–8. doi:10.1161/CIRCIMAGING.114.002546 .[18] T. M. Fernandes, D. S. Poch, W. R. Auger, Treatment of chronic thromboembolic pulmonary hypertension: The roleof medical therapy and balloon pulmonary angioplasty, Methodist Debakey Cardiovasc J 12 (4) (2016) 205–212. doi:10.14797/mdcj-12-4-205 .[19] S. Quadery, A. Swift, C. Billings, A. Thompson, C. Elliot, J. Hurdman, A. Charalampopoulos, I. Sabroe, I. Armstrong,N. Hamilton, P. Sephton, S. Garrod, J. Pepke-Zaba, D. Jenkins, N. Screaton, A. Rothman, A. Lawrie, T. Cleveland,S. Thomas, S. Rajaram, C. Hill, C. Davies, C. Johns, J. Wild, R. Condliffe, D. Kiely, The impact of patient choiceon survival in chronic thromboembolic pulmonary hypertension, European Respiratory Journal (2018). doi:10.1183/13993003.00589-2018 .[20] A. Kandathil, M. Chamarthy, Pulmonary vascular anatomy & anatomical variants, Cardiovascular Diagnosis and Therapy8 (3) (2018) 201–207. doi:10.21037/cdt.2018.01.04 .
21] M. M. Madani, Surgical treatment of chronic thromboembolic pulmonary hypertension: Pulmonary thromboendarterec-tomy, Methodist Debakey Cardiovasc J 12 (4) (2016) 213–218. doi:10.14797/mdcj-12-4-213 .[22] K. Stam, R. W. van Duin, A. Uitterdijk, I. Krabbendam-Peters, O. Sorop, A. H. Jan Danser, D. J. Duncker, D. Merkus,Pulmonary microvascular remodeling in chronic thrombo-embolic pulmonary hypertension, American Journal of Physiol-ogy - Lung Cellular and Molecular Physiology 315 (6) (2018) L951–L964. doi:10.1152/ajplung.00043.2018 .[23] T. Jujo, S. Sakao, H. Ishibashi-Ueda, K. Ishida, A. Naito, T. Sugiura, A. Shigeta, N. Tanabe, M. Masuda, K. Tatsumi,Evaluation of the microcirculation in chronic thromboembolic pulmonary hypertension patients: The impact of pulmonaryarterial remodeling on postoperative and follow-up pulmonary arterial pressure and vascular resistance, PLoS ONE 10 (8)(2015) 1–18. doi:10.1371/journal.pone.0133167 .[24] C. Gerges, M. Gerges, R. Friewald, P. Fesler, P. Dorfm¨uller, S. Sharma, K. Karlocai, N. Skoro-Sajer, J. Jakowitsch,B. Moser, S. Taghavi, W. Klepetko, I. M. Lang, Microvascular Disease in Chronic Thromboembolic Pulmonary Hy-pertension: Hemodynamic Phenotyping and Histomorphometric Assessment, Circulation (2020) 376–386 doi:10.1161/CIRCULATIONAHA.119.041515 .[25] M. Spazzapan, P. Sastry, J. Dunning, D. Nordsletten, A. de Vecchi, The use of biophysical flow models in the surgicalmanagement of patients affected by Chronic Thromboembolic Pulmonary Hypertension, Frontiers in Physiology 9 (MAR)(2018) 1–13. doi:10.3389/fphys.2018.00223 .[26] W. Yang, M. Dong, M. Rabinovitch, F. P. Chan, A. L. Marsden, J. A. Feinstein, Evolution of hemodynamic forces inthe pulmonary tree with progressively worsening pulmonary arterial hypertension in pediatric patients, Biomechanics andModeling in Mechanobiology 18 (3) (2019) 779–796. doi:10.1007/s10237-018-01114-0 .[27] M. Yin, A. Yazdani, G. E. Karniadakis, One-dimensional modeling of fractional flow reserve in coronary artery disease:Uncertainty quantification and Bayesian optimization, Computer Methods in Applied Mechanics and Engineering 353(2019) 66–85. doi:10.1016/j.cma.2019.05.005 .[28] C. G. Bezerra, P. A. Lemos, R. A. Feij´oo, L. O. M¨uller, G. D. M. Talou, P. J. Blanco, C. A. Bulant, Comparison of1D and 3D Models for the Estimation of Fractional Flow Reserve, Scientific Reports 8 (1) (2018) 1–12. doi:10.1038/s41598-018-35344-0 .[29] R. L. Spilker, J. A. Feinstein, D. W. Parker, V. M. Reddy, C. A. Taylor, Morphometry-based impedance boundaryconditions for patient-specific modeling of blood flow in pulmonary arteries, Annals of Biomedical Engineering 35 (4)(2007) 546–559. doi:10.1007/s10439-006-9240-3 .[30] R. Clipp, B. Steele, An evaluation of dynamic outlet boundary conditions in a 1D fluid dynamics model, MathematicalBiosciences and Engineering 9 (1) (2012) 61–74. doi:10.3934/mbe.2012.9.61 .[31] W. Huang, R. T. Yen, M. McLaurine, G. Bledsoe, Morphometry of the human pulmonary vasculature., Journal of appliedphysiology (Bethesda, Md. : 1985) 81 (5) (1996) 2123–2133. doi:10.1002/ar.1092200410 .[32] M. S. Olufsen, Structured tree outflow condition for blood flow in larger systemic arteries, American Journal of Physiology-Heart and Circulatory Physiology 276 (1) (1999) H257–H268. doi:10.1152/ajpheart.1999.276.1.H257 .[33] A. R. Clark, M. H. Tawhai, Temporal and spatial heterogeneity in pulmonary perfusion: a mathematical model topredict interactions between macro- and micro-vessels in health and disease, The ANZIAM Journal 59 (4) (2018) 562–580. doi:10.1017/S1446181118000111 .[34] M. J. Chambers, M. J. Colebank, M. U. Qureshi, R. Clipp, M. S. Olufsen, Structural and hemodynamic properties of murinepulmonary arterial networks under hypoxia-induced pulmonary hypertension, Proceedings of the Institution of MechanicalEngineers, Part H: Journal of Engineering in Medicine (2020) 095441192094411 doi:10.1177/0954411920944110 .[35] M. S. Olufsen, N. A. Hill, G. D. A. Vaughan, C. Sainsbury, M. Johnson, Rarefaction and blood pressure in systemic andpulmonary arteries, Journal of Fluid Mechanics 705 (2012) 280–305. doi:10.1017/jfm.2012.220 .[36] A. Birolleau, G. Po¨ette, D. Lucor, Adaptive bayesian inference for discontinuous inverse problems, application to hyperbolicconservation laws, Communications in Computational Physics 16 (1) (2014) 1–34. doi:10.4208/cicp.240113.071113a .[37] A. Fedorov, R. Beichel, J. Kalpathy-Cramer, J. Finet, J.-C. Fillion-Robin, S. Pujol, C. Bauer, D. Jennings, F. Fennessy,M. Sonka, et al., 3d slicer as an image computing platform for the quantitative imaging network, Magnetic resonanceimaging 30 (9) (2012) 1323–1341.[38] L. Antiga, M. Piccinelli, L. Botti, B. Ene-Iordache, A. Remuzzi, D. A. Steinman, An image-based modeling frameworkfor patient-specific computational hemodynamics, Medical & biological engineering & computing 46 (11) (2008) 1097.[39] M. Colebank, L. M. Paun, M. U. Qureshi, N. Chesler, D. Husmeier, M. S. Olufsen, L. E. Fix, Influence of image seg-mentation on one-dimensional fluid dynamics predictions in the mouse pulmonary arteries, Journal of The Royal SocietyInterface 16 (159) (2019) 20190284. doi:10.1098/rsif.2019.0284 .[40] M. U. Qureshi, G. D. A. Vaughan, C. Sainsbury, M. Johnson, C. S. Peskin, M. S. Olufsen, N. A. Hill, Numericalsimulation of blood flow and pressure drop in the pulmonary arterial and venous circulation, Biomechanics and Modelingin Mechanobiology 13 (5) (2014) 1137–1154. doi:10.1007/s10237-014-0563-y .[41] M. I. Townsley, Structure and composition of pulmonary arteries, capillaries, and veins, Comprehensive Physiology 2 (1)(2012) 675–709. doi:10.1002/cphy.c100081 .[42] M. U. Qureshi, M. Colebank, L. M. Paun, L. Ellwein Fix, N. C. Chesler, M. A. Haider, N. A. Hill, D. Husmeier, M. S.Olufsen, Hemodynamic assessment of pulmonary hypertension in mice: a model-based analysis of the disease mechanism,Biomechanics and Modeling in Mechanobiology 18 (1) (2019) 219–243. doi:10.1007/s10237-018-1078-8 .[43] M. Colebank, M. Umar Qureshi, M. S. Olufsen, Sensitivity analysis and uncertainty quantification of 1-D models ofpulmonary hemodynamics in mice under control and hypertensive conditions, International Journal for Numerical Methodsin Biomedical Engineering (March) (2019) 1–29. doi:10.1002/cnm.3242 .[44] S. J. Lighthill, Mathematical biofluiddynamics, SIAM, 1975.[45] S. R. Lammers, P. H. Kao, H. J. Qi, K. Hunter, C. Lanning, J. Albietz, S. Hofmeister, R. Mecham, K. R. Stenmark, . Shandas, Changes in the structure-function relationship of elastin and its impact on the proximal pulmonary arterialmechanics of hypertensive calves, American Journal of Physiology - Heart and Circulatory Physiology 295 (4) (2008)1451–1459. doi:10.1152/ajpheart.00127.2008 .[46] D. F. Young, F. Y. Tsai, Flow Characteristics in Models of Arterial Stenoses - II. Unsteady Flow, Journal of Biomechanics6 (1973) 547–559.[47] G. Porenta, D. F. Young, T. R. Rogge, A finite-element model of blood flow in arteries including taper, branches, andobstructions, Journal of Biomechanical Engineering 108 (1984) 161–167.[48] W. Zhong, X. Ji, C. Li, J. Fang, F. Liu, Determination of permeability and inertial coefficients of sintered metal porousmedia using an isothermal chamber, Applied Sciences (Switzerland) 8 (9) (2018). doi:10.3390/app8091670 .[49] J. L. Lage, B. V. Antohe, Darcy ’ s Experiments and the Deviation to Nonlinear Flow, Journal of Fluids Engineering 122(2000) 619–625.[50] M. J. Bom, S. P. Schumacher, R. S. Driessen, P. A. van Diemen, H. Everaars, R. W. de Winter, P. M. van de Ven,A. C. van Rossum, R. W. Sprengers, N. J. Verouden, A. Nap, M. P. Opolski, J. A. Leipsic, I. Danad, C. A. Taylor,P. Knaapen, Non-invasive procedural planning using computed tomography-derived fractional flow reserve, Catheterizationand Cardiovascular Interventions (2020). doi:10.1002/ccd.29210 .[51] A. R. Pries, D. Neuhaus, P. Gaehtgens, Blood viscosity in tube flow: dependence on diameter and hematocrit., TheAmerican journal of physiology 263 (6 Pt 2) (1992) H1770–H1778. doi:10.1152/ajpheart.1992.263.6.H1770 .[52] C. T. Kelley, Iterative methods for optimization, SIAM, 1999.[53] J. Su, C. C. Logan, A. D. Hughes, K. H. Parker, N. M. Dhutia, C. C. Danielsen, U. Simonsen, Impact of chronic hypoxiaon proximal pulmonary artery wave propagation and mechanical properties in rats, American Journal of Physiology-Heartand Circulatory Physiology 314 (6) (2018) H1264–H1278. doi:10.1152/ajpheart.00695.2017 .[54] K. H. Parker, An introduction to wave intensity analysis, Medical and Biological Engineering and Computing 47 (2) (2009)175–188. doi:10.1007/s11517-009-0439-y .[55] J. Mynard, D. J. Penny, J. J. Smolich, Wave intensity amplification and attenuation in non-linear flow: Implications for thecalculation of local reflection coefficients, Journal of Biomechanics 41 (16) (2008) 3314–3321. doi:10.1016/j.jbiomech.2008.10.002 .[56] G. Simonneau, A. Torbicki, P. Dorfm¨uller, N. Kim, The pathophysiology of chronic thromboembolic pulmonary hyper-tension, European Respiratory Review 26 (143) (2017) 160112.[57] K. S. Burrowes, A. R. Clark, A. Marcinkowski, M. L. Wilsher, D. G. Milne, M. H. Tawhai, Pulmonary embolism:Predicting disease severity, Philosophical Transactions of the Royal Society A: Mathematical, Physical and EngineeringSciences 369 (1954) (2011) 4255–4277. doi:10.1098/rsta.2011.0129 .[58] B. T. Tang, S. S. Pickard, F. P. Chan, P. S. Tsao, C. A. Taylor, J. A. Feinstein, Wall shear stress is decreased in thepulmonary arteries of patients with pulmonary arterial hypertension: an image-based, computational fluid dynamics study,Pulmonary circulation 2 (4) (2012) 470–476.[59] Y. Yang, Y. Xiang, Y. Li, G. Chu, H. Zou, M. Arowo, J. Chen, 3d cfd modelling and optimization of single-phase flow inrotating packed beds, The Canadian Journal of Chemical Engineering 93 (6) (2015) 1138–1148.[60] M. Vel´azquez Mart´ın, A. Albarr´an Gonz´alez-Trevilla, P. Escribano Sub´ıas, Fractional Flow Reserve-guided PulmonaryAngioplasty in Chronic Thromboembolic Pulmonary Hypertension, Revista Espa˜nola de Cardiolog´ıa (English Edition)69 (9) (2016) 863. doi:10.1016/j.rec.2015.12.029 .[61] I. M. Lang, P. Dorfm¨uller, A. V. Noordegraaf, The pathobiology of chronic thromboembolic pulmonary hypertension,Annals of the American Thoracic Society 13 (2016) S215–S221. doi:10.1513/AnnalsATS.201509-620AS .[62] K. R. Chaudhary, M. Taha, V. J. Cadete, R. S. Godoy, D. J. Stewart, Proliferative versus degenerative paradigms inpulmonary arterial hypertension: have we put the cart before the horse?, Circulation Research 120 (8) (2017) 1237–1239.[63] E. M. Lau, D. Abelson, N. Dwyer, Y. Yu, M. K. Ng, D. S. Celermajer, Assessment of ventriculo-arterial interaction inpulmonary arterial hypertension using wave intensity analysis, European Respiratory Journal 43 (6) (2014) 1801–1804. doi:10.1183/09031936.00214113 .[64] M. U. Qureshi, M. J. Colebank, D. Schreier, D. Tabima, M. Haider, N. C. Chesler, M. S. Olufsen, Characteristic impedance:Frequency or time domain approach?, Physiological Measurement 39 (1) (2018). doi:10.1088/1361-6579/aa9d60 .[65] M. Sch¨afer, N. Wilson, D. D. Ivy, R. Ing, S. Abman, L. P. Browne, G. Morgan, M. Ross, D. McLennan, A. J. Barker,B. Fonseca, M. Di Maria, K. S. Hunter, U. Truong, Noninvasive wave intensity analysis predicts functional worseningin children with pulmonary arterial hypertension, American Journal of Physiology - Heart and Circulatory Physiology315 (4) (2018) H968–H977. doi:10.1152/ajpheart.00227.2018 .[66] E. Mahmud, M. M. Madani, N. H. Kim, D. Poch, L. Ang, O. Behnamfar, M. P. Patel, W. R. Auger, Chronic Throm-boembolic Pulmonary Hypertension: Evolving Therapeutic Approaches for Operable and Inoperable Disease, Journal ofthe American College of Cardiology 71 (21) (2018) 2468–2486. doi:10.1016/j.jacc.2018.04.009 .[67] D. Boulate, F. Perros, P. Dorfmuller, J. Arthur-Ataam, J. Guihaire, L. Lamrani, B. Decante, M. Humbert, S. Eddahibi,P. Dartevelle, E. Fadel, O. Mercier, Pulmonary microvascular lesions regress in reperfused chronic thromboembolic pul-monary hypertension, Journal of Heart and Lung Transplantation 34 (3) (2015) 457–467. doi:10.1016/j.healun.2014.07.005 .[68] M. Madani, T. Ogo, G. Simonneau, The changing landscape of chronic thromboembolic pulmonary hypertension manage-ment, European Respiratory Review 26 (146) (2017).[69] A. Valent´ın, L. Cardamone, S. Baek, J. D. Humphrey, Complementary vasoactivity and matrix remodelling in arterialadaptations to altered flow and pressure, Journal of the Royal Society Interface 6 (32) (2009) 293–306. doi:10.1098/rsif.2008.0254 ..