TPMS technology to infer biomarkers of macular degeneration prognosis in in silico simulated prototype-patients under the study of heart failure treatment with sacubitril and valsartan

Unveiling the mechanism of action of a drug is key to understand the benefits and adverse reactions of drug(s) in an organism. However, in complex diseases such as heart diseases there is not a unique mechanism of action but a wide range of different responses depending on the patient. Exploring this collection of mechanisms is one of the clues for a future personalised medicine. The Therapeutic Performance Mapping System (TPMS) is a Systems Biology approach that generates multiple models of the mechanism of action of a drug. This is achieved by (1) modelling the responses in human with an accurate description of the protein networks and (2) applying a Multilayer Perceptron-like and sampling method strategy to find all plausible solutions. In the present study, TPMS is applied to explore the diversity of mechanisms of action of the drug combination sacubitril/valsartan. We use TPMS to generate a range of mechanism of action models explaining the relationship between sacubitril/valsartan and heart failure (the indication), as well as evaluating their relationship with macular degeneration (a common/recurrent adverse effect). We found that a lower response in terms of heart failure treatment is more associated to macular degeneration development, although good response mechanisms can also associate to the adverse effect. A set of 30 potential biomarkers are proposed to identify mechanisms (or patients) more prone to suffering macular degeneration when presenting good heart failure response. As each molecular mechanism can be particular not only of cells but also individuals, we conclude that the study of the collection of models generated using TPMS technology can be used to detect adverse effects personalized to patients.


Introduction
The algorithm of TPMS uses as input signals the values of activation (+1) and inactivation (-1) of the drug protein targets. The outputs are then the values of activation and inactivation of the proteins defining the indication's phenotype (retrieved from the BED), named effectors. Each node of the protein network receives as input the output of the connected nodes and each link receives a weight ( # ). The sum of inputs is transformed by a hyperbolic tangent function to generate the score of the node, which becomes the "output signal" towards the connected nodes.
The # parameters are obtained by optimization, using a Stochastic Optimization Method based on Simulated Annealing. 22 The models are trained by using the restrictions defined by the BED and the specific conditions set by the user. We ranked all solutions by the number of restrictions

Measures to compare sets of MoAs
To understand the relationships between all potential mechanisms we defined measures of comparison between different sets of solutions. We expect that a drug will revert the conditions of a disease phenotype. Consequently, a drug should inactivate the active protein effectors of a pathology-phenotype and activate the inactive ones. Here, we defined several measures in order to study and compare sets of MoAs from different views (see further details in supplementary material).

TSignal
To quantify the intensity of the response of a MoA and compare it with others, we created a measure called TSignal. The TSignal is the average of the output signals of the protein effectors (equation in supplementary material).

Distance between two sets of MoAs
We calculated the distance between two or more sets of MoAs in order to determine their similarity. For that, we used a modified Hausdorff distance (MHD) introduced by Dubuisson and Jain. 23 Details of equations are in supplementary material.

Potential biomarkers extracted from MoAs
For HF, MoAs are ranked by their TSignal and split in four quartiles: the first quartile (top 25%) contains MoAs with higher intensity of the response (TSignal), which in turn reduces the values of the effectors associated with a disease phenotype (we named them as Low-disease MoAs).
The fourth quartile (bottom 25%) collects MoAs with lower intensity of response (thus, we named as High-disease MoAs). On the other hand, for MD, the first quartile (top 25%) contains MoAs with higher intensity, which in this case, as an adverse event, it increases the values of the effectors associated to the comorbidity (we name them as High-adverseEvent MoAs). The fourth quartile (bottom 25%) collects MoAs with lower intensity of response (thus, we named as Low-adverseEvent MoAs).
We use the comparison between both groups of High-and Low-MoAs to identify the bestclassifier proteins. Best-classifier proteins are those that classify the best between High-and Lowgroups and are determined by a Data-Science strategy (see supplementary material). Bestclassifier proteins are strongly related to the intensity of a response and are differently distributed between High-and Low-MoAs. We only select the 200 proteins (or pair of proteins) reproducing with higher accuracy the classification. Assuming the hypothesis that the selected MoAs are representative of individual prototype patients, these proteins can be used as biomarkers to classify a cohort of patients by the activity or inhibition of the protein.
Each best-classifier protein has different output signals in the Low-and High-group MoAs and the distribution in both sets can be compared. We use a Mann-Whitney U test to compare the two distributions of output signals and select those proteins for which the difference is significant (p-value< 0.01), having an average output signal in Low-HF with opposite sign to the average output signal in High-HF (i.e. positive vs. negative or vice versa). We name these as differential bestclassifier proteins. By following this strategy, we can identify two groups of differential bestclassifier proteins: those active in Low-disease (positive output signal in average) and inactive in High-disease (negative output signal in average), and those active in High-disease and inactive in Low-disease.

