Multiscale modelling of Potts shunt as a potential palliative treatment for suprasystemic idiopathic pulmonary artery hypertension: a paediatric case study.

Potts shunt (PS) was suggested as palliation for patients with suprasystemic pulmonary arterial hypertension (PAH) and right ventricular (RV) failure. PS, however, can result in poorly understood mortality. Here, a patient-specific geometrical multiscale model of PAH physiology and PS is developed for a paediatric PAH patient with stent-based PS. In the model, 7.6mm-diameter PS produces near-equalisation of the aortic and PA pressures and [Formula: see text] (oxygenated vs deoxygenated blood flow) ratio of 0.72 associated with a 16% decrease of left ventricular (LV) output and 18% increase of RV output. The flow from LV to aortic arch branches increases by 16%, while LV contribution to the lower body flow decreases by 29%. Total flow in the descending aorta (DAo) increases by 18% due to RV contribution through the PS with flow into the distal PA branches decreasing. PS induces 18% increase of RV work due to its larger stroke volume pumped against lower afterload. Nonetheless, larger RV work does not lead to increased RV end-diastolic volume. Three-dimensional flow assessment demonstrates the PS jet impinging with a high velocity and wall shear stress on the opposite DAo wall with the most of the shunt flow being diverted to the DAo. Increasing the PS diameter from 5mm up to 10mm results in a nearly linear increase in post-operative shunt flow and a nearly linear decrease in shunt pressure-drop. In conclusion, this model reasonably represents patient-specific haemodynamics pre- and post-creation of the PS, providing insights into physiology of this complex condition, and presents a predictive tool that could be useful for clinical decision-making regarding suitability for PS in PAH patients with drug-resistant suprasystemic PAH.

Introduction

Pulmonary arterial hypertension (PAH) is a rare disease in paediatric patients that is associated with significant morbidity and mortality. In the majority of paediatric patients, PAH is idiopathic or associated with congenital heart disease (Ivy et al. 2013). Idiopathic pulmonary artery hypertension (iPAH) is a rare, chronic disorder of the pulmonary vasculature characterised by cellular changes in the vascular walls, which cause progressive constriction, obstruction or obliteration of the small pulmonary vessels in the lungs, thereby increasing the resistance to pulmonary blood flow. In response, the right ventricle (RV) progressively adapts to pump blood through the high-pressure pulmonary vasculature by developing hypertrophy and dilatation, eventually resulting in cardiac failure and death.

iPAH in children, probably due to the highly adaptive RV myocardium, is characterised by the ability to sustain very high pulmonary arterial pressures above systemic levels for a long time (Barst et al. 2011). Furthermore, due to often non-specific symptoms, iPAH in children is typically diagnosed relatively late (Hoeper et al. 2013; Ivy et al. 2013). Although the implementation of the so-called triple therapy strategy in paediatric iPAH has lead to significant improvement in prognosis of these patients (Shu et al. 2021), there is a subgroup showing refractoriness of vascular resistance to medical treatment. In the settings of donor organ shortage, for children with progressive suprasystemic PAH presenting with an inadequate response to drug therapy and progressive RV failure, there is a need for alternative approaches to avoid further RV deterioration.

Based on the superior long-term survival and lower RV failure incidence in PAH patients with Eisenmenger syndrome due to untreated congenital heart disease compared to iPAH (Baruteau et al. 2014; Grady and Eghtesady 2016; Sizarov et al. 2016), surgical creation of a post-ventricular right-to-left shunt was suggested as a promising palliative treatment for severely ill children with drug-refractory suprasystemic iPAH (Blanc et al. 2004; Baruteau et al. 2014; Grady and Eghtesady 2016). The reversed flow version of the so-called Potts shunt (PS), an anastomosis between the left pulmonary artery (LPA) and the descending aorta (DAo) (Potts et al. 1946), in cases where pulmonary arterial pressures are above systemic ones, allows blood to flow from the LPA to the DAo, thereby nearly equalising pressures and partially decompressing the RV. Furthermore, the blood flow through the PS allows to limit desaturation and the risk of paradoxical embolism to the lower body, thus sparing the brain (Hansmann 2017; Boudjemline et al. 2013). The vicinity of the DAo to LPA in humans creates an attractive possibility to place a covered stent between the lumens of these two vessels and create an anastomosis percutaneously, thus allowing substantial reduction of treatment invasiveness in these patients (Sizarov et al. 2016). Recently, feasibility of percutaneous stent-based PS creation in adults and children with iPAH has been reported (Boudjemline et al. 2017; Esch et al. 2013).

Whilst small surgical and interventional series of PS procedures have demonstrated sustained improvement in functional capacities and prolonged survival in the majority of patients, the clinical response to PS creation is mixed with substantial mortality (Baruteau et al. 2014; Grady and Eghtesady 2016; Boudjemline et al. 2017), the mechanisms of which are poorly understood. Although cardiovascular magnetic resonance imaging can provide useful insights (Schäfer et al. 2019), the full-spectrum of haemodynamic changes due to the PS creation is not accessible. In contrast, computational models that can accurately and efficiently predict the patient-specific post-procedural haemodynamic changes based upon pre-operative characteristics can prove to be a useful tool in understanding this complex condition. Furthermore, such models could present a substantial aid in clinical decision-making to assess the suitability of PS creation in individual patients with particular haemodynamic characteristics. With the exception of a recent non patient-specific study using the well-known CircAdapt lumped parameter model (LPM) (Delhaas et al. 2018)—which simulated haemodynamic changes due to increased pulmonary arterial pressures, and subsequent creation of PS with varying shunt diameters in a reference patient, till now there is no comprehensive computational investigation of this subject. The goals of our study are: (i) to present a geometric multiscale model that is capable of representing relevant 3D anatomy and associated flow features in the case of suprasystemic iPAH palliated with PS; (ii) to demonstrate that this model can be tuned to a patient-specific case while reproducing pre- and post-operative clinically measured haemodynamic parameters with acceptable errors; (iii) to shed light on pre- and post-operative local and global haemodynamics dictated by model physics; (iv) to assess how PS diameter and length affect global and local haemodynamics; and (v) to build a ‘standalone LPM’ model in the process and assess the extent to which it can reproduce the GMM output.

The development and output of two models—a geometric multiscale model (GMM) model and a standalone LPM—is reported. The developed models are then applied to a paediatric iPAH case with percutaneously created PS to investigate the accuracy of representation of the clinically measured pre- and post-operative haemodynamics. The fully-tuned GMM is validated against patient-specific data and provides insights into blood flow features through and around the PS. The results from the computationally expensive GMM are compared to the less demanding LPM exploring the validity of the latter in providing solutions for patient-specific parameter estimation and assessment of global haemodynamics in such a complex condition. Furthermore, for the first time, 3D flow features are comprehensively assessed along with changes in global haemodynamics in response to stent-based PS creation of varying diameters and lengths. This validated computational model is adaptable and presents a predictive tool that could be useful for clinical decision-making regarding suitability for PS in PAH patients with drug-resistant suprasystemic PAH.

Materials and methods

Patient characteristics and measurements

To go beyond the generalized assumptions of haemodynamic PAH parameters, data from clinical measurements were obtained for a 13-year-old patient with morphologically normal heart and suprasystemic iPAH complicated by the RV failure despite the triple vasodilator therapy, who received the PS using a covered stent implanted through percutaneous approach at the Necker University Hospital for Sick Children, Paris, France, as published previously (Boudjemline et al. 2017). In this particular patient, despite the substantial technical difficulties during the procedure, there was no acute circulatory deterioration after PS creation with near equalization of the systemic and pulmonary artery (PA) pressures. Pre-operative measurements included an electrocardiogram, pressure tracings from heart catheterisation, Doppler flow velocity tracings from echocardiography, and CT-angiography imaging. Body surface area (BSA) of the patient was .

Pre-operatively measured pressures were available for the right atrium (RA), RV, DAo, main PA (mPA), and the pulmonary capillary wedge pressure, which was used as an indirect estimate of the left atrial (LA) pressure. Post-operative haemodynamic measurements were limited to pressure tracings for the DAo and mPA. Depending upon the completeness of the measurements, the pressure traces are averaged over 3-5 cardiac cycles. The measured pre-operative pressures used in this study are summarised in Table 1.

Table 1

Measured pressures in the pre- and post-operative states

	Pre-operative				Post-operative
	P_DAo	P_mPA	P_LA	P_RA	P_DAo	P_mPA
Systolic [mmHg]	94	112	–	–	97	102
Diastolic [mmHg]	53	67	–	–	51	54
Mean [mmHg]	69	85	4	6	71	76

The average pre-operative cardiac output (CO) for the patient was 3.4 L/min, based upon three measurements of thermodilution between the RA and RV under sedation and further three measurements within the same catheterisation but under nitric oxide inhalation and 100% oxygen. Durations of pre-operative atrial and ventricular systole were obtained from Doppler echocardiography and electrocardiography, from which a cardiac cycle time-period of 0.9s was extracted. Pre-operative end-diastolic ventricular volumes (EDV), corresponding end-systolic atrial volumes (ESV), and the myocardial wall volumes for each cardiac cavity, , were determined by 3D reconstruction of the contrast medium-stained cavities and their surrounding walls as visible on the cardiac CT imaging. The measured pre-operative CO and cardiac cycle time-period provided stroke volumes (SV) for the LV and the RV, which yielded corresponding ejection fractions (EF) when combined with the CT-reconstructed EDVs. These measurements are summarised in Table 2.

