Deep Learning Analyses to Delineate the Molecular Remodeling Process after Myocardial Infarction

Specific proteins and processes have been identified in post-myocardial infarction (MI) pathological remodeling, but a comprehensive understanding of the complete molecular evolution is lacking. We generated microarray data from swine heart biopsies at baseline and 6, 30, and 45 days after infarction to feed machine-learning algorithms. We cross-validated the results using available clinical and experimental information. MI progression was accompanied by the regulation of adipogenesis, fatty acid metabolism, and epithelial–mesenchymal transition. The infarct core region was enriched in processes related to muscle contraction and membrane depolarization. Angiogenesis was among the first morphogenic responses detected as being sustained over time, but other processes suggesting post-ischemic recapitulation of embryogenic processes were also observed. Finally, protein-triggering analysis established the key genes mediating each process at each time point, as well as the complete adverse remodeling response. We modeled the behaviors of these genes, generating a description of the integrative mechanism of action for MI progression. This mechanistic analysis overlapped at different time points; the common pathways between the source proteins and cardiac remodeling involved IGF1R, RAF1, KPCA, JUN, and PTN11 as modulators. Thus, our data delineate a structured and comprehensive picture of the molecular remodeling process, identify new potential biomarkers or therapeutic targets, and establish therapeutic windows during disease progression.


Introduction
Myocardial infarction (MI) occurs when blood flow suddenly stops due to occlusion of a coronary artery, leading to local ischemia in the heart [1][2][3]. This condition triggers an adverse myocardial remodeling response involving a wide variety of signaling pathways, potentially including extracellular matrix (ECM) dysregulation, cardiomyocyte apoptosis, and cardiogenic processes, such as myocyte concentric and eccentric hypertrophy, slippage, accumulation of interstitial tissue, a molecular shift towards an embryonic pattern, and loss of cardiac energy reserves [4,5]. Cardiac remodeling secondary to MI is a well-defined process that may ultimately lead to heart failure (HF) [6][7][8], but the specific dynamic molecular mechanisms underlying this progression have not been fully characterized. This process has been examined in several studies using small experimental animal models, mainly providing data about individual genes or proteins [9][10][11][12]. Although these investigations have provided valuable information regarding cardiac remodeling, they have been insufficient for capturing this complex process as a whole.
Increasing evidence supports the idea that specific biological processes (e.g., proteinprotein interactions or epigenetic regulation) are likely influenced by the biological context, such as a specific tissue or a certain disease [13][14][15][16]. Compared to small animal models, large experimental animal models more closely resemble human cardiac physiology, function, and anatomy; therefore, their use constitutes a key step in translating experimentally obtained information into clinical use. In addition, vast amounts of data are constantly being generated, and compiling, analyzing, and interpreting this information as a whole constitutes an overwhelming task. For this purpose, new technologies are rapidly emerging that combine different engineering approaches and bioinformatics. Within this context, systems biology arose as an interdisciplinary field of study with the aim of unravelling the key interactions within complex biological networks following a holistic approach based on computational and mathematical models.
Using transcriptomic information from infarcted and remote in vivo swine heart tissues, we extensively analyzed the temporal and region-specific myocardial gene expression patterns in response to MI [17]. To determine the underlying molecular cause of the observed cardiac remodeling, we incorporated this transcriptomic data into a deep learning model, the Therapeutic Performance Mapping System (TPMS ® ) ( Figure 1) [18]. TPMS ® applies artificial intelligence and pattern recognition techniques to combine interactomic data with the molecular and clinical responses observed in patients. First, the interactomic information is used to generate a skeleton for computational models, which act as a network of potential mechanistic interactions. Second, this network is fitted into a deep learning model and trained using clinical and molecular responses observed in patients. This ultimately results in the generation of a mathematical model capable of both reproducing existing knowledge and discerning the mechanisms of action hidden under thousands of molecular interactions that are otherwise inaccessible.
In the present study, we incorporated transcriptomic data obtained from infarcted porcine hearts at different time points into their corresponding mathematical configurations to gain insights into the molecular evolution of MI over time and its development into HF. (1) First, we characterized myocardial infarction (MI) at the molecular level via manual curation of the literature and using a compendium of massive public databases describing the molecular interactions of interest. (2) In parallel, we used experimental transcriptomic data to define the molecular behavior of MI at each time point. (3) We used the information to generate a map of proteins regulated by the obtained experimental data. (4) We fed the mathematical models to identify patterns in the data, (5) identify key regulatory proteins, and (6) predict new mechanisms of action. (7) The final result included a scientific justification of every prediction. ANN: artificial neural network.
In the present study, we incorporated transcriptomic data obtained from infarcted porcine hearts at different time points into their corresponding mathematical configurations to gain insights into the molecular evolution of MI over time and its development into HF.

