In this study, we used an image-based nonlinear one-dimensional CFD model to predict arterial hemodynamic features that would be difficult, if not impossible, to measure experimentally and which may be clinically relevant to disease progression. Based on an existing set of imaging data pre- and post-CTEPH, we developed subject-specific 1D arterial network models. We used an efficient multi target optimization scheme for calibration that incorporated experimentally measured MPA systolic and diastolic pressures and areas, as well as dynamic flow data at the LPA and RPA. Our approach yielded simulation results with excellent agreement with measured data, and predicted pulmonary arterial stiffness, TAWSS and OSI distribution pre- and post-CTEPH. The increase in pulmonary arterial stiffness post-CTEPH will drive upstream (i.e., right ventricular (Liu et al. 2022)) and downstream (i.e., pulmonary capillary (Wang and Chesler 2011)) remodeling to accelerate disease progression. Moreover, changes in WSS are relevant to the growth or dissolution of thrombi (Sakariassen et al. 2015). Thus, these model results allow for a more comprehensive understanding of the effects of CTEPH, enhancing our ability to study this condition and potentially improve patient care and treatment strategies.
Modeling and Optimization
Hemodynamics modeling is a powerful tool for investigating hypothesized mechanisms of pulmonary vascular disease development and progression. There has been limited application of these tools to CTEPH in comparison to other forms of pulmonary hypertension, such as pulmonary arterial hypertension (PAH) (Clipp and Steele 2009; Zambrano et al. 2018; Dong et al. 2021). The use of subject-specific modeling has the potential to provide clinicians with a better understanding of patient specific disease states, optimizing clinical decision-making and ultimately enabling personalized medicine. Several computational studies calibrated their models using existing experimental data. Zambrano et al. (2018) developed subject-specific computational models of pulmonary vascular hemodynamics for both a PAH patient and a healthy individual. The parameters of the model, including Windkessel parameters (\({R}_{p},{R}_{d},\) and \(C\)) for each terminal vessel in the outflow branches, were adjusted using in vivo measurements of pressure from RHC. Additionally, the linearized stiffness E of the large arteries was adjusted based on the relationship between pressure and diameter from MRI. Clipp and Steele (2009) generated 1D CFD models of the pulmonary arterial networks of lambs, incorporating considerations such as geometry, compliance, structured tree boundary conditions, and respiration. Similar to Zambrano et al. (2018), parameters describing the boundary conditions in the Clipp et al. study were inferred by comparing the computed inlet pressure waveform to experimental data. In this study, we not only used pressure calibration data (systolic and diastolic) to infer boundary conditions but also used area (systolic and diastolic) and flow (time-series) calibration data, enhancing model accuracy and reliability. As shown in Fig. 4 for a representative animal and consistent with all animals, model predictions quantitatively matched measured pressures and flows (high 𝑅2 values).
As mentioned earlier, there are few studies that used computational modeling to predict hemodynamics in CTEPH. Colebank et al. (2021) developed 1D CFD models of pulmonary arterial networks based on computed tomography imaging in CTEPH patients. They simulated four different CTEPH disease cases to test physiological hypotheses related to remodeling. Tsubata et al. (2023) performed 3D CFD simulations to model subject-specific data from healthy subjects, CTEPH patients and CTEPH patients treated with balloon pulmonary angioplasty. However, neither paper had an abundance of pre- and post-CTEPH data from individual subjects. A novel aspect of our study is that we used paired data sets pre- and post-CTEPH, which enables a subject-specific interpretation of CTEPH progression.
Pulmonary Arterial Parameters
Pulmonary hypertension is known to increase pulmonary arterial stiffness, which is measured in clinical and preclinical studies using a wide range of direct and indirect methods. Su et al(2017) used wave intensity analysis to indirectly measure pulmonary arterial stiffness in 10 healthy subjects and 21 subjects with pulmonary hypertension (11 with PAH and 10 with CTEPH). They reported that the wave reflection index in these subjects was significantly higher than in healthy controls. We analytically determined the mechanical properties of the PAs using systolic and diastolic pressure and a linear wall mechanics model. In contrast, previous studies(Qureshi et al. 2019) required parameter estimation to determine pulmonary arterial stiffness. As shown in Fig. 5a and consistent with the literature, we observed vessel stiffening after CTEPH induction, which may contribute to disease progression (Sun and Chan 2018).
Pulmonary arterial stiffening in CTEPH is likely a combination of load-dependent stiffening and structural stiffening (Wang and Chesler 2011; Pewowaruk et al. 2021). That is, the increase in pulmonary vascular resistance due to CTEPH increases pressure and thus increases the operating point of arteries on the non-linear pressure-diameter curve (Wang and Chesler 2011), called load-dependent stiffening. In addition, proliferative and fibrotic changes in the pulmonary arterial wall can cause structural remodeling that increases vessel stiffness (Kim et al.). The combination of increased resistance and stiffness elevates right ventricular afterload and contributes to the development of right heart failure(Wang and Chesler 2011). Qureshi et al. (2019) used a model similar to ours to study hemodynamic changes in normoxic and hypoxic mice. Their Windkessel resistance increased in hypoxa while Windkessel compliance decreased. Figure 5b-c illustrate an increase in Windkessel proximal and distal resistance with CTEPH, primarily attributable to the obstruction and narrowing of pulmonary arteries and arterioles, respectively. Additionally, as depicted in Fig. 5d, a decrease in Windkessel compliance occurs with CTEPH, likely due to a combination of load-dependent and structural stiffening in the distal arteries and arterioles.
Flow and shear stress
Changes in the flow of blood within the lungs are believed to play a crucial role in the onset and worsening of pulmonary vascular diseases (Clark and Tawhai 2018). Here, changes in flow magnitude and distribution were observed in both lungs for all subjects from baseline to CTEPH. The change in distribution with CTEPH, in particular the left-right asymmetry, is due to uneven microspheres distribution (Mulchrone et al. 2019), which may not reflect the clinical disease state.
Due to these flow changes as well as artery dilation, our 1D model predicts WSS decreases throughout the arterial network with CTEPH. Clinical studies that use image-based methods to compute WSS in individuals have shown that those with PH exhibit lower WSS than normotensive subjects (Barker et al. 2015; Terada et al. 2016). Modeling studies have found the same result. Zambrano et al. (2021) performed 3D image-based, computational fluid dynamics simulations, in which interactions between blood flow and wall deformation were included and predicted that subjects with PAH have lower TAWSS than controls. Tang et al. (2011) used the same approach as Zambrano et al. (2021) and similarly predicted lower mean WSS in the proximal arteries in the PAH group. From patient-specific 3D computational fluid dynamics simulations in CTEPH, Tsubata et al. (2023) observed lower systolic WSS in CTEPH patients in comparison to controls as well as to patients after balloon pulmonary angioplasty treatment. Our results are consistent with the results of these image-based and computational studies, as shown in Fig. 6.
At the MPA, we found that TAWSS significantly decreased following CTEPH in four of the five subjects; in the one exception, that we attribute this to poor image quality at baseline that likely increased the estimated arterial network area. A significant reduction in LPA TAWSS was observed for all subjects in our study. The observed decrease in TAWSS can be attributed to a dramatic drop in flow to the left lung due to the severity of the obstruction. This finding is consistent with the suggestion that WSS reduction in the proximal PAs is an indicator of PAH disease severity (Yang et al. 2019). At the RPA, TAWSS decreased in all subjects except one (not the same subject for which MPA TAWSS increased). Overall, the decrease in TAWSS in the RPA was less dramatic and more variable than in the LPA (Fig. 6).
In the intralobar arteries, TAWSS is higher compared to the extralobar arteries in both baseline and CTEPH conditions. This observation is consistent with Pillalamarri et al(2021) who found that TAWSS increased with increasing distance from the heart based on a computational investigation of patient-specific pulmonary hemodynamics across 32 subjects in distinct WHO groups of PH. The Pillalamarri et al. (Pillalamarri et al. 2021)study included one subject with CTEPH, in whom TAWSS decreased in the distal part of the network, consistent with our results.
The biological consequences of decreased WSS in the pulmonary circulation include altered endothelial cell signaling that drives vasoconstriction and the production of pro-inflammatory factors (Allen et al. 2023). Decreased TAWSS may also contribute to the synthesis and accumulation of collagen in the arterial wall and proliferation of smooth muscle cells (Ryu et al. 2022). This pulmonary arterial remodeling can result in a stiffer and less compliant pulmonary artery network (Friesen et al. 2019), which further exacerbates the severity of PH and can lead to right heart failure.
Oscillatory shear index
Like TAWSS, OSI can drive changes in endothelial cell structure and function that affect disease progression. In endothelial cells, high OSI is associated with increased reactive oxygen species, which can drive vasoconstriction and inflammation (Peiffer et al. 2013). In the large extrapulmonary arteries, especially the MPA and RPA, OSI magnitude was quite small at baseline and only became somewhat larger with CTEPH. In contrast, Tsubata et al. (2023) found a significant increase in proximal artery OSI in CTEPH subjects compared to controls. Their OSI magnitudes were on average much larger than ours (0.05 verses 0.35). Bartolo et al. (2022) used a 1D framework similar to ours to compute OSI in a human pulmonary arteriovenous network and found OSI magnitudes similar to ours, which suggests that the 1D framework and/or the assumed velocity profile affects OSI calculations.
In the intralobar arteries, there was no clear pattern regarding the impact of CTEPH on OSI (Fig. 8d). To further investigate temporal abnormalities in intralobar flow dynamics, we considered the combination of TAWSS and OSI (Fig. 9) represented by the dimensionless parameter \(\phi\). Tsubata et al. (2023) suggested that the combination of decreased WSS and increased OSI can contribute to thrombogenicity. Our results indicate that the percentage of vessels that had low TAWSS (< 5dyne/cm2) and high OSI (> 0.05) increased from baseline to CTEPH in all subjects. When the vascular endothelium experiences a decrease in wall shear stress and an increase in OSI, it inhibits the production and release of nitric oxide, a potent vasodilator (Terada et al. 2016) As a result, the abnormal flow dynamics as seen in our study may be linked to impaired vasodilation and vasoconstriction in the pulmonary arterioles.
Limitations
Several limitations exist in our study, notably the small sample size, which contributed to heterogeneity. Limited image resolution in some subjects, especially one subject at baseline and a different subject after CTEPH, resulted in challenges extracting accurate subject-specific arterial networks and changes in those arterial networks with CTEPH. Since time-series pressure data were not recorded for either baseline or CTEPH conditions, we calibrated our model based on the scalar values of systolic and diastolic pressures (as well as scalar values of systolic and diastolic MPA area and time series LPA and RPA flow). The inclusion of time-series MPA pressure data would increase model fidelity. Finally, we assumed Blunt-like flow with a power law radial distribution for all arteries and held the power fixed at g = 9. However, the range of Womersely number in the models was 2 to 40. As a consequence, in the small arteries, which will have lower Womersely numbers, TAWSS is overestimated. The impact on OSI is difficult to predict a priori; comparison to Tsubata et al. (2023) suggests OSI is underestimated in our model. Adapting g to local diameter or using an explicit Womersley profile in in the 1D CFD model would increase model complexity and solution time but increase accuracy.