Table 2

Pre-operative measurements for heart chamber output, end-diastolic volumes, ejection fractions, myocardial volumes, and systole durations

Pre-operative	LA	RA	LV	RV
Cardiac output(CO) [L/min]	–	–	–	3.4
Activation duration (t_max) [s]	0.16	0.14	0.4	0.38
End-diastolic volume (EDV) [ml]	–	–	66.9	91.8
End-systolic volume (ESV) [ml]	22.7	51.9	–	–
Ejection fraction (EF) [-]	–	–	0.76	0.55
Myocardial volume (V_w) [ml]	6.4	10.3	51.7	69.0

Models

Two closed-loop computational models have been developed: a stand-alone LPM; and a GMM model consisting of a reconstructed patient-specific three-dimensional (3D) flow domain, coupled to an LPM. The fast-to-compute LPM provides quick solutions (in the order of seconds) for patient-specific parameter estimation and global haemodynamics assessment, and the GMM, while the computationally expensive (run-time in the order of days), provides detailed local flow information in the vessels and the PS.

Geometric multiscale model (GMM)

A schematic of the GMM is shown in Fig. 1. The 3D domain is shown in blue and comprises sections of the large systemic and pulmonary vasculature relevant to the PAH physiology. The 3D anatomical geometry is prepared for computational modelling using Amira software (Thermo Fisher Scientific, Darmstadt, Germany) and Mimics Innovation Suite (Mimics Innovation Suite kernel description 2019). The systemic vessels represented in the 3D model include the ascending aorta (AAo), aortic arch with its branches—the innominate artery (IN), the left carotid artery (LC), and left subclavian artery (LS)—and the thoracic portion of the descending aorta (DAo). The pulmonary vasculature is represented by the main pulmonary artery (mPA), the right pulmonary artery (RPA), and the left pulmonary artery (LPA) with its two branches (denoted as LPA-1 and LPA-2). To reduce the complexity of simulations and make the results applicable also to surgical PS cases, post-operative models are generated by modelling the PS as a cylinder of constant diameter placed in perpendicular fashion between the LPA and DAo adjacent walls (Fig. 2a).

Fig. 1

GMM for the iPAH and PS physiology: the 3D geometry (see Fig. 2a) is represented in blue and the remaining circulation is represented as a lumped parameter model (LPM) with en electrical analogy to blood-flow. represents volumetric flow-rate; represents pressure, represents volume

Fig. 2

3D reconstruction of the patient anatomy and CT angiogram. Labels are defined in Fig. 1

Initially, a stent-based PS with 7.6 mm diameter PS is generated, corresponding to the clinical case under consideration (see Sect. 2.1). To determine the effect of the stent diameter on flow patterns and global haemodynamics, additional post-operative geometries are created using cylinders with diameters of 5, 6, 7, 8, 9, and 10 mm. For simplicity, the simulations of the shunts with varying diameters were performed with the length of modelled stents constrained by the distance between the LPA and DAo, thus, not protruding into the vascular lumens, a situation corresponding to either surgically created side-to-side anastomosis or a spool-shaped covered stent with complete apposition of its flaring ends to the vessel walls (Chigogidze et al. 2006). The currently available stents used to create the PS in a clinical setting, however, invariably result in their ends protrusion into the vessel lumens, as was the also the case in the patient considered in this study (Fig. 2b). Therefore, additional simulations of the PS creation were run with maintaining diameter at 7 mm, while varying the lengths of the modelled stent: four further post-operative geometries were created to investigate the effect of stent protrusion into the LPA and DAo, with the stent lengths of 10, 15, 20, and 25 mm, while keeping the shunt flow orientation perpendicular to the walls of anastomosed vessels. All models were meshed using TetGen (Si 2015; Sahni et al. 2008) available within the SimVascular Suite (Updegrove et al. 2017; Lan et al. 2018), with additional local and regional mesh refinement ensuring preservation of key geometrical features.

The 3D flow domain is coupled to an LPM, which represents the global circulatory system with a hydraulic-electric analogy (Shi et al. 2011). The LPM partitions the circulatory system into multiple blocks or compartments, each of which is described by a series of elements representing the linear viscous losses (R), nonlinear viscous losses (K), blood inertia (L), and/or vessel compliance (C), depending upon the salient characteristics of the local vasculature. The resulting instantaneous pressure-flow relationships are shown in Table 3. The four heart chambers are each described by a single fibre model (see Section 2.2.3) and the atrioventricular and semilunar valves are modelled dynamically (see Sect. 2.2.4), with the contraction of the heart and valve function represented by time-varying capacitors (compliances) and diodes, respectively.

Table 3

Pressure (P) and flow rate (Q) relationships applying the hydraulic-electrical analogy. The last two rows show the heart chambers and heart valves, respectively, where V represents volume of the chamber

Component	Pressure-flow relationship
	P_1 – P_2 = RQ
	P_1 – P_2 = KQ|Q|
	P_1 – P_2 = LdQdt
	Q_1 – Q_2 = CdPdt
	Q_1 – Q_2 = dVdt
	P_1 – P_2 = BQ|Q|+MdQdt

In the GMM, Fig. 1, downstream of each 3D outlet (see for example the highlighted path in orange downstream of the DAo), the arterial structure is represented by a three-element RCR Windkessel model, while the subsequent venous structure is represented by an additional two-element RC Windkessel model. Specifically, the three-element RCR Windkessel model comprises a proximal resistance of the large vessels, , a distal resistance of the microvasculature, , and the total arterial compliance, C, representing the elastance of the large vessels. The two-element RC Windkessel model comprises a venous resistance, , and the total venous compliance, .

Deriving the stand-alone LPM from the GMM

For creating the stand-alone LPM, the 3D flow-domain of the GMM is substituted by lumped parameter components while minimising the error between the GMM solution and the stand-alone LPM (Pant et al. 2013, 2014). Based on the GMM solution, it was found that the while the pressure drops in the PAs and the DAo section downstream the PS could be ignored, the pressure drop in the AAo, supraoartic branches, and DAo section upstream the Potts shunt required an appropriate lumped representation. Each of these segments is described by a linear resistance , a nonlinear quadratic resistance , and an inductance . The values of these parameters are found through a regression analysis on the GMM solution (Pant et al. 2013, 2014). The resulting stand-alone LPM, excluding the components that yielded near-zero values in the regression analysis, is shown in Fig. 3. The post-operative stand-alone LPM includes a shunt law model, shown in magenta colour in Fig. 3, which is described in further detail in Sect. 2.2.5.

Fig. 3

Stand-alone LPM for the iPAH and PS physiology. The lumped components that represent the 3D regions (see Figure 1) are shown in blue. In the post-operative state, the lumped representation of the Potts shunt is shown with components in magenta. represents volumetric flow-rate; represents pressure, represents volume

Single fibre model for heart chamber

Each heart chamber is described by a single fibre wrapped around the cavity (Arts et al. 1991; Bovendeerd et al. 2006; Pant et al. 2016, 2017). This model describes the following relationship between the cavity pressure , cavity volume , the wall volume , and the stress in the fibre, :The dependence on volume alone, as opposed to shape of the cavity, makes the single fibre model particularly suitable for 0D modelling, as it has been shown that under rotational symmetry, the shape of the chamber has little effect upon the cavity pressure and fibre stress (Arts et al. 1991). The fibre stress comprises an active component and a passive component, , whereby:The active component of stress is described bywhere c represents contractility, represents the maximum active sarcomere stress, f(l) represents the force-length relationship in a sarcomere of length l, represents the time-dependent variation of active stress with representing time elapsed since activation, and represents the active viscous stress dependent on the sarcomere shortening velocity. For a sarcomere length l, f(l) is described aswhere the length parameters , , and are sarcomere material constants. With and representing the cavity volume and sarcomere lengths at zero transmural pressure, the sarcromere stretch for a generic sarcomere length l and cavity volume is given byThe time-dependent activation is described aswhere is the time since the beginning of cavity activation, is the total systole duration of the specific cardiac chamber, and controls the shape of the activation curve. The active viscous stress is described as:where is the initial sarcomere shortening velocity and governs the shape of the velocity-stress relation, and the shortening velocity isLastly, the passive stress component, , is modelled by:where and are sarcomere material constants. Thus, Eqs. (1)–(9) define the relationship between pressure, volume and flow rate for each cardiac chamber.

Heart valve model

