Computer models for studying interventions used in treating refractory pulmonary hypertension
Seong Woo Han, Charles Puelz, Craig G. Rusin, Daniel J. Penny, Ryan Coleman, Charles S. Peskin
CComputer models for studying interventionsused in treating refractory pulmonary hypertension
Seong Woo Han BA , Charles Puelz PhD , Craig G. Rusin PhD ,Daniel J. Penny MD, PhD, MHA , Ryan Coleman MD , and Charles S. Peskin PhD Courant Institute of Mathematical Sciences, New York University Department of Pediatrics-Cardiology, Baylor College of MedicineJanuary 29, 2021
Abstract
This paper describes computer models for three palliations used in treating refractory pul-monary hypertension. These procedures effectively create an atrial septal defect, a ventricularseptal defect, or a patent ductus arteriosus, also known as a Potts shunt. The aim in all threecases is to generate a right-to-left shunt, allowing for either pressure or volume unloading ofthe right side of the heart in the setting of right ventricular failure, while allowing for main-tenance of cardiac output. However, these shunts are created at the expense of introducingde-oxygenated blood into the systemic circulation, thereby lowering the systemic arterial oxy-gen saturation. The models developed in this paper are based on compartmental descriptions ofhuman hemodynamics and oxygen transport. An important parameter included in our modelsis the cross-sectional area of the surgically created defect. Numerical simulations are performedto compare different interventions and various shunt sizes and to assess their impact on hemo-dynamic variables and oxygen saturations.
Pulmonary hypertension refers to a spectrum of cardiovascular and/or pulmonary diseases thatinvolve elevations in a person’s pulmonary vascular resistance (PVR). Over time, this elevatedPVR causes pathologic remodeling of the right ventricle and the pulmonary vasculature, ultimatelyresulting in right ventricular failure and death. In this paper, our focus is on interventions forrefractory pulmonary hypertension, which corresponds to disease that is unresponsive to standardmedical treatments [5, 11]. Even with aggressive pharmacotherapy, patients often will ultimatelyrequire either lung transplantation or palliative surgical or catheter-based procedures, with the goalsof either approach being aimed at providing relief to the ailing right ventricle and extending thepatient’s life for an unknown period of time.Palliative shunts used in the setting of a failing right ventricle due to elevated PVR can beclassified into two main categories: (1) pre-tricuspid shunts, meaning the shunt occurs prior toblood crossing the tricuspid valve, which is what occurs in the setting of an atrial septal defect(ASD), and (2) post-tricuspid shunts, meaning the shunt occurs after the blood passes the tricuspidvalve, which are where the ventricular septal defect (VSD) and patent ductus arteriosis (PDA)/Pottsshunt occur. This classification scheme is important, as pre-tricuspid shunts are viewed as volume-unloading shunts for the right ventricle, meaning they can unload excess volume from a failing1 a r X i v : . [ q - b i o . T O ] J a n ight ventricle, but do not directly affect the pressure the right ventricle has to pump against.Post-tricuspid shunts, on the other hand, are pressure-unloading shunts, meaning they provide alower-resistance pathway for blood to traverse, decreasing the resistance the failing right ventriclehas to pump against. Pressure-unloading shunts are often preferred, as it is the pressure load theright ventricle has to pump against that results in its failure and patient demise. However, becausethese shunts, regardless of location, allow for right-to-left shunting, they result in a decrease insystemic oxygen saturations of varying degrees and severity. Shunt effectiveness is determined notonly by location of the shunt, but by the size of the shunt as well. Small septal defects may quicklybecome restrictive over time, reducing their effectiveness at either pressure- or volume-unloadingthe right ventricle. The same is true for narrow shunts (either PDA/Potts shunt), particularly astheir length increases.The effects of various shunts on pressures, flows, and oxygen saturations are often not clear inpractice. Furthermore, shunt flows are highly sensitive to shunt size, a parameter that can be variedwithin the modeling framework developed herein. The goal of this paper is to use computationalmodels to study the impact of several possible shunts, used for treating refractory pulmonary hyper-tension, on important hemodynamic variables and oxygen saturations. The three shunts consideredhere are either within the (1) atrial septum, (2) ventricular septum, or (3) between the main pul-monary artery and aorta [10, 2]. We refer to these surgically created defects respectively by thefollowing names: (1) an atrial septal defect (ASD), (2) a ventricular septal defect (VSD), or (3) aPotts shunt. For each intervention, we develop and apply computational models to study both thebenefit in terms of reduced pulmonary artery pressure and also the detriment in terms of systemicarterial oxygen desaturation, as functions of the cross-sectional area of the shunt.The complexities associated with pulmonary hypertension have motivated the use of compu-tational models for “non-invasively” studying disease progression, diagnosis, and the performanceof possible treatments. We recall several important contributions of physics-based models for pul-monary hypertension. Delhass et al. constructed compartmental models for two possible interven-tions considered in this paper, the ASD and the Potts shunt [4]. Our paper extends their resultsto a comparison of the ASD and Potts shunt with the VSD. Gerringer et al. built and calibratedcompartmental models with animal data to study the effect of progressing pulmonary hypertensionon resistance and compliance [6]. Tewari et al. also used calibrated compartmental models to inves-tigate changes in important hemodynamic parameters after the onset of disease [26]. Vessel networkmodels, which incorporate spatial variations of blood flow and pressure, were used by Acosta et al.to derive early diagnostic indicators of disease [1]. Qureshi et al. also used vessel network models tostudy several classes of pulmonary hypertension as well as to simulate control and diseased animalmodels [18, 17]. There have also been modeling efforts to understand the impact of pulmonaryhypertension on remodeling of heart tissue. Raush et al. constructed three-dimensional solid me-chanics models of the ventricular chambers that were coupled to a mathematical description oftissue remodeling under the high pressure loads associated with hypertension [19].The rest of this paper is organized as follows. Section 2 describes the mathematical modelsfor blood flow and oxygen transport and the numerical methods used to approximate the resultingequations. That section also includes a discussion of parameter selection and the shunt modelderived from the Gorlin equation. Section 3 details results from our models for each of the threepossible surgical interventions. Our models are used to study the dependence of several importanthemodynamic variables on shunt size, and we also investigate the time-dependent details of theshunt flow waveform. Conclusions are provided in section 5.2 Circulation models and numerical methods
In this section, we present the models used in this work and the numerical methods by which themodel equations are solved. The following subsections describe a hemodynamic model, a cardiacchamber model that specifies the time-varying compliance of each heart chamber in the hemody-namic model, a shunt model based on the Gorlin equation that makes it possible to include shuntsof specified cross-sectional area in the hemodynamic model, and finally an oxygen transport modelthat calculates oxygen saturations throughout the model circulation based on the flows obtainedfrom the hemodynamic model.
The circulation is represented by a collection of compartments corresponding to compliance cham-bers. These chambers are connected by resistors that are equipped with valves [16]. We use thefollowing compliance relation for each of the N compliance chambers, numbered i = 1 , , ..., N : V i = ( V d ) i + C i P i , i = 1 , ..., N. (1)The parameter C i is the compliance of chamber i , which is assumed to be constant for arteries andveins but time-varying for the heart chambers. The variable V i is the volume of compliance chamber i , the variable P i is the pressure of that chamber, and the parameter ( V d ) i is the dead volume, thatis, the volume of the chamber when the pressure is zero. We assume the flow from chamber i tochamber j is governed by a pressure-flow relationship of the following form: Q ij = S ij R ij ( P i − P j ) = S ij G ij ( P i − P j ) , i, j = 1 , ..., N, (2)where S ij = (cid:40) , P i > P j , , P i ≤ P j . (3)Equation (2) describes the flow through a resistance that is equipped with a valve. The conductance G ij is the reciprocal of the resistance R ij . Conductance is convenient because it can be set equalto zero to represent the absence of a connection between two chambers. The variable S ij , which isdetermined by P i and P j according to equation (3), denotes the state of the valve, with S ij = 1 whenthe valve is open, and S ij = 0 when the valve is closed. Note that the words “open” and “closed”have the opposite meaning here from their use in electricity, where a closed switch is conductingand an open switch is non-conducting.Equipping every connection with a valve does not involve any loss of generality. Between anypair of chambers i and j , our framework allows for two connections of the type described above,one with a valve that allows flow only from i to j , and another with a valve that allows flow onlyfrom j to i . To model a situation in which there is no valve in a connection between chambers i and j , we need only set G ij equal to G ji . To model a leaky valve we may set G ij and G ji topositive but unequal values. Lastly, to model the situation in which there is no connection at allbetween chambers i and j , we set G ij = G ji = 0. Thus, our framework allows for a great variety ofconnection types and patterns merely by specifying the (non-symmetric) N by N matrix G .Upon differentiating equation (1) with respect to time and using the principle that the rateof change of volume is equal to inflow minus outflow, together with equation (2), one obtains the3ollowing system of ordinary differential equations for the pressures as functions of time: ddt ( C i P i ) = N (cid:88) j =1 ( S ji G ji ( P j − P i ) − S ij G ij ( P i − P j ))= N (cid:88) j =1 ( S ij G ij + S ji G ji )( P j − P i ) . (4)We assume here that all of the dead volumes are constant but allow for the possibility that someof the compliances, specifically those of the heart chambers, are functions of time. How thesecompliances are specified will be described in subsection 2.2. Equation (4) will be modified later toinclude shunt flows modeled by the Gorlin equation; see subsection 2.3.Our numerical scheme for equation (4) is the backward Euler method: C ni P ni − C n − i P n − i ∆ t = N (cid:88) j =1 ( S nij G nij + S nji G nji )( P nj − P ni ) . (5)This is a system of equations for the unknown pressures at time step n . It is a nonlinear system,because S ij is a function of P i and P j , see equation (3). The reason for using the backward Eulermethod here is its unconditional stability. If two compliance chambers are connected by a verylarge conductance (that is, by a very small resistance), their pressures will equilibrate on a veryfast time scale, and we do not want to be required to use a small enough time step to resolve thedetails of that rapid equilibration. This situation actually arises in the circulation whenever thereare two chambers with an open heart valve between them, since an open valve (at least when it isnon-stenotic) has a very high conductance.The procedure that we use to solve the nonlinear system (5) is based on the following observation:given the valve states, equation (5) is reduces to a linear system that is easy to solve for the pressures.Also, given the pressures, it is easy to evaluate the valve states from equation (3). Thus, theprocedure starts with a guess for the valve states (a good guess is the valve states that were foundon the previous time step), solves equation (5) for the pressures, resets the valve states according tothe pressures via equation (3), and so on. The process stops when the valve states (and thereforethe pressures) stop changing, and in practice this happens very quickly. On most time steps, theinitial guess, that the valves states are the same as they were at the previous time step, turns outto be correct. When the valve states stop changing, the problem stated in equation (5) is actuallysolved (except, of course, for round-off error), not merely solved to within some tolerance. This isbecause the valve states are discrete.For further discussion of the methodology described here, see [8]. As in that reference, ourmodels use six compliance chambers corresponding to the left and right ventricles and the systemicand pulmonary arteries and veins. We do not separately model the atria, but instead treat eachatrium as part of the venous system to which it is connected, and moreover we do not take intoaccount the time dependence of the atrial compliances. The ventricular compliances are, of course,time dependent in our model, but not in the same way as in [8], see subsection 2.2.In order to model severe pulmonary hypertension, the pulmonary resistance is taken to be 1.5times greater than the systemic resistance [21]. This is very different from the normal case in whichthe pulmonary resistance is approximately 10 times smaller than the systemic resistance [9, 14, 28].A possible consequence of pulmonary hypertension is right-heart hypertrophy, making the normallythin-walled right ventricle into a chamber that more closely resembles the normal left ventricle[22, 15, 3]. To model this, we use typical left-ventricular parameters for both ventricles; see the4arameters Resistance ( R ) Dead Volume ( V d ) Compliance ( C )Units mmHg/(L/min) L L/mmHgS 17.5 - -P 26.25 - -Mi 0.01 - -Ao 0.01 - -Tr 0.01 - -Pu 0.01 - -SA - 0.825 0.00175PA - 0.0382 0.00175SV - 0 1.75PV - 0 0.08Table 1: Parameters for the circulation model. Pulmonary resistance is chosen to be 1.5 timesgreater than the systemic resistance to simulate refractory pulmonary hypertension.Abbreviations: S, systemic organs; P, lungs; Mi, mitral valve; Ao, aortic valve; Tr, tricuspid valve;Pu, pulmonic valve; SA, systemic arteries; PA, pulmonary arteries; SV, systemic veins; PV,pulmonary veins.next subsection. Another consequence of pulmonary hypertension is right heart failure that leadsto increased blood volume [25]. Accordingly, we use a total blood volume of 8.68 L, instead of anormal blood volume, which would be about 5 L. This change is needed to elevate the systemicvenous pressure sufficiently to fill the hypertrophied right heart and produce a viable cardiac output.With these exceptions, we are using the same parameters as in [8]. The hemodynamic parametersused in the present model, other than the cardiac chamber parameters, are stated in Table 1. This subsection details the time-varying elastance model used for the left and right ventricles,adapted from [13]. Note that the elastance, denoted by E , is the reciprocal of the compliance C .For a cardiac chamber, the elastance, and therefore the compliance, is a given function of time.In this work, we use the same elastance function E ( t ) , with the same parameters, for the left andright ventricles. This choice is reasonable because the large right-sided pressures associated withpulmonary hypertension lead to remodeling and thickening of the right ventricular wall [7]. In severepulmonary hypertension, these changes result in a right ventricular pressure/volume characteristicssimilar to that of the left ventricle [23]. Maximum and minimum ventricular elastances are denoted E max and E min . E max is the end-systolic elastance and E min is the end-diastolic elastance. Thefunctional form of the elastance E ( t ) is given during the time interval [0 , T ] as follows: E ( t ) = k (cid:18) g ( t )1 + g ( t ) (cid:19) (cid:18)
11 + g ( t ) −
11 + g ( T ) (cid:19) + E min , (6)where g ( t ) = (cid:18) tτ (cid:19) m , g ( t ) = (cid:18) tτ (cid:19) m . (7)Here T is the period of the heartbeat. Note that equation (6) makes E (0) = E ( T ) . (We have madea slight modification of the formula used in [13] to ensure this.) Outside of the interval [0 , T ] , we5arameters Symbol Units Left Ventricle Right VentricleMinimal elastance E min mmHg/L 0.08 × × Maximal elastance E max mmHg/L 30.00 × × Contraction exponent m - 1.32 1.32Relaxation exponent m - 27.4 27.4Systolic time constant τ min 0.269T 0.269TDiastolic time constant τ min 0.452T 0.452TDead Volume V d L 0.010 0.010Period of heartbeat T min 0.0125 0.0125Table 2: Parameters for the time varying compliances in the heart model.define E ( t ) as a periodic function with period T , so that E ( t ) = E ( t + T ) for all t . The parameter k is chosen so that the maximum value of E ( t ) is E max . The formula for k to achieve this is k = E max − E min max t ∈ [0 ,T ] [( g ( t )1+ g ( t ) )( g ( t ) − g ( T ) )] (8)The maximum in the denominator of the formula for k is computed by evaluating the expressionthat needs to be maximized at a collection of equally spaced points within the interval [0 , T ] , andthen choosing the largest of the values of that expression that are found. Although this proceduredoes not yield the exact maximum value, it comes close enough for practical purposes, especiallysince the goal is to find the maximum value, rather than the time at which it occurs. Parametervalues used for the heart chambers are provided in Table 2. The constant τ controls the timescaleof contraction, τ controls the duration of systole, and m and m govern the speed of contractionand relaxation respectively. Note that τ and m are estimated from previously employed values [24],and the values for E min and E max are similar to those used by [12]. The Gorlin equation is used to calculate flows through surgically created shunts (ASD, VSD, orPotts shunt) in our model [20]. This allows us to specify the cross-sectional area of the connectionthat the surgeon creates, and to study how the shunt size affects hemodynamic variables and thetransport of oxygen.To derive the shunt model, consider two chambers, denoted by the indices 1 and 2, separatedby a wall with a hole in it that corresponds to the shunt. Let A be the cross-sectional area of thehole. We assume that the velocity of the blood as it goes through the hole is much larger than thevelocity in the two chambers, so that we may consider the fluid in each of the two chambers as ifit were at rest. Let Q denote the volume of blood flow per unit time through the hole, with thedirection from chamber 1 to chamber 2 considered positive. Then, the spatially averaged velocityof blood flow in the hole itself is given by v = Q/A . (9)Let P and P be the pressures in the two chambers, and let P be the pressure within the hole.Suppose, for example, that Q >
0. By Bernoulli’s equation in the upstream chamber up to the holeitself, one has P = P − ρv = P − ρ A Q . (10)6n the region downstream of the hole, Bernoulli’s equation does not apply because the flow thereis dominated by turbulent eddies that dissipate energy. The result is that the pressure is relativelyconstant in the downstream region, in particular that P = P . It follows that P − P = ρ A Q , Q > . (11)By the same reasoning, for flow in the other direction P − P = ρ A Q , Q < . (12)Equations (11) and (12) can be combined as follows: P − P = ρ A | Q | Q, (13)This shows that the hydraulic resistance of the hole is given by R shunt = ρ A | Q | . (14)Note that the above formula for R shunt is independent of viscosity. In reality, there is a very smallviscous resistance as well, so we modify the above formula to read R shunt = R visc + ρ A | Q | . (15)In most situations R visc is negligible, but we should include it to prevent R shunt from being zero,which would otherwise happen in principle every time that Q changes sign. Since R visc is includedonly for this reason, we choose the very small value R visc = 0.1 mmHg/(L/min). Evaluating theconductance of the hole from equation (15), one obtains G shunt = 1 R visc + ρ A | Q | (16)In making use of equation (16), one must be careful about units. In this paper, we use what maybe called physiological units, in which volume is measured in liters, pressure in mmHg, and time inminutes. The constants A and ρ need to be expressed in these units. The use of liters for volumeimplies that our unit of length is the decimeter (dm), which is equal to 10 cm. Thus, A needs tobe expressed internally in terms of dm , although in the presentation of our results, we use cm since these units have more meaning to the reader. To express density in physiological units, oneneeds the units of mass. The units of force are mmHg · dm , and the units of acceleration aredm/min , so the units of mass are mmHg · dm · min . Dividing this by dm , one obtains the unitsof density as mmHg · (min/dm) . After taking units carefully into account in this way, the densityin physiological units is ρ = 0 . · mmHg · min dm . (17)Another complication in the use of equation (16) is that the shunt conductance G shunt is flow-dependent. A simple idea here would be to use the shunt flow on the previous time step to setthe shunt conductance for the present time step, but instead of this, we use a fixed-point iteration,with the shunt flow on the previous time step as the initial guess. At each step of the fixed-point7teration, equation (16) is used to set the shunt conductance based on the latest guess for the shuntflow. Then, the shunt conductance is inserted into the appropriate two places in the conductancematrix G (one entry for each flow direction, since there is no valve involved in the shunt). Finally,all of the pressures and flows for the circulation are computed, including the shunt flow. The benefitof doing the fixed-point iteration can be seen in Figure 1 since it removes the numerical oscillationsseen in the blood flow waveform. In practice, 10 fixed-point iterations are used for each time step,and this achieves good enough agreement between the flow that is used to set the shunt conductanceand the flow that is calculated on the basis of that shunt conductance.Figure 1: A comparison of shunt flow results with the ASD model, obtained with and withoutfixed-point iterations. The shunt area is 1 cm . The left panel shows flow computed without thefixed-point iteration. The right panel shows flow computed with the fixed-point iteration. An important consequence of the surgical interventions considered in this paper is the mixing ofoxygenated and doxygenated blood. Our approach to the modeling of oxygen transport follows Tuand Peskin [27]. Time-varying oxygen concentrations for each compliance chamber are describedby the following system of differential equations: ddt ([ O ] i V i ) = N (cid:88) j =1 j (cid:54) = i ([ O ] j Q ji − [ O ] i Q ij + M ji ) . (18)The variable [ O ] i is the oxygen concentration in compliance chamber i , the variable Q ji is the bloodflow from j to i , and the parameter M ji is the rate at which oxygen is added to the stream of bloodthat is flowing from chamber j to chamber i . Note that M ji is positive if oxygen is being addedto the blood stream, and negative if oxygen is being removed. The correctness of equation (18)relies on the fact that all of the flows that appear in it are positive or zero. This is a benefit of ourformulation in which every connection is equipped with a valve, as described above in subsection2.1. These equations describe conservation of oxygen during transport between chambers, metabolicconsumption of oxygen within systemic organs, and replenishment of oxygen within the lungs. Aftercomputing the flows at time step n , those flow values are used to numerically update the oxygenconcentrations from time step n − to n as follows: [ O ] ni V ni − [ O ] n − i V n − i ∆ t = N (cid:88) j =1 j (cid:54) = i ([ O ] n − j Q nji − [ O ] n − i Q nij + M ji ) . (19)8ote that this is the forward Euler method insofar as the oxygen concentrations are concerned,although it differs from the forward Euler method by using the flows at time step n . The mannerin which M is determined for use in these equations is described below.We use the millimole (mmol) as the unit for the amount of oxygen. It then follows from ourother choices of units that the units of oxygen concentration are mmol/L and the units of the rate ofoxygen consumption by the body are mmol/min. A standard concentration of hemoglobin in bloodis 2.5 mmol/L, and since each hemoglobin molecule can carry four oxygen molecules, the oxygenconcentration when hemoglobin is fully saturated is 10 mmol/L.There are only two places in our model where the variable M that appears in equation (18)is nonzero. One of these is in the connection from the pulmonary arteries (pa) to the pulmonaryveins (pv). We assume that M pa,pv is such that the stream of blood flowing from the pulmonaryarteries to the pulmonary veins becomes fully saturated with oxygen during its passage through thepulmonary capillaries. This gives the equation M pa,pv = (10 mmol/L − [ O ] pa ) Q pa,pv . (20)Equation (20) is used to set M pa,pv at every time step. Note that this is not the same as setting [ O ] pv = 10 mmol/L. The reason for this is that there may be other streams of blood entering thepulmonary venous compartment besides the one coming from the pulmonary arteries. In particular,since we regard the left atrium as being part of the pulmonary venous compartment, this will bethe case when we are simulating a surgically created atrial septal defect.The other nonzero value of M in our model is M sa,sv , which is negative, since it representsoxygen consumption by the tissues. This oxygen is extracted from the stream of blood that flowsfrom the systemic arteries (sa) to the systemic veins (sv). In the simulations reported here, we keep M sa,sv constant, and we calculate its value by noting that a normal cardiac output is 5.6 L/min anda normal amount of oxygen extraction by the systemic tissue is 30%. These values result in thefollowing: − M sa,sv = 0 . · (10 mmol/L ) · (5 . L/min ) = 16 . mmol/min (21)The initial value for the oxygen concentration is set to 10 mmol/L in all compartments. Simulationsare run until all hemodynamic variables and oxygen concentrations have reached a periodic steadystate. The values reported in the following section are averages over the last five of the computedcardiac cycles. In this section, we examine the simulated effects for each of the three interventions on several hemo-dynamic variables. First, we investigate changes in pressure and oxygen saturation in the systemicand pulmonary arteries as the shunt size is varied. Second, we examine shunt flow direction to de-termine whether the shunt is indeed right-to-left, as anticipated, or is perhaps bidirectional. In thesesimulations, 100 time steps are used for each cardiac cycle. The heart rate is 80 beats/minute, andeach computer experiment is run for 12.5 minutes of simulated time, thus covering 1,000 simulatedheartbeats. This duration is sufficient for the simulation to achieve a periodic steady state for allvariables in all cases. When we report a single value for any quantity as the result of a simulation,it is the average of that quantity over the last five cardiac cycles of that simulation.9 .1 Pressures and oxygen saturations
First, we consider the pressures and oxygen saturations for each intervention. Figure 2 shows thesystemic and pulmonary arterial blood pressures (mean values) as functions of the shunt area.Results for the atrial septal defect (ASD) are in the left panel, results for the ventricular septaldefect (VSD) are in the middle panel, and results for the Potts shunt are in the right panel. TheASD results show that this intervention is not successful in lowering the pulmonary arterial pressure,which is perhaps consistent with the fact that an ASD is a volume-unloading shunt. The insets inthe figure confirm that there is actually a very small effect, but it is certainly not one that wouldbe therapeutic. The VSD intervention lowers the mean pulmonary artery pressure from about 105mmHg to about 85 mmHg, which could be beneficial. This result is achieved by a shunt withcross-sectional area equal to 0.6 cm . It is interesting to note that beyond this value for the area,the pulmonary arterial pressure slightly increases as the shunt size increases. The Potts shunt mostsubstantially lowers the mean pulmonary arterial pressure, from about 105 mmHg to about 75mmHg. Unlike in the VSD case, the pulmonary arterial pressure with the Potts shunt decreasesmonotonically with increasing shunt size, but most of the benefit has already occurred with a shuntsize of 0.5 cm . Our model suggests that there is little benefit in using a larger Potts shunt sizethan this value. Figure 3 shows the pressures in pulmonary artery and systemic artery for the threeinterventions, all on the same plot, as functions of the shunt area.Figure 2: A comparison of pressures in the pulmonary artery (blue) and systemic artery (red) asthe shunt area is varied in the ASD, VSD, and Potts shunt. The inset figure shows a zoomed inportion of the results.We further investigate oxygen transport for each intervention in Figure 4. This figure depictssystemic flow, oxygen saturation, and the rate of oxygen delivery for the three interventions. Sys-temic flow is relevant here because it is used in the computation of oxygen delivery to the systemictissues. Oxygen saturation is the oxygen concentration (in mmol/L) expressed as a percentage of10 mmol/L, which is the maximum possible oxygen concentration in our model. The rate at whichoxygen is delivered to the systemic tissues is calculated by multiplying the systemic flow by thesystemic arterial oxygen concentration.All three interventions increase systemic flow. The increase is substantial, however, only in thecase of VSD, where the increase in systemic flow has a substantial effect on oxygen delivery tothe systemic tissues, as discussed below. All three interventions decrease systemic arterial oxygensaturation. This is inevitable, since the interventions by design are allowing deoxygenated bloodto bypass the lungs. The effect is smallest in the case of the ASD, but since the ASD interventionhad essentially no benefit, the fact that it also does the least harm is not really of interest. TheVSD and Potts shunt produce similar decreases in systemic arterial oxygen saturation, but thesetwo interventions look quite different from the point of view of oxygen delivery to the systemic10igure 3: A comparison of pressures in pulmonary artery and systemic artery for the threeinterventions: the left panel shows the pulmonary artery pressures for the three interventions asthe shunt area is varied. The right panel shows the systemic artery pressures for the threeinterventions as the shunt area is varied. Note different pressure scales in the two panels.tissues. The increase in systemic flow in the VSD case seems to compensate nicely for the drop insystemic arterial oxygen saturation. Even at the shunt size of 0.4 cm , where the oxygen deliveryis the smallest in the VSD case, it is only about 4% smaller than in the pre-intervention state(corresponding to a shunt area of zero). Recall, however, that the optimal reduction in pulmonaryarterial blood pressure occurs in the case of VSD at a shunt size of 0.6 cm , and at this shunt size,the reduction of oxygen delivery is even smaller – only about 2.5%.Figure 4: A comparison of systemic flow, oxygen saturation, and oxygen delivery rate for the threeinterventions: the left panel shows the systemic flow, the middle panel shows the oxygensaturation, and the right panel shows the rate of oxygen delivery to the systemic tissues. Figures 5, 6, and 7 depict shunt flow waveforms for each of the three interventions: the ASD, VSD,and Potts shunt respectively. Two shunt sizes, 0.1 cm and 1 cm , are considered. The last 5cardiac cycles of a 1000-cycle computer simulation are shown in each case. In all cases, our signconvention is that right-to-left flow is considered positive. Note that the ASD shunt flow is stronglybidirectional, whereas the VSD and Potts shunt flows are exclusively right-to-left. The VSD flowis essentially zero, however, during part of each cardiac cycle, whereas the Potts shunt flow nevercomes close to zero, but instead is sustained right-to-left flow throughout the cardiac cycle.Although bidirectional flow does not play an important role in this study, since the ASD inter-vention does not appear to be useful, it is an important feature of our methodology that it evaluates11he shunt flow as a function of time, and not merely the mean flow. This is because bidirectionalflow can exchange oxygen between two compartments even when there is no mean shunt flow at all,and this can have a substantial impact on oxygen transport. Indeed, in the congenital heart diseasecalled transposition of the great arteries, the pulmonary and systemic circulations form parallelloops, and there cannot be any mean flow from one to the other. Survival of the patient after birthis then completely dependent on the existence of a bidirectional shunt [27].Figure 5: Shunt flow waveforms in the ASD. The left panel shows flow when the shunt area is 0.1cm , and the right panel flow when the shunt area is 1 cm .Figure 6: Shunt flow waveforms in the VSD. The left panel shows flow when the shunt area is 0.1cm , and the right panel flow when the shunt area is 1 cm .Figure 7: Shunt flow waveforms in the Potts shunt. The left panel shows flow when the shunt areais 0.1 cm , and the right panel shows flow when the shunt area is 1 cm .12 Limitations
The specific results reported above may certainly depend upon the specific parameters chosen. Inorder to apply the methodology of this paper with confidence to any particular patient, it will benecessary to identify the relevant cardiovascular parameters of that patient. Making the modelpatient-specific is also a way to test the validity of the model, since the pre-operative state of apatient can be used to identify patient-specific parameters, and then the model can be used topredict what the immediately post-operative state of the patient will be. Comparison with theactual post-operative state will then be a strong test of the model. Future work will thereforebe directed toward the development of a methodology for identifying the model parameters thatcorrespond best to the state of a particular patient, so that the model can be made useful in clinicalpractice.
In this paper, we have presented a methodology that can be used to study surgical interventionsthat are designed to alleviate the detrimental effects of refractory pulmonary hypertension. Wehave illustrated the use of this methodology by comparing three such interventions, all of which aredesigned to allow some blood flow to bypass the lungs: an atrial septal defect, a ventricular septaldefect, and a Potts shunt. For each intervention, we have simulated a range of defect sizes from 0to 1 cm . Our results are that the ASD is ineffective at lowering blood pressure in the pulmonaryartery, but that the VSD and Potts shunt are both effective, with a greater effect being produced bythe Potts shunt. These results are consistent with the fact that an ASD is volume-unloading whilea VSD or Potts shunt is pressure-unloading. Both the VSD and Potts shunt lower the systemicarterial oxygen saturation in our study, but this is partially compensated by an increase in systemicflow, so that oxygen delivery to the systemic tissues is lowered to a lesser degree than the systemicarterial oxygen saturation. The increase in systemic flow is substantially greater for the VSD thanfor the Potts shunt, with the result that oxygen delivery is only slightly reduced in the VSD case,but more substantially reduced in the case of the Potts shunt. Thus, when judged by reduction ofpulmonary arterial pressure alone, the Potts shunt appears to perform better. When consideringoxygen delivery, the VSD appears to be the best choice, but it is important to note that this shuntresults in blood mixing in the ventricles, and consequently, desaturated blood is delivered to thebrain, coronaries arteries, and lower body. In contrast, the Potts shunt delivers desaturated bloodonly to the lower body, while fully saturated blood is delivered to the brain and heart. This isby virtue of the fact that mixing occurs downstream from the carotid arteries. All of the aboveeffects are quantified in our study as functions of the size of the defect in each case, and this kindof information could be useful to a clinician who needs to decide how large a defect to create. Charles Puelz was supported in part by the Research Training Group in Modeling and Simulationfunded by the National Science Foundation via grant RTG/DMS-1646339.
References [1] S Acosta, C Puelz, B Rivière, DJ Penny, KM Brady, and CG Rusin. Cardiovascular mechanicsin the early stages of pulmonary hypertension: a computational study.
Biomechanics and odeling in Mechanobiology , 16(6):2093–2112, 2017.[2] AE Baruteau, E Belli, Y Boudjemline, D Laux, M Lévy, G Simonneau, A Carotti, M Hum-bert, and D Bonnet. Palliative Potts shunt for the treatment of children with drug-refractorypulmonary arterial hypertension: updated data from the first 24 patients. European Journalof Cardio-Thoracic Surgery , 47(3):e105–e110, 2015.[3] D Chemla, V Castelain, P Herve, Y Lecarpentier, and S Brimioulle. Haemodynamic evaluationof pulmonary hypertension.
European Respiratory Journal , 20(5):1314–1331, 2002.[4] T Delhaas, Y Koeken, H Latus, C Apitz, and D Schranz. Potts shunt to be preferred aboveatrial septostomy in pediatric pulmonary arterial hypertension patients: a modeling study.
Frontiers in Physiology , 9:1252, 2018.[5] CD Etz, HA Welp, TDT Tjan, A Hoffmeier, E Weigang, HH Scheld, and C Schmid. Medicallyrefractory pulmonary hypertension: treatment with nonpulsatile left ventricular assist devices.
The Annals of Thoracic Surgery , 83(5):1697–1705, 2007.[6] JW Gerringer, JC Wagner, D Vélez-Rendón, and D Valdez-Jasso. Lumped-parameter mod-els of the pulmonary vasculature during the progression of pulmonary arterial hypertension.
Physiological Reports , 6(3):e13586, 2018.[7] S Giusca, E Popa, MS Amzulescu, I Ghiorghiu, IM Coman, BA Popescu, M Delcroix, JU Voigt,C Ginghina, and R Jurcut. Is right ventricular remodeling in pulmonary hypertension depen-dent on etiology? an echocardiographic study.
Echocardiography , 33(4):546–554, 2016.[8] FC Hoppensteadt and CS Peskin.
Modeling and simulation in medicine and the life sciences .Springer Science & Business Media, 2012.[9] G Kovacs, A Olschewski, A Berghold, and H Olschewski. Pulmonary vascular resistances duringexercise in normal subjects: a systematic review.
European Respiratory Journal , 39(2):319–328,2012.[10] JA Leopold. Catheter-based therapies for patients with medication-refractory pulmonary ar-terial hypertension.
Circulation: Cardiovascular Interventions , 8(11):e003332, 2015.[11] M Levy, DS Celermajer, E Bourges-Petit, MJ Del Cerro, F Bajolle, and D Bonnet. Add-ontherapy with subcutaneous treprostinil for refractory pediatric pulmonary hypertension.
TheJournal of Pediatrics , 158(4):584–588, 2011.[12] F Liang, S Takagi, R Himeno, and H Liu. Multi-scale modeling of the human cardiovascu-lar system with applications to aortic valvular and arterial stenoses.
Medical and BiologicalEngineering and Computing , 47(7):743–755, 2009.[13] JP Mynard, MR Davidson, DJ Penny, and JJ Smolich. A simple, versatile valve model for usein lumped parameter and one-dimensional cardiovascular models.
International Journal forNumerical Methods in Biomedical Engineering , 28(6-7):626–641, 2012.[14] N Naderi. Chapter 11 - Hemodynamic Study. In Majid Maleki, Azin Alizadehasl, and MajidHaghjoo, editors,
Practical Cardiology , pages 183 – 191. Elsevier, 2018.1415] AV Noordegraaf, KM Chin, F Haddad, PM Hassoun, AR Hemnes, SR Hopkins, SM Kawut,D Langleben, J Lumens, and R Naeije. Pathophysiology of the right ventricle and of thepulmonary circulation in pulmonary hypertension: an update.
European Respiratory Journal ,53(1), 2019.[16] CS Peskin and C Tu. Hemodynamics in congenital heart disease.
Computers in Biology andMedicine , 16(5):331–359, 1986.[17] MU Qureshi, MJ Colebank, LM Paun, LE Fix, N Chesler, MA Haider, NA Hill, D Husmeier,and MS Olufsen. Hemodynamic assessment of pulmonary hypertension in mice: a model-basedanalysis of the disease mechanism.
Biomechanics and Modeling in Mechanobiology , 18(1):219–243, 2019.[18] MU Qureshi, GDA Vaughan, C Sainsbury, M Johnson, CS Peskin, MS Olufsen, and NA Hill.Numerical simulation of blood flow and pressure drop in the pulmonary arterial and venouscirculation.
Biomechanics and Modeling in Mechanobiology , 13(5):1137–1154, 2014.[19] MK Rausch, A Dam, S Göktepe, OJ Abilez, and E Kuhl. Computational modeling of growth:systemic and pulmonary hypertension in the heart.
Biomechanics and Modeling in Mechanobi-ology , 10(6):799–811, 2011.[20] T Reynolds. The determination of aortic valve area by the gorlin formula: what the cardiacsonographer should know.
Journal of the American Society of Echocardiography , 3(4):331–335,1990.[21] AK Roy, SP Gaine, and KP Walsh. Percutaneous atrial septostomy with modified butterflystent and intracardiac echocardiographic guidance in a patient with syncope and refractorypulmonary arterial hypertension.
Heart, Lung and Circulation , 22(8):668–671, 2013.[22] JJ Ryan and SL Archer. The right ventricle in pulmonary arterial hypertension: disordersof metabolism, angiogenesis and adrenergic signaling in right ventricular failure.
CirculationResearch , 115(1):176–188, 2014.[23] K Ryo, A Goda, T Onishi, A Delgado-Montero, B Tayal, HC Champion, MA Simon, MA Math-ier, MT Gladwin, and J Gorcsan III. Characterization of right ventricular remodeling in pul-monary hypertension associated with patient outcomes by 3-dimensional wall motion trackingechocardiography.
Circulation: Cardiovascular Imaging , 8(6):e003176, 2015.[24] N Stergiopulos, JJ Meister, and N Westerhof. Determinants of stroke volume and systolic anddiastolic aortic pressure.
American Journal of Physiology-Heart and Circulatory Physiology ,270(6):H2050–H2059, 1996.[25] S Stickel, W Gin-Sing, M Wagenaar, and JSR Gibbs. The practical management of fluidretention in adults with right heart failure due to pulmonary arterial hypertension.
EuropeanHeart Journal Supplements , 21(Supplement_K):K46–K53, 2019.[26] S Tewari, S Bugenhagen, Z Wang, D Schreier, B Carlson, N Chesler, and D Beard. Analysisof cardiovascular dynamics in pulmonary hypertensive c57bl6/j mice.
Frontiers in Physiology ,4:355, 2013.[27] C Tu and CS Peskin. Hemodynamics in transposition of the great arteries with comparison toventricular septal defect.
Computers in Biology and Medicine , 19(2):95–128, 1989.1528] J Widrich and M Shetty. Physiology, pulmonary vascular resistance. In