Protein signaling networks are static views of dynamic processes where proteins go through many biochemical modifications such as ubiquitination and phosphorylation to propagate signals that regulate cells and can act as feed-back systems. Understanding the precise mechanisms underlying protein interactions can elucidate how signaling and cell cycle progression occur within cells in different diseases such as cancer. Large-scale protein signaling networks contain an important number of experimentally verified protein relations but lack the capability to predict the outcomes of the system, and therefore to be trained with respect to experimental measurements. Boolean Networks (BNs) are a simple yet powerful framework to study and model the dynamics of the protein signaling networks. While many BN approaches exist to model biological systems, they focus mainly on system properties, and few exist to integrate experimental data in them. In this work, we show an application of a method conceived to integrate time series phosphoproteomic data into protein signaling networks. We use a large-scale real case study from the HPN-DREAM Breast Cancer challenge. Our efficient and parameter-free method combines logic programming and model-checking to infer a family of BNs from multiple perturbation time series data of four breast cancer cell lines given a prior protein signaling network. Because each predicted BN family is cell line specific, our method highlights commonalities and discrepancies between the four cell lines. Our models have a Root Mean Square Error (RMSE) of 0. 31 with respect to the testing data, while the best performant method of this HPN-DREAM challenge had a RMSE of 0. 47. To further validate our results, BNs are compared with the canonical mTOR pathway showing a comparable AUROC score (0. 77) to the top performing HPN-DREAM teams. In addition, our approach can also be used as a complementary method to identify erroneous experiments. These results prove our methodology as an efficient dynamic model discovery method in multiple perturbation time course experimental data of large-scale signaling networks. The software and data are publicly available at https: //github. com/misbahch6/caspo-ts.
Regarding the training of BNs with respect to multiple perturbation datasets, CellNOpT (CNO) [16] assembles BNs from a Prior Knowledge Network (PKN) and phosphoproteomic datasets. Their tool has been implemented using stochastic search algorithms (more precisely, a genetic algorithm), to suggest multiple BNs explaining the data [17]. However, stochastic search methods cannot generate a complete set of solutions, hence they cannot guarantee a global optimal solution. In [11,12], the authors overcome this problem by proposing caspo, an approach based on ASP to infer BNs explaining the underlying protein signaling network. This approach can generate all possible optimal Boolean models as compared to the CellNOpt approach. The authors in [14], presented a framework based on integer linear programming (ILP) to learn the subset of interactions best fitted to the experimental data. Recently, another approach based on ILP has been proposed to reconstruct BNs from experimental data. Their learning approach do not require the information about the activation/repression properties of the network’s edges [13]. The methods mentioned above are very useful but restrain themselves to learn from only two time points, assuming the system has reached an early steady-state when the measurements are performed. This assumption prevents us from capturing interesting characteristics like loops [3]. To overcome this issue, the caspo time series (caspo-ts) method was proposed in [8]. This method learns BNs from multiple perturbation phosphoproteomic time series data given a PKN. The proposed method is based on ASP and a model-checking step is needed to detect true positive BNs. They tested their approach on synthetic data for a small PKN (≈17 nodes and ≈50 edges) [8]. More recently, an approach based on genetic algorithms was proposed to learn context specific networks given a PKN and experimental information about stable states and their transitions but it does not scale well with large networks and finding a global optimum is not guaranteed [18].
The structure of the protein signaling network was generated by mapping the experimentally measured phosphorylated proteins (HPN-DREAM dataset) to their equivalents from literature-curated databases and connecting them together within one network (see Materials and methods). The reference network (Fig 2) was built using the ReactomeFIViz app (also called the ReactomeFIPlugIn or Reactome FI Cytoscape app) [25], which accesses the interactions existing in the Reactome and other databases [25,26]. The PKN shown in Fig 2 consists of 64 nodes (7 stimuli, 3 inhibitors, and 23 readouts) and 178 edges. The learning and testing datasets used in this study were extracted from the HPN-DREAM challenge and correspond to time series protein measurements upon different perturbations of four breast cancer cell lines—UACC812, BT20, BT549, and MCF7 [20,21] (see Materials and methods). Since readout signals are measured on variable ranges depending on the protein, a normalization step was necessary. The learning dataset had a few noisy, inconsistent and incomplete time series data points. The caspo-ts system identified these inconsistencies existing in the time series data. The recurrent experimental inconsistency observed was an oscillation in the protein signal upon experimental inhibition of the same protein. To resolve the above mentioned issues, we performed the following data processing steps on the learning dataset: The experimental errors pointed in steps 2-5 were raised as warning or exceptions by caspo-ts. Steps 1 to 5 were applied on the learning dataset. Only step 1 was applied on the testing dataset. In this section, we show the generated BNs for each cell line. For this, we used caspo-ts to learn the BNs from the PKN (Fig 2) and the phosphoproteomic data of four breast cancer cell lines—BT20, BT549, MCF7, and UACC812. We inferred a family of cell line specific BNs for each cancer cell line and they are shown in the Supplementary Figures (S1, S2, S3 and S4 Figs). As explained in the caspo-ts modeling framework section, the caspo-ts method produces BNs fulfilling two criteria, (i) satisfaction of the over-approximation criteria (see Materials and methods) and (ii) optimality with respect to the RMSE objective function. ASP-optimal solutions were fast to collect, their computation time ranged from 36 seconds to 3 minutes depending on the cell line (S1 Table). Afterwards, these ASP-optimal BNs were given to the model-checker for further verification. This second step is more complex and we put a restriction for the computation time of 7 days for each cell line. The number of verified BNs varies from one cell line to another, depending on a number of factors such as the number of perturbations, the order of answer sets in the solutions space, and the perturbation order. The total number of verified ASP-optimal BNs within the 7 days time-frame were 231,52,188 and 150 for the BT549, MCF7, BT20 and UACC812 cell lines respectively. We obtained 191,21, and 72 true positive BNs for BT549, MCF7, and BT20 cell lines respectively with an optimal fit to the data. For the UACC812 cell line, we were unable to obtain true positive BNs within the 7 day time limit for verification. Hence, we kept the first 20 BNs from the 150 ASP-optimal BNs for the UACC812 cell line. The caspo-ts method uses the ASP solver (clingo), which is able to exhaustively enumerate all solutions. The clingo solver by default uses an enumeration scheme, in which, once a solution is found, it backtracks to the first point from where the next solution can be found. This typically leads to the situation where successive solutions only change in a small part. As a result, caspo-ts may enter a solution space where BNs are clustered together. We have observed that given the size of the PKN and the small number of perturbations in the experimental data, the solution space of the caspo-ts can be rather very large containing billions of BNs making it difficult to enumerate true positive BNs (because of the model checking overhead) in a reasonable time if it gets stuck in a cluster of false positive BNs. An aggregated network was built (Fig 3) by combining the BN families (with 191,21,72, and 20 BNs for BT549, MCF7, BT20, and UACC812 cell lines respectively) obtained for the four cell lines by keeping the hyper-edges (Boolean functions) having a frequency higher than 0. 3 within each BN family. The frequency is calculated by counting the number of common Boolean functions and dividing it by the total number of Boolean functions within the BN family of each cell line. This aggregated network contains 34 nodes and 74 Boolean functions involving 36 AND gates. As compared to the PKN (Fig 2), the inferred networks are highly specific to each cell line. In Fig 3, all cell lines share only 4% of Boolean functions which are shown in thick black colored edges. This shows that the inferred BNs of these four breast cancer cell lines are very diverse and different from each other. To measure cell lines similarity, we calculated the similarity score by applying the Graph Similarity Measure (see Materials and methods) on the family of BNs (with 191,21,72, and 20 BNs for BT549, MCF7, BT20, and UACC812 cell lines, respectively). This algorithm receives two parameters as input: (1) one gold standard BN and (2) a family of BNs. It outputs a score in [0; 1], measuring the average of the similarity scores between each BN in the family and the gold standard BN. In our case, the gold standard BN is the aggregation of one family of BNs. The similarity scores between all pairs of breast cancer cell lines are shown in Table 1. Fig 3 agrees with the results presented in Table 1 as we can see the clear discrepancies among the four cell lines. It can be seen that 23% of the Boolean functions are shared among BT549 and MCF7, and also between BT20 and UACC812. BT20 shares the least number of Boolean functions (15%) with BT549. This table revealed pronounced differences among different cell lines of breast cancer. We also analyzed the diversity of Boolean functions among the family of BNs within the same cell line. The similarity among Boolean functions from BT20 (0. 73) and MCF7 (0. 63) is higher than the ones from BT549 (0. 43) and UACC812 (0. 46) cell lines. There are a total of 69 distinct Boolean functions shown in Fig 4 along with their respective frequencies. It is interesting to note that the B549 and UACC812 cell lines have more distinct models among their family of BNs with a variable frequency range. This shows that these cell lines have different mechanisms agreeing with the results obtained through graph similarity measure given in Table 1. Fig 5 shows the common Boolean functions along with their frequency in all BNs. Interestingly, only 4% of the Boolean functions are shared in all cell lines and 88% of these shared functions have the same frequency. In this figure, there is only one Boolean function which is frequent in 3 cell lines and has a lower frequency in BT20. The union of the BNs learned for each cell line is displayed in the Supplementary Figures (S1, S2, S3 and S4 Figs). The caspo-ts method revealed that cell line specific reactions are clustered around the AKT, MAPK3, and PIK3R1 proteins. PI3K is an important factor for cancer development in HER2 amplified cancers (UACC812) as compared to non-HER2 amplified (BT20, BT549 and MCF7) cancer cell lines. We can see from the Supplementary Figures (S1, S2, S3 and S4 Figs) that PIK3R1 exists in all cell lines but is rather more connected in the UACC812 cell line with 10 incoming edges while in others with only 1 incoming edge. The PIK3R1 node in UACC812 (S4 Fig) has a centrality measure of 0. 37 while in the other three cell lines the centrality measure is less than 0. 11. The centrality measure is used to quantify the most important node within a network i. e. , the number of times a node has been used as a bridge (along the shortest path) to connect to other nodes in the network [27]. It has been established that P1K3R1 (the regulatory unit of PI3K) plays an important role in suppressing tumors [28,29]. Recently, it has been found that PIK3R1 is mutated in 3% of breast cancer cell lines[30]. Nonetheless, it is worth studying the impact of the PIK3R1 regulatory unit in breast cancer. The performance of the caspo-ts method is evaluated using three criteria: 1) RMSE calculation using a typical learning and testing data approach, 2) random data comparison, 3) AUROC (Area Under the Operating Curve) score. The BNs are learned using the learning dataset (see Materials and methods) only. The prediction accuracy is evaluated by comparing the RMSE of trajectories in the testing dataset with those recovered by the learned networks (see Eq 1). There are two types of RMSE—discrete and model. The discrete RMSE is imposed by the discretization of the method. Since we use a discrete learning approach, our recovered traces will be in {0,1} and this introduces an error with respect to continuous measurements in [0; 1]. The model RMSE refers to the learned BN error with respect to the normalized time series data; that is, the model RMSE is at least as large as the discrete RMSE. When the difference between these two is zero then the inferred BNs are able to recover the discrete trajectories without any error. If the model RMSE is greater than the discrete RMSE then the inferred BNs have some errors in the recoverability of the discrete time series data. To check how our method performs in case of random time series, we have calculated the RMSE score for random data and compared it with learning and testing data. Next, the validity of these networks is verified by comparing them with the canonical MTOR signaling pathway using two parameters, i. e. , true positive rate (TPR) and false positive rate (FPR).
As a future direction, we are planning to investigate several aspects of the caspo-ts method, such as (i) the order of the solution space of over-approximated Boolean models; (ii) the computational time for checking reachability; (iii) designing an efficient experimental design strategy and applying it prior to selecting the most informative experiments. Because caspo-ts uses an ASP solver to enumerate BNs, in the resulting sequence of solutions similar BNs are typically clustered together. This can be problematic for large scale problems where we cannot explore the whole solution space in reasonable time. We are currently working on sampling to randomly select BNs from the solution space. Further, we are also studying another technique, which allows for shuffling the order in which solutions are enumerated [32]. We are planning to implement this by dynamically modifying the heuristic of the ASP solver at execution time. Finally, to reduce the false positive rate, we are planning to use multi-shot ASP solving [33] allowing us to customize the search and modify the underlying ASP program at runtime. In our case, we can call the model checker during solving to learn and add constraints to prune wrong BNs early on.
The DREAM portal provides unrestricted access to complex, pre-tested data to encourage the development of computational methods. In this study, we are focused on the HPN-DREAM challenge, which was motivated by the fact that the same perturbation may lead to different signaling behaviors in different backgrounds, making it necessary to build a model which can perform unseen predictions (absent from the learning data). The main goal of the HPN-DREAM challenge is to learn signaling networks efficiently and effectively to predict the dynamics of breast cancer [19]. PKNs are available in different databases such as Reactome, PID, and KEGG among others [26,34–45]. We can construct a PKN through a tool such as ReactomeFIViz [25] which is available as a Cytoscape [46] plugin. A PKN alone cannot be used to build reliable dynamical models or to explain underlying biological behaviors [16,47], especially in the case of multiple perturbations data because of the need of specificity. In order to overcome this issue, methods have been proposed which take into account both literature based knowledge and experimental data to build logic models [3,7, 11,12,16]. In the caspo-ts modeling framework section, we have given the formal definitions of the inputs and the output (BN) of the caspo-ts. Here, we formally describe the over-approximation criteria. Finally. we give pseudo encodings of the input of the ASP part of the caspo-ts. This work introduces the study of a graph similarity measure in order to check the variability among the families of BNs generated by caspo-ts. We compare the reactions existing in the gold standard network (A) with the family of BNs (B) and is based on the Jaccard similarity coefficient which measures the similarity of these models.
Traditional canonical signaling pathways help to understand overall signaling processes inside the cell. Large scale phosphoproteomic data provide insight into alterations among different proteins under different experimental settings. Our goal is to combine the traditional signaling networks with complex phosphoproteomic time-series data in order to unravel cell specific signaling networks. In this study, we have applied the caspo time series (caspo-ts) approach which is a combination of logic programming and model checking, over the time series phosphoproteomic dataset of the HPN-DREAM challenge to learn cell specific BNs. The learned BNs can be used to identify the cell specific topology. Our analysis suggests that caspo-ts scales to real datasets, outputting networks that are not random with a lower fitness error than the models used by the 178 methods which participated in the HPN-DREAM challenge. On the biological side, we identified the cell specific and common mechanisms (logical gates) of the cell lines.
phosphorylation
learning
medicine and health sciences
breast tumors
ontologies
protein interaction networks
signaling networks
social sciences
cancers and neoplasms
neuroscience
learning and memory
oncology
data management
cognitive psychology
network analysis
information technology
data processing
research and analysis methods
computer and information sciences
mathematical functions
proteins
mathematical and statistical techniques
breast cancer
proteomics
biochemistry
psychology
post-translational modification
biology and life sciences
cognitive science
Computational discovery of dynamic cell line specific Boolean networks from multiplex time-course data