The adopted model (Mynard et al. 2012; Pant et al. 2016, 2017) determines valve dynamics based upon the instantaneous pressure difference across the valve. Neglecting viscous losses, the relationship between the pressure drop across the valve and the volumetric flow-rate Q through the valve is described by the Bernoulli relation as:where B and M represent the valve nonlinear resistance and inertance of the valve, respectively, and depend on the geometrical parameters of the valve and the fluid density aswhere and are the effective area and effective length of the valve, respectively. The time varying dynamics of is described as:where is a variable describing the valve state, and and are the minimum and maximum effective areas, respectively. Based on Mynard et al. (2012), the valve state depends on the pressure difference across the valve and is described by:where and are parameters that govern the valve opening and closing rates, respectively. The valve starts to open when and is fully open when is equal to 1. Similarly, the valve starts to close when and is fully closed when is equal to 0.

LPM for shunt

Here, the lumped model for shunts proposed in the literature (Migliavacca et al. 2000, 2001) is adopted for modelling the PS. The pressure drop across the PS is described as:where is the pressure in the LPA, is the pressure in the DAo, is the instantaneous volume flow-rate in the shunt and , and are functions of the shunt/vessel geometry and the fluid properties. Equation (14) is often expressed in terms of the shunt diameter, , such that:where , and are proportionality constants. In lieu of experimental data, and are derived computationally using regression analysis on the GMM solutions for all the post-operative geometries with shunt diameters varying from 5 mm to 10 mm. This analysis also showed that the last term in Eq. (15) accounts for only a 0.2% difference in pressure drop across the PS. Thus,this term is neglected and the pressure drop across the shunt is described by:

Solution methodology

For the 3D part of the GMM, blood is assumed to be an incompressible Newtonian fluid with density 1.06 gcm and dynamic viscosity 0.04 gcms, governed by the Navier–Stokes equations and a rigid wall assumption. A no-slip condition is imposed on the walls (i.e. surfaces that are neither inlets nor outlets) and a backflow stabilisation coefficient of 0.5 is employed to control backflow at the outlets (Moghadam et al. 2011). The inlets and outlets of the 3D model are coupled to the LPM by interface conditions such that time-dependent pressures yielded by the LPM are uniformly applied to the boundaries of the 3D domain and the volumetric flow rates at the 3D inlets and outlets are imposed as forcing terms to the LPM (Moghadam et al. 2013). The GMM was developed within the SimVascular Suite (Updegrove et al. 2017; Lan et al. 2018), which is an open-source software utilising the finite element method for solving the Navier–Stokes fluid flow equations, and widely utilised for patient-specific modelling, medical device design, and surgical planning. A constant time-step of seconds is specified, which was selected following a comparison of results from three different GMM simulations with fixed time increments of 1.0, 2.0, and 5.0 10 seconds. The simulations were run for multiple cardiac cycles until the relative changes in all the variables were less than 0.5% of those in the previous cardiac cycle.

The stand-alone LPM network, Figure 3, leads to a set of ordinary differential equations which are solved using a fourth-order Runge–Kutta scheme with constant time-step of 1.0 10 seconds. Within 25 cardiac cycles, the relative changes in the state variables were less than 0.1% when compared to those in the previous cardiac cycle. These converged (stabilised) values for the pre-operative state variables are then used to define the initial conditions for the post-operative state in order to expedite convergence.

Patient-specific parameter specification and estimation

For a reference set of parameters, the lumped parameters representing the 3D sections are first estimated through a regression analysis on the GMM solutions, see Sect. 2.2.2 and the references therein. With this, the stand-alone LPM includes correct 3D representation, but lacks patient-specific parameters. Subsequently, patient-specific parameter estimation is performed on the stand-alone LPM, Fig. 3, as it is significantly cheaper in terms of computation time relative to the GMM.

For parameter specification and estimation, the measurements acquired in the patient, Tables 1 and 2, are utilised. In the absence of direct pressure measurements in the ascending aorta, the pre-operative aortic pressure is assumed comparable to the pressure in thoracic portion of the descending aorta . This surrogate measurement of is used for estimation of model parameters as it better constrains the estimates compared to , and allows for the computation of systemic vascular resistance for further parameter estimation (see Sect. 2.4.2). In what follows, where manual tuning is used the parameters are adjusted to reproduce the following target measurements: surrogate pre-operative aortic and PA pressures (systolic, diastolic, and mean), the LV and RV EDVs, LV and RV stroke volumes, and the CO.

Heart and valve models

From the pre-operative patient measurements (see Sect. 2.1 and Table 2), many heart model parameters can be directly specified. In particular, the cardiac cycle time-period of 0.9s, the myocardium wall volumes (), and the activation durations () for each heart chamber are directly specified. For the geometrical parameters of the heart model, the unknown parameters are the chamber volumes at zero transmural pressure, , which are manually tuned, and found to be 20.6 ml, 39.1 ml, 28.9 ml, and 43.8 ml for the LA, RA, LV, and RV, respectively. Table 4 displays the parameters pertaining to the sarcomere active and passive stress laws in the single fibre model, which are based upon the experimental values provided in (Bovendeerd et al. 2006) and then manually tuned.

Table 4

Patient-specific pre-operative parameters for the active and passive stress laws

Parameters describing the sarcomere active and passive stresses
T_a0[kPa]	71.0	T_p0[kPa]	0.9	c_v[-]	0.0	c_p[-]	12.0
v_0[μm/s]	10.0	l_0[μm]	2.0	l_a0[μm]	1.5	l_ar[μm]	2.0
EaLA[-]	1.0	EaRA[-]	1.0	EaLV[-]	1.05	EaRV[-]	0.48

Table 5 shows the opening and closing rates of the valves, provided in (Mynard et al. 2012), together with the effective valve areas, , which are initially approximated from literature (Mynard et al. 2012) and then manually tuned to ensure that the flow rates for the atrioventricular (AV) valves are consistent with the pre-operative echocardiographic velocity traces.

Table 5

Patient-specific pre-operative parameters for the valve model: mitral valve (MiV), aortic valve (AoV), tricuspid valve (TrV), and pulmonary valve (PuV)

Valve model parameters
	MiV	AoV	TrV	PuV
K_vo[cm s/g]	0.03	0.012	0.03	0.02
K_vc[cm s/g]	0.04	0.012	0.04	0.02
A_eff[cm^2]	2.1	2.2	2.1	2.3

Global circulation parameters

To specify the remaining circulation parameters, first the total vascular resistances and compliances are computed, which are then appropriately distributed to individual branches and components. Pulmonary vascular resistance (PVR) is derived from the measured CO, the measured mean , and the measured mean . Similarly, the total systemic vascular resistance (SVR) is derived from the measured CO, the measured mean and the measured mean . For the systemic circulation, a mean flow-split between the upper and lower body is assumed to be 30:70, thus providing total resistances for the upper and lower body compartments. The flow within the aortic arch branches is assumed to be proportional to the outlet areas. In the absence of flow-split measurments, this method of flow-distribution based on the area of the outlet branches has been previously used in other studies (Troianowski et al. 2011; Arbia et al. 2015). The flow-split between the right and left pulmonary circulations (and further within the two LPA branches) is also assumed to be proportional to the outlet areas leading to the corresponding compartments. The total resistances subsequent to each outlet are computed based on the these mean flow rates and assuming no pressure drop in the 3D region. This assumption of flow-splits (based on outlet areas) to compute total resistances may be alleviated in future studies if reliable measurements or estimates of the flow-splits are available. The outlet areas are specified in Appendix A.

The total resistances subsequent to each outlet are split into individual branch components (for example, such a branch corresponding to the DAo outlet is highlighted in Fig. 1 in brown colour) as follows. The resistance is assumed to be 10 times larger than (LaDisa et al. 2011) and the venous resistance is assumed to be 10% of the total branch resistance. Thus, in each branch these resistances , , and are assumed to be 8.2%, 81.8%, and 10% of the total branch resistance , respectively.

With the resistances defined, the compliances within the three-element Windkessel models are related to the time constant of the exponential pressure decay such that C = / (Spilker and Taylor 2010). For the PAs, is defined by:where SV is the stroke volume and PP is the pulse pressure (Ross et al. 2013; Chemla et al. 2016), enabling the compliances for the pulmonary artery Windkessels to be calculated. Similarly, using Eq. (17), the compliances for the systemic branch arteries are calculated using SVR instead of PVR. The venous compliances are calculated for each branch from the relationship = , with is set equal to 1s for the lungs and equal to 2s for the upper and lower body branches based on the literature (Baretta 2014; Kilner et al. 2009; Presson et al. 1998; Spilker et al. 2007). Such assumptions based on the literature, which are here applied to all branches identically, are necessary in the absence of time-varying and simultaneous measurements of pressure and flow-rates in the major arteries and veins. When such measurements are available, the individual compliances may be estimated directly from this data. The AAo compliance and mPA compliance were manually adjusted to align with the target pre-operative measurements.

For the post-operative stand-alone LPM, the proportionality constants in Eq. (16) describing the behaviour of flow across the shunt are derived from regression analysis on the GMM solutions for all non-protruding shunt geometries with diameters of 5, 7, and 9 mm (see Appendix B.1). This analysis ensures that 3D behaviour of flow across the shunt is well captured in the stand-alone LPM. The set of all patient-specific lumped parameters describing the global circulation are summarised in Table 6.