Myocardial Infarction Model
Female Landrace × Large White pigs (n = 9) were premedicated with intramuscular (IM) azaperone (10 mg/kg) (Stresnil ® , Laboratorios Esteve, Barcelona, Spain), and then administered intravenous (IV) pentobarbital sodium (15 mg/kg) (Tiobarbital ® 1 g, B. Braun, Melsungen, Germany). These pigs underwent endotracheal intubation with 2% inhaled isoflurane as the anesthesia. During the procedure, IV fentanyl (0.75 mg/kg/45 min) (Fentanest ® , Kern Pharma, Madrid, Spain) was administered as an analgesic, and IV atracurium besylate (1.5 mg/kg bolus) (Tracrium ® , GlaxoSmithKline, Brentford, UK) was used (1) First, we characterized myocardial infarction (MI) at the molecular level via manual curation of the literature and using a compendium of massive public databases describing the molecular interactions of interest. (2) In parallel, we used experimental transcriptomic data to define the molecular behavior of MI at each time point. (3) We used the information to generate a map of proteins regulated by the obtained experimental data. (4) We fed the mathematical models to identify patterns in the data, (5) identify key regulatory proteins, and (6) predict new mechanisms of action. (7) The final result included a scientific justification of every prediction. ANN: artificial neural network.

Myocardial Infarction Model
Female Landrace × Large White pigs (n = 9) were premedicated with intramuscular (IM) azaperone (10 mg/kg) (Stresnil ® , Laboratorios Esteve, Barcelona, Spain), and then administered intravenous (IV) pentobarbital sodium (15 mg/kg) (Tiobarbital ® 1 g, B. Braun, Melsungen, Germany). These pigs underwent endotracheal intubation with 2% inhaled isoflurane as the anesthesia. During the procedure, IV fentanyl (0.75 mg/kg/45 min) (Fentanest ® , Kern Pharma, Madrid, Spain) was administered as an analgesic, and IV atracurium besylate (1.5 mg/kg bolus) (Tracrium ® , GlaxoSmithKline, Brentford, UK) was used to induce muscular relaxation. After a left lateral thoracotomy, MI was induced via permanent ligation of the circumflex artery as described previously [19]. IM tulathromycin (2.5 mg/kg) (Draxxin ® , Pfizer Animal Health, New York, NY, USA) was administered as an antibiotic prophylactic, and a transdermal fentanyl patch was applied to facilitate postoperative analgesic care. These surgical interventions were monitored by electrocardiogram (ECG), capnography, and pulse oximetry, and with noninvasive measurement of arterial blood pressure and temperature. Six days (n = 3), 30 days (n = 3), or 45 days (n = 3) after MI, animals were randomly sacrificed by an IV overdose of potassium chloride solution. At each temporal stage, we analyzed three paired myocardial samples from the infarct core and non-infarcted remote myocardium. As a physiological control, we analyzed myocardial samples from healthy animals (n = 3).

Tissue Collection and RNA Extraction
Following sternotomy, hearts were washed in ice-cold buffered saline solution to remove blood residue. Biopsies were obtained from the infarct core (center of the scar on the left ventricle), remote myocardium (non-infarcted interventricular septum), and control myocardium (healthy animals). To ensure RNA stabilization, biopsies were preserved in Allprotect Tissue Reagent (Qiagen, Hilden, Germany) at room temperature. Total RNA was isolated using the RNeasy Fibrous Tissue Mini Kit (Qiagen, Hilden, Germany). RNA purity and integrity were assessed by spectrophotometry (NanoDrop ND-1000, NanoDrop Technologies, Thermo Fisher Scientific, Waltham, MA, USA) and nanoelectrophoresis (2100 Bioanalyzer, Agilent Technologies, Santa Clara, CA, USA).

Microarray Gene Expression Analysis
We obtained microarray expression profiles using the GeneChip ® Porcine Genome Array (Affymetrix). We processed 200 ng of total RNA from each sample, labeling, fragmenting, and hybridizing it to the GeneChip ® following the manufacturer's instructions. Arrays were scanned using an Affymetrix GeneChip ® Scanner 3000 7 G. The raw expression data were preprocessed using the robust multichip average (RMA) normalization method [20], yielding 24,123 probe sets on a log2 basis.

