A vascular biology network model focused on inflammatory processes to investigate atherogenesis and plaque instability

Selventa Knowledgebase (SK)

RCR-based HYP generation process

The RCR methodology utilized for network augmentation has been described previously[15] and a detailed description may be found in the supplementary methods (Additional file2). Briefly, RCR analysis identifies potential HYPs for the statistically significant mRNA State Changes observed in the transcriptomics datasets. These upstream controllers are termed HYPs, as they represent statistically significant hypotheses that are potential explanations for the observed mRNA State Changes (Figure 1C). Detailed descriptions of the probabilistic scoring metrics (richness and concordance) can be found in Catlett et al.[15], whereas the use of causal assertions in the construction of the V-IPN are further described in the Additional file2: Supplementary Methods.

V-IPN construction: model structure and boundaries

The workflow for the creation of the V-IPN is illustrated in Figure 1A. The initial literature-based network scaffold was defined by specific cell, tissue, species and disease contexts (e.g., ECs, aorta, human and atherosclerosis) known to be implicated in vascular pathobiology. The V-IPN nodes and edges comprising the scaffold were assembled in a sequential process by first using causal connections derived from knowledge published in the scientific literature and captured by the SK (Figure 1B).

The literature-derived framework was further augmented with nodes derived from the RCR analysis of vascular inflammation transcriptomics datasets (referred to as “model building” datasets, Table 1). RCR analysis yielded several dozen additional HYPs that were vetted for biological relevance and incorporated into the network as new nodes. Such nodes were connected to the literature scaffold using causal relationships captured by the SK. The resulting integrated network was manually reviewed by scientists with expertise in vascular biology and inflammation. The modular framework consists of six subnetworks that accompany this manuscript in XGMML and .XLS formats (Additional file3). The network architecture may be viewed from the XGMML files using freely available network visualization software such as Cytoscape (http://www.cytoscape.org/).

Gene expression datasets used for V-IPN construction and evaluation

Four transcriptomics datasets from isolated human ECs (GSE13139 [Hs_EC_GFP_oxLDL_vs_ct])[24] and atherosclerotic human coronary arteries (GSE40231 [Hs_athCA_vs_ctIMA])[25], as well as murine aortas (E-MTAB-1696 [Mm_Ao_16w_ApoE_CS_vs_sham]), were analysed by RCR. RCR results (HYPs) were then used to evaluate network performance by determining HYP-level coverage and odds ratios (OR) across the six subnetworks constituting the V-IPN. Names describing the species and experimental settings for each dataset were created and they were used throughout the results and discussion section to facilitate comparative analyses.

Calculation of coverage and odds ratio

Evaluating the RCR results for each dataset in the context of each of the six V-IPN subnetworks was estimated by calculating coverage and odds ratio (OR). Figure 2 depicts a schematic representation of coverage (sensitivity) and OR, as well as the equations involved in their calculation. The dataset’s coverage was calculated as the fraction of possible HYPs in each subnetwork that are significant HYPs in a given dataset (Figure 2A, B and C). The OR was calculated as the odds of having significant HYPs in the network divided by the odds of having non-significant HYPs in the network (Figure 2B and C). Thus, sensitivity is an estimate of subnetwork coverage (overlap), whereas OR estimates the odds of significant HYP enrichment for a specific dataset-subnetwork pair. An OR higher than one implies that the odds of having significant HYPs in a given subnetwork are higher than the odds of having not significant HYPs in that subnetwork (Figure 2D). The larger the OR of a given dataset, the better the network encompasses the biology embedded in the dataset.

CS generation and ApoE^-/- mice exposure

We set up a study in ApoE-deficient mice in which we investigated the effects of CS on cardiovascular endpoints including plasma lipid profiles, and transcriptomics of aortas. The E-MTAB-1696 [Mm_Ao_16w_ApoE_CS_vs_sham] transcriptomics dataset was generated from aortas displaying evidence of atherosclerotic plaques in ApoE^-/- mice exposed to cigarette smoke (CS). All animal experimental procedures and CS exposure were approved by an Institutional Animal Care and Use Committee (IUCAC) and are described in detail in the Additional file2: Supplementary methods. Total cholesterol measurements in plasma, atherosclerotic plaque measurements in the aortic arch and immunohistochemical stainings in vascular tissues of ApoE^-/- mice were conducted according to methods detailed in the Addditional file2: Supplementary methods

V-IPN construction and biological integration: description of modular framework and boundaries

To capture the diverse array of biological processes involved in the development of atherosclerotic plaques, the V-IPN network model was constructed using a modular approach that represents key processes related to vascular inflammation and atherogenesis. The biology modelled started with a literature-derived network scaffold followed by a RCR analysis of two murine and one human transcriptomics datasets (“Model Building” Datasets, Table 2) to enhance the representation of biological disease mechanisms from both species. This RCR-based enhancement of the networks helped uncover additional disease-relevant mechanisms not readily identified during the literature-based component of model building. This model building step strengthened the network’s capability to interpret datasets from multiple species. The RCR predictions (HYPs) were included in the network if they had been reported to be mechanistically linked to the process of interest. Biological processes represented in the V-IPN were integrated into six distinct subnetworks that captured key pathobiological events in vascular disease: Endothelial Cell Activation, Endothelial Cell-Monocyte Interaction, Foam Cell Formation, Platelet Activation, Smooth Muscle Cell Activation and Plaque Destabilization (Figure 3). The first five subnetworks describe fundamental atherogenic mechanisms underlying vascular inflammatory responses in discrete cellular populations, whereas the Plaque Destabilization subnetwork represents a collection of molecular events that occur in advanced, unstable atherosclerotic lesions. While each subnetwork was constructed to reflect the biological relationships that are involved in specific processes, the subnetworks contain shared elements. For example, since the transcription factor NF-kB is a pleiotropic protein involved in multiple inflammatory pathways, the node “transcriptional activity of NFkB” exists in multiple V-IPN subnetworks, which display varying network connectivity depending on the focus of the subnetwork. Twenty five RCR-predicted HYPs representing mostly pro-inflammatory biological signaling were present in four of the six subnetworks and Ox-LDL was the only HYP common to all subnetworks (Additional file4: Figure S1). Absolute and relative numbers for all overlapping nodes between subnetworks are shown in Additional file5: Figure S2. The provided XGMML encoding of the subnetworks allows for the assembly of the V-IPN as a single, agglomerated network using freely available network visualization software such as Cytoscape (http://www.cytoscape.org/).

Assessment of V-IPN subnetwork-level HYP coverage

Test datasets

To test the ability of the V-IPN to identify biological processes modulated by a variety of experimental conditions and models related to atherogenesis, we evaluated a series of in vitro and in vivo transcriptomics datasets by RCR, the details of which are summarized in Tables 2 and3. Lists of HYPs for each dataset meeting significance criteria (richness and concordance p-values < 0.05) were evaluated for subnetwork-level coverage, i.e. assessment of HYPs identified as significant in the datasets and also present in the V-IPN. In order to statistically validate the processes being represented in each of the transcriptomics datasets, ORs integrating sensitivity and specificity metrics into a single number were calculated. The overall V-IPN coverage ranged from 9 to 24%, with the human dataset from atherosclerotic coronary arteries (Hs_athCA_vs_ctIMA) exhibiting the highest degree of coverage (24%) and the highest OR (2.47) (Figure 2D). Datasets from Jurkat cells (Hs_JurkT_ars_vs_ct) and an immortalized EC line (Hs_EC_GFP_oxLDL_vs_ct) showed the lowest HYP coverage (9%) (Figure 2D). A summary of the absolute HYP coverage across all six V-IPN subnetworks for each dataset utilized in the network evaluation is presented in Table 4. A significant proportion of transcriptomics data-derived HYPs were captured across the six subnetworks. HYPs predicted from the Hs_athCA_vs_ctIMA dataset displayed the highest level of coverage across four V-IPN subnetworks, particularly within the Plaque Destabilization subnetwork where coverage approached 40%. The high range of HYP coverage utilizing a dataset derived from human coronary arteries displaying advanced atherosclerotic lesions underscores the validity of the V-IPN in capturing vascular disease-related processes.

	Dataset name	Dataset ID	N° of state changes (SC)	N°of HYPs in dataset	NoHYPs overlapping withEC activation	N° HYPs overlapping withplatelet activation	N° HYPs overlapping withEC-monocyte interaction	N° HYPs overlapping withfoam cell formation	N° HYPs O/L with SMC Activation	N° HYPs O/L with Plaque Destabilization	N° HYPs O/L with V-IPN
Total N° of possible HYPs in each subnetwork for human/mouse				2410/ 2354	155/149	59/55	30/27	118/113	88/84	104/102	340*/334
Mouse	Mm_Ao_78w_ApoE_vs_wt	GSE10000	3872	449	59	11	13	42	32	43	110
	Mm_Ao_16w_ApoE_CS_vs_sham	E-MTAB-1696	1928	245	24	8	4	23	19	18	52
Human	Hs_EC_GFP_oxLDL_vs_ct	GSE13139	398	152	16	4	4	16	9	14	31
	Hs_EC_LOX1_oxLDL_vs_ct	GSE13139	752	224	30	8	8	21	16	21	52
	Hs_EC_oxPAP_vs_ct	GSE20060	609	309	45	12	8	21	23	27	77
	Hs_athCA_vs_ctIMA	GSE40231	3643	309	42	17	11	36	20	37	80
	Hs_NHBE_CDKinh_rel_vs_blk_8h	E-MTAB-1272	1655	296	27	5	7	19	18	14	53
	Hs_MVEC-C_hpx_vs_ct_24h	GSE11341	459	192	27	8	4	16	16	18	48
	Hs_MVEC-L_hpx_vs_ct_24h	GSE11341	504	198	22	4	4	15	13	14	42
	Hs_JurkT_ars_vs_ct	GSE46909	2381	157	15	5	3	3	9	8	31

Reverse causal reasoning (RCR) analysis was conducted using richness and concordance p values <0.05. HYPs that met both criteria were considered to be statistically significant. *Total number of unique human HYPs among the six subnetworks. Total is not 551 HYPs, the sum obtained from Figure 3, due to node/HYP overlap between the subnetworks. Mm_Ao_78w_ApoE_vs_wt (GSE10000) is included in this table as a building dataset reference. O/L: overlapping. Dataset Name Dataset ID N° of state changes (SC) N° of HYPs in dataset N° HYPs O/L with EC Activation N° HYPs O/L with Platelet Activation N° HYPs O/L with EC-Monocyte Interaction N° HYPs O/L with Foam Cell Formation N° HYPs O/L with SMC Activation N° HYPs O/L with Plaque Destabilization N° HYPs O/L with V-IPN.

Figure 4 depicts the coverage and OR of the overlap for each dataset onto each V-IPN subnetwork. Overall, datasets from immortalized cell lines in culture (Hs_EC_GFP_oxLDL_vs_ct; Hs_JurkT_ars_vs_ct) exhibited a lower coverage and OR compared to datasets from intact murine (Mm_Ao_78w_ApoE_vs_wt; Mm_Ao_16w_ApoE_CS_vs_sham) and human (Hs_athCA_vs_ctIMA) vascular tissues. (Mm_Ao_78w_ApoE_vs_wt), a model construction dataset, is included here only as a reference. Within the EC datasets, primary cells treated with Ox-PAPC (Hs_EC_oxPAP_vs_ct) showed the largest coverage compared to both Ox-LDL-treated HAEC datasets (Hs_EC_GFP_oxLDL_vs_ct; Hs_EC_LOX1_oxLDL_vs_ct) and MVECs (Hs_MVEC-C_hpx_vs_ct_24h; Hs_MVEC-L_hpx_vs_ct_24h). Within the HAEC study, a larger HYP coverage was shown by LOX-1-transfected ECs (Hs_EC_LOX1_oxLDL_vs_ct) compared to cells transfected with GFP (Hs_EC_GFP_oxLDL_vs_ct), suggesting that vascular inflammatory processes are indeed initiated by overexpressing LOX-1.

Hs_athCA_vs_ctIMA exhibited statistically significant ORs for all of the six subnetworks and the highest coverage of all datasets (except Mm_Ao_78w_ApoE_vs_wt, a model construction dataset) in four of the six subnetworks: Plaque Destabilization, Platelet Activation, EC-Monocyte Interaction and Foam Cell formation. This remarkable finding underlines the value of contrasting gene expression datasets from atherosclerotic arteries (e.g. coronary arteries) with normal vessels from the same subjects (e.g., internal mammary arteries) as performed for this dataset[25]. Interestingly, the murine dataset derived from ApoE^-/- aortas of old mice (Mm_Ao_78w_ApoE_vs_wt) displayed a very similar OR pattern to the human dataset across all subnetworks, with the exception of the Platelet Activation subnetwork, which indicates that largely similar biological pathways underlie atherosclerosis in both species (Figure 4). This result also suggests that the process of platelet activation may play a larger role in the development of atherosclerotic plaques in humans compared to advanced-age murine models of atherosclerosis.

HYP scoring and HYP directionality

The predicted directionality of all significant HYPs for all subnetworks, i.e. increased or decreased predictions based on the downstream gene expression is included in a supplementary file for each dataset examined ( Additional file7: Datasets_Analysis_Dashboards). A representative mapping of the Plaque Destabilization subnetwork depicted in Figure 5A shows bar plots for all possible HYPs that were predicted as significant in the Mm_Ao_16w_ApoE_CS_vs_sham, Mm_Ao_78w_ApoE_vs_wt, and Hs_athCA_vs_ctIMA datasets. This visualization highlights Hs_athCA_vs_ctIMA as the dataset exhibiting the largest HYP coverage and enrichment values compared to two murine datasets (Mm_Ao_16w_ApoE_CS_vs_sham and Mm_Ao_78w_ApoE_vs_wt). The subnetwork depiction also shows biological components being largely up-regulated as they relate to the process of plaque destabilization in human disease. A representative example of the gene expression data downstream of HYP taof(Stat1) scored for the Hs_athCA_vs_ctIMA and the Mm_Ao_78w_ApoE_vs_wt datasets is depicted in Figure 5B and C, respectively. The HYP contains 94 and 82 measured RNA abundance nodes and a total of 31 and 50 differentially expressed RNAs mapped to the network from the human and murine datasets, respectively. Twenty-five and 46 genes supporting upregulated activity and 6 and 4 supporting downregulated activity from the human and murine datasets, respectively. The proportion of RNA downstream nodes supporting an increased activity of STAT1 over those supporting a downregulated activity may be related to the key cell signaling role of STAT1 in cells embedded in advanced murine atherosclerotic lesions. Four additional examples of HYP scoring displaying the underlying gene expression data of HYPs for the Mm_Ao_78w_ApoE_vs_wt dataset were included in the supplement. Selected HYPs included “macrophage activation”, “Ccl5”, “monocyte adherence” and “kaof(Chuk)” (Additional file8: Figure S3a-d).

The HYP coverage of Hs_athCA_vs_ctIMA was examined not only across all V-IPN subnetworks, but also across all subnetworks and models we have previously published[16–20]. A bar plot visualization of all overlapping networks is depicted in Additional file9: Figure S4 as coverage (bar length) and OR (grey color intensity). Subnetworks displaying the largest coverage and highest ORs are within the IPN (Inflammatory Process Network), TRAG (Tissue Repair and Angiogenesis) and DACS (DNA damage, Autophagy, Cell death, and Senescence) models and include Dendritic Cell Migration, Neutrophil Chemotaxis, Natural Killer (NK) Cell Activation, Epithelial Cell Barrier Defense, Macrophage activation, Immune Regulation of Angiogenesis and MAP kinases (Mapk). Many of these subnetworks constitute biological processes that have also been implicated in the development of atherosclerotic lesions.

V-IPN evaluation of preclinical data translatability

V-IPN coverage of predicted HYPs from human and murine in vivo datasets

To test the power of the V-IPN at capturing mechanisms and biological pathways implicated in advanced vascular lesion development in human arteries, we evaluated coverage across the V-IPN subnetworks with a gene expression dataset from human coronary arteries isolated from CAD patients undergoing bypass surgery (Hs_athCA_vs_ctIMA)[25]. Lesion stage- and species-specific mechanistic differences were assessed by comparing the HYP coverage of each V-IPN subnetwork by the human dataset with two murine aortic datasets from early (Mm_Ao_16w_ApoE_CS_vs_sham) and advanced (Mm_Ao_78w_ApoE_vs_wt) atherosclerosis. Aortic tissue was collected from sexually mature 13–16 week-old ApoE^-/- mice following 30 days of CS exposure (Mm_Ao_16w_ApoE_CS_vs_sham) and from 78 week-old, unexposed ApoE^-/- mice (Mm_Ao_78w_ApoE_vs_wt). The Mm_Ao_78w_ApoE_vs_wt dataset was used in this evaluation as a reference for stage- and species-specific comparisons with the Mm_Ao_16w_ApoE_CS_vs_sham and human Hs_athCA_vs_ctIMA datasets, respectively.

V-IPN coverage of HYPs common to the murine and human datasets revealed that more HYPs were common between Mm_Ao_78w_ApoE_vs_wt and Hs_athCA_vs_ctIMA than between any other dataset pair comparison (Mm_Ao_16w_ApoE_CS_vs_sham // Hs_athCA_vs_ctIMA; Mm_Ao_16w_ApoE_CS_vs_sham // Mm_Ao_78w_ApoE_vs_wt) or the combination of the three datasets. The number of HYPs shared between Mm_Ao_78w_ApoE_vs_wt and Hs_athCA_vs_ctIMA was 2–5 times higher than any other possible comparison within the three datasets (Table 5). This observation prompted us to conduct murine-murine and murine-human dataset comparisons to further examine the molecules and pathways shared by these datasets. This approach allowed us to evaluate whether distinct species- and lesion stage-specific mechanisms may play a role in vascular lesion development.

Mm_Ao_16w_ApoE_CS_vs_sham	Mm_Ao_78w_ApoE_vs_wt	Hs_athCA_vs_ctIMA	No of HYPs	EC activation	Platelet activation	EC-monocyte interaction	Foam cell formation	SMC activation	Plaque destabilization
X	X	X	18	11 (7)	3 (5)	2 (7)	8 (7)	9 (10)	9 (9)
X	X	-	8	2 (1)	0 (0)	0 (0)	4 (3)	3 (3)	1 (1)
X	-	X	6	2 (1)	0 (0)	2 (7)	3 (3)	2 (2)	4 (4)
-	X	X	32	18 (12)	6 (10)	5 (17)	15 (13)	4 (5)	16 (16)
X	-	-	20	9 (6)	5 (8)	0 (0)	8 (7)	5 (6)	4 (4)
-	X	-	51	28 (18)	2 (3)	6 (21)	15 (13)	16 (18)	17 (17)
-	-	X	24	11 (7)	8 (14)	2 (7)	10 (8)	5 (6)	8 (8)
-	-	-	583	73 (47)	36 (61)	12 (41)	55 (47)	43 (49)	44 (43)

Numbers between parentheses are percentage of total number of possible HYPs.

Early vs advanced murine atherosclerosis datasets

The predicted HYPs from the murine datasets showed a low degree of overlapping HYP coverage (0-3%, Table 5), suggesting that distinct molecular pathways are active at various stages of atherogenesis in the same species. The Mm_Ao_16w_ApoE_CS_vs_sham dataset alone covered 0-7%. We sought to evaluate the coverage of this dataset across the V-IPN subnetworks to determine the mechanistic similarities and differences underlying lesion formation following CS-exposure when compared to advanced-age lesions in the same species. To determine the extent of vascular lesion development in the CS-exposure model, aortic histomorphometry and plasma lipid profiles were performed at the end of the study. Total cholesterol, LDL and VLDL from animals exposed to CS for 30 days were significantly increased, as was the size of atherosclerotic plaques in the aortic arch (Figure 6, panels A to D). H&E staining of cross-sections of the aortic roots showed typical intimal thickenings in areas close to the aortic valve leaflets; immunohistochemical staining of these lesions with a MAC3 antibody indicated that they were infiltrated by numerous macrophages (Figure 6, panels E to G).

Predicted HYPs that were common to the mouse datasets from young and old mice (Mm_Ao_16w_ApoE_CS_vs_sham and Mm_Ao_78w_ApoE_vs_wt) indicated a low degree of overlapping across the V-IPN subnetworks (0-3%, Table 5). Coverage analysis showed some mechanisms and molecules being activated solely in the Mm_Ao_16w_ApoE_CS_vs_sham dataset including PPARA, CD40LG, RAC1, PGE2 and SREBF2. Two predicted HYPs were found in four out of the six subnetworks (decreased PPARA and increased CD40LG) including the Plaque Destabilization subnetwork (Figure 5). AGTR1A (angiotensin II receptor 1A) was also a predicted HYP shared by both datasets in three subnetworks. These results indicate that in addition to a small set of mechanisms shared by early and advanced murine lesions, distinct biological pathways may also contribute to lesion formation in CS exposure-induced early vascular lesions compared to those implicated in advanced, older lesions.

Advanced murine vs. advanced human atherosclerosis datasets

In contrast to the murine dataset comparison results, the degree of overlapping HYPs between Mm_Ao_78w_ApoE_vs_wt and Hs_athCA_vs_ctIMA ranged from 10-16% (Table 5), except for the SMC Activation subnetwork (5%). Many significant HYPs were shared between these two datasets within the Plaque Destabilization subnetwork, (Figure 5A). Significant HYPs related to vascular pathobiology were mapped to the V-IPN; they included the processes of inflammation, angiogenesis, and monocyte and macrophage differentiation. Common HYPs included the canonical transcriptional regulators AP1, IRF3, REL and SPI1, nuclear receptors and signal transducers (e.g., NCOR1, STAT1), growth factors (e.g., VEGFA), and chemoattractants (e.g., CCL5). Cytokines involved in the differentiation and function of macrophages and lymphocytes CSF1, CSF2, IFNG, IL1B and IL6, were also shared between the two datasets. A gene functional clustering analysis of common HYPs using DAVID (http://david.abcc.ncifcrf.gov/home.jsp) revealed inflammation, cytokine activity, chemotaxis and the toll-like receptor signaling pathways as the top ranked functional categories (Additional file11: Table S3). This HYP coverage analysis indicates that a shared repertoire of biological mechanisms underlies the development of advanced vascular lesions in both humans and mice.

Early murine vs. advanced human atherosclerosis

A comparison between Mm_Ao_16w_ApoE_CS_vs_sham and Hs_athCA_vs_ctIMA revealed a low degree of overlapping HYPs, ranging from 0-7% spanning all subnetworks (Table 5). Furthermore, coverage analysis within the Plaque Destabilization subnetwork indicated only a few HYPs being shared between these two datasets (increased CCL2, increased “response to hypoxia”, and increased TGFB1) (Figure 5A). Only two additional HYPs, HIF1A and PPARD, were shared by both datasets in other subnetworks. This analysis suggests that despite a comparable atherosclerotic plaque morphology, the molecular events leading to its development in a CS exposure murine model are quite distinct from the molecular pathways leading to plaque formation and destabilization in advanced human lesions.

Pathobiological content of the V-IPN

Early mechanisms in vascular disease development such as arterial cell dysfunction are amenable to controlled experimental perturbations in animal models or in vitro settings. In contrast, advanced atherosclerotic lesions are challenging to recreate experimentally, which has led to a paucity of sound data on the precise cellular and molecular mechanisms leading to plaque instability and eventual rupture. Disease modelling approaches, such as the implementation of the V-IPN reported here, overcome these barriers by integrating current knowledge and large gene expression datasets into networks reflecting the pathobiology of interest. Each new dataset mapping on the network has the potential to expand our mechanistic knowledge and further refine the network’s structure and content.

Impaired endothelial production of prostacyclin (PGI2) and nitric oxide (NO) in early arterial lesions facilitate vasoconstriction, inflammation and oxidative stress at a time when no morphological changes in the vessel wall have occurred[29]. Cell adhesion molecules (CAMs) are well represented in the V-IPN as they have been extensively documented in the migration of inflammatory cells from the vascular lumen to the subendothelial space[30]. In the presence of hypercholesterolemia and a local excess of reactive oxygen species, oxidation of lipids and lipoproteins leads to activation of phagocytic cells that set in motion a cascade of inflammatory events including the release of growth and chemotactic factors as well as the migration and proliferation of SMCs[31]. Excess lipid uptake by macrophages promotes their differentiation into foam cells, which is the hallmark of fatty streaks observed in early atherosclerosis. We captured the signaling mechanisms responsible for these phenomena in the SMC Activation and Foam Cell Formation subnetworks, respectively. In humans, the final stage of atheroma formation is reached after years of continuous exposure to environmental noxious stimuli, such as cigarette smoke and a diet rich in saturated fats and poor in antioxidants[29, 32]. The natural disease progression results in plaque growth and positive vessel remodeling to maintain a functional lumen size. Although atherosclerotic plaques remain clinically silent for decades, they may evolve to become advanced lesions that are prone to calcification, cap thinning, hemorrhage and rupture. Platelets play a major role in the advanced stages of plaque development where neovascularization, thrombosis, plaque erosion and rupture constitute fatal complications[33]. Interestingly, platelets also participate in early atherogenic events by promoting EC activation and forming microthrombi in fatty streaks[34]. These platelet-related pathways have been captured in the Plaque Destabilization and Platelet Activation subnetworks. Thus, the pathobiology incorporated in the V-IPN represents a comprehensive implementation of our current knowledge of atherogenesis, which includes the primary cellular players, as well as the array of biological processes involved, ranging from EC dysfunction to the formation of fatty streaks, atheromas, and subsequent plaque instability.

Comparisons of RCR-based analyses of transcriptomics datasets from cells in culture and intact murine and human tissues rendered a series of powerful insights. The underlying molecular findings as well as the pathobiological relevance are described in detail below and summarized in Additional file12: Table S4.

The V-IPN captured predicted HYPs from primary HAECs and an immortalized HAEC line that may account for the divergent phenotypes exhibited by the two cell types

ECs uptake oxidized lipids and Ox-LDL through the scavenger receptor LOX-1[8]. Oxidized lipids induce inflammatory responses mediated by NFkB, including upregulation of cytokines, chemokines and CAMs that result in substantial leukocyte recruitment[35, 36]. Predicted HYPs from HAECs overexpressing LOX-1 or GFP (Hs_EC_LOX1_oxLDL_vs_ct; Hs_EC_GFP_oxLDL_vs_ct) differed significantly from the HYP profile of HAECs stimulated by proinflammatory oxidized phospholipid, Ox-PAPC (Hs_EC_oxPAP_vs_ct). A distinction between the three datasets is indicated by the relatively low number of predicted HYPs observed in Ox-LDL-treated HAECs; 5 and 14 HYPs for GFP- and LOX-1-transfected cells, respectively, compared to 40 significant HYPs predicted from Hs_EC_oxPAP_vs_ct, a dataset generated from HAECs treated with Ox-PAPC (Additional file10: Table S2). The marked differences in the number of significant predicted HYPs from each dataset could reflect differing magnitudes of signaling events elicited upon stimulating ECs with atherogenic lipids. Indeed, Ox-PAPC, a purified component of Ox-LDL, may have wider and more potent effects on EC biology compared to Ox-LDL[37]. The analysis presented herein could indicate that Ox-PAPC treatment is a more efficient process by which to induce the molecular features that most closely resemble chronic development of atherosclerotic plaques. Alternatively, the differences could be explained by the phenotypic differences between the cell types used and the number of human donors represented in each dataset. Cells from dataset Hs_EC_LOX1_oxLDL_vs_ct were obtained from an immortalized (SV40-induced) human aortic EC line[38], whereas primary HAEC cultures from 96 donors were used to generate Hs_EC_oxPAP_vs_ct expression data[39]. Cellular transformation of immortalized cell lines is associated with phenotypic changes at multiple levels including gene expression, biochemical, metabolic and proliferative capacity[40, 41]. Therefore, transformed cell lines may have more limited abilities to respond to experimental induction by atherogenic lipids.

Mechanisms identified by the V-IPN discriminate between early and late vascular lesions in ApoE^-/- mice and highlight the commonalities between advanced murine and human atherosclerotic lesions

We have conducted RCR analysis on gene expression datasets from a study utilizing the ApoE^-/- mouse strain, a well-established model of atherosclerosis[42]. The minimal overlap of common HYPs (6 total, Table 5) between aortas of young (Mm_Ao_16w_ApoE_CS_vs_sham) and old (Mm_Ao_78w_ApoE_vs_wt) adult ApoE^-/- mice demonstrates a substantial divergence of atherogenic processes taking place in young ApoE^-/- mice (16-week old) exposed to CS for 30 days compared to advanced atherosclerotic disease observed in older ApoE^-/- mice (78 weeks). Strikingly, predicted HYPs from human (Hs_athCA_vs_ctIMA) and old murine (Mm_Ao_78w_ApoE_vs_wt) datasets generated from intact arteries harboring advanced atherosclerotic lesions shared a significantly larger set of causal mechanisms (32 total HYPs, Table 5) indicating that similar mechanisms underlie the development of advanced atherosclerotic lesions in both species. The pathobiological picture that emerged when overlaying predicted HYPs from both datasets onto the V-IPN is one where a complex series of proliferative, apoptotic and inflammatory events driven by intracellular transducers and transcription regulators are all taking place simultaneously. A functional clustering analysis of the common HYPs from both species using DAVID revealed functional categories consistent with the mechanisms described above (Additional file11: Table S3). The vast majority of this set of commonly mapped HYPs, which were predicted increased in the human (Hs_athCA_vs_ctIMA) and old murine (Mm_Ao_78w_ApoE_vs_wt) datasets, should serve as an initial reference for future simulations using murine and human vascular datasets. Our results highlight the power of the V-IPN to assess, at the molecular level, the degree of translatability from murine morphologic[43] and tomographic data[44] on plaque instability[45] to the human clinical setting using transcriptomics data.

V-IPN evaluation reveals that distinct molecular pathways contribute to atherogenesis in a CS-exposure murine model of disease

Among the predicted HYPs from the mouse Mm_Ao_16w_ApoE_CS_vs_sham and Mm_Ao_78w_ApoE_vs_wt datasets, coverage analysis across the V-IPN revealed a remarkably low degree of overlapping mechanisms across the V-IPN subnetworks (0-3%, Table 5). This result suggests that distinct biological pathways contribute to early lesion development in ApoE^-/- mice exposed to CS. A comparison between Mm_Ao_16w_ApoE_CS_vs_sham and Hs_athCA_vs_ctIMA also revealed a low extent of overlapping HYPs, ranging from 0-7% spanning all subnetworks (Table 5). Furthermore, broad coverage analysis of the predicted HYPs indicated that the Mm_Ao_16w_ApoE_CS_vs_sham dataset contained 52 HYPs found across the V-IPN, whereas the Hs_athCA_vs_ctIMA dataset contained 80 HYPs (Table 4). This result suggests that a more diverse array of biological mechanisms underlie advanced atherosclerotic lesion development in the human disease. Indeed, coverage analysis specifically within the Plaque Destabilization subnetwork indicated only a few unique mechanisms being activated in the Mm_Ao_16w_ApoE_CS_vs_sham samples (e.g. decreased PPARA and increased CD40LG) (Figure 5). A full characterization of discrete cellular events distinguishing acute atherogenesis in a murine model from advanced-stage murine and human disease further demonstrates the utility of the V-IPN in evaluating novel datasets to better understand comparability between species and/or disease model systems.

V-IPN-generated causal paths were consistent with advanced, unstable lesions, in coronary atherosclerotic arteries of CAD patients

Unlike all other datasets used for model construction and evaluation, Hs_athCA_vs_ctIMA represents paired gene expression profiles from 37 human atherosclerotic coronary arteries and control internal mammary arteries (IMA); each pair was obtained from the same subject (STAGE study)[25]. A principal component analysis (PCA) of the gene expression profiles identified two groups of control and diseased arteries (Additional file13: Figure S5). All pairs were clearly clustered, thus providing a sense of how distinct the gene profile of advanced coronary atherosclerosis was, compared to the reference IMA. Pathological evidence indicates that unstable plaques are complex lesions with common morphological features including a thin fibrous cap, a large lipid core, a network of vasa vasorum and focal inflammation[46]. HYP coverage analysis of the human dataset (Hs_athCA_vs_ctIMA) was distinctly consistent with mechanisms driving the morphological features of unstable plaques. Remarkably, the Plaque Destabilization and Platelet Activation subnetworks exhibited the largest ORs for HYP overlap across all datasets examined and through the six V-IPN subnetworks (Figure 4). The ORs for the Hs_athCA_vs_ctIMA dataset were even higher than those obtained for Mm_Ao_78w_ApoE_vs_wt; in the Platelet Activation (2.9 vs. 1.1) and Plaque Destabilization subnetworks (4.1 vs. 3.4); Mm_Ao_78w_ApoE_vs_wt is a dataset from advanced murine aortic lesions used for network construction. This analysis highlights model-level coverage of additional mechanisms unique to the human condition and validates the strength of the V-IPN to differentiate expression data derived from advanced human and murine lesions. Consistent with a decreased activity of SMCs leading to fibrous cap thinning in lesions that are prone to rupture, the SMC Activation network exhibited lower coverage and ORs compared to both murine datasets, Mm_Ao_16w_ApoE_CS_vs_sham and Mm_Ao_78w_ApoE_vs_wt (2.1 vs. 2.7 and 2.8) (Figure 4). Significant HYPs categorized within pathobiological functions linked to plaque destabilization are described below.

Model limitations

RCR-based models do not operate with integrated feedbacks and non-linear elements that contribute to regulating a dynamic output. In order to make accurate predictions, mathematical models incorporate feedback elements tuned to match phenotypic constrains (e.g., blood pressure values). In RCR, all biological feedbacks are implicitly integrated in the datasets. RCR-based models do not dynamically model the regulatory processes controlling biological pathways. RCR-based models are tools to extract biological processes embedded in large sets of molecular data driven by specific experimental perturbations, and to contextualize those findings within a body of knowledge.

Glossary

Reverse Causal Reasoning (RCR): A computational methodology for identifying potential upstream controllers leading to differential molecular profiles.

Selventa Knowledgebase (SK): A network representing a working set of knowledge fit for a specified use. The SK is used as a substrate for RCR. It encodes prior scientific knowledge as a network of nodes that are connected by edges.

Biological Expression Language (BEL): The knowledge representation language used to build the SK.

Node: A biological entity or process in the SK.

Edge A causal relationship (e.g., increase, decrease, subset) connecting two nodes in the knowledgebase.

State Change (SC): A differential measurement across a sample group (e.g., treated and control) that is converted to a discrete value of increase, decrease, or no change, based on two statistical metrics: richness and concordance.

Hypothesis (HYP): A small, directed causal network containing an upstream node representing a biological entity or process connected by a causal increase, decrease or ambiguous edges to downstream nodes representing measured entities.

HYP upstream node: A controller of downstream nodes in a HYP and a potential explanation for state changes (SC) mapped to the downstream nodes.

HYP downstream nodes: Nodes in a HYP mapped to quantities measured in the dataset.

HYP causal edges: The causal relationships (i.e., increases or decreases) connecting the HYP upstream node to each downstream node.

Richness: A measure of the relevance of a HYP to the changes observed in an experimental dataset.

Concordance: A measure of the consistency of the direction of the changes observed in an experimental dataset.

Coverage (sensitivity): An estimate of the fraction of possible HYPs in a subnetwork that are significant in a dataset. Coverage is a measure of HYP enrichment.

Odds ratio (OR): The probability of having significant dataset HYPs within a network. The higher the OR, the better the network encompasses the biology embedded in a given dataset.

X: Protein abundance of X

taof(X): Transcriptional activity of X

exp(X): RNA expression of X

gtpof(X): GTP-bound activity of X

kaof(X): Kinase activity of X

paof(X): Phosphatase activity of X

catof(X): Catalytic activity of X

X P@Y: Abundance of X phosphorylated at Y