An explainable artificial intelligence approach for predicting cardiovascular outcomes using electronic health records

Understanding the conditionally-dependent clinical variables that drive cardiovascular health outcomes is a major challenge for precision medicine. Here, we deploy a recently developed massively scalable comorbidity discovery method called Poisson Binomial based Comorbidity discovery (PBC), to analyze Electronic Health Records (EHRs) from the University of Utah and Primary Children’s Hospital (over 1.6 million patients and 77 million visits) for comorbid diagnoses, procedures, and medications. Using explainable Artificial Intelligence (AI) methodologies, we then tease apart the intertwined, conditionally-dependent impacts of comorbid conditions and demography upon cardiovascular health, focusing on the key areas of heart transplant, sinoatrial node dysfunction and various forms of congenital heart disease. The resulting multimorbidity networks make possible wide-ranging explorations of the comorbid and demographic landscapes surrounding these cardiovascular outcomes, and can be distributed as web-based tools for further community-based outcomes research. The ability to transform enormous collections of EHRs into compact, portable tools devoid of Protected Health Information solves many of the legal, technological, and data-scientific challenges associated with large-scale EHR analyses.


Introduction
The application of data-science methods to electronic health record (EHR) databases promises a new, global perspective on human health, with widespread applications for outcomes a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Children's Hospital under an IRB that waived consent (see ethics statement). We refer to this cross-institution extract as the Utah Data Resource. Because the aggregate is comprised of exact dates and other protected patient information, the data cannot be made publicly available. Information regarding how qualified researchers might apply for data access can be found here https://irb.utah.edu/about/contact/. However, All Probabilistic Graphical Models described in this paper are available through the web using the following link: https://pbc.genetics. utah.edu/lemmon2021/bayes/.
Results for three different starting cohort sizes are shown. Note how stratification lowers the strength of p-values as a function of the size of the stratum. This effect is exacerbated when more than a few potentially confounding variables are controlled for, and stratification quickly results in cohorts that are too small for discovery activities, as the comorbidities fail to achieve statistical significance. For example, using a starting corpus of 9,525 records, stratification followed by χ 2 analysis fails to detect the well-known comorbid relationship between AF and Stroke for female patients aged 50-59 when white ancestry is included in the stratum description. By contrast, the PBC approach maintains power across all comparisons. For more on these points, see [10].

Comorbidities of heart transplant
We evaluated every pairwise combination of diagnoses, procedures, and medications mentioned in our EHR corpus for comorbid associations, using PBC [10] to adjust on a patient-by-patient basis for the potentially confounding demographic variables shown in Fig 1. Fig 2A summarizes the results of this computation as a patient disease network. The network provides a visual overview of the entire EHR corpus, wherein every node (state) is a diagnosis, procedure, or medication, and edges denote Bonferroni significant comorbid relations between terms. Given a node of interest, heart transplant, for example, its comorbid diagnoses and associated procedures and medications can be recovered by following edges to that node back to their terms.
The transition probabilities associated with each edge provide means to calculate the pairwise contributions of each term to the outcome's observed (marginal) frequency in the EHR corpus. This provides a way to intuit an outcome's comorbidity landscape, and calculate the expected flux of patients through that region of the network. These patient 'trajectories' provide a framework for cost prediction and service allocation activities. For example, the trajectory for adult heart transplant (2B) tracks the time course of diagnoses, procedures and medication use preceding and following heart transplantation. Thus, one can follow the trajectory of ischemic heart disease, flowing through the diagnosis of heart failure, cardiogenic shock, administration of the vasoactive medication milrinone, and culminating in heart transplantation with subsequent downstream complications. Crucially, this methodology provides precise measures of patient flux between these nodes.