Compilation of Transcriptomic Data
First, the data were filtered to discard all entries with contradictory information (i.e., two entries for the same gene name with negative and positive ratio values) and to identify the number of uniquely altered genes. Next, the swine transcriptomics were translated into their human equivalents via reciprocal best hits (RBHs) with BLAST and Gene Name Correspondence. Pig-to-human RBHs were identified using the InParanoid database [21]. We only included genes with an associated protein product. Gene information was mapped one-to-one to its protein product for its introduction into the protein network. For each gene, the correlation between RNA and protein levels was assessed using an RNA-to-protein ratio, as described previously [22]. If RBHs were not found for a protein, we used the reviewed UniProt entry for the human protein with a matching gene name (Supplementary  Table S1). Next, the proteins were labeled according to whether a protein is activated or inhibited under physiological conditions, and we used this information as a reference for restricting detection to the variability from these values, as defined by p-values adjusted using the Benjamini-Hochberg procedure to control the false discovery rate (FDR) at 0.01 [23,24]. Finally, proteins with human UniProt IDs within each cohort were used as molecular restrictions for our models.

Molecular Characterization of Pathology
We integrated published data with our results to define a set of molecular profiles characterizing MI and adverse myocardial remodeling, which were used to build the protein network and mathematical model.
We carefully and extensively reviewed articles in the PubMed database (abstract or fulltext depending on the inclusion of information on the molecular definition of disease) that included the following search strings: myocardial infarction-"myocardial infarction" "molecular"[All Fields] AND English [lang]). If the evidence for candidate proteins was not enough to include them on the effector list, specific searches for the protein were performed to assess its validity. We further characterized the pathophysiological processes at the protein level using a set of 202 unique proteins (121 related to MI and 136 to cardiac remodeling, with some overlap) that centered our analysis on the pathological conditions of interest in the human biological network.

Therapeutic Performance Mapping System (TPMS) Generation of Mathematical Models
We studied the pathology at 6, 30, and 45 days of evolution in samples from the infarct core (C6, C30, and C45, respectively), which included 4737, 4730, and 4203 differentially expressed proteins, respectively. We also studied samples from the remote area surrounding the infarct localization at the same time points (R6, R30, and R45, respectively), which included 122, 117, and 21 differentially expressed proteins, respectively. Using these data, we identified the main pathophysiological processes that were altered in MI and cardiac remodeling.
We trained the mathematical models using a large collection of well-established pathophysiological signals and clinical information relevant to the examined pathologies. Molecular descriptions of these input-output signals were obtained primarily via manual literature mining and from a compendium of massive databases that accumulate biological and clinical knowledge (e.g., KEGG, MINT, REACTOME, BIOGRID, and DrugBank). The model inputs included information on drugs, pathologies, and protein/gene relationships that could inhibit or activate one or more nodes of the protein network (their targets), triggering a perturbation through the system. The model outputs were experimental microarray data regarding upregulated or downregulated genes/proteins and clinical information. Using the model's input and output data, mathematical algorithms were generated to trace a change or perturbation from one to another, elucidating the mechanism of action explaining this connection.
The generated collection of known input-output physiological signals can be envisioned as a list of physiological principles that are representative of all humans or of particular pathophysiological conditions. These sets of rules are collated to form a "truth table" that every constructed mathematical model must satisfy. Transcriptomic information was compiled, analyzed, and included in the models when it met the reliability requirements described above.
The constructed models were used to determine the weight of the relative value of each protein (node) relationship. However, the very high number of links exponentially increased the number of parameters to solve. To optimize the system, we used two different approaches: one based on randomized systems (Monte Carlo-based system) [25], and the other based on information derived from the network topology [26].

Analyzing and Solving the Mathematical Models
The TPMS technology employs two different but complementary strategies to solve mathematical models. The artificial neural network (ANN) strategy can identify relationships among regions of the network (generalization). This strategy provides a predictive value that infers the probability that a specific relationship exists between two or more sets of proteins. In this case, we tested each differential protein against the described cardiac remodeling signature. Next, the model's predictive capacity is cross-validated using different sets of data towards what is described in the literature and databases.
The second strategy is the sampling method, which allows the observed effects to be traced back to specific molecules or drugs. Sampling methods are only applied once a key region of the protein map has been identified using ANNs or is suggested by experimental work. Once a response (indication, adverse effects, specific clinical improvements, etc.) has been identified and linked to a specific stimulus (compound, pathology, protein expression alteration, etc.) using ANNs, sampling methods enable analysis of the mechanism of action and elucidation of the hidden relationships between them.

Biomarker Identification
Using the characterization of MI, we performed multivariate analysis to evaluate combinations of two, three, or five proteins and identify the combinations that best classified either the solutions of the models or the microarray experiments. To assess the probability of correctly predicting and classifying the set of samples of high-throughput data, we used a "leave one out" strategy, which involved analyzing a subset of samples that were not previously included in the models.
The combinations were obtained via four strategies. To classify the solutions of the models, we used a model-based strategy and a model/HT-based strategy. Both strategies identified the protein combinations that better classified the solutions of the models to their corresponding cohort. The model/HT-based strategy additionally filtered the protein combinations using the high-throughput data (i.e., what was measured in the microarray), reducing the number of proteins and facilitating the identification of a combination with a higher generalization capability.
To classify the microarray experiments, we used an HT-based strategy and an HT/ model-based strategy. Both strategies identified the combinations of proteins that better classified the microarray experiments to their corresponding cohort. The HT-based strategy did not use the information from the disease models, whereas the HT/model-based strategy filtered the combinations by the proteins found to be relevant in the models.
Using these strategies, we generated a list of 106 individual proteins that represent potential biomarkers of the description of the pathology at a given time point (Supplementary  Table S2).

Generating A Molecular Model of Cardiac Remodeling
Starting from the differential proteins for each cohort, we performed a relationship analysis with ANNs to identify proteins mechanistically related to cardiac remodeling. Proteins having a predictive value > 50% (corresponding to~20% probability of accepting false positives) were included for further evaluation (Table 1). Upon identifying the differential proteins more closely related to cardiac remodeling, two analyses were performed. A gene set enrichment analysis (GSEA) based on functional annotation of the differentially expressed genes was used to identify relevant biological processes (KEGG pathways, Gene Ontology (GO) function, and GO process terms; Supplementary Table S3) [27,28]. We also applied triggering-protein analysis to identify the key proteins in a network that strongly affect the rest of the dataset (i.e., the proteins triggering the observed output). These proteins were further investigated and sampling methods used to explore their exact relationship to cardiac remodeling.

Data Integration
We integrated the in silico analysis of the microarray gene expression data to identify relationships across signature datasets, providing an independent assessment and functional validation. To cross-validate our findings, we used several complementary methods for complex data integration. The Metascape bioinformatics tool [29] was used for pathway analyses, the Perseus software platform (V1.6.5.0) [30] for high-dimensional omics data analysis, and the STRING [31] online tool to explore relevant protein-protein interactions (PPI) at 0.9 confidence.

High-Dimensional Data Integration of Microarray Gene Expression during 6 Weeks of MI Progression
Each distinct tissue dataset (C6, C30, C45, R6, R30, and R45) revealed substantial changes over 6 weeks of MI progression. To determine whether these changes were consistent across data types and to validate their time-dependent evolution in functional pathways, we adopted three independent strategies.
First, we sought to assess the internal structure of the data in a way that best explained its variance. Due to the high complexity of the datasets, we used the Perseus platform to perform principal component analysis (PCA) for the reduction of dimensionality. This approach clearly separated the core infarcted area from the remote area and the control group, which clustered together ( Figure 2).

Data Integration
We integrated the in silico analysis of the microarray gene expression data to identi relationships across signature datasets, providing an independent assessment and fun tional validation. To cross-validate our findings, we used several complementary met ods for complex data integration. The Metascape bioinformatics tool [29] was used f pathway analyses, the Perseus software platform (V1.6.5.0) [30] for high-dimensional om ics data analysis, and the STRING [31] online tool to explore relevant protein-protein i teractions (PPI) at 0.9 confidence.

High-Dimensional Data Integration of Microarray Gene Expression during 6 Weeks of MI Progression
Each distinct tissue dataset (C6, C30, C45, R6, R30, and R45) revealed substanti changes over 6 weeks of MI progression. To determine whether these changes were co sistent across data types and to validate their time-dependent evolution in function pathways, we adopted three independent strategies.
First, we sought to assess the internal structure of the data in a way that best e plained its variance. Due to the high complexity of the datasets, we used the Perseus pla form to perform principal component analysis (PCA) for the reduction of dimensionalit This approach clearly separated the core infarcted area from the remote area and the co trol group, which clustered together (Figure 2). Second, we employed Metascape algorithms to sort out the most affected pathwa in the gene set according to their GO terms ( Figure 3A). We were able to further clust them according to their biological significance and associated p-values ( Figure 3B,C These data suggested that MI progression was accompanied by a great alteration in mit chondrial metabolism, as well as in the regulation of adipogenesis, fatty acid metabolism and epithelial-mesenchymal transition processes. Second, we employed Metascape algorithms to sort out the most affected pathways in the gene set according to their GO terms ( Figure 3A). We were able to further cluster them according to their biological significance and associated p-values ( Figure 3B,C). These data suggested that MI progression was accompanied by a great alteration in mitochondrial metabolism, as well as in the regulation of adipogenesis, fatty acid metabolism, and epithelial-mesenchymal transition processes.
Finally, GSEA and overlapping tests allowed the contextualization of this information and distribution of the GO terms across core vs. control, and remote vs. control datasets. This revealed intrinsic differences between the core and remote regions of the heart after MI (Figure 4A-C). Specifically, the core region exhibited greater deterioration of cardiac contraction, depolarization during action potentials, and of heart rate regulation ( Figure 4A, upper panels), whereas the remote region exhibited a heavy enrichment of cellular respiration processes ( Figure 4A, lower panels).
Notably, seven proteins overlapped between the different time points regardless of the myocardial region (Table 2). Of the 105 proteins, 52 (49.5%) were recognized by both the models and the transcriptomic data analysis (model/HT or HT/model identification method). As these proteins were identified in models of both the pathology and the experimental data, they were considered more robust. Moreover, we found two distinct combinations of proteins from the infarct core data at each time point, which had generalization capabilities and accuracies of 100% (p ≤ 0.05; Table 3). No reliable combinations were found using remote myocardium data. Finally, GSEA and overlapping tests allowed the contextualization of this information and distribution of the GO terms across core vs. control, and remote vs. control datasets. This revealed intrinsic differences between the core and remote regions of the heart after MI ( Figure 4A-C). Specifically, the core region exhibited greater deterioration of cardiac contraction, depolarization during action potentials, and of heart rate regulation ( Figure 4A, upper panels), whereas the remote region exhibited a heavy enrichment of cellular respiration processes ( Figure 4A, lower panels).

Time-Dependent Identification of MI-Derived Biomarkers Strongly Related to Cardiac Remodeling
To specifically investigate the adverse cardiac remodeling processes that occur after MI, we integrated the high-throughput microarray data with molecular information available in the literature. Analysis of all tissues (infarcted core and remote) and time points (6, 30, and 45 days) identified 105 proteins ( Figure 5B; Supplementary Table S2). Of these,

O75676
KS6A4 Models: the combinations of proteins that better classify the solution of the model to the corresponding cohort; models/HT: the combinations of proteins that better classify the solution of the model to the corresponding cohort filtered by the proteins acting according to the high-throughput (HT) data; HT/models: the combinations of proteins that better classify the microarray experiments to the corresponding cohort filtered by the proteins relevant in the models; 1: detected/positive; 0: negative; −: not detected.
Models: the combinations of proteins that better classify the solution of the model to the corresponding cohort; models/HT: the combinations of proteins that better classify the solution of the model to the corresponding cohort filtered by the proteins acting according to the highthroughput (HT) data; HT/models: the combinations of proteins that better classify the microarray experiments to the corresponding cohort filtered by the proteins relevant in the models; 1: detected/positive; 0: negative; −: not detected. a Indicates whether the protein is secreted and/or has been detected in plasma. b Indicates if there is a previously known relationship between the protein and MI. Models/HT: the combinations of proteins that better classify the solution of the model to the corresponding cohort filtered by the proteins according to the high-throughput (HT) data; MI: myocardial infarction; C6: infarct core area 6 days after infarction; C30: infarct core area 30 days after infarction; C45: infarct core area 45 days after infarction.

Description of the Molecular Mechanistic Relationships Defining MI Evolution
Using the differentially expressed proteins found to be related to cardiac remodeling in the ANN analysis, we performed enrichment analysis to identify and understand the mechanism of action of cardiac remodeling in the infarct core at each time point. We identified a total of 355 processes that were altered in the infarct core region for at least one time point (221 upregulated and 134 downregulated). As shown in Figure 5A and Supplementary Table S3, there was a substantial overlap between the processes altered in the infarct core region at all three time points (C6, C30, and C45), including 8 KEGG pathways, 21 GO function terms, and 72 GO process terms. Other processes were specifically altered at one time point or two time points-mainly during the C6-C30 and C30-C45 transitions, although six GO process terms overlapped between C6 and C45.
Importantly, we found that 101 processes were altered in the infarct core throughout the progression to 45 days postinfarction. Most of the upregulated processes were related to ECM formation (organization, assembly, and adhesion) and to ECM components (hyaluronic, heparin, collagen, and glycosaminoglycan-related processes). We also detected alterations in processes related to muscle contraction, with an upregulation of actin-and calcium-related processes and downregulation of cardiac muscle contraction. Most of the downregulated processes were related to metabolism.
Some of the altered processes exhibited time-dependent behavior. Inflammation clearly appeared to be upregulated 6 days post-MI and decreased in later stages. Processes related to angiogenesis, cell proliferation, and morphogenesis started to increase 30 days post-MI and become more relevant by 45 days post-MI. Apoptosis also appeared to gain importance in the late stages. The upregulated pathways at 45 days suggested some normalization of metabolism, though several processes remained altered ( Figure 6). Cells 2021, 10, x FOR PEER REVIEW 14 of 20 importance in the late stages. The upregulated pathways at 45 days suggested some normalization of metabolism, though several processes remained altered ( Figure 6).

Identification of Key Proteins driving Post-MI Alterations
After the mechanistically affected pathways were identified, triggering-protein analysis allowed us to determine the mechanistic relationship between paired sets of proteins. Each analyzed time point presented a list of particular source proteins (Table 4) with some degree of overlap between the different cohorts. Overall, the transcriptomic experiments in the infarct core region identified 20 proteins (primarily involved in growth factor signaling) as the main drivers, or source proteins, of post-MI alterations. KPCA, PTN11, JAK2, 1433B, ERBB2, VEGFA, JUN, RAF1, and IGF1R were identified as source proteins

Identification of Key Proteins Driving Post-MI Alterations
After the mechanistically affected pathways were identified, triggering-protein analysis allowed us to determine the mechanistic relationship between paired sets of proteins. Each analyzed time point presented a list of particular source proteins (Table 4) with some degree of overlap between the different cohorts. Overall, the transcriptomic experiments in the infarct core region identified 20 proteins (primarily involved in growth factor signaling) as the main drivers, or source proteins, of post-MI alterations. KPCA, PTN11, JAK2, 1433B, ERBB2, VEGFA, JUN, RAF1, and IGF1R were identified as source proteins at all three time points. MAPK03, SRC, STAT1, NGF, RAC1, and TF65 were identified as common source proteins in C6 and C30. CBLB was common to C6 and C45, though its activity changed from day 6 (activated) to day 45 (inhibited). GCR, CDC42, and GRB2 were identified as source proteins only in C30, whereas MAPK01 was found specifically in C45. Table 4. Identification of time-dependent source proteins that explain the molecular mechanism of action in myocardial infarction.

Time-Points of Protein Expression in the Core MI Area C6 C30 C45
Proteins available for analysis a 1182 1150 1033 Maximum % of proteins explained by other differential proteins b 49% 48% 48% Number of source proteins c 16 18 11 % of explainable proteins explained by triggering proteins d 92% 93% 89% Proteins were identified by triggering-protein analysis. a The number of differential proteins related to cardiac remodeling available to be analyzed for each time point. b The percentage of evaluable proteins that were linked to another of the evaluated proteins (i.e., explainable by other proteins). c The number of source proteins selected. d The percentage of explainable proteins (i.e., proteins in the second row) that are explained by or linked to the selected source proteins. C: infarcted core region.

Unraveling the Molecular Mechanism of Action Relating the Source Proteins to Cardiac Remodeling
Mechanistic analysis of the relationship between the identified source proteins and cardiac remodeling revealed an overlap between the different time points. The common pathways between the source proteins and cardiac remodeling (Figure 7) involved IGF1R, RAF1, KPCA, JUN, and PTN11. Analysis of the downstream signaling revealed a mechanism of action potentially involving 18 proteins, 15 of which were differentially expressed during at least one time point (Table 5). The ERK kinases MAPK01 and MAPK03 were apparently activated by PTPN11 and PKCA, respectively, and both were upregulated at all three time points. The other upstream proteins (IGF1R and RAF1) were inhibited. In addition, PTPN11 was functionally inhibited by MAPK03, resulting in a reduction in MAPK01 activity.
With regard to cellular proliferation, survival, and differentiation, we found that TSC2 was upregulated in C45. We also detected heavy upregulation of the ECM remodeling protein MMP14 and of BAD and P53, which are key components of apoptosis and survival. Accordingly, at all three time points, we detected downregulation of MTOR, SL9A1 (hypertrophy-related), and THB (impaired myocyte contractility).

Discussion
Our first goal was to assess the data structure of the transcriptomic information obtained in vivo. We integrated this highly dimensional data and analyzed its composition, revealing that MI progression was accompanied by a regulation of adipogenesis, fatty acid metabolism, and epithelial-mesenchymal transition processes. Moreover, this analysis uncovered robust segregation between the infarct core samples versus the remote myocardium and controls, indicating major differences between these areas at the regulatory level. We specifically investigated the biological significance underlying these structural variations, and found that the infarct core region of the heart was enriched in affected processes relating to heart muscle contraction and membrane depolarization. In contrast, the remote myocardium was heavily affected by respiratory metabolism deficiencies. These results helped contextualize the MI.
Adverse ventricular remodeling is a hallmark of MI evolution and progression towards HF [32][33][34]; therefore, we investigated which of the identified mechanisms were related to adverse cardiac remodeling. Using machine learning techniques to analyze all tissues and time points, we identified 105 altered genes that are most likely related to cardiac remodeling in MI. This provided a list of potential biomarkers that may indicate the pathological evolution according to each time point. Further analysis revealed two distinct combinations of genes from the infarct core data at each time point, that exhibited generalization capabilities and accuracies of 100%, suggesting that these genes are of great importance in describing the whole progression.
With this information, we performed enrichment analyses to examine the mechanism of action of cardiac remodeling in the infarct core at each time point. A total of 355 processes were altered in the infarct core for at least one time point, and 101 were altered at all time points (i.e., 6, 30, and 45 days after infarction). The high number of upregulated ECMrelated processes was consistent with previous reports that glycosaminoglycan mediates cardiac remodeling by facilitating the inflammatory response [35]. In addition to its role in remodeling, the ECM also plays an important role in cell communication [36] and contributes to helping damaged cardiac tissue. To avoid ischemic conditions, the organism will attempt to revascularize and repair the affected area, promoting angiogenesis. Accordingly, angiogenesis was among the first morphogenic responses that we detected as being sustained over time. We also detected the apparent upregulation of other processes related to the morphogenesis of different organs, suggesting that postischemic cells present in the heart could be trying to recapitulate embryogenic processes and bolster regeneration.
Upon identification of the altered processes, we performed a protein-triggering analysis to establish the key genes mediating the adverse remodeling response. We then modeled the behaviors of these genes, generating the first description of an integrative mechanism of action for MI progression. This mechanistic analysis overlapped at different time points, and the pathways that were common between the source proteins and cardiac remodeling involved IGF1R, RAF1, KPCA, JUN, and PTN11 as modulators. The downstream signaling of these proteins produces a mechanism of action potentially involving 18 proteins, 15 of which were differentially expressed in at least one cohort.
The present study is the first to perform an in-depth analysis of the molecular changes characterizing the progression of MI. Our results contain a large amount of intricate information. As expected, MI evolution encompasses a myriad of altered processes and specific proteins that regulate time-dependent stages and determine the final extent of the pathology.

Conclusions
We have integrated dynamic transcriptomic regulation and all available information on MI and cardiac remodeling to delineate a clear, structured, and simplified picture of the complete post-MI remodeling process. Furthermore, our study also elucidates new potential targets and therapeutic windows for the treatment of MI complications and novel avenues for MI research, with the ultimate goal of fully unravelling the whole pathology.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/cells10123268/s1, Table S1: Initial protein filtering process to generate the truth table, Table S2: List of potential biomarkers to describe MI progression, Table S3

Informed Consent Statement: Not applicable.
Data Availability Statement: Transcriptomic data is publicly available in the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/, accessed on 8 February 2019) database with series accession number GSE34569.

Conflicts of Interest:
The authors declare no conflict of interest.