The Right Heart Network and Risk Stratification in Pulmonary Arterial Hypertension.

BACKGROUND: Prognosis in pulmonary arterial hypertension (PAH) is closely related to indexes of right ventricular function. A better understanding of their relationship may provide important implications for risk stratification in PAH. RESEARCH QUESTION: Can clinical network graphs inform risk stratification in PAH? STUDY DESIGN AND METHODS: The study cohort consisted of 231 patients with PAH followed up for a median of 7.1 years. An undirected, correlation network was used to visualize the relationship between clinical features in PAH. This network was enriched for right heart parameters and included N-terminal pro-hormone B-type natriuretic peptide (NT-proBNP), comprehensive echocardiographic parameters, and hemodynamics, as well as 6-min walk distance (6MWD), vital signs, laboratory data, and diffusing capacity for carbon monoxide (Dlco). Connectivity was assessed by using eigenvector and betweenness centrality to reflect global and regional connectivity, respectively. Cox proportional hazards regression was used to model event-free survival for the combined end point of death or lung transplantation.
RESULTS: A network of closely intertwined features centered around NT-proBNP with 6MWD emerging as a secondary hub were identified. Less connected nodes included Dlco, systolic BP, albumin, and sodium. Over the follow-up period, death or transplantation occurred in 92 patients (39.8%). A strong prognostic model was achieved with a Harrell's C-index of 0.81 (0.77-0.85) when combining central right heart features (NT-proBNP and right ventricular end-systolic remodeling index) with 6MWD and less connected nodes (Dlco, systolic BP, albumin, sodium, sex, connective tissue disease etiology, and prostanoid therapy). When added to the baseline risk model, serial change in NT-proBNP significantly improved outcome prediction at 5 years (increase in C-statistic of 0.071 ± 0.024; P = .003).
INTERPRETATION: NT-proBNP emerged as a central hub in the intertwined PAH network. Connectivity analysis provides explainability for feature selection and combination in outcome models.

FOR EDITORIAL COMMENT, SEE PAGE 1138

Study Question: Can network graphs inform risk stratification models in PAH?

Results: NT-proBNP emerged as a central hub of a highly intertwined network in PAH; network graphs and centrality analysis provide explainability for risk stratification models in PAH.

Interpretation: Risk prediction can best be accomplished by combining a small number of central right heart features with functional status and markers of end-organ dysfunction; moreover, serial changes of the central hub further calibrate risk prediction.

Survival in pulmonary arterial hypertension (PAH) is closely associated with right heart failure.1, 2, 3, 4, 5, 6 In the multicenter National Institutes of Health registry of D’Alonzo et al, right heart failure accounted for more than two-thirds of deaths, and both right atrial pressure and cardiac index emerged as strong predictors of outcome. Since then, large registry studies have refined risk prediction in PAH by incorporating B-type natriuretic peptide (BNP) or its N-terminal pro-hormone (NT-proBNP), functional class, and other selected clinical or laboratory data.2, 3, 4, 5, 6, 7 These registries include the French Pulmonary Hypertension Network (FPHN), the Registry to Evaluate Early and Long-term Pulmonary Arterial Hypertension Disease Management (REVEAL), and the Pulmonary Hypertension Outcomes Risk Assessment (PHORA) prediction models.2, 3, 4, 5, 6

In parallel to these registry-based scores, much attention has been focused on identified imaging predictors of outcome. Several right heart metrics have been considered, including tricuspid annular plane systolic excursion (TAPSE), right ventricular (RV) longitudinal strain (RVLS), RV ejection fraction, RV end-systolic dimension or remodeling indexes, and, more recently, right atrial volume or right atrial reservoir function.9, 10, 11, 12, 13 Interpreting imaging studies in PAH is sometimes challenging in view of the multiplicity of indexes and variable coverage among studies (Fig 1A). In fact, many studies focus on a specific right heart metric with little consideration given to the strong associations that exist between them. A better understanding of these relationships may provide valuable insights for risk stratification models. In recent years, graph theory has emerged as a powerful tool to model pairwise relationships between nodes or features. Network graphs are commonly used in social sciences or molecular biology and increasingly used in clinical medicine., They provide a method for visualizing data architecture and identifying central features as well as potential complementarity between parameters (Fig 1B). For example, a study by Oldham et al used network graphs to identify key features of invasive cardiopulmonary exercise testing, which were then used to identify clinically relevant clusters.

Figure 1

Context of the study. A, Coverage of right sided heart metrics varies among imaging studies. Network graphs potentially help identify key features for actionable prediction models or cluster analysis. B, Network graphs visualize and quantify data connectivity through centrality analysis and network topology.

The main objective of the current study was to model data architecture in PAH by using network graphs. We used a dataset that was enriched for right heart parameters as well as comprehensive clinical and laboratory data. Our first objective was to use network graphs to visualize and quantify feature connectivity in PAH. Our second objective was to determine by what extent connectivity is related to all-cause mortality or lung transplantation in PAH. The final objective was to determine how connectivity in PAH informs risk stratification models.

Study Design and Methods

Study Population