Multimorbidity network for heart transplant supports conditional outcome risk calculations
Although trajectories provide intuitive and useful overviews of the comorbidity landscape, effective outcomes research requires calculating the joint contributions of conditionally dependent multimorbid terms on an outcome. We leverage Probabilistic Graphical Models as an explainable AI solution for this computationally intensive task. Fig 3A illustrates a multimorbidity network derived from a temporalized Probabilistic Graphical Model for the predisposing comorbidities of adult heart transplant presented in Fig 2B. Because the edges in a multimorbidity network denote conditional dependencies between terms, rather than transition probabilities, the multimorbidity network's topology is necessarily different from the trajectory topology shown in Fig 2B. The PGM provides easy means to calculate outcomes risk for any combination of variables in it. For example, a prior diagnosis of cardiomyopathy (nonischemic) increases the risk of heart transplantation 86±35 fold, whereas a diagnosis of viral myocarditis confers a 59±21 fold increase in risk. The strongest single variable for heart transplant risk is the use of the vasoactive medication milrinone, which increases risk 175±30 fold. Note that we are not suggesting milrinone causes heart transplant-rather that the prescription of milrinone in a patient's medical record is a powerful predictor of future heart transplant. The utility of PGMs for outcomes research is best illustrated by their application to problems of complex multimorbid outcomes analyses, where conditional dependencies of these variables interact to further modulate risk for the outcome under study. For example, we can explore the role of heart disease etiology on transplant risk in the context of milrinone infusion. Thus, a cardiomyopathy patient requiring milrinone has a 407±101 fold increased risk for heart transplant. Likewise, a patient with viral myocarditis requiring milrinone therapy has a 346±93 fold increased risk for heart transplant; while milrinone use in a patient with ischemic heart disease confers a 205±28 fold increased risk of heart transplant. Moreover, while both cardiomyopathy and ischemic heart disease have similar increased risks for heart transplant in isolation (86±35 fold and 64±14 fold, respectively), cardiomyopathy patients who require milrinone therapy are at far greater risk for heart transplant than patients with ischemic heart disease requiring milrinone. Additional conditional queries conducted with the PGM are presented in Fig 3A. This list is by no means exhaustive-the PGM is capable of answering an astonishing number of queries-3 25 to be precise. We encourage the reader to explore these by following the link to the corresponding web application https://pbc.genetics. utah.edu/lemmon2021/bayes. In this context, the explainable nature of PGMs lays the foundation for massively parallel testing of novel hypotheses between multiple, complex clinical variables of interest.
The comorbidity landscape for pediatric heart transplant is dramatically different from that of adults, as it includes a large contribution from congenital heart defects (CHD) and palliative surgical procedures. Fig 3B presents a multimorbidity network for 13 common CHD terms defined by echocardiogram and identified by PBC as comorbid with pediatric heart transplant. A prior diagnosis of dilated cardiomyopathy (DCM), defined as genetic/idiopathic DCM, increases a child's risk for heart transplant 102.2±33.6-fold, over the marginal probability of transplant. Among single ventricle forms of CHD, patients with hypoplastic left heart syndrome (HLHS) are at the greatest risk for heart transplant (56.8±17.8-fold), as compared to tricuspid atresia (17.1±11.8-fold) or laterality defects (25.8-fold ± 8.5). Again, the utility of PGMs for complex multimorbid analyses is highlighted by the ability to calculate the additional risk for heart transplant in a child with a laterality defect, if that child also requires the Norwood surgery (51.3±10.5-fold).

Multimorbidity network for sinoatrial node dysfunction supports multimorbidity risk calculations for a range of clinical and demographic health predictors
Fig 4A extends the investigations to include the impacts of these same pediatric heart surgeries in the context of various CHD phenotypes on a different clinical outcome, sinoatrial node dysfunction (SND). The Fontan surgery dominates the landscape of pediatric SND, increasing the risk 19.6±6.4-fold over the marginal probability of SND. Moreover, Fontan surgery is the only clinical variable with a direct connection to SND; the other clinical variables connect indirectly to SND via the Fontan node. Thus, the relative risk of SND for specific forms of single ventricle CHD (HLHS, tricuspid atresia, unbalanced AVSD) following the Fontan surgery are similar (Fig 4), indicating that the Fontan surgery itself is the primary indicator of future SND, rather than the underlying form of CHD that required the procedure. Collectively, the preceding analyses demonstrate how multiple nets can be used in tandem to address complex multimorbidity outcomes questions.
Multimorbidity networks also provide powerful means to investigate the impacts of various demographic factors upon outcomes. The net in Fig 4B models the multimorbid landscape surrounding SND in adult patients. As SND and AF are both risk factors for each other [31], we temporalized the network (see Methods) to analyze clinical variables that precede SND. The ancestry and ethnicity nodes enable explorations of demographic impacts upon SND and its comorbidities. Thus, in the University of Utah Hospital system, a Hispanic patient with AF has a 61±6 fold increased risk of SND, compared to 30±1 fold risk for white ancestry and 40±7 fold risk for African Americans. These results underscore the potential of our approach to inform ethnic/racial health care disparities with precise, quantitative results, and in the context of a specific health care system. Moreover, these findings illustrate how our approach can empower these discussions despite demographic skews in the underlying EHR corpus (see S2 and S3 Tables); an important finding for the Utah health system.