Results and discussion
We applied TPMS to the HPN using as input signals the drug targets of sacubitril/valsartan (NEP / AT1R) and as output signals the proteins associated with Heart Failure (HF) extracted from the BED. As described in the methodology, out of all MoAs found by TPMS, we selected 200 satisfying the largest number of restrictions (and at least 80% of them) to perform further analysis.
To rank the MoAs according to the intensity of the signal arriving from the drug, we calculated the TSignal of every MoA associated with HF and MD, i.e. the average output signal arriving to the protein effectors of both pathologies. According to the TSignal of HF and following the procedure described in Materials and Methods, we defined two groups of MoAs: Low-HF, containing the MoAs with a higher intensity of the response and therefore a healthier phenotype, and High-HF, with the MoAs of lower intensity of the response and therefore an increased HF disease phenotype. We also calculated the TSignal of MD and define two groups of MoAs: Low-MD, producing a reduced adverse effect, and High-MD, producing an increased adverse effect.
Note that TPMS was only executed once, optimising the results to satisfy the restrictions on HF data. The values of MD are obtained by measuring the signal arriving at the MD effectors, which are part of the HPN and also receive signal. This procedure was chosen because we treat HF as the indication of the drug (sacubitril/valsartan), while MD is a potential adverse effect.
In the following sections, we analysed and compared the four groups of MoAs and searched for biomarkers that can potentially identify MD as an adverse effect of the drug (and consequently classify patients). We used the web server GUILDify v2.0 15 to analyse the comorbidity of both phenotypes and analyse the relevance of the proposed biomarkers in a different context.  Figure 1a).

Analysis of MoAs of high/low intensity response associated to HF
We selected the 200 best-classifier proteins after defining the two groups of MoAs. These are defined as the proteins (or pairs of proteins) that can allocate better the MoAs in High-HF and Low-HF. Among these proteins, we identified the differential best-classifier proteins. These are proteins that have output signals significantly different between Low-HF and High-HF (Mann-Whitney U test, adjusted p-value<0.01) and for which the average has opposite sign between the cohorts (i.e. active in Low and inactive in High or vice versa). We identified two groups of We calculated the Gene Ontology (GO) enriched functions of the differential best-classifier proteins that are active in Low-HF MoAs and inactive in High-HF MoAs. The enrichment was calculated using the software FuncAssociate. 24 Among the enriched functions we find processes associated with the SCAR complex, the positive regulation of actin nucleation, the regulation of neurotrophin TRK receptor and dendrite extension. We did the same for GO functions of the differential best-classifier proteins that are inactive in Low-HF and active in High-HF. We find functions such as phosphatidylinositol kinase activity, MAP kinase activity, DNA damage induced protein phosphorylation and superoxide anion generation. Although some functions are enriched in both sets, such as Fc gamma receptor signalling, the majority of functions enriched in different groups of differential best-classifier proteins are different (see Supplementary Table 2).
Some of the proteins and functions highlighted in the current analysis have been related to myocardial function. On the one hand, our differential best-classifier proteins active in Low-HF MoAs and inactive in High-HF MoAs point towards an important role for actin nucleation and polymerization mechanisms in drug response (reflected by the functions regulation of actin nucleation, regulation of Arp2/3 complex-mediated actin nucleation, SCAR complex, filopodium tip, or dendrite extension); in fact, the alteration of actin nucleation and polymerization mechanisms has been reported in heart failure. 25-27 Interestingly, a role for the activation of another differential best-classifier candidate, ATGR2, has been proposed to mediate some of the beneficial effects of angiotensin II receptor type 1 antagonists, such as valsartan. 28,29 On the other hand, results for the differential best-classifier proteins inactive in Low-HF and active in High-HF point towards the importance of phosphatidylinositol kinase mediated pathways (phosphatidylinositol-3,4-bisphosphate 5-kinase activity) and MAP kinase mediated pathways (MAP kinase kinase activity, best classifier proteins MAPK1, MAPK3, MAPK11, MAPK12 or MAPK13) in response to sacubitril/valsartan; furthermore, both signalling pathways have been associated to cardiac hypertrophy and subsequent heart failure. 30,31

Limitations
Although TPMS returns the amount of signal from the drug arriving to the rest of the proteins in the HPN, this signal is only a qualitative measure. We are not using data about the dosage of the drug or the quantity of expression of the proteins. However, we are already working to make TPMS move towards the growing tendency of Quantitative Systems Pharmacology. The quantification of the availability of drugs in the target tissue for each patient opens the opportunity to have an accurate patient simulation to do in silico clinical trials.

Conclusions
There is a need of systems biology-based methods to simulate the different responses of a drug in patients. The specific case of sacubitril/valsartan stands out because of the amount of resources invested in the safety of the drug and the concern on the risk of MD. In this study, we apply TPMS to uncover the different MoAs of sacubitril/valsartan over HF and reveal its relationship with MD. We hypothesize that all MoAs coexist in cells, but in a population of patients some cells may have more proclivity to certain MoAs than others. We define a prototype-patient as one with a single MoA for a drug and study the in silico trial of HF treatment with sacubitril/valsartan and its potential adverse effect MD. TPMS achieves this by modelling an accurate representation of the HPN and applying a Multilayer Perceptron-like and sampling method strategy to find all plausible solutions. We found that HF low responder MoAs are more associated to MD development at the same time, although good responders are also associated to MD. Different sets of proteins have been found to classify the mechanisms according to HF and MD response, which include functions such as PI3K and MAPK kinase signalling pathways, involved in HF-related cardiac hypertrophy, or fibrinolysis and coagulation processes (e.g. FGB, SERPINE1 or SERPING1) and growth factors (e.g. FGF1 or PDGF) related to MD induction. We propose 30 biomarkers to identify patients potentially developing MD under the successful treatment with sacubitril/valsartan. Out of this 30, we propose 10 biomarkers that have been found to be involved in the comorbidity between HF and MD predicted by a different approach (GUILDify), including well-known HF and MD effectors also related to the mechanisms of sacubitril/valsartan and/or HF, such as AGER, NRG1, ITGB5 or IL1A.