Consecutive patients with PAH followed up at Stanford University between November 2008 and March 2014 were considered for inclusion in the study. The diagnosis of PAH was based on a resting mean pulmonary arterial pressure (MPAP) ≥ 25 mm Hg, pulmonary arterial wedge pressure ≤ 15 mm Hg, and pulmonary vascular resistance (PVR) > 240 dynes⋅s/cm5 (3 Wood units) based on the definition at the time of enrollment in the registry. An echocardiogram was required for inclusion in the study as well as an NT-proBNP laboratory value. Patients with congenital heart disease (n = 14), porto-pulmonary hypertension (n = 18), or with a clinical diagnosis of moderate to severe chronic obstructive or restrictive lung disease (n = 12) were excluded. The final cohort consisted of 231 patients, with the majority of patients having prevalent PAH (n = 175 [76%]). The study was approved by the Stanford Institutional Review Board (IRB #14083) and conducted in agreement with the Declaration of Helsinki seventh revision (2013).

Clinical Data Collection

Clinical data included demographic characteristics, anthropomorphic measures, etiology of PAH, right heart catheterization data, echocardiography data, NT-proBNP, comprehensive metabolic panel and complete blood count, 6-min walk distance (6MWD), and diffusing capacity of carbon monoxide (Dlco). The primary outcome consisted of all-cause mortality or lung transplantation; outcomes were verified by using chart review and the United States Social Security Death Index. Because the current study was enriched for echocardiographic parameters, the date of echocardiography was chosen as the reference date for outcome analysis. We also calculated the REVEAL 1.0 score, the French Pulmonary Hypertension Network (FPHN), and the World Symposium of Pulmonary Hypertension score.3, 4, 5

As part of the clinical registry, the majority of the data points were available in proximity to the echocardiogram. In 90% of patients, right heart catheterization was obtained within 4 months of the index echocardiography. For NT-proBNP, 61% of patients had an NT-proBNP value within 1 week and 78% within 1 month; of the remaining patients, 36 had stable longitudinal NT-proBNP values < 300 ng/L; the other 15 patients with elevated values had an NT-proBNP value that remained in the same risk category. The majority of laboratory data were obtained within 1 month of the echocardiography. Only a few variables required imputations (eg, red blood cell distribution [n = 2], 6MWD [n = 4], and Dlco [n = 27]).

Quantitative Echocardiography

Studies were acquired by using a Philips iE33 ultrasound system and analyzed according to the American Society of Echocardiography recommendations by a cardiologist blinded to the clinical outcome (M. A.). A focused apical four-chamber right ventricular view was used to measure parameters of right atrial and right ventricular size and function. Measures included RV end-systolic remodeling index (RVESRI), RV longitudinal strain (RVLS), RV fractional area change, TAPSE by a two-dimensional method, RV areas including relative RV to left ventricular area, right atrial volume, left ventricular eccentricity index, left ventricular ejection fraction, estimated RV systolic pressure (RVSP), and estimated right atrial pressure (eRAP). RVESRI was measured as previously described as the ratio of RV lateral length to septal length, and it represents an internally scaled end-systolic dimension parameter. The eRAP was determined according to the American Society of Echocardiography recommendations but we also assigned an eRAP of 20 mmHg in the presence of a dilated inferior vena cava with minimal collapse index (< 20%).

Right Heart Catheterization

Using an internal jugular or femoral vein approach, RAP, pulmonary arterial pressures (systolic, mean [MPAP], and diastolic pressure) and pulmonary capillary wedge pressure were measured. Cardiac output was derived by using the indirect Fick method and assumed oxygen consumption. PVR was indexed to body surface area, and relative PVR (PVRr) was measured as the ratio of PVR and systemic vascular resistance.

Statistical Methods

Baseline characteristics are presented as median and interquartile range for continuous variables, and categorical variables are presented as number and percentage. To construct the correlation network, pairwise Spearman rank correlations were calculated by using the R package Hmisc (version 4.4-0), and undirected networks were plotted with igraph (version 1.2.5). Only pairwise correlations with Bonferroni-adjusted P values < .05 were included and displayed via the Fruchterman-Reingold method (supplementary document additional statistical explanation section):. We selected the Spearman correlation rather than the Pearson correlation because it is better suited to model monotonic nonlinear relationships. Missing values were imputed for Dlco by using the k-nearest neighbors’ method (impute.knn function) in the R package impute (version 1.60.0) and verified that imputation did not change centrality measures (e-Fig 1).

Centrality measures the importance of a node in a network in terms of connectivity. We used two commonly used measures of centrality that reflect both global and regional connectivity; that is, eigenvector and betweenness centrality, both calculated by using igraph., Eigenvector centrality measures a node’s importance by considering both the direct connection to the node as well as the number of connections to its neighbors (connected nodes); as such, it measures the “reach” of the overall connections in a network. Because the network is enriched for cardiopulmonary unit metrics, eigenvector centrality will measure overall cardiopulmonary connections. It has the advantage of representing an overall cardiopulmonary centrality measure without preferentially selecting a specific metric. Betweenness centrality, conversely, quantifies the relative number of shortest paths passing through a specific node; as such, it will highlight nodes that connect directly to other nodes, similar to a hub between communication towers. A node that has both the highest eigenvector and betweenness centrality will usually occupy a central topology in a simple network. Nodes were color-coded by clinical domains, with size reflecting betweenness centrality. e-Table 1 summarizes the parameters included in the correlation network. Focused correlation matrices of highly connected variables are presented.

Cox proportional hazards regression was used for univariable and multivariable outcome analyses. A stepwise approach was used to develop a multivariable model. We first determined independent predictors from each cluster grouped according to eigenvector centrality tertiles. Parameters that were highly co-linear were excluded; for highly co-linear parameters, the choice of parameter was guided by its hazard ratio as well as by literature review (e-Table 2). We then assessed their incremental value in a multivariable analysis. Age, sex, connective tissue disease (CTD), and prostanoid therapy, as well as sodium, were also included for their known association with outcome in PAH. We used categorical or ordinal data because these are easier to translate in clinical practice; thresholds were based on previously studies or accepted laboratory reference values.,, Kaplan-Meier analysis was based on the Cox prognostic index. Alternate models were considered to account for potentially missing data in an imaging-focused model (without NT-proBNP or Dlco) and a biomarker-focused model (without imaging or Dlco). Although predefined thresholds were used, we also performed receiver-operating characteristic analysis with thresholds identified based on the Youden index as well as fixed 80% sensitivity and specificity.

For longitudinal analysis, we focused on the most central node of the network and evaluated its complementarity value to the baseline model. P values < .05 were considered statistically significant. Analyses were performed by using R packages Hmisc (version 4.4-0) and igraph (version 1.2.5) and MedCalc Statistical Software version 19.6 (MedCalc Software Ltd.).

Results

Of the 231 patients, 182 (79%) were female, and 73 (32%) had CTD. Median age was 48 years (38.6-57.3), MPAP was 50 ± 16 mm Hg, PVR indexed on body surface area was 19.2 (12.5-26.8) Wood units ∗ m2, and the median NT-proBNP was 496 ng/L (120-1,507) (Table 1). Of the patients with CTD, 41 (56%) had scleroderma, 16 (22%) had systemic lupus erythematosus, 12 (16%) patients had mixed CTD or vasculitis, and five patients had rheumatoid arthritis (7%). No patients were in cardiogenic shock at the time of their echocardiography, although six patients were hypotensive (systolic BP < 90 mm Hg and cardiac index < 2 L/min/m2 at the time of their right heart catheterization).

Table 1

Baseline Patient Characteristics

Characteristic	Value
Sample size	231
Demographic characteristics
Age, y	48.0 [38.6 to 57.3]
Female sex	182 (78.8%)
BMI, kg/m2	26.7 [29.9 to 31.6]
Connective tissue disease–associated PAH	73 (31.6%)
Vital signs
Systolic BP, mm Hg	113 [106 to 126]
Diastolic BP, mm Hg	70 [63 to 77]
Heart rate, beats/min	81.0 [71.0 to 88.7]
Right heart catheterization
Mean pulmonary arterial pressure, mm Hg	50.1 ± 16.1
Right atrial pressure, mm Hg	7.0 [5.0 to 11.8]
Pulmonary artery wedge pressure, mm Hg	10.0 [8.0 to 13.0]
Cardiac index, L/min/m2	2.0 [1.67 to 2.35]
Pulmonary vascular resistance index, WU ∗ m2	19.2 [12.5 to 26.8]
Relative pulmonary vascular resistance, %	48.5 [33.3 to 0.61.5]
Echocardiographic data (selected)
Relative RV end-systolic area	1.06 [0.67 to 1.51]
RV end-systolic remodeling index	1.47 [1.34 to 1.60]
RV lateral longitudinal strain, %	–16.3% [–20.1 to –12.9]
Severe tricuspid regurgitation (3 to ≥ 4)	51 (22%)
Estimated RAP (eRAP)	3 mm Hg (15%)8 mm Hg (54%)≥ 15 mm Hg (31%)20 mm Hg (9%)
Pericardial effusion (> 0.5 cm)	57 (25%)
Left ventricular ejection fraction (%)	60.7 [56.0 to 66.2]
Functional status
NYHA functional class
I	13 (5.6%)
II	91 (39.4%)
III	107 (46.3%)
IV	20 (8.7%)
6MWD, m	429.8 [315.6 to 516.0]
Diffusing capacity of the lungs for carbon monoxide corrected for hemoglobin (% predicted)	74.0% [60.0% to 87.0%]
Laboratory data
NT-proBNP, ng/L	407 [100 to 1,300]
eGFR, mL/1.73 m2	62.1 [49.1 to 74.1]
Albumin, g/L	3.8 [3.4 to 4.0]
Sodium, mEq/L	138 [136 to 140]
Hemoglobin, g/L	14.0 [12.3 to 15.1]
RBC distribution width, %	14.9 [13.5 to 16.7]
HbA1C > 6.5%, %	14 (6.1%)
Comorbidities
Obesity (BMI > 30 kg/m2)	78 (33.8%)
Diabetes mellitus	14 (6.1%)
Systemic hypertension	13 (5.6%)
Chronic kidney disease (eGFR < 45 mL/min/1.73 m2)	33 (14.3%)
Therapeutic consideration
Incident cases (no PAH-specific therapy initiated)	55 (23.8%)
On prostanoid therapy	82 (35.5%)
Endothelin receptor antagonist	68 (29.4%)
Phosphodiesterase-5 inhibitors	114 (49.4%)
Spironolactone	63 (27.2%)
Anticoagulation	44 (19.0%)

Data are presented as mean ± SD if normally distributed and as median [interquartile range] if not. 6MWD = 6-min walk distance; eGFR = estimated glomerular filtration rate; HbA1C = glycosylated hemoglobin; NT-proBNP = N-terminal pro-B type natriuretic peptide; NYHA = New York Heart Association; PAH = pulmonary arterial hypertension; RAP = right atrial pressure; RV = right ventricular; WU = Wood unit.

Clinical Network in PAH

We constructed an undirected network of clinical features in PAH (Fig 2). The network was enriched for echocardiographic parameters and included a total of 43 features (e-Table 1). NT-proBNP emerged as a central hub based on both eigenvector and betweenness centrality. Several echocardiographic parameters, including RVLS, RVESRI, and right atrial volume index, also had high eigenvector centrality. Among features with intermediate eigenvector centrality (second tertile), most notable were RAP, 6MWD, New York Heart Association functional class, and total bilirubin. Several other features had lower eigenvector centrality, including vital signs, markers of end-organ dysfunction (glomerular filtration rate [GFR], sodium, albumin, total bilirubin, hemoglobin, RBC distribution width), Dlco, and CTD etiology, as well as age and sex. In terms of topology, a shorter path length was observed between 6MWD, estimated GFR, and RBC distribution width and eRAP than with RV structural or functional indexes.

Figure 2

A-C, Right heart failure domains and correlation network in PAH. A, Clinical features in PAH organized according to physiological domains. B, Correlation network in PAH. Size corresponds to betweenness centrality, color to clinical domains, and color of the edges according to the positive (red) or negative (blue) correlation. C, eigenvector and betweenness centrality for the different clinical features showing that NT-proBNP is the most central feature by both metrics. 6MWD = 6-min walk distance; Alb = albumin; ALP = alkaline phosphatase; ALT = alanine transferase; CTD = connective tissue disease etiology; Dlco = diffusing lung capacity of carbon monoxide; EI = left ventricular eccentricity index in end-systole; eGFR = estimated glomerular filtration rate (according to Modification of Diet in Renal Disease formula); eRAP = estimated right atrial pressure; Hb = hemoglobin; HbA1C = glycosylated hemoglobin; HR = heart rate; LAAI = left atrial area indexed; LVEF = left ventricular ejection fraction; MPAP = mean pulmonary arterial pressure; Na = sodium; NT-proBNP (or NT-BNP) = N-terminal pro-B type natriuretic peptide; NYHA = New York Heart Association functional class; PAH = pulmonary arterial hypertension; Peff = pericardial effusion; PFT = pulmonary function test; PLT = platelet count; PVRI = pulmonary vascular resistance indexed to body surface area; PVRr = pulmonary vascular resistance divided by the systemic vascular resistance; RA = right atrial; RAAr = right atrial area divided by left atrial area; RAP = right atrial pressure; RAVI = right atrial volume index; RDW = RBC distribution width; RV = right ventricular; RVEDAI = right ventricular end-diastolic area index; RVESAI = right ventricular end-systolic area index; RVESRI = right ventricular end-systolic remodeling index; RVESAr = right ventricular end-systolic area divided by left ventricular end-systolic area; RVFAC = right ventricular fractional area change; RVLS = right ventricular longitudinal strain (absolute value); RVSP = right ventricular systolic pressure based on the tricuspid regurgitation Doppler signal; SBP = systolic BP; SVRI = systemic vascular resistance indexed on body surface area; TAPSE = tricuspid annular plane systolic excursion; Tbili = total bilirubin; TR = tricuspid regurgitation.

As expected, the network region surrounding NT-proBNP was densely populated with RV echocardiography parameters. A strong relationship was found between parameters of right heart size or function (Figs 3A, 3B). The relationship between two-dimensional measured TAPSE and other functional metrics was weaker, however. A nonlinear relationship was noted between NT-proBNP (log scale) and RVESRI (rho = 0.59; P < .001) as illustrated in Figure 3C. Among patients with eRAP ≥ 15 mm Hg (n = 70), only three (4.2%) had normal RVESRI and NT-proBNP. On multivariable regression analysis, NT-proBNP was independently associated with RVESRI (r partial = 0.38), eRAP (r partial = 0.16), tricuspid regurgitation (TR) severity grade (r partial = 0.35) and estimated GFR (r partial = –0.19) with an R2 of 0.46 (P < .001).

Figure 3

A-C, Association between parameters of PAH. A, Spearman correlation between right heart size and remodeling metrics. B, Spearman correlation between functional metrics as well as RVESRI. C, Nonlinear monotonic relationship between NT-proBNP (log scale) and RVESRI with overlay of elevated RAP estimated by echocardiography. eRAP = estimated right atrial pressure; NT-proBNP (or NT-BNP) = N-terminal pro-B type natriuretic peptide; PAH = pulmonary arterial hypertension; RAAr = right atrial area divided by left atrial area; RAP = right atrial pressure; RAVI = right atrial volume index; RVEDAI = right ventricular end-diastolic area index; RVESAI = right ventricular end-systolic area index; RVESRI = right ventricular end-systolic remodeling index; RVESAr = right ventricular end-systolic area divided by LV end-systolic area; RVFAC = right ventricular fractional area change; RVLS = right ventricular longitudinal strain (∗ in absolute value); TAPSE = tricuspid annular plane systolic excursion.

Feature Connectivity and Outcome

During a median follow-up time of 7.1 years (2.9-9.6 years), the combined outcome of death or lung transplantation occurred in 92 patients (39.8%) and included 66 deaths and 26 lung transplantations. The 1- and 5-year event-free survival was 89% and 72%, respectively (Fig 4A). Hazard ratios were stable over the course of the study, with a strong association (r = 0.99) between 5-year and longer term outcome (Fig 4B). Of the 45 potential predictors, 33 were significantly associated with event-free survival on univariable analysis (e-Table 2). Among the echocardiographic parameters, RVESRI had the nominally highest absolute log2 transformed hazard ratio followed by TAPSE (two-dimensional)/RVSP ratio and right heart structural and functional parameters. Hazard ratios (absolute log2 transformed) were related to eigenvector centrality (rho = 0.42; P = .006) but not betweenness centrality. This relationship was, however, far from linear (Fig 4C). For example, some features with high eigenvector centrality such as RVSP or PVRr had lower hazard ratios, whereas other features with low eigenvector centrality (eg, markers of end-organ dysfunction) had a higher hazard ratio.

Figure 4

A-C, Relationship between network connectivity and outcome in PAH. A, Transplant-free survival during the entire follow-up period: every patient enrolled in the registry was followed up for up to 5 years. B, Stability in hazard ratio over time with strong correlation of hazard ratio at 5 years and during long-term follow-up. C, Eigenvector centrality (decreasing order) and heat map of hazard ratio normalized to the highest standardized value of NT-proBNP of 2.1. 6MWD = 6-min walk distance; Alb = albumin; ALP = alkaline phosphatase; ALT = alanine transferase; CTD = connective tissue disease etiology; Dlco = diffusing lung capacity of carbon monoxide; EI = left ventricular eccentricity index in end-systole; eGFR = estimated glomerular filtration rate (according to Modification of Diet in Renal Disease formula); eRAP = estimated right atrial pressure; Hb = hemoglobin; HbA1C = glycosylated hemoglobin; HR = heart rate; LAAI = left atrial area indexed; LVEF = left ventricular ejection fraction; MPAP = mean pulmonary arterial pressure; Na = sodium; NT-proBNP (or NT-BNP) = N-terminal pro-B type natriuretic peptide; NYHA = New York Heart Association functional class; PAH = pulmonary arterial hypertension; Peff = pericardial effusion; PLT = platelet count; PVRI = pulmonary vascular resistance indexed on body surface area; PVRr = pulmonary vascular resistance divided by the systemic vascular resistance; RAAr = right atrial area divided by left atrial area; RAP = right atrial pressure; RAVI = right atrial volume index; RDW = RBC distribution width; RVEDAI = right ventricular end-diastolic area indexed; RVESAI = right ventricular end-systolic area indexed; RVESRI = right ventricular end-systolic remodeling indexed; RVESAr = right ventricular end-systolic area divided by LV end-systolic area; RVFAC = right ventricular fractional area change; RVLS = right ventricular longitudinal strain; RVSP = right ventricular systolic pressure based on the TR Doppler signal; SBP = systolic BP; SVRI = systemic vascular resistance indexed on body surface area; TAPSE = tricuspid annular plane systolic excursion; Tbili = total bilirubin; TR = tricuspid regurgitation.

Multivariable Cox Regression Model

To develop a multivariable outcome model, we used a stepwise approach according to centrality-based clusters (Table 2). Factors that emerged as independent in each cluster were then included in the final model. In selecting the parameters to consider, we avoided highly co-linear variables. For RV dimension, RVESRI was chosen because it is internally scaled, minimizing potential pseudo-normalization with obesity or fluid overload states (e-Table 2). The final model (Fig 5A) combined features with high (NT-proBNP and RVESRI), intermediate (6MWD), and lower (Dlco, systolic BP, albumin, sodium, sex, and CTD etiology) eigenvector centrality. This model provided good discrimination of outcome, with a Harrell’s C-index of 0.81 (0.77-0.85). Five risk strata were generated based on the normally distributed prognostic index (weighted sum of coefficients) and provided good separation of outcome as shown by the Kaplan-Meier curves (Fig 5B). e-Figure 2 shows Kaplan-Meier curves of the validated risk score. Sensitivity analysis was performed in patients with prevalent PAH and patients not on prostanoid therapy, which produced a Harrell’s C-index of 0.84 (0.80-0.88) and 0.75 (0.68-0.81), respectively. We also considered alternate models depending on available parameters; that is, an imaging-focused model (no NT-proBNP or Dlco) and a biomarker-focused model (no echocardiographic parameters or Dlco) (Table 3). Although the χ2 of these models was lower, the C-statistic for the 5-year outcome was not statistically different (P = .10 vs imaging-focused model and P = .29 vs laboratory-focused model).

Table 2

Multivariable Cox Regression Model Based on Eigenvector Centrality

Cluster with	Parameter		Cluster-Based Analysis	Multivariable Model
Higher centrality	NT-proBNPRVESRIPVR/SVR ratio	< 300 ng/L[300-1,400] ng/dL> 1,400 ng/L[1.35-1.70]> 1.70> 2/3	2.26 [1.23-4.15]2.94 [1.54-5.56]2.51 [1.14-5.56]6.13 [2.51-14.94]1.59 [1.00-2.52]	0.41 [0.22-0.75]2.78 [1.26-6.12]4.04 [1.67-9.77]
Intermediate centrality	eRAPNYHA functional class6MWD	> 15 mm HgIII-IV< 165 m> 445 m	1.93 [1.24-3.0]1.79 [1.10-2.92]0.54 [0.34-0.85]	2.09 [1.05-4.19]
Lower centrality	Resting HRSBPeGFRSodiumAlbuminRDWDlco	> 92 beats/min< 100 mm Hg< 45 mL/min/1.73 m2< 136 MEq/L< 3.5 g/L> 15%< 60%: < 80%	1.65 [0.99-2.74]3.86 [2.26-6.62]1.78 [1.10-2.87]1.62 [1.01-2.62]2.08 [1.28-3.36]1.80 [1.14-2.85]2.13 [1.36-3.34]	3.09 [1.72-5.56]1.96 [1.22-3.15]1.89 [1.17-3.03]1.77 [1.09-2.85]
Demographic and etiology	Male sexCTDAge (per 10 y)		…	1.74 [1.02-2.96]1.65 [0.98-2.78]
Therapy	Prostanoid therapyIncident		2.21 [1.46-3.33]…	3.08 [1.94-4.89]
C-index				0.81 [0.77-0.85]
χ2				116

The variables that emerged significant in each cluster were considered in the final multivariable model; age sex and CTD etiology were also included. 6MWD = 6-min walk distance; CTD = connective tissue disease; Dlco = diffusing lung capacity of carbon monoxide; eGFR = estimated glomerular filtration rate; eRAP = estimated right atrial pressure; HR = heart rate; NT-proBNP = N-terminal pro-B type natriuretic peptide; NYHA = New York Heart Association; PVR = pulmonary vascular resistance; RDW = RBC distribution width; RVESRI = right ventricular end-systolic remodeling indexed; SBP = systolic BP; SVR = systemic vascular resistance.

Figure 5

A-D, Cox regression model in PAH. A, Independent parameters retained in the multivariable model. B, Kaplan-Meier transplant-free survival according to risk strata. C, Longitudinal changes in NT-proBNP within 1 year with 38% exhibiting a change in their NT-proBNP risk category (n = 160). D, Receiver-operating characteristic curves for the incremental value of serial changes in NT-proBNP when added to the baseline model (∗ indicates significant change in ROC). 6MWD = 6-min walk distance; Alb = albumin; CTD = connective tissue disease etiology; BNP = B-type natriuretic peptide; Dlco = diffusing lung capacity of carbon monoxide; Na = sodium; NT-proBNP = N-terminal pro-B type natriuretic peptide; PAH = pulmonary arterial hypertension; Peff = pericardial effusion; RVESRI = right ventricular end-systolic remodeling indexed; SBP = systolic BP; Tx = therapy.

Table 3

Alternative Models Considering Different Data Availability

Cluster	Parameter		Imaging Focused	Laboratory Focused
High centrality	NT-proBNPRVESRI	< 300 ng/L[1.35-1.70]> 1.70	…3.97 [1.86-8.47]7.12 [3.12-16.27]	0.27 [0.15-0.47]…
Intermediate centrality	eRAP	> 15 mm Hg	1.92 [1.02-3.62]	…
	6MWD	< 165 m	2.16 [1.08-4.34]	2.13 [1.07-4.26]
Lower centrality	SBPHR	< 100 mm Hg> 92 beats/min	4.17 [2.33-7.47]1.84 [1.12-3.03]	3.18 [1.77-5.68]
	AlbuminSodium	< 3.5 g/L<136 mEq/L	2.23 [1.39-3.57]	2.14 [1.34-3.43]2.19 [1.39-3.47]
	Dlco	< 80%	…	…
	Male sexCTD		2.09 [1.23-3.55]2.02 [1.21-3.37]	1.94 [1.16-3.25]1.87 [1.13-3.09]
Therapy	Prostanoid therapy		2.49 [1.59-3.92]	3.34 [2.12-5.25]
C-index	Full		0.78 [0.73-0.83]	0.79 [0.75-0.83]
χ2			107	98

The “…” indicates parameters not considered in the model. Imaging-focused models did not include NT-proBNP or Dlco; laboratory-focused models did not include imaging or Dlco. 6MWD = 6-min walk distance; CTD = connective tissue disease; Dlco = diffusing lung capacity of carbon monoxide; eRAP = estimated right atrial pressure; HR = heart rate; NT-proBNP = N-terminal pro-B type natriuretic peptide; RVESRI = right ventricular end-systolic remodeling indexed; SBP = systolic BP.

Thresholds for Discrimination of Outcome

Using receiver-operating characteristic analysis, we assessed thresholds for the main outcome based on the Youden index and the fixed 80% sensitivity or specificity (Table 4). For the majority of parameters, including NT-proBNP, Dlco, and PVRr, values were consistent with commonly used thresholds in PAH or laboratory reference limits. For total bilirubin, however, a lower threshold was observed

Table 4

Receiver-Operating Curve Analysis and Associated Optimal and Fixed Sensitivity and Specificity Criteria

Parameter	AUC	Youden Index		80% Sensitivity	80% Specificity
		Value	(se, sp)	Value (sp)	Value (se)
NT-proBNP	0.72 (0.66-0.78)	544	(74, 65)	398 (55)	1,363 (46)
RVESRI	0.72 (0.66-0.78)	1.47	(74, 61)	1.42 (48)	1.59 (48)
RVESA ratio	0.72 (0.66-0.78)	1.26	(62, 76)	0.89 (50)	1.4 (55)
RVLS	0.64 (0.57-0.70)	–17.4	(77, 47)	–18.3(41)	–12.8 (34)
PVRr	0.60 (0.53-0.66)	0.55	(50, 67)	0.34 (29)	0.62 (35)
Dlco	0.66 (0.59-0.72)	60	(44, 81)	82 (41)	62 (44)
SBP	0.61 (0.55-0.68)	99	(20, 97)	123 (30)	106 (37)
HR	0.59 (0.53-0.66)	87	(41, 77)	71 (31)	90 (30)
6MWD	0.66 (0.59-0.72)	409	(64, 62)	490 (35)	320 (41)
eGFR	0.62 (0.55-0.68)	62	(73, 47)	65 (36)	46 (36)
Sodium	0.64 (0.58-0.70)	136	(45, 75)	139 (37)	136 (40)
Albumin	0.63 (0.56-0.69)	3.5	(42, 81)	4.0 (32)	3.5 (43)
Total bilirubin	0.64 (0.57-0.70)	0.7	(51, 73)	0.2 (31)	0.9 (39)
Hemoglobin
Female	0.61 (0.54-0.69)	11.2	(31, 90)	14.6 (28)	11.8 (31)
Male	0.71 (0.57-0.83)	15.2	(67, 74)	16.5 (35)	14.5 (50)
RDW	0.69 (0.63-0.75)	14.9	(72. 63)	14.6 (53)	16.9 (33)

se = sensitivity; sp= specificty; 6MWD = 6-min walk distance; AUC = area under the curve; Dlco = diffusing lung capacity of carbon monoxide; eGFR = estimated glomerular filtration rate; HR = heart rate; NT-proBNP = N-terminal pro-B type natriuretic peptide; PVRr = relative pulmonary vascular resistance; RDW = RBC distribution width; RVESA = right ventricular end-systolic area; RVESRI = right ventricular end-systolic remodeling indexed; RVLS = right ventricular longitudinal strain; SBP = systolic BP.

Serial Changes in NT-proBNP (Central Feature)

In 160 patients (69%), follow-up NT-proBNP values were available within 1 year. Thirty-eight percent of patients transitioned to a different NT-proBNP risk category at follow-up (Fig 5C). The low risk NT-proBNP group remains the most stable over time, with 88% remaining in the same risk category. Serial changes in relative NT-proBNP (natural logarithmic scale) were associated with event-free survival (C-index of 0.64 [0.58-0.70]). When added to the baseline model, the C-statistic for 5-year outcome improved (from 0.78 [0.71-0.84] to 0.85 [0.78-0.90] corresponding to an increase in C-statistic of 0.071 ± 0.024;P = .0030) (Fig 5D), with Harrell’s C-index for long-term outcome improving from 0.75 (0.70-0.80) to 0.80 (0.75-0.84). Factors associated with relative change in NT-proBNP included CTD etiology (r partial = 0.18), hypoalbuminemia (r partial = 0.16), hyponatremia (r partial = 0.18), and NT-proBNP risk strata (r = –0.21).

Discussion

In the current study, network graph analysis not only identified NT-proBNP as a central hub in PAH but also provided explainability for risk stratification models. We found that risk prediction in PAH was best achieved by combining a small number of central right heart features with less connected parameters reflecting functional capacity, end-organ dysfunction, and etiology. Moreover, risk was not static as shown by the serial NT-proBNP analysis.

One of the cornerstones of network analysis is the concept of centrality, which quantifies a node’s importance. Because this network was enriched with right atrial and right ventricular parameters, eigenvector centrality represented a global measure of right heart connectivity. It is therefore not surprising that eigenvector centrality was more strongly related to outcome than betweenness centrality. As an approach original to our study, we not only used centrality to identify a central hub but also to cluster parameters according to their connectivity.

NT-proBNP emerged as the central hub of our clinical network, with independent associations with RV dimension, right atrial pressure, and renal function. Nagaya et al were the first to highlight the prognostic value of BNP at baseline and follow-up in PAH. More recently, using the REVEAL Registry, Frantz et al confirmed the value of baseline and follow-up BNP as a predictor of 5-year outcome. Consistent with these studies, we also found that baseline and serial changes in NT-proBNP were strongly related to event-free survival. Its reproducibility, availability, and strong prognostic value explains why BNP or NT-proBNP are included in the majority of risk scores in PAH.

In recent years, much attention has been focused on identifying prognostic imaging markers in PAH. Many indexes have been considered, including TAPSE, RV longitudinal strain, RV end-systolic dimension, RVEF, or, more recently, right atrial size or function or pulmonary artery pulsatility. Although most of these parameters have been strongly associated with event-free survival, recent studies highlight the value of RV end-systolic dimension. In fact, RV end-systolic dimension emerged as a strong marker of ventriculo-arterial coupling and adaptation and is sensitive to increases in afterload; in addition, RV dimensions are less likely to pseudo-normalize in the presence of severe TR.,23, 24, 25 Few studies, however, have directly compared RV end-systolic dimension with RVLS or TAPSE., In the study of Ryo et al, RV end-systolic volume index quantified by using three-dimensional echocardiography was more predictive than RV global area strain. In addition, we have found that RVESRI was also more predictive than RVLS. Risk prediction does nor, howvever, rely on a single parameter in isolation. It is therefore not surprising, given the intertwined nature of connections in PAH, that other combinations may emerge across studies. For example, Fine et al has previously shown the incremental prognostic value of NT-proBNP and RVLS in a large study from the Mayo Clinic., More recently, Ghio et al also reported, using a pooled analysis, that TAPSE, TR severity, and inferior vena cava size were key features in defining risk. These studies show that different combination of right heart metrics may yield equivalent discrimination of outcome the choice of metric should be guide by center expertise and reproducibility of measures.

Our study also highlights the importance of markers of end-organ dysfunction. These markers are particularly useful for differentiating patients at intermediate and high risk of clinical worsening.,, For example, Dlco, which depends on both pulmonary blood flow and microvascular lung architecture, has been consistently associated with worse prognosis in PAH., Low Dlco could additionally be found in patients with pulmonary venoocclusive disease, which also carries a worse prognosis. Similarly, as highlighted by the REVEAL and PHORA studies, impairment of BP regulation or renal, hepatic, or hematologic function is associated with increased mortality in PAH.,, The pathophysiology of end-organ dysfunction in PAH is complex, and several mechanisms such as neurohormonal activation, oxidative stress, and nutritional status may affect BP regulation, serum sodium, RBC distribution width, uric acid, or albumin.27, 28, 29, 30 Disease trajectories are also important to consider in risk models. This has been incorporated in the REVEAL 2.0 by including recent all-cause hospitalization (≤ 6 months) as one of the predictive factors. As shown in the current study and others, serial changes in NT-proBNP (this study’s central feature) may also refine risk prediction.

Network graphs provide explainability for alternate models in PAH. For example, despite having few or no imaging parameters, multicenter risk scores achieve strong prediction of outcome. Consideration of the highly connected network is also important when evaluating novel biomarkers. For example, novel biomarkers need to be placed in the context of other routinely available biomarkers; in fact, co-linear biomarkers are unlikely to improve risk prediction (although they can provide invaluable mechanistic insights).31, 32, 33

Taken together, clinicians often favor models in which complexity matches actionability. Feature selection in models should also be guided foremost by physiology but also informed by data architecture, its predictive ability, reliability, cost, risk, and availability for serial measurements. The simplest model would include demographic characteristics, etiology, functional class, health status, vital signs, laboratory data, and recent trajectories. Moreover, ongoing iteration of the PHORA models will continue to incorporate tolerance to missing data (or alternate feature selection) and use Bayesian network modeling approaches, thus providing a nuanced and pragmatic approach to risk stratification.

The current study has several limitations. First, this was a single-center study; the analysis, however, includes a comprehensive number of features and long-term follow-up. Second, although echocardiography, laboratory, and right heart catheterization data were not collected simultaneously in all participants, they were collected as part of a registry with representative nature of the data. We also verified that centrality was not affected by imputation of missing data (especially relevant for Dlco). Third, because oxygen consumption was assumed rather than directly measured, this likely affected the prognostic value of cardiac index. Finally, TAPSE was measured by using two-dimensional methods, which in our laboratory have excellent correlations with M-mode measurements.

Interpretation

Network graph analysis identified NT-proBNP as a central feature in PAH and provided explainability for risk stratification models in PAH.