Multimorbidities of congenital malformations augmented by maternal health data
The impact of maternal health on health outcomes in the child is an area of intense investigation. The Multimorbidity network shown in Fig 5A places a child's risk for congenital malformations in the context of a maternal diagnosis of pregnancy-induced hypertension (HTN-PREG) during that pregnancy, leveraging outcomes data for over 130,000 births at the University of Utah Hospital system over the last 15 years. HTN-PREG elevates the risk of cardiac and circulatory congenital anomalies 1.83±0.03-fold, an effect not due to maternal age differences (S1 Fig). The multimorbidity network also illuminates the strong dependencies between clinical variables and allows for quantitative assessments of risk. For example, a diagnosis of Down Syndrome is associated with a 25.9±0.8-fold increased risk for a congenital cardiac anomaly (S4A Table). Moreover, a child with a congenital cardiac anomaly is a priori 9.2±0.9-fold more likely to have a nervous system anomaly than baseline (S4B Table). The impact of maternal health on a child's The clinical variables were chosen based on Bonferroni-corrected ICD10 and RXnorm billing codes significantly associated (preceding) with heart transplant. Each node represents a diagnosis, procedure, or medication code and each edge represents a conditional dependence between nodes. For detailed description of the clinical variables, please refer to S5 risk of CHD is further explored in Fig 5B. Our ability to seamlessly combine and compute upon maternal/child EHR data highlights the extensibility of our approach to study health outcomes across generations in order to define the impacts of maternal health on childhood outcomes.

Web-based outcomes calculators
We repackaged the multimorbidity networks as stand-alone web-based outcomes calculators. This allows users to interact with a multimorbidity network as an 'app', whereby they can use slider buttons to toggle values of its states and to select an outcome of interest. These web-apps are available here: https://pbc.genetics.utah.edu/lemmon2021/bayes/bayes.

Ethics statement
Human subjects approval for this study was obtained following review by the University of Utah Institutional Review Board, IRB_00095807 under a waiver of consent and authorization. Patient data was not anonymized prior to the start of the study. All authors completed Human Subjects research requirements.

Utah data resource
The University of Utah maintains an Enterprise Data Warehouse (EDW)-a central storage and search facility for all clinical data collected from all affiliated University hospitals and clinics across the Intermountain West. SQL queries were used to aggregate data from various tables and collect the following information: (1) gender, ancestry, ethnicity, and age for each patient; (2) list of patient visits, along with visit dates, and medical terms associated with each visit, including diagnostic codes, procedure codes, and medications ordered. ICD9 and ICD10 diagnosis codes consist of 18,000 and 142,000 codes respectively, while procedural codes (CPT) include around 10,000 codes. In all, we collected records for 1.6 million patients, 21 million visits and 166 million diagnosis (DX), procedure (PX) and medication (RX) codes. See S1-S5 Tables for additional details.
We combined these data with the Primary Children's Hospital's database of echocardiographic variables (diagnoses, ventricular function, valve gradients, chamber/vessel sizes, etc.) dating back to 2006 for 65,618 probands, 44,254 of which also appear longitudinally in the University's EDW. These data contain 529,317 mother-child pairs with EHR data, 14,155 of which include a child with echo data, allowing us to study maternal contributions to congenital heart disease (CHD). Collectively, these data comprise the Utah Data Resource (UDR). For the purposes of computation, custom encryption is applied to the UDR to produce data free of protected health information (PHI) and unintelligible without its cyphers. We can then generate statistics on this PHI free data in a variety of compute environments, decrypting the results on PHI approved machines.  In this analysis, a patient's diagnoses are inferred via billing codes. Thus, the investigations and risk calculations presented herein reflect medical practice within the University of Utah Hospital network and Primary Children's Healthcare. How closely they approximate underlying universal ('true') risks is still unknown. Moving forward, we note that the methods described below provide powerful means for large-scale cross institutional comparisons aimed at discovering differences in medical practice and billing trends.

Patient disease network
We used a Poisson Binomial based methodology called PBC [10] to discover comorbidities within our EHR corpus. Standard methods such as stratification seek to control for confounding variables through 'stratifying' by age and gender (for instance) and calculating comorbidity statistics for each strata, under the restrictive assumption a every patient in a stratum has the same probability of manifesting each morbidity. However this approach fails to scale, since the use of many confounding variables leads to strata too small to detect a statistical significance comorbidity. In contrast, PBC models the effects of age, gender, race, ethnicity, insurance type, and the length and density of each patient's medical record. These input features are used to determine per-patient probabilities for each medical term, using a Poisson binomial test. The result is much greater statistical power [10].
PBC was used to find significant connections among every possible combination of ICD diagnoses, procedures, and RxNorm medication terms, thereby creating a patient disease network [10]. Patient disease network is a term borrowed from Capobianco et. al 3 and denotes a network comprising all significant connections among diagnoses, procedures and medications (Bonferroni p-value cutoff 10E-9.48). We only considered terms appearing in at least 15 patients. This filter reduced the number of unique terms to 39,055 ICD10 diagnosis codes, 5,716 CPT procedure codes, and 1,764 RxNorm medication codes. We used Minimum Description Length clustering [32] to visualize the data, so that nodes with similar combinations of edges would lay near one another in the network. We also determined the patient flux between every pair of nodes. The result is shown in Fig 2A, which provides a visual representation of our patient disease network for the entire EHR corpus.
In keeping with previous work [13,[33][34][35][36] on patient disease networks, we refer to a sub-portion of the network, focused on a single outcome as a trajectory, or term trajectory. Fig 2B shows a trajectory for adult heart transplant. Trajectories provide means to display additional features of the network, such as transition probabilities (which correspond to patient flux between nodes), and the marginal frequencies of outcomes and comorbid terms within the EHR corpus. Collectively, this information allows for better intuition of the disease landscape surrounding an outcome. The trajectory is also a useful starting point for cost and service allocation calculations.

Multimorbidity networks
While trajectories describe transition probabilities between two comorbid terms, they provide no means to determine the combined effects of multiple comorbid diagnoses, and associated clinical procedures and medications upon an outcome. We have employed Probabilistic Congenital abnormalities of the gastrointestinal tract; Nervous: Nervous system congenital abnormalities; Eye: Congenital abnormalities of the Eye; CleftLip: Cleft lip. Panel B. Multimorbidity landscape for child's risk of congenital heart defects in the context of pregnancy-induced hypertension. N = 125,014 mothers. ASD, atrial septal defect; VSD, ventricular septal defect, HLHS, hypoplastic left heart syndrome; Coarctation, coarctation of the aorta; TOF, tetralogy of fallot; BAV, bicuspid aortic valve. For detailed description of the clinical variables, please refer to S5 Table. The target node (HTN-PREG) is colored red and nodes with direct connections to the target (ie, within the Markov blanket) are circled red. Values in Tables represent mean ± standard deviation. https://doi.org/10.1371/journal.pdig.0000004.g005 Graphical Models (PGMs) to overcome this limitation. We learned the structures of the PGMs using the python3 package "pomegranate" [28], which provides a Bayesian Information Criterion (BIC)-based DP-A � exact structure search algorithm [37,38,46]. The exact search algorithm explores the entire applicable space of conditional dependencies in order to discover the optimal network structure for the data. Parameter learning for this optimal network is accomplished using the loopy belief propagation algorithm [39]. We use the same package for our inference and multimorbidity risk calculations. The visual interpretation was designed using the graph_tool [40] Python3 package and D3.js Java library.
For each Probabilistic Graphical Model, a maximum of 25 comorbid features were selected using PBC and validated by experts in the medical field (TAM, DW, MDP, BEB, RUS, MTF). Features that were judged to be of clinical relevance, importance or interest for the field under study were selected and used as inputs to learn the PGM structure and infer risk. These selected features became the inputs used to learn the PGM structure and infer risk. The patient's features were described in a categorical data format, (e.g. indicating the ancestry, ethnicity, or insurance type) or "present/absent" binary variables in case of medical diagnoses and procedures. A continuous feature (e.g. age, BMI, blood pressure) were discretized based on established clinical thresholds. Because the PGMs only present the facts about the data, PGMs themselves cannot discover or infer the temporal order of the events (unless specified as a Dynamic PGM). To overcome this issue, for our temporalized PGMs we have imposed the order (discovered using PBC; see [10] for additional details.) on the EHR extraction process prior to learning the Probabilistic Graphical Model structure. When trained on temporalized data, PGMs are forced to learn temporal conditional probabilities. Missing data are handled inherently by the Probabilistic Graphical Model structure learning process. That is, no patients were excluded due to missing data and no missing data was imputed. The resulting temporalized structures we call multimorbidity networks.
Probabilistic Graphical Models represent conditional dependencies in the dataset as a directed acyclic graph (DAG); however, it is important not to confuse directionality with causality or temporal ordering. In keeping with best practice, the multimorbidity networks are visualized in their undirected, moralized form, in which every node is connected to its Markov blanket. A single constructed multimorbidity Network provides an inference engine capable of answering O(3 n ) personalized conditional risk queries, where n denotes the number of features describing a patient's condition, and the base of the exponent is 3, because in case of binary health records data there are three states for each node that can be specified: present, absent, or status unknown.

Confidence values
Risk estimates derived from Probabilistic Graphical Models are maximum likelihood estimates given the optimal structure under the BIC and an assumed uniform prior probability of any distinct EHR. To obtain standard deviation values for these estimates, we created 100 nets in parallel [41] from bootstrap replicates of the same data used to create Figs 3, 4 and 5. We then queried the resulting replicate nets, and calculated standard deviations of risks of outcomes of interest.

Discussion
The ability to model dependencies among multiple risk factors is crucial for meaningful outcomes research. Unfortunately, traditional techniques, such as logistic regression, have limited ability to capture so-called 'conditional dependencies' between variables, which are the heart and soul of multimorbid analyses. Although mixture and generalized linear models with mixed effects can (in principle) overcome this weakness, these techniques are limited because a new model must be designed for every question. Neural nets provide one possible alternative. Although they can account for non-linear interactions in the data and are scalable [7], Neural nets are often referred to as 'black boxes' (i.e., lacking explainability) [14,15,20,21,[42][43][44][45][46] due to the difficulties in determining precisely how and why different input variables were used to produce the outputs.
Because we sought not merely to predict outcomes, but also to understand the relationships between multiple clinical variables and outcomes, we selected an 'explainable' AI solution, rather than a black box approach. Probabilistic Graphical Model-based [23][24][25]46] multimorbidity networks offer best-practice solutions to this problem. Moreover, they effectively model data without recourse to a fixed decision protocol (e.g decision trees), and are resilient to missing/unknown data. Crucially, the contributions of different combinations of variables to an outcome can be precisely and easily determined.
Explainability comes at a cost; unlike Neural nets, which are incredibly scalable, multimorbidity networks can model a maximum of only 30 or so variables at once [28,37,38]. It is therefore necessary to pre-identify high impact variables when modeling an outcome, a need fulfilled by PBC [10]. We argue that the ability to rigorously investigate interrelations among 30 or so primary determinants represents a giant step toward understanding cardiovascular disease.
Our results illustrate how multimorbidity networks provide explainable solutions for understanding the joint impacts of diagnoses, medications, and medical procedures on cardiovascular health outcomes. We emphasize that the necessarily brief results reported here hardly exhaust the contents of these machineries. Consider that a multimorbidity network with n nodes supports~3 n possible queries. The net shown in Fig 4B, for example, supports~3 14 different queries-a number that gives some indication both of the complexity of the data being extracted from the EHR corpus by our approach, and the value of these multimorbidity networks to further outcomes research.

Conclusion
The analyses presented here provide a first step toward a global description of heart disease and associated comorbidities across the USA intermountain west. However, the map we seek resides not so much in the results reported here, as it does in the products of our analyses: the PGM multimorbidity networks. As we have explained, these networks support multitudes of queries, and when used in combination, support both wide-ranging and focused explorations of a disease landscape. Given the right datasets, we have shown that the approach can provide new insights, such as the mother-child cross-generational cardiovascular multimorbidities we described. However, our approach also has limitations. Our exact approach allows us to model at most~30 health conditions at a time. In future work we would like to relax this limiting factor by allowing approximate solutions that enable us to scale up the complexity of the multimorbidity networks to thousands of health conditions. Another area for innovation regards incorporation of continuous variables, as current software packages do not allow us to incorporate such variables at scale, however there is no theoretical limitation preventing their use in a PGM framework.
A major strength of our approach is that these outcomes machineries can be redistributed as web-based tools. Indeed, the multimorbidity Networks described here have been made available online [pbc.genetics.utah.edu/lemmon2021/bayes], with the hope that the wider scientific community will find them useful for their own outcomes research. The ability to transform enormous collections of EHR data into compact, portable machines for outcomes research, with no exchange of PHI, solves many of the legal, technological, and data-scientific challenges associated with large-scale EHR analyses.
Supporting information S1 Fig. Distribution density plot of mother's age at pregnancy, with and without hypertension complicating pregnancy. Blue line: mothers with diagnosis of hypertension complicating pregnancy (N = 11,523 mothers). Red line: mothers without diagnosis of hypertension complicating pregnancy (N = 113,491 mothers). (TIF) S1  Fig 5 of the main text. For example, a child with a known diagnosis of Down Syndrome has a 25.9-fold increased risk of a cardiac congenital anomaly over the marginal risk of cardiac anomaly. HTN-PREG, hypertension complicating pregnancy (AKA pregnancy-induced hypertension). For detailed description of the clinical variables please refer to S5 Table. (TIF) S5