Table 6

Patient-specific pre-operative parameters for Windkessel models and the lumped parameters representing 3D regions: R in units of g/cm4s, K in units of g/cm7, L in units of g/cm4, and C in units of cm4s2/g

Pre-operative parameters: windkessel and lumped representation of 3D regions
RpIN	701.0	RdIN	7273	C^IN	1.90e-04	RpLC	1218.4	RdLC	12375	C^LC	1.12e-04
RpLS	1965.9	RdLS	27674	C^LS	5.00e-05	RpDAo	35.2	RdDAo	1669	C^DAo	8.28e-04
RpLPA-1	1001.7	RdLPA-1	10017	C^LPA-1	1.62e-04	RpLPA-2	915.2	RdLPA-2	9152	C^LPA-2	1.62e-04
RpRPA	231.3	RdRPA	2313	C^RPA	7.00e-04	C^AAo	4.4e-04	C^mPA	1.5e-04	RvsUB	480.4
CvsUB	4.16e-03	RvsLB	205.9	CvsLB	9.71e-03	RvsLL	584.5	CvsLL	1.71e-03	RvsRL	282.8
CvsRL	3.54e-03	R3DAAo	10.1	K3DAAo	0.167	K3DIN	0.839	L3DIN	8.04	K3DLC	1.32
L3DLC	17.97	K3DLS	88.7	L3DLC	38.58	R3DDAo	21.2	K3DDAo	0.98	L3DDAo	8.41

Heart parameters governing pulmonary artery hypertension

In the model, the generation of higher RV pressures due to the PAH relative to a healthy state is mechanistically governed by the two phenomena: first, through a combination of RV wall volume and , as varying these parameters changes the amount of myocardial muscle and its stretch, see Eq. (5), which affects the fibre stresses and consequently the pressures; and second, through the active stress produced by the myocardial sarcomeres, see Eq. (3). Since our study is concerned with immediate effects of the PS creation, the myocardial wall volume, as it was measured in the patient, was fixed in the model. Thus, to generate higher RV pressures corresponding with the clinical measurements, the active stress produced by the sarcomeres in the RV wall must be higher in comparison to a healthy state. To achieve this, the shape of the time-dependent contraction , which is governed by the exponent in Equation (6), is altered. The alternative approach to produce higher RV pressures is to alter the contractility coefficient c in Eq. 3. This is akin to effectively varying the maximum stress as it is the product of c and that appears in Eq. 3. This approach is also tested, and it is found that varying the contractility coefficient does not influence global haemodynamics (see Appendix C), when compared to varying in Equation 3.

Results

All results are presented from the last cardiac cycle in the GMM simulation after the solution has converged. In the main manuscript, the results for the GMM are presented, while the stand-alone LPM results are shown in Appendix B.

Comparison of simulation output with pre- and post-operative measurements

The comparison of simulation output with pre- and post-operative measurements is performed for the 7.6mm diameter PS corresponding to the procedural outcome in the patient. While the model does not include protruding ends for the stent for this comparison, in Sects. 3.5 and 4.3 it is shown that protrusion has minimal effect on global haemodynamics. Tables 7-I shows the GMM output against the measurements for the quantities that are used for patient-specific parameter estimation. For the pre-operative state, the GMM solution is generally within 2.0% of the clinical measurements and only differs by a maximum of 5% for the mean aortic pressure. This agreement between the pre-operative computational solution and the clinical measurements is corroborated by the pressure waveforms presented in Fig. 4.

Table 7

Pre- and post-operative , , EDV, SV, EF, and CO from the GMM against clinical measurements

	P_AAo [mmHg]			P_mPA [mmHg]			EDV [ml]		SV [ml]		EF		CO [L/min]
	Systolic	Diastolic	Mean	Systolic	Diastolic	Mean	LV	RV	LV	RV	LV	RV	LV	RV
I. Pre-operative measurements and model output (parameter estimation)
Measurement	94.0	53.0	69.0	112.0	67.0	85.0	66.9	91.8	51.0	51.0	0.76	0.55	3.4	3.4
GMM	95.4	52.5	65.4	111.1	68.3	84.6	67.3	92.0	50.7	50.7	0.75	0.55	3.38	3.38
II. Post-operative measurements and model output (validation)
Measurement	97.0	51.0	71.0	102.0	54.0	76.0	–	–	–	–	–	–	–	–
GMM	104.1	60.8	74.3	106.9	61.3	76.8	68.9	83.4	42.7	59.6	0.62	0.71	2.85	3.97

Fig. 4

Comparison of pre-operative pressure waveforms generated by the GMM against the measurements over one cardiac cycle

For qualitative validation, the simulation output for the pre-operative state is compared against the Doppler velocity tracings for flow across the AV valves, the only clinical data that was not directly used for model parameter estimation (Fig. 5). As can be observed, the model accurately captures the flow velocity profiles across the AV valves. In particular, the flow over mitral valve (MiV) reflects normal LV diastolic function (E>A) in contrast to the flow profile across the tricuspid valve (TrV) demonstrating impaired filling and diastolic dysfunction of the pressure-loaded RV (E

Fig. 5

Comparison of GMM solution for pre-operative volumetric flow rates through the atrioventricular valves with Doppler measured velocity tracings

Further validation of the model and its predictive capabilities are assessed by comparing the post-operative GMM simulation output against the respective clinical measurements (Table 7-II). The GMM output is in reasonable agreement with the measured pressures, reproducing the mean mPA pressure and mean AAo pressure within 1% (absolute error = 0.8 mmHg) and 5% (absolute error = 3.3 mmHg) of the measurements, respectively. The predicted systolic pressures in the GMM exceed the measurements by 5% (absolute error = 4.9 mmHg) for the mPA and 7% (absolute error = 7.1 mmHg) for the aorta. The diastolic pressures exceed the measurements by 13% (absolute error = 7.3 mmHg) and 19% (absolute error = 9.8 mmHg) for the mPA and the AAo, respectively.

Global haemodynamic changes due to the PS creation

The key advantage of modelling is that detailed haemodynamic changes due to PS creation can be assessed at all locations within circulation. Figures 6, 7, and 8 depict such changes occurring in one cardiac cycle in pressures, flow-rates and absolute volumes passing through various compartments, respectively, for the scenario with the non-protruding PS of 7.6mm diameter. Major changes observed (Figs. 6, 7, and 8 and Table 7) in haemodynamics due to the PS creation as predicted by the GMM are as follows:

Fig. 6

GMM results for the 7.6 mm diameter PS showing pre- to post-operative changes in pressure at key locations in the arterial network. PV loops are additionally included. In the PV loop plots, the x-axis represents volume [ml] and y-axis represents pressure [mmHg]. In all other plots the x-axis represents time [s] and the y-axis represents pressure [mmHg]. The valve parameters are dimensionless. For a key to symbol nomenclature, please see Fig. 1. represents the pressure gradient across the PS

Fig. 7

GMM results for the 7.6 mm diameter PS showing pre- to post-operative changes in flow-rate and volumes at key locations in the arterial network. In all the plots the x-axis represents time, and y-axis for volumes, , is in [ml], while for the flow-rates, , is in [ml/s]. For a key to symbol nomenclature, please see Fig. 1

Fig. 8

GMM results for the 7.6 mm diameter PS showing pre- to post-operative changes in volume of blood flowing in one cardiac cycle through the circulatory system. All numerical values are for volumes in ml, and values in red represent pre-operative state while those in blue represent post-operative state

The mean flow-rate through the shunt is ml/s (corresponding to a volume displacement of ml in one cardiac cycle) with a mean pressure gradient of mmHg across the PS.

The pressure in the PAs decreases while the pressure in the aorta increases, becoming nearly equal to each other.

The LV output decreases by 16% while the RV output increases by 18%.

The LV ejection fraction decreases (from 0.75 pre-PS to 0.62 post-PS) while that of the RV improves (from 0.55 pre-PS to 0.71 post-PS).

The LV EDV remains nearly unchanged (increase by 2%), while LV ESV increases substantially (by 58 %), resulting in a net decrease in LV SV (by 16%), and hence the LV output (by 16%).

Even though the LV output decreases, the flow from the LV into the aortic arch branches, and hence to the upper body and the superior vena cava (SVC), increases by 16%. In contrast, the LV contribution to the flow into the DAo decreases by 29%.

The total flow in the DAo downstream the PS, however, increases due to RV contribution to the flow through the PS, corresponding to an increase of the flow in the inferior vena cava (IVC) (by 18%).

The work done by the LV (assessed by the area within the ventricular PV-loop) decreases from mmHg-ml pre-operatively to mmHg-ml post-operatively.

The RV EDV decreases by 9% and its end-systolic volume decreases by an even larger amount (42%), resulting in a net increase in RV stroke volume (by 18%), and hence the RV output (by 18%).

Even though the RV output increases, the net flow into the RPA and the LPA downstream the PS decreases (by 16% and 15%, respectively) due to the post-ventricular right-to-left shunt. Correspondingly, flow to the left and right pulmonary veins also decreases (by 15% and 16%, respectively), reflecting the blood flow diversion through the proximal LPA into the shunt, the volume of which exceeds the volume corresponding to the increase in RV output.

The work done by the RV (the area within the PV-loop) increases from mmHg-ml pre-operatively to mmHg-ml post-operatively.

Defining as the sum of flow-rates in the left and right pulmonary veins and as the sum of flow-rates in the inferior and superior venae cavae, the ratio changes from 1.0 pre-operatively to 0.72 in the post-operative state.

The valve function remains largely unaffected for all the valves.

The operating volumes (minimum and maximum volumes during a cardiac cycle) of the LA increase, which is accompanied by a corresponding increase in the LA pressures; while the operating volumes of the RA decrease, with a corresponding decrease in the RA pressures.

The pressure changes in the inferior and superior venae cavae are less than 2 mmHg, while those in the left and right pulmonary veins are less than 3mmHg.

Spatio-temporal haemodynamics around the PS

The GMM presents an advantage that local flow features in the 3D domain can be assessed. Figure 9 shows streamlines that originate from the mPA (thus showing only the RV contribution to the blood flow) at different points in the cardiac cycle. The streamlines are coloured by the velocity magnitude and highlight three phenomena in the simulated case of shunt flow going in direction perpendicular to the walls of anastomosed vessels:: During the cardiac cycle, the pressure gradient across the shunt increases in systole reaching its maximum at approximately peak-systole and then reduces to a slightly negative value at end-systole, quickly recovering to approximately zero pressure gradient for most of diastole (Fig. 6). The corresponding flow-rate through the shunt follows a pattern close to that of the pressure gradient (Fig. 7), given the pressure gradient is driving the flow.

Fig. 9

Velocity streamlines that originate from the mPA for the 7.6 mm diameter PS at different times of the cardiac cycle

the flow velocities in the PS are significantly higher in comparison to those in the PAs (during peak systole the velocity magnitudes within the shunt average to approximately 225 cm/s);

the jet from the PS impinges with high velocity on to the opposite wall of the DAo, particularly in the systolic phase; and

while most of the flow from the shunt is diverted to the distal DAo, the streamlines (with pathlines showing a similar pattern) indicate that a small amount of flow is going to the aortic arch branches.

Figure 10 shows the pressure variation in the 3D flow domain at various points in the systolic phase, and shows the high pressure created in the region of the DAo opposite to the shunt by the impinging jet. Similarly, Fig. 11 shows wall shear stress (WSS) variation in the 3D region at various points in the systolic phase, and highlights high WSS regions within the DAo wall created by the shunt jet.

Fig. 10

Wall pressure for the 7.6 mm diameter PS at different times within the systolic phase

Fig. 11

Wall shear stress (WSS) for the 7.6 mm diameter PS at different times within the systolic phase

Varying shunt diameter

The results for simulations of the pre-operative state and post-operative state with PS diameters varying between 5mm and 10mm are presented in Tables 8 and 9. The length of each PS in these simulations is constrained by the distance between the LPA and DAo such that the PS does not protrude into either vessel. The aortic pressure increases and the PA pressure decreases with increasing shunt diameter (Table 8). The systolic pressures are affected significantly more than the diastolic pressures. For example, when comparing the extreme diameters of 5mm to 10mm, the systolic increases by nearly 6% and the systolic decreases by 7%. In contrast, the diastolic pressures are largely unaffected—the diastolic and both vary by less than 1%.

Table 8

GMM: , , EDV, SV, EF and CO for different PS diameters

Diameter [mm]	P_AAo [mmHg]			P_mPA [mmHg]			EDV [ml]		SV [ml]		EF		CO [L/min]
	Systolic	Diastolic	Mean	Systolic	Diastolic	Mean	LV	RV	LV	RV	LV	RV	LV	RV
pre-op	95.4	52.5	65.4	111.1	68.3	84.6	67.3	92.0	50.7	49.9	0.75	0.54	3.38	3.38
5	101.2	60.2	72.9	111.0	61.7	78.9	69.1	86.5	44.9	58.2	0.65	0.67	2.99	3.88
6	102.3	60.7	73.6	109.5	61.3	77.9	69.2	85.2	43.8	59.0	0.63	0.69	2.92	3.93
7	103.3	60.8	74.0	107.9	61.3	77.2	69.1	84.0	43.0	59.4	0.62	0.71	2.87	3.96
7.6	104.1	60.8	74.3	106.9	61.3	76.8	68.9	83.4	42.7	59.6	0.62	0.71	2.85	3.97
8	104.6	60.8	74.4	106.4	61.3	76.7	68.7	83.0	42.4	59.7	0.62	0.72	2.83	3.98
9	105.9	60.8	74.7	104.8	61.3	76.2	68.2	82.2	41.9	60.0	0.61	0.73	2.79	4.00
10	106.9	60.8	75.0	103.3	61.2	75.9	67.8	81.5	41.4	60.2	0.61	0.74	2.76	4.02

Table 9

GMM: work done (area under the PV-loop) by the ventricles and mean volumetric flow rate, , total volume displaced in a cardiac cycle, , and pressure drop, , across the shunt for different PS diameters

D_SH	W_LV	W_RV	Q¯_SH	V¯_SH	ΔP¯_SH	Q_p/Q_s
[mm]	×10^3[mmHg-ml]	×10^3[mmHg-ml]	[ml/s]	[ml]	[mmHg]	[-]
pre	4.13	4.95	–	–	–	1
5	3.83	5.80	15.1	13.6	8.18	0.78
6	3.75	5.86	17.0	15.3	6.28	0.75
7	3.69	5.87	18.0	16.2	5.20	0.73
7.6	3.66	5.87	18.7	16.9	4.58	0.72
8	3.65	5.87	19.5	17.6	4.11	0.72
9	3.64	5.85	20.5	18.5	3.51	0.70
10	3.63	5.82	21.3	19.2	2.88	0.69

With the smaller shunt diameters of 5 mm and 6 mm, the PS is not able to reduce the PA pressures to near systemic levels and the pulmonary-to-systemic systolic arterial pressure ratio only reduces from a pre-operative value of 1.16 to post-operative values of 1.1 and 1.07 for the 5 mm and 6 mm diameter shunts, respectively (Fig. 12a,b). Increasing the shunt diameter to 7 mm and 8 mm, resulted in reduction of the PA pressures to near systemic levels, with a post-operative pulmonary-to-systemic systolic arterial pressure ratio of 1.04 and 1.02, respectively (Fig. 12c,d). With further increments in shunt diameter to 9 mm and 10 mm, more blood is pumped from the pulmonary circulation to the DAo resulting in the exceeding the in some part of systole (Fig. 12e,f). The mean pressure gradient across the shunt decreases and the flow-rate through the shunt increases with increasing shunt diameters (Table 9). Both of these show a near-linear variation in the range of diameters tested (Fig. 13a). The volumetric flow-rate profiles through the shunt differ for varying diameters (Figure 13b). For a diameter of 7 mm, the flow rate is positive throughout the cardiac cycle. However, increasing the shunt diameter to 8 mm results in a small period of reversed flow, which is amplified with further increments in diameter.

Fig. 12

GMM: post-operative and for different shunt diameters

Fig. 13

GMM: effect of varying shunt diameter on flow-rate and pressure gradient across the shunt

The ventricular PV loops show that the LV EDV remains almost constant, while its ESV increases as the shunt diameter increases (Fig. 14a). As a result, the LV output progressively decreases with increasing shunt diameters (Table 8). This reduction in LV cardiac output leads to successive worsening of LV ejection fraction. For the RV, both EDV and ESV decrease with increasing shunt diameters (Fig. 14b). Decrease in ESV is, however, larger than that in EDV resulting in larger stroke volumes, and hence higher RV outputs and improving RV ejection fractions (Table 8). An important observation from the PV loops is the change in work performed by the ventricles. The area enclosed by the PV loop can be seen as a measure of the ventricular stroke work, which is a product of the stroke volume and the mean ventricular pressure. The work done by the LV after PS creation decreases with increasing shunt diameters (consistent with decreasing stroke volume pumped against higher aortic pressures), while the work done by the RV increases with increasing shunt diameters up to a certain point and slightly decreases for the larger diameters (consistent with increasing stroke volume pumped against lower PA pressures) (Table 9). The latter is an important observation from the point of view of RV dynamics after PS creation: together with the reduction in RV pressures due to increasing shunt diameters, the work done by the RV nevertheless increases due to increased volumes pumped by the ventricle. However, in contrast to this increase in the RV stroke volume, its EDV decreases after PS creation, and this change in EDV is more pronounced for larger shunt diameters.

Fig. 14

GMM: effect of shunt diameter on pressure-volume loops for a left ventricle and b right ventricle

The 3D assessment of flow with varying PS diameters is performed for the extreme diameters of 5mm and 10mm. Figure 15 shows velocity streamlines originating from the mPA for these two cases at two time points in systole. The corresponding surface pressures and wall shear stresses are similarly shown in Fig. 16. One can see that flow complexity is higher with higher diameter shunts. The larger PS diameter results in significantly higher flow through the shunt, but this higher flow keeps the maximum velocities comparatively less affected due to the diameter increase (Fig. 15). Higher shunt diameters lead to higher pressures in the aorta, and for high diameters—9mm and 10mm shunts as seen in Fig. 12e,f— the AAo pressure can become equal or even exceed the mPA pressure during periods of systole (also see Fig. 16c). Smaller shunt diameters lead to smaller localised regions of high pressure on the DAo wall, due to the directed and more focused shunt flow jet of smaller size (Fig. 15b). Finally, changing the shunt diameter does not have a significant effect in reducing regions of abnormally high WSS (Fig. 16).

Fig. 15

Velocity streamlines for 5 mm and 10 mm diameter shunts

Fig. 16

Surface pressure and wall shear stress (WSS) for 5 mm and 10 mm diameter shunts at 0.09 s, 0.15 s, and 0.22 s

Varying length of the communicating stent (protrusion)

The results for changing length of the stent used to create the PS leading to protrusion into the LPA and the DAo are presented in Table 10, which shows the PA and aortic pressures for shunts with lengths of 10, 15, 20, and 25 mm against the non-protruding PS with length of approximately 4.2 mm. The shunt diameter for all the cases of varying lengths is kept constant at 7mm. Figure 17 shows the corresponding pressure traces in the AAo and mPA. Finally, Fig. 18 illustrates the flow-rate across the shunt during the cardiac cycle with varying shunt lengths. Interestingly, increasing the length of the communicating stent from 10mm up to 25mm does not result in substantial increase of the pressure gradient across the PS. Similarly, volumetric flow-rate through the shunt with increasing length is also not substantially affected with the only minor changes observed in the diastole when the shunt flow-rates are small. These observations imply that at sufficiently large shunt diameter, the global haemodynamics is not significantly affected by the shunt length as it does not affect the resistance to fluid flow in the principal directions.

Table 10

and for 7 mm diameter PS with different lengths

Length	P_AAo			P_mPA
	Systolic	Diastolic	Mean	Systolic	Diastolic	Mean
[mm]	[mmHg]	[mmHg]	[mmHg]	[mmHg]	[mmHg]	[mmHg]
≈ 4.2 (no protrusion)	103.3	60.8	74.0	107.9	61.3	77.2
10	103.6	60.7	74.1	108.4	61.2	77.2
15	104.1	60.6	74.3	108.4	61.2	77.2
20	104.7	60.5	74.5	108.5	61.2	77.3
25	104.8	60.5	74.6	108.8	61.2	77.5

Fig. 17

GMM: effect of shunt length on and

Fig. 18

GMM: effect of shunt length on flow across the shunt

For the 3D assessment of flow, velocity vectors in the PA in a plane perpendicular to the PS are considered. At peak systole (s) and peak deceleration phase (s), these velocity vectors for all the four shunt lengths are shown in Fig. 19. It shows that as the protruded areas of the stent interfere with the flow in the LPA: recirculating zones of flow are formed in the downstream region of the protrusion. Blood velocity in the eye of such von Kármán-type of vortices is low and these recirculating regions are larger and higher in number as the stent’s protrusion length increases. Thus, while the global flow features are not affected significantly, the local flow features, as demonstrated by the 3D assessment, are substantially altered by the protrusion length.

Fig. 19

Flow field visualisation in the PA for different shunt lengths at 0.3–0.45 s

Stand-alone LPM output

The stand-alone LPM also reliably predicts the global haemodynamics and its changes due to PS creation. While the GMM results are of higher fidelity, and are therefore used in the main text to assess physiology, the stand-alone LPM model results are of comparable, albeit of slightly lower, fidelity. In order to avoid overcrowding of the main text, the results from the stand-alone LPM model are presented separately in Appendix B. Importantly, with the exception of the local flow-features and assessment of varying shunt lengths, the stand-alone LPM accurately predicts all other haemodynamic parameters, including those that vary due to changes in shunt diameter.

Discussion

Model validation

LPMs, known interchangeably as a zero-dimensional (0D) models, are frequently employed to analyse global haemodynamics in the circulatory system over a range of physiological and pathological conditions (Shi et al. 2011). However, analysing complex local flow fields, typically associated with surgical or interventional procedures, often necessitates a more detailed geometric model in the region of interest. Hence, two closed-loop computational models have been developed: a stand-alone LPM; and a GMM consisting of a reconstructed patient-specific three-dimensional (3D) flow domain, coupled to an LPM. In this study, due to a lack of the complete set of time-varying traces of pressure and flow-rate measurements in the major arteries and veins, we adopt a largely manual approach to parameter estimation. Alternatively, if such measurements are available, it would be possible to employ automated parameter estimation methods, such as those based on data-assimilation methods (Pant et al. 2016, 2017, 2014), which can also account for measurement uncertainty.

The results of Sect. 3.1 show that when comparing the post-operative predictions with clinical measurements of the arterial pressures, the largest relative errors are in the diastolic phase, but the absolute differences do not exceed 10 mmHg. The larger disparities in the post-operative state may be attributed to insufficiently robust parameter estimation, largely due to an incomplete set of pre-operative measurements available for the patient, the inherent uncertainty in the used clinical measurements (measurement device errors, procedural errors, and the fact that the measurements do not correspond to steady-state haemodynamics simulated in the model, where the haemodynamic outputs do not vary from one cardiac cycle to another), potentially non-unique set of model parameters, and auto-regulatory/adaptive mechanisms, which may be present in the post-operative state and not accounted by the model. Divergences between the model output and clinical measurements may also be present due to lumped representations of the heart and valve models, with inherent approximations. Notably, a previous study for PS analysis in a non patient-specific manner (Delhaas et al. 2018) utilised the three-wall segment model for the heart (Lumens et al. 2009), which explicitly accounts for ventricular interactions through the inter-ventricular septum. These interactions are ignored in this study as septal bulging has been shown to have little effect on global haemodynamics (Lumens et al. 2010), but their inclusion may play a role in achieving better agreement of model output with post-operative measurements. Lastly, explicit modelling of 3D fluid-structure interaction may also help in reducing the discrepancies between model output and measurements.

Although many of the aforementioned limitations can be addressed with more sophisticated models, it will also magnify the problem of estimating additional parameters for these models, while clinical measurements remain limited and uncertain. Overall, given that the GMM reproduces all the pre-operative clinical measurements within 5% differences, and that differences in post-operative predictions versus the measurements are less than 20%, the GMM presented here as well as the stand-alone LPM depicted in Appendix B can be considered as a reasonably good model for assessing the effects of PS creation: apart from the diastolic aortic pressure which in the clinical situation surprisingly decreased after PS creation, both models predict variations in the same direction and similar magnitude as can be seen in the measurements.

Haemodynamic changes induced by the Potts shunt

The results of Section 4.2 show that the connection of the high-pressure PAs with the low-pressure aorta naturally results in reduction of PA pressures, increase of aortic pressures, and a commensurate flow from the LPA to the DAo through the shunt. The changes induced in global haemodynamics, see Fig. 8, are best explained by consideration of the effective resistances faced by the LV and RV in the post-operative state. The RV, due to the bypass created by the PS, experiences a lower effective resistance (alternatively, a lower afterload), which results in a larger stroke volume generated by the RV, explaining its higher output. Between the left and the right pulmonary circulations, the reduction in resistance occurs in the left branch, explaining the increased flow in the proximal LPA section and a reduced flow in the distal LPA and the RPA. The flow volume through the shunt is, however, larger than the increase in the RV output with the difference occurring at the cost of flow reduction in the RPA. Effectively, the PS steals all the increase in the RV output and additionally some flow from RPA, thus resulting in decreased effective pulmonary blood flow despite higher RV output when compared to the pre-operative state. This apparently adverse effect of the PS, in terms of reduced return of oxygenated blood from the lungs to the LV, may be beneficial in the long-term as reduced flow through the pathological pulmonary vascular bed can potentially result in regression of the arterial wall media hypertrophy similar to situation seen after correction of certain congenital heart defects (Hsu et al. 2016).

The LV, on the other hand, largely due to the increased aortic pressures created by the PS and, possibly additionally, because of the flow competition between the streams coming from the shunt and aortic arch (also potentially due to interference with the protruding end of the covered stent), experiences a higher resistance (alternatively, a higher afterload) in the post-operative state. This results in a lower volume of blood ejected by the LV, explaining its lower output. The increase in the effective resistance for the LV occurs, however, in aortic segment downstream of the PS level, explaining (i) an increased flow to the aortic arch branches, and hence to the SVC, and (ii) the reduced LV contribution to flow in the DAo and to the lower body compartment. The flow in the DAo downstream the PS is the sum of oxygenated blood flow from the LV and the deoxygenated blood flow from the RV coming through the shunt, resulting in an increased IVC flow compared to the pre-operative state.

In addition to the increased afterload for the LV after PS creation, there is an issue of reduced pulmonary venous return to the LV due to the post-ventricular right-to-left shunt. Together with reduced LV stroke volume after the PS creation, it could result in decreased perfusion of the upper body and coronary circulation with the oxygenated blood. The simulations for this particular patient demonstrate increase of the blood flow originating from the AAo to the aortic arch branches due to the higher aortic pressures after the PS creation. This results in a better perfusion of the upper body with oxygenated blood. An extrapolation of this reasoning suggests that the higher AAo pressures should also result in increase of the coronary arterial perfusion and eventually sustained (or increased) oxygen delivery to the ventricular myocardium despite the decreased LV output. In the clinical series of PS, there are, however, some patients who demonstrated acute haemodynamic deterioration with cardiac arrest after PS creation, presumably due to a catastrophic drop in the LV output and coronary arterial perfusion. Thus, while ‘compensatory’ increase of coronary perfusion may be seen as a positive course of events, degree of reduction in the LV output remains an overall cause for concern as it determines the level of blood oxygenation going to the coronary circulation. While physics dictates that coronary perfusion will be higher due to higher aortic pressures after PS creation, in the scenario of extreme right-to-left shunting that induces severe reduction in LV output, the increased flow to the coronaries (and also the aortic arch branches) will occur largely retrogradely with the deoxygenated blood coming from the RV through the PS. Furthermore, in patients with long-standing severe iPAH, the myocardium of the LV has been shown to become atrophic with contractile cardiomyocyte dysfunction (Manders et al. 2014). The substantial reduction in the LV output after PS creation in combination with pre-existent LV myocardial atrophy is also of potential concern, as in such iPAH patients the degree of LV stroke volume decrease will be so severe that the LV output may not compete with the deoxygenated blood flow coming into the ascending aorta from the high-pressure pulmonary circulation through the PS. This, in turn, can result in perfusion of the upper body and the coronary circulation largely by deoxygenated blood leading rapidly to myocardial ischaemia and circulatory arrest, as was seen in the patients with acute deterioration upon PS creation (Boudjemline et al. 2017). Interestingly, simulations for the patient in this study, despite the clinical history of well tolerated creation of the PS, have demonstrated small amount of the blood flow going retrogradely from the shunt into aortic arch. Patients from the published series of percutaneously created PS, who suffered from acute circulatory deterioration during the procedure, had retrograde flow in aortic arch clearly visible on Doppler echocardiography (Boudjemline et al. 2017), corroborating the hypothesis of circulatory arrest due to myocardial ischaemia because of catastrophic decrease of LV stroke volume. Further development of the model should address the factors potentially influencing the myocardial oxygen delivery and ventricular performance in the setting of the right-to-left shunt through the PS and arterial pressure changes in the patients with severe suprasystemic iPAH and ventricular dysfunction.

As can be seen from the areas within the PV-loops, creation of the PS causing a post-ventricular right-to-left shunt results in approximately 11% decrease in the work done by the LV while operating at higher pressures. This shows that the contribution of the smaller LV stroke volume towards decreasing the workload overpowers the contribution of higher LV afterload towards increasing it. For the RV, however, even though a reduction in afterload is achieved by the PS creation, the increase in volume returning to the ventricle overpowers and results in a net 18% increase in the RV workload, accompanied presumably by increase in myocardial oxygen consumption. These findings are consistent with the main conclusion of the previous LPM study on PS for severe PAH, where the authors also reported that the shunt does not result in unloading of the RV work despite its partial decompression (Delhaas et al. 2018). The increase of the work performed by the already dilated RV after PS creation is, thus, due to the increased stroke volume accompanied, however, by the significant decrease in afterload. Given the confirmed clinical improvement of the patients with severe suprasystemic iPAH who survived the PS creation on intermediate term (Baruteau et al. 2015), the clinical relevance of such an increased RV workload remains yet unclear: it is accompanied by an enhanced EF and lower absolute RV volumes, so of less stressed fibers, which suggest enhanced contractile function and potential for favourable remodelling. Nonetheless, to avoid circulatory collapse, increased workload for the RV myocardium after the PS creation should be met with increased oxygen delivery through the coronary circulation, sustained perfusion of which depends on LV performance and oxygenation level of the blood entering the coronary arteries. Oxygen delivery to the myocardium as yet, is not accounted for in the model. Inclusion of this dynamic interplay between myocardial perfusion with variable oxygen concentration, associated changes in cardiac mechanics, and the ability of RV to do more work within these constraints, should be further analysed to determine criteria for suitability of individual patients for the PS creation.

The reduction in LV output is accompanied by an increase in the operating volumes for the LA (Fig. 7). Effectively, both the maximum and minimum LA volumes are higher when compared to the pre-operative state. However, their difference (maximum minus minimum) is smaller, as the volume leaving LA in each cycle is smaller due to the lower volume returning from the pulmonary veins. On the contrary, for the RA, both the maximum and minimum RA volumes are lower, but their difference is larger due to the higher volume returning from the systemic veins. With an understanding of the volume changes, as described above, the changes in pressure follow an expected and predictable pattern (Figure 6).

3D flow features and effects of varying shunt diameter and length

The jet of flow coming from the shunt into the DAo impinges on its wall opposite to the shunt and creates high impulse forces and localised high pressure regions (Figs. 9 and 10), which potentially may be damaging to the tissue and over time lead to adaptive changes of the aortic wall. The relatively high velocity values in the shunt and the DAo downstream the PS lead to high wall shear stresses (Fig. 11). Such abnormally high wall shear stress may affect endothelial cell function, vascular biology, and promote development of aneurysms in the long-term (Dolan et al. 2013).

Finally, as discussed above, the streamlines in Fig. 9 (with pathlines showing similar patterns) point to some deoxygenated blood originating from the mPA and traversing through the PS to the aortic arch branches. This proportion is insignificant in the current case, where post-operative flow of oxygenated blood from the LV into aortic arch branches is higher relative to pre-operative state, and, thus it is unlikely that oxygenation of the brain is compromised. However, as mentioned above in Sect. 4.2, excessive right-to-left shunting may result in catastrophic decrease of LV contribution to the upper body perfusion allowing deoxygenated blood from the shunt traversing retrogradely to the aortic arch branches and to the coronaries in significant amounts. This is also evidenced by the streamlines in Fig. 15, where a higher right-to-left shunting achieved in a PS with larger diameters PS shows a larger proportion of streamlines traversing to the aortic arch branches. This finding of increasing right-to-left volume shunting with increasing PS diameters in this study (Table 9) is in some conflict with the previous study using the CircAdapt LPM, where authors report ‘hardly and change change in ratio for PS diameters above 7mm’ (Delhaas et al. 2018). However, the patient-specific case assessed in this study (with mmHg) aligns with the upper region (mPAP>80 mmHg) of the top-left subplot of Figure 2 in Delhaas et al. (2018), where the changes are aligned with our findings within the range of diameters considered (Table 9).

The results of the model reported here show a clear near-linear relationship between the shunt diameters and both cycle-averaged pressure-gradient across the shunt and flow-rate through the shunt (Fig. 13a). Thus, there is a clear potential and need to determine an optimum shunt diameter for an individual patient that balances pressure equalisation achieved by the PS and the extent of the right-to-left shunting determining the post-operative LV output and RV workload. This balance will depend in every particular iPAH patient on different factors, which are needed to be addressed in a further study before any recommendations regarding a ‘universal’ shunt diameter, such as a percentage of the DAo diameter, can be provided.

Figure 13b highlights the importance of shunt diameter selection to ensure that blood is pumped from the pulmonary to systemic circulations across the entire cardiac cycle, where shunt diameters larger than 8mm caused some back-flow from the DAo to the LPA during diastole. The general trend is that the amount of back-flow increases with increasing shunt diameters, although the direction of the flow across the shunt also depends on the vascular resistances, which may vary in response to different factors, such as physical or emotional exertion or inflammation. To avoid this undesirable left-to-right shunting, a valved PS may be considered (Baruteau et al. 2015).

The results of varying shunt lengths in Sect. 3.5 show that while global haemodynamics is largely unaltered with changing shunt lengths, the local flow-features are significantly more disturbed due to PS protrusion. In particular, the protrusions result in von Kármán type of vortices, within which the blood velocities are small, and hence the risk of thrombus formation may be elevated, especially with longer protrusions.

With the presented GMM, patient-specific assessment of haemodynamics is facilitated, which may be used for assessing PS suitability and choice of PS dimensions of diameter and length. However, for general conclusions about the effects of PS, a larger cohort must be modelled in future studies. Lastly, while the GMM allows for assessment of the PS’s orientation—for example, the angles between the PS and DA/LPA centrelines, PS tapering/flaring at the ends, PS curvature, etc.—and subsequent optimisation of PS geometry for optimal haemodynamics, these aspects are not evaluated in this study.

Role and utility of the stand-alone LPM

The results from the stand-alone LPM are presented and discussed in Appendix B. Here, it is highlighted that the LPM models display the same trends for the assessment of global haemodynamics as the GMM. This is not surprising as the global haemodynamics, also in the GMM, is determined primarily by the LPM components of the model and the shunt diameter, both of which are replicated in the stand-alone LPM. Indeed, the stand-alone LPM does not account for shunt length protrusion, but the results of Sect. 3.5 show that this has minimal effect on global haemodynamics.

More importantly, it should be noted that since the stand-alone LPM runs in seconds, in contrast to the days needed to obtain results from the GMM, it presents an ideal tool for parameter estimation, where large series of multiple simulations are necessary to comprehensively assess and predict the haemodynamic changes. This approach has been adopted in this study, where the estimation of parameters was largely done manually with the fast run times of the stand-alone LPM presenting, thus, an immense utility. Furthermore, since the LPM captures the trends in global haemodynamics correctly, it can be used for quick assessment of how different clinically measurable pre-operative parameters would affect global response of the circulation to PS creation, and thus, determine the suitability of an individual iPAH patient for the procedure. Also, such a predictive stand-alone LPM would provide a useful tool during clinical decision-making about optimal shunt diameter for the patient with particular set of pre-operative parameters balancing the RV partial decompression with the tolerable decrease in LV output achieved by the PS.

The GMM was useful to create the LPM components for the shunt and for the aortic part above the shunt. The later may be established based on geometrical information, while the current shunt law might be general enough for other iPAH patients. Overall, it is only when local flow features are deemed of high importance, that the GMM presents unparalleled advantages. This however, may be of less importance in the global assessment of the suitability of an individual patient for the PS creation.

It should be noted that the standalone LPM may not be very suitable for the assessment of oxygen delivery. This is because any standalone LPM is incapable of determining the split of shunted flow between the lower body (PS flow traversing downwards in DAo) and the upper body (PS flow traversing upwards to the aortic arch branches, and potentially the coronaries). To assess oxygen delivery, the standalone LPMs must assume this split, which may not be realistic. The GMM on the other hand is capable of quantifying this split and hence providing a more accurate assessment of oxygen delivery.

Conclusions

This work presents the development of a patient-specific geometrical multiscale computational model of blood flow without and with PS, and its stand-alone equivalent lumped-parameter model of the whole circulation. Model parametrization is based on clinical data obtained in a paediatric patient with suprasystemic idiopathic PAH accompanied by RV failure who received the stent-based PS. These computational models are, respectively, representative and predictive of the patient-specific haemodynamics pre- and post-creation of the PS. The results show sensitivity of the local and global haemodynamics to the PS diameter, for which an optimum value should be probably based on patient-specific factors, other than just the DAo size, affecting local and global haemodynamics. In particular, in the range of diameters considered, a near-linear relationship between the shunt diameter and pressure-drop/shunt-flow is observed. The assessment of 3D flow features shows that the jet of shunt flow impinges onto the DAo wall creating regions of high pressure and high wall shear stress, which may be detrimental to endothelial cell function and promote the development of aneurysms. PS protrusion only affects local haemodynamics, with thrombogenic flow features. Overall, this study provides insights into physiology of this complex condition and hints towards potential causes of the circulatory collapse in some patients, such as critical LV output decline and ensuing retrograde aortic flow from the RV through PS inducing myocardial ischaemia, and higher post-operative RV workload. Complemented with oxygen delivery assessment, this model may become a useful predictive tool in clinical decision-making regarding suitability for this type of palliative treatment in an individual PAH patient with drug-resistant suprasystemic PAH.

Limitations and future work

The primary limitations of this work are that only one patient has been analysed, oxygen delivery was not included in the analysis, and that regulatory/adaptive mechanisms have been ignored in modelling the haemodynamic changes. Considering the first limitation, the goal of this study was, however, to develop and validate the computational model accurately capturing the haemodynamics specific to suprasystemic PAH and predicting changes due to the PS creation. To verify whether computed representation of haemodynamics reasonably corresponds to clinical measurements, the model was fine-tuned to one representative patient-specific case. To further confirm the validity of presented approach, a larger study is required where complete pre- and post-operative data for multiple patients are acquired, with a defined research and measurement protocol, so that parameters can be automatically and uniquely estimated for each case. Less uncertain and reliable sets of measurements are key to data-driven modelling, robust calibration of the existing models, and the development of new models that account for all relevant physics. Instructively, the availability of such measurements can also be useful in refining the model physics (Yang et al. 2016). The second limitation can be addressed by incorporating models predicting local oxygen delivery beyond the generalised assumptions of upper vs. lower body saturations as reported elsewhere (Delhaas et al. 2018), into the GMM. These models, including a detailed representation of the coronary circulation, may be immensely beneficial for getting deeper insights into the PS-induced changes in the ventricular workload and function. In addition, the septal shift is not included in the model, which could be addressed using the three-wall segment model describing mechanics of ventricular interaction (Lumens et al. 2009). Although equalization of ventricular pressures will lead to a better shape of the LV post-operatively, its effect on global haemodynamics, when compared to the much stronger influence of the shifts in volumes and pressures after the PS creation is probably much smaller. The final limitation is harder to address as the effectiveness of endogenous biological compensatory mechanisms are highly individual and depends strongly on large number of cardiovascular and non-cardiac factors. In particular, a key missing component now is modelling the adaptation of RV and LV ventricular function (e.g. contractility) to changes in their output, oxygen delivery, and work performed. With such an addition to the model, it will be possible to develop second order, time-dependent logic for assessing post-operative haemodynamics. The development of such a model and its relation to patient-specific characteristics, such as level of pre-operative ventricular dysfunction, PVR, SVR, etc. forms the primary area of future investigation.

Table 11

Outlet areas in the 3D domain (see Figure 2a for labels)

Outlet	Area [mm^2]
IN	43.4
LC	25.5
LS	11.4
DAo	78.8
LPA-1	119.5
LPA-2	130.9
RPA	517.7

Table 12

Pre- and post-operative , , EDV, SV, EF and CO from the stand-alone LPM against clinical measurements

	P_AAo [mmHg]			P_mPA [mmHg]			EDV [ml]		SV [ml]		EF		CO [L/min]
	Systolic	Diastolic	Mean	Systolic	Diastolic	Mean	LV	RV	LV	RV	LV	RV	LV	RV
Pre-operative measurements and model output (parameter estimation)
Measurement	94.0	53.0	69.0	112.0	67.0	85.0	66.9	91.8	51.0	51.0	0.76	0.55	3.4	3.4
LPM	93.4	52.7	65.7	111.7	68.7	85.1	67.0	92.0	51.2	51.2	0.76	0.56	3.41	3.41
Post-operative measurements and model output (validation)
Measurement	97.0	51.0	71.0	102.0	54.0	76.0	–	–	–	–	–	–	–	–
LPM	100.3	61.6	75.2	105.8	62.3	77.6	68.5	86.7	45.1	61.8	0.66	0.71	3.00	4.12

Table 13

Stand-alone LPM: , , EDV, SV, EF and CO for different PS diameters

Diameter [mm]	P_AAo [mmHg]			P_mPA [mmHg]			EDV [ml]		SV [ml]		EF		CO [L/min]
	Systolic	Diastolic	Mean	Systolic	Diastolic	Mean	LV	RV	LV	RV	LV	RV	LV	RV
pre-op	93.4	52.7	65.7	111.7	68.7	85.1	67.0	92.0	51.2	50.5	0.76	0.55	3.42	3.42
5	98.7	60.7	73.5	110.1	62.1	79.0	68.0	88.6	46.5	59.3	0.68	0.67	3.10	3.95
6	99.4	61.3	74.4	108.5	62.1	78.2	68.1	88.1	45.9	60.3	0.67	0.69	3.06	4.02
7	100.0	61.5	75.0	106.9	62.2	77.7	68.3	87.2	45.4	61.1	0.66	0.70	3.02	4.07
7.6	100.3	61.6	75.2	105.8	62.3	77.6	68.5	86.7	45.1	61.8	0.66	0.71	3.00	4.12
8	100.5	61.7	75.3	105.1	62.3	77.4	68.6	86.1	44.8	61.7	0.65	0.72	2.99	4.12
9	101.1	61.8	75.6	103.1	62.4	77.1	68.8	85.0	44.3	62.2	0.64	0.73	2.95	4.15
10	101.6	61.9	75.8	101.1	62.4	76.8	69.0	84.1	43.8	62.6	0.64	0.74	2.92	4.17

Table 14

Stand-alone LPM: work done (area under the PV-loop) by the left ventricle and the right ventricle for varying shunt diameters

D_SH	W_LV	W_RV	Q¯_SH	V¯_SH	ΔP¯_SH	Q_p/Q_s
[mm]	×10^3[mmHg-ml]	×10^3[mmHg-ml]	[ml/s]	[ml]	[mmHg]	[-]
Pre	4.13	5.01	–	–	–	1
5	3.94	5.84	14.7	13.2	9.48	0.78
6	3.91	5.91	16.6	15.0	7.48	0.76
7	3.87	5.95	18.0	16.2	6.26	0.74
7.6	3.85	5.97	18.6	16.7	5.73	0.73
8	3.83	5.97	19.2	17.3	5.40	0.72
9	3.82	5.97	20.3	18.3	4.63	0.71
10	3.78	5.95	21.3	19.2	3.94	0.70