Addressing the Needle in a Haystack Problem in Time Behavioural Study of 3rd Order Gene Combinations in WNT3A Stimulated HEK 293 Cells

1. Integration, Innovation and Insight

At any particular time point, it is often the case that genes are working synergistically in combinations, even though their expression measurements are singular in nature. The problem explodes combinatorially with even a small set of recorded genes, when one steps to explore higher (here 3rd) order combinations. With the huge total number of combinations, it becomes nearly impossible for any biologist to study the system wide dynamics of any pathway as well as locate combinations of genuine interest. This study demonstrates how biologists can use a machine learning based search engine to address the needle in a haystack problem of discovering meaningful combinations of higher order in a vast search forest, while cutting down the time required to search the same. Further wet lab test might assist in intervening the pathway at a combinatorial level, in time.

2. Significance

Sinha [2] recently demonstrated the use of machine learning based search engine to rank/reveal gene combinations at 2nd order for the time series data by Gujral and MacBeath [1] and showed how it is possible to locate combinations of priority that might be working synergistically, using sensitivity methods and powerful support vector ranking algorithm. However, the problem explodes combinatorially with even a small set of 71 recorded genes in the study by Gujral and MacBeath [1], when one steps to explore 3rd order combinations. With the total number of ⁷¹$C_3$ (= 57155) combinations, it becomes nearly impossible for any biologist to study the system wide dynamics of any pathway. Also, the amount of time usually needed to search for and test a combination is far more than the search down by the machine learning based search engine. Here, I extend the research work by Sinha [2] to conduct a behavioral study of 3rd order gene combinations using individual gene expressions measured in time, in WNT3A stimulated HEK 293 cells. 1

3. Introduction

The details of the machine learning based search engine has been recently published in Sinha [2] and deployed to explore the 2nd order combinations of genes in the data set provided by Gujral and MacBeath [1]. Nevertheless, here, I point to the fundamentals of the published work for completeness.

3.1. A Combinatorial Problem

Sensitivity analysis plays a major role in computing the strength of the influence of involved factors in any phenomena under investigation. When applied to expression profiles of various intra/extracellular factors that form an integral part of a signaling pathway, the variance and density based analysis yields a range of sensitivity indices for individual as well as various combinations of factors. These combinations denote the higher order interactions among the involved factors. Computation of higher order interactions is often time consuming but it gives a chance to explore the various combinations that might be of interest in the working mechanism of the pathway. For example, in a range of fourth order combinations among the various factors of the Wnt pathway, it would be easy to assess the influence of the destruction complex formed by APC, AXIN, CSKI and GSK3 interaction. But the effect of these combinations vary over time as measurements of fold changes and deviations in fold changes vary. So it is imperative to know how an interaction or a combination of the involved factors behave in time and Sinha [2] develops a procedure to track the behaviour by exploiting the influences of these involved factors.

3.2. A Possible Solution

In this work, after estimating the individual effects of factors for a higher order combination, the individual indices are considered as discriminative features. A combination, then, is a feature set in higher order (≥2 ,i.e multivariate). With an excessively large number of factors involved in the pathway, it is difficult to search for important combinations in a wide search space over different orders. Exploiting the analogy with the issues of prioritizing webpages using ranking algorithms, for a particular order, a full set of combinations of interactions can then be prioritized based on these features using a powerful ranking algorithm via support vectors Joachims [3]. Recording the changing rankings of the combinations over time reveals how higher order interactions behave within the pathway and when an intervention might be necessary to influence the interaction within the pathway.

3.3. Wnt Signaling and Secretion

Sharma [4]’s accidental discovery of the Wingless played a pioneering role in the emergence of a widely expanding research field of the Wnt signaling pathway. A majority of the work has focused on issues related to • the discovery of genetic and epigenetic factors affecting the pathway Thorstensen et al. [5] &Baron and Kneissel [6], • implications of mutations in the pathway and its dominant role on cancer and other diseases Clevers [7], • investigation into the pathway’s contribution towards embryo development Sokol [8], homeostasis Pinto et al. [9] & Zhong et al. [10] and apoptosis Pećina-Šlaus [11] and • safety and feasibility of drug design for the Wnt pathway Kahn [12], Garber [13], Voronkov and Krauss [14], Blagodatski et al. [15] & Curtin and Lorenzi [16].

The Wnt phenomena can be roughly segregated into signaling and secretion part. The Wnt signaling pathway works when the WNT ligand gets attached to the Frizzled(FZD)/LRP coreceptor complex. FZD may interact with the Dishevelled (DVL) causing phosphorylation. It is also thought that Wnts cause phosphorylation of the LRP via casein kinase 1 (CK1) and kinase GSK3. These developments further lead to attraction of Axin which causes inhibition of the formation of the degradation complex. The degradation complex constitutes of AXIN, the $\beta$-catenin transportation complex APC, CK1 and GSK3. When the pathway is active the dissolution of the degradation complex leads to stabilization in the concentration of $\beta$-catenin in the cytoplasm. As $\beta$-catenin enters into the nucleus it displaces the GROUCHO and binds with transcription cell factor TCF thus instigating transcription of Wnt target genes. GROUCHO acts as lock on TCF and prevents the transcription of target genes which may induce cancer. In cases when the Wnt ligands are not captured by the coreceptor at the cell membrane, AXIN helps in formation of the degradation complex. The degradation complex phosphorylates $\beta$-catenin which is then recognised by F BOX/WD repeat protein $\beta$-TRCP. $\beta$-TRCP is a component of ubiquitin ligase complex that helps in ubiquitination of $\beta$-catenin thus marking it for degradation via the proteasome.

Contrary to the signaling phenomena, the secretion phenomena is about the release and transportation of the WNT protein/ligand in and out of the cell, respectively. Briefly, the WNT proteins that are synthesized with the endoplasmic reticulum (ER), are known to be palmitoyleated via the Porcupine (PORCN) to form the WNT ligand, which is then ready for transportation. It is believed that these ligands are then transported via the EVI/WNTLESS transmembrane complex out of the cell Bänziger et al. [17] & Bartscherer et al. [18]. The EVI/WNTLESS themselves are known to reside in the Golgi bodies and interaction with the WNT ligands for the later’s glycosylation Kurayoshi et al. [19] & Gao and Hannoush [20]. Once outside the cell, the WNTs then interact with the cell receptors, as explained in the foregoing paragraph, to induce the Wnt signaling. Of importance is the fact that the EVI/WNTLESS also need a transporter in the from of a complex termed as Retromer.

4. Methods

Please refer to sections of Sinha [2] for methods, design of study and analysis of data for 2nd order combinations. The same method and design of study is used to generate results for 3rd order combinations presented in this study.

5. Time Series Data

Gujral and MacBeath [1] present a set of 71 WNT-related gene expression values for 6 different times points over a range of 24-hour period using qPCR. The changes represent the fold-change in the expression levels of genes in 200 ng/mL WNT3A-stimulated HEK 293 cells in time relative to their levels in unstimulated, serum-starved cells at 0-hour. Gujral and MacBeath [1] state that qPCR data are the means of three biological replicates. Only genes whose mean transcript levels changed by more than two-fold at one or more time points during the 24-hour time course were considered significant. Positive (negative) numbers represent up (down) -regulation. We have already covered the issues related to these data sets in detail in Sinha [21]. Readers are requested to go through them in the pointed reference. The tools of study which are used here have been published in another foundational work in Sinha [21].

6. Design of Experiment

6.1. Pipeline for Time Series Data

For the case of time series data, interactions among the contributing factors are studied by comparing triplets of fold-changes at single time points. The prodecure begins with the generation of distribution around measurements at single time points with added noise is done to estimate the indices. A distribution is generated for the fold changes at single time points. Then for every gene, there is a vector of values representing fold changes as well as deviations in fold changes for different time points and durations between time points, respectively. Next a listing of all $C_k^n$ combinations for k number of genes from a total of n genes is generated. k is $\ge 2$ and $\le (n-1)$. Each of the combination of order k represents a unique set of interaction between the involved genetic factors. After this, the datasets are combined in a specifed format which go as input as per the requirement of a particular sensitivity analysis method. Thus for each $p^{th}$ combination in $C_k^n$ combinations, the dataset is prepared in the required format from the distributions for two separate cases which have been discussed above. (See .R code in mainScript-1-1.R). After the data has been transformed, vectorized programming is employed for density based sensitivity analysis and looping is employed for variance based sensitivity analysis to compute the required sensitivity indices for each of the p combinations. This procedure is done for different kinds of sensitivity analysis methods.

After the above sensitivity indices have been stored for each of the $p^{th}$ combination, the next step in the design of experiment is conducted. Since there is only one recording of sensitivity index per combination, each combination forms a training example which is alloted a training index and the sensitivity indices of the individual genetic factors form the training example. Thus there are $C_k^n$ training examples for $k^{th}$ order interaction. Using this training set $SV{M}_{learn}^{Rank}$ Joachims [3] is used to generate a model on default value C value of 20. In the current experiment on toy model C value has not been tunned. The training set helps in the generation of the model as the different gene combinations are numbered in order which are used as rank indices. The model is then used to generate score on the observations in the testing set using the $SV{M}_{classify}^{Rank}$ Joachims [3]. Note that due to availability of only one example per combination, after the model has been built, the same training data is used as test data to generates the scores. This procedure is executed for each and every sensitivity analysis method. This is followed by sorting of these scores along with the rank indices (i.e the training indices) already assigned to the gene combinations. The end result is a sorted order of the gene combinations based on the ranking score learned by the $SV{M}^{Rank}$ algorithm. Finally, this entire procedure is computed for sensitivity indices generated for each and every fold change at time point and deviations in fold change at different durations. Observing the changing rank of a particular combination at different times and different time periods will reveal how a combination is behaving.

Note that the following is the order in which the files should be executed in R, in order, for obtaining the desired results (Note that the code will not be explained here) - • use source("mainScript-1-1.R") with arguments for Dynamic data • source("SVMRank-Results-D.R"), to rank the interactions (again this needs to be done separately for different kinds of SA methods), • use source("Combine-Time-files.R"), if computing indices separately via previous file, • source("Sort-n-Plot-D.R") to sort the interactions. Note that the sorting is chages the interaction ranking in time. Thus • use source("Interaction-Priority-Intime.R") to find the prioritized ranking of each and every interaction over the different time points and finally • use source("Print-Ranking-AND-Interaction-Rank.R") to print individual ranking of the required input factor with other interaction factors. Table 1 shows the ranking scores in discending order for 3rd order combinations using four different sensitivity methods (i.e rows) and at five different time points (i.e columns).

7. Results & Discussion

7.1. Time Series Data by Gujral and MacBeath [1]

NOTE - Ranking was assigned on scores that were sorted in DECREASING values. So, 1 was assigned to highest score and vice versa.

Results for the 3^rd order interactions are presented here. The results first discuss the behaviour of interactions across the snapshots of time using the computed sensitivities on fold change measurements per time snapshot. The analysis was done using 4 different sensitivity indices. Out of the $C_3^{71}$ combinations, I consider/present only those combinations that show a ranking within first 10,000 out of 57,155. This choice is liberal and biologists/oncologists can have a more stricter choice as per need. Two observations are made, • the ranking of a particular combination is conserved (i.e within the 10,000 range) in a particular time point or in the early phase or late phase of WNT3A stimulation, across the majority of the four sensitivity methods, which is a strict criteria of assessment or • the ranking of a particular combination is conserved across time points/phase (i.e they are within the 10,000 range) and the majority of the four sensitivity methods, which is relaxed criteria of assessment. Applying this filter helps reveal important combinations of interest that might be working synergistically at a higher order level in the cell.

Regarding technical points of implementation, the rankings were generated without scaling/normalizing the time series data provided by Gujral and MacBeath [1]. For estimating the sensitivity indices, a small gaussian distribution using the function rnorm that generates a vector of normally distributed random variables given a vector length n (here 9, the 10th one is the mean/recorded gene regulation itself), a population mean $\mu$ and population standard deviation $\sigma$. The syntax for using rnorm is as follows: rnorm(n, mean, sd). Further, I use the jitter funtion to add a little bit of noise to the data. This helps to see if the generated rankings are robust or not.

7.2. Enumeration and Ranking of $C_3^{71}=57155$ Combinations from Gujral and MacBeath [1]

In the supplementary section, I present four files, each containing the rankings of 3rd order combinations, that wary in time (shown for 5 time points). Each file represents the rankings computed using a particular sensitivity method. A particular row represents a particular rank given by a number. Following this, are combinations taking up that particular rank at different time points. The methods used are Hilbert Schmidt Independence Criterion indices (HSIC) indices (with rbf and linear kernel in Da Veiga [22]) and Sobol indicies (with 2002 implementation in Saltelli [23] and martinez implementation in Martinez [24] and Baudin et al. [25]). Of importance to note is that using these files, one can see a which priority/ranking one wants to investigate and see what combinations are taking up that priority in time. Changing combinations for a particular rank indicate that at a particular point in time, a particular combination might be predominant among others.

7.3. Conserved Machine Learning Rankings for Tested PORCN-WNT-X Combinations

The Drosophila segment polarity gene product Porcupine (Porc) was first identified as being necessary for processing Wingless (Wg), a Drosophila Wnt (Wnt) family member. Tanaka et al. [26] identified Mouse (Mporc) and Xenopus (Xporc) homologs of porc and found that they encode endoplasmic reticulum (ER) proteins with multiple transmembrane domains. Further, Mporc mRNA was differentially expressed during embryogenesis and in various adult tissues, demonstrating that the alternative splicing is regulated to synthesize the specific types of Mporc. In transfected mammalian cells, they found all types of Mporc affected the processing of mouse WNT1, WNT3A, WNT4, WNT6, and WNT7B but not WNT5A. Lastly, they also found that all types of Mporc co-immunoprecipitated with various WNT proteins. Their results suggested that Mporc may function as a chaperone-like molecule for WNT.

Liu et al. [27] indicate that post-translational modification of WNTs includes lipid modification and glycosylation. The former is performed by PORCN. PORCN is a membrane-bound O-acyltransferase located in the endoplasmic reticulum and can add palmitoleate groups to WNT proteins that is necessary for WNT ligand secretion, and it is a member of the membrane-bound O-acyltransferases (MBOATs). Lipid modification is necessary for Wnt activity, and the opposite is true for glycosylation as observed by Willert et al. [28]. Liu et al. [29] developed a screen for small molecules that blocked WNT secretion and discovered LGK974, a potent and specific small-molecule PORCN inhibitor. They show that LGK974 inhibits WNT signaling, including reduction of the WNT-dependent LRP6 phosphorylation and the expression of WNT target genes, like as AXIN2. The inhibitor is effective in multiple tumor models at well-tolerated doses. Together, their findings provide a strategy and and a tool for targeting WNT-driven cancers through the inhibition of PORCN. Further down the line, Madan et al. [30] developed a novel potent, orally available PORCN inhibitor, ETC-1922159 that blocked the secretion and activity of all WNTs. ETC-1922159 is remarkably effective in treating RSPO-translocation bearing colorectal cancer (CRC) patient-derived xenografts. This is the first example of effective targeted therapy for this subset of CRC. By this demonstration they show that inhibition of WNT signaling by PORCN inhibition holds promise as differentiation therapy in genetically defined human cancers.

Based on these experimental tests and documented literature, the synergy of PORCN-WNT can be used to see if the above machine learning based engine gives appropriate ranking to 3rd order combinations of PORCN-WNT-X (X, a particular gene/protein). If the rankings are appropriate, then we can infer that the search engine indeed points to combinatorial synergies, whether tested or unexplored, at biological level. Gujral and MacBeath [1] recorded the regulations of PORCN along with WNT1, WNT2B, WNT3A, WNT4 and WNT5A.

Here, I present and demonstrate the conservation of rankings of PORCN-WNT-X combinations across different sensitivity methods. Using the linear kernel and HSIC sensitivity analysis method, Table 2 shows rankings of combinations within the first 10,000 range (with low numerical value meaning a very high priority/role) mostly during the first phase (or after t = 1 hour of WNT3A stimulation). These point to the possible role of combinations during the early phase of WNT3A stimulation. As time passes, the rankings of these combinations get lower ranks (i.e higer numerical values) pointing to their down play of role when the effect of WNT3A stimulation has subsided in the late phase. These 3rd order synergies indicate the efficacy of the machine learning based search engine in finding meaningful combinations that might be of interest to (developmental)biologists, molecular biologists and oncologists.

A total of 2415, 3rd order combinations involving PROCN were obtained from a full set of $C_3^{71}$ = 57155 combinations. Out of these 2415 combinations, those related to PORCN-WNT synergy are selected. Further, from this selected set, using the above criteria for conserved rankings, I report/tabulate the meaningful combinations that might be working synergistically. Table 3, Table 4 and Table 5 show the rankings for the same combinations as in Table 2, but using rbf kernel for HSIC, 2002 implementation for SOBOL and martinez implementation for SOBOL, respectively. As on tallies the rankings of across these tables for a particular combination, one finds that the role of the combination of interest is conserved. This conservation points to the existence of the biological synergy, whether the combination has been tested or unexplored/untested. At least at the 2nd order, considering the combinations of PORCN-WNT which have already been established in wet lab experiments (in above literature), the tabulated combinations with their appropriate ranks show the promise of the machine learning search engine in effectively locating the PORCN-WNT combinations. Further, the presented rankings point to combinations of PORCN-WNT-X, i.e the 3rd order combinations. So, considering all of the $C_3^{71}$ combinations, the machine learning search engine quickly ranks them and helps in tackling the needle in a haystack problem of finding combinations in a vast search forest, in a very short period of time. Even $C_3^{71}$ is a big range of combinations to deal with, for any biologist/oncologist in a wet lab setting. The manhine learning search engine is a tool that will assist many biologists/oncologists cut the time of search and also zoom in for particular combination of interest.

7.3.1. Examining the Behaviour of CTNNB1-PORCN-WNT3 Combination

Here we take up the case of CTNNB1-PORCN-WNT3 to examine its behaviour in time with respect to the recordings in Gujral and MacBeath [1] and the rankings of the combination across the four sensitivity methods.

Armadillo repeat is a repetitive amino acid sequence found in $\beta$-catenin. $\beta$-catenin is a protein that in humans is encoded by the CTNNB1 gene. $\beta$-Catenin was first discovered in McCrea et al. [31] as a component of a mammalian cell adhesion complex which is responsible for cytoplasmatic anchoring of cadherins. It acts as an intracellular signal transducer in the Wnt signaling pathway as shown by Peifer et al. [32] and Kemler [33]. It is known that the WNTs affect the downstream $\beta$-catenins. Also, PORCN help in the secretion of the WNTs and affects their activity. So a 3rd order axis is known to exist. Gujral and MacBeath [1] in their study recorded the activity of CTNNB1 also. Figure 1 shows the recordings of CTNNB1, PORCN and WNT3, each measured individually in Gujral and MacBeath [1]. Graphically, Figure 2 shows the changing rankings of CTNNB1-PORCN-WNT3 combination in time, across different sensitivity analysis methods (i.e in Table 2, Table 3, Table 4 and Table 5). Using the criteria used for considering a combination as showing conserved rankings, Figure 2 shows rankings in the range of 1 to 10,000, which lie above the threshold line (see figure). For HSIC linear, HSIC rbf and SOBOL martinez, it was found that the rankings were within the first 10,000 range, thus showing a majority out of four chosen sensitivity methods. Though passing a relaxed criteria across time, the machine learning search engine does point to the existence of the biological synergy between CTNNB1, PORCN and WNT3. Similar interpretations can be made for rankings of other combinations from the above tables.

7.4. Enumeration of top 10, 3rd Order Combinations for WNT3A Stimulated Response Genes

I now present the last section of the manuscript by tabulating 3rd order combinations of some of the WNT3A stimulated response genes, recorded by Gujral and MacBeath [1]. Out of the 71 response genes, I take a few genes with their family members for consideration. Further, I only present here combinations of interest, that were filtered using the machine learning search engine and those that passes the criteria for being termed conserved in rankings. Note that these combinations have either been established or require wet lab testing.

Table 6, shows the top 10 gene combinations involving a particular family of gene, that show conserved rankings across all four sensitivity methods. There are certain patterns that emerge from these rankings. I only cover these patters as experimental validations of them already exists. Nevertheless, on the basis of these experimentally established patterns, one can test remaining combinations that have yet to be explored and have been tabulated here.

7.4.1. Adenomatosis Polyposis Coli (APC)

For APC, there are 5 combinations involving paired like homeodomain 2 (PITX2) with APC, thus depicting possible synergy between APC-PITX2 that need to explored. Kuraguchi et al. [34] showed that genetic deletion of APC in embryonic mouse oral epithelium (K14-Cre;${\mathrm{APC}}^{cko/cko}$) resulted in supernumerary tooth formation, thus suggesting that WNT signaling and the levels of APC are crucial determinants of tooth initiation. Following this, Wang et al. [35] observed that in ${\mathrm{APC}}^{cko/cko}$ mice, which express Cre recombinase uniformly throughout skin ectoderm and oral and dental epithelium, died at birth. Although their tooth germs appeared normal at E13.5, by E14.5 the mutant teeth were severely disrupted, with numerous irregular epithelial buds protruding from the oral epithelium into jaw mesenchyme, and intense expression of FGF8, SHH, PITX2, p21 and FGF4 transcripts and elevated levels of $\beta$-catenin. In Gujral and MacBeath [1], APC was found to be down regulated ($-ive$ numbers), while PITX2 was upregulated ($+ive$ numbers). The search engine confirms the existence of this biological synergy by pointing out 3rd order combinations involving APC-PITX2 interaction.

7.4.2. V-myc Avian Myelocytomatosis Viral Oncogene Homolog (MYC)

For MYC, there are 5 combinations involving casein kinase 2 $\alpha$ 1 (CSNK2A1) with MYC, and 5 combinations involving SUMO specific peptidase 2 (SENP2) with MYC. This depicts the possible synergy between MYC-CSNK2A1 and MYC-SENP2 that need to explored. Yang et al. [36] observed that CSNK2A1-mediated MAX phosphorylation increased C-MYC and $\beta$-catenin binding and regulated HMGB1 promoter activity through E-BOX. These further lead to promotion of cell growth, migration, and invasion and progression of cholangiocarcinogenesis. SENP1 is frequently overexpressed and correlates with the high expression of c-MYC, in breast cancer tissues. Sun et al. [37] found that SENP1, deSUMOylates c-MYC, resulting in its stabilization and activation. In Gujral and MacBeath [1], MYC was found to be up regulated ($+ive$ numbers) along with SENP2, and CSNK2A1 was found to show a transition from down regulation to up regulation and then down regulation. The search engine confirms the existence of this biological synergy by pointing out 3rd order combinations involving MYC-SENP2 and MYC-CSNK2A1 interaction.

7.4.3. fRizzled Class Receptor (FZD)

For FZD6, there are 5 combinations involving porcupine (PORCN) with FZD6. This depicts the possible synergy between FZD6-PORCN that need to explored. WNT ligands require palmitoylation by PORCN for their secretion and interaction with FZD receptors. Ghimire and Deans [38] in their analysis of PORCN CKOs suggest that the contribution of WNT signaling may be to establish the asymmetric distributions of FZD3, FZD6, and VANGL2 at the basolateral junctions between cochlear-supporting cells rather than as a diffusible attractant. In Gujral and MacBeath [1], FZD6 was found to be up regulated ($+ive$ numbers) along with PORCN. The search engine confirms the existence of this biological synergy by pointing out 3rd order combinations involving FZD6-PORCN interaction.

7.4.4. Glycogen Synthase Kinase 3 (GSK3)

For GSK3A (or $\alpha$), there are 3 combinations involving frizzled related protein (FRZB) with GSK3A. This depicts the possible synergy between GSK3A-FRZB that need to explored. Jaka et al. [39] found that after silencing the FRZB gene, in the case of GSK3B (or $\beta$), lower level of phosphorylation was observed, and a lower P-GSK3B/GSK3B ratio, which suggested an increase in the activity of this kinase. Gujral and MacBeath [1], FRZB was found to be down regulated ($-ive$ numbers), while GSK3A was up regulated ($+ive$ numbers). The search engine confirms the existence of this biological synergy by pointing out 3rd order combinations involving FRZB-GSK3A interaction.

7.4.5. Dishevelled Segment Polarity Protein (DVL)

For DVL1, there are 5 combinations involving histone acetyltransferase E1A binding protein p300 (EP300) with DVL1. For DVL2, there are 4 combinations involving Jun proto-oncogene, AP-1 transcription factor subunit (JUN) with DVL2. This depicts the possible synergy between DVL1-EP300 and DVL2-JUN that need to explored. Zhong et al. [40] show that EP300 mutation and loss of GATA6 function bypassed the antidifferentiation activity of WNT signaling, rendering these cancer cells resistant to WNT inhibition. They point that consistent with the WNT-dependent nature of pancreatic cancer, many components of the WNT-signaling cascade were essential for cell growth. However, the extent of the dependencies was variable, which in some cases might be due to functional redundancy (e.g., DVL1 and DVL3). So there is a glimpse of possible synergy between DVL1 and EP300, that has not been explored. Gan et al. [41] found that DVL and c-JUN form a complex with $\beta$-catenin–T-cell factor 4 (TCF-4) on the promoter of WNT target genes and regulate gene transcription. The complex forms via two interactions of nuclear DVL with c-JUN and $\beta$-catenin, respectively, both of which bind to TCF. Here, the interaction between DVL and JUN is already established. Gujral and MacBeath [1], DVL-1/2 was found to be up regulated ($+ive$ numbers), while EP300 and JUN were also up regulated ($+ive$ numbers). The search engine confirms the existence of the biological synergy between the above components by pointing out 3rd order combinations involving DVL1-EP300 and DVL2-JUN interaction.

7.4.6. Low Density Lipoprotein Receptor-Related Protein (LRP)

For LRP6, there are 4 combinations involving F-box and WD repeat domain containing 11 (FBXW11) with LRP6. This depicts the possible synergy between LRP6-FBXW11 that need to explored. Wang et al. [42] state that upon activation of the pathway by the binding of Wnt ligand to Frizzled and LRP5–LRP6 receptors, the axin complex is inhibited and results in the accumulation of soluble $\beta$-catenin that can enter the nucleus, where it interacts with transcription factors of the TCF/LEF1 family to regulate a series of target genes. It is assumed that APC helps phosphorylated $\beta$-catenin to dissociate from AXIN, creating a catalytic cycle of binding and release of the substrate. Others have suggested that APC acts either upstream of the phosphorylation reactions, by transporting $\beta$-catenin to the complex or downstream of the phosphorylation reactions, by recruiting the ubiquitin ligase bTrCP (FBXW11) to the complex. Holt et al. [43] observe that FBXW11 targets include $\beta$-catenin, key mediator of WNT signaling, critical to digital, neurological, and eye development. There might be an indirect synergy between LRP6-FBXW11. Gujral and MacBeath [1], LRP6 was found to be up regulated ($+ive$ numbers), along with FBXW11. The search engine confirms the possible existence of the biological synergy between the above components by pointing out 3rd order combinations involving LRP6-FBXW11 interaction.

7.4.7. C-Terminal Binding Protein (CTBP)

For CTBP1, there are 6 combinations involving fibroblast growth factor 4 (FGF4) with CTBP1. For CTBP2, there are 5 combinations involving catenin $\beta$ 1 (CTNNB1) with CTBP2. This depicts the possible synergy between CTBP1-FGF4 and CTBP2-CTNNB1 that need to explored. Wang et al. [42] observe that during the epithelial–mesenchymal transition (EMT) process, TGF$\beta$ induced isoform switching of FGF receptors, causing the cells to become sensitive to FGF2. Addition of FGF2 to TGF$\beta$-treated cells perturbed EMyoT by reactivating the MEK-Erk pathway and subsequently enhanced EMT through the formation of MEK-Erk-dependent complexes of the transcription factor $\delta$EF1/ZEB1 with the transcriptional corepressor CTBP1. Kim et al. [44] demonstrate that CTBP2 associates with major components of the $\beta$-catenin (CTNNB1) destruction complex and limits the accessibility of $\beta$-catenin to core transcription factors in undifferentiated embryonic stem cells (ESCs). Thus the synergies between the components have been established. Gujral and MacBeath [1], CTBP-1/2 and CTNNB1 were found to be up regulated ($+ive$ numbers), while FGF4 was down regulated for a major period of time (apart from being upregulated). The search engine confirms the possible existence of the biological synergy between the above components by pointing out 3rd order combinations involving CTBP1-FGF4 and CTBP2-CTNNB1 interactions.

7.4.8. Cyclin D (CCND)

For CCND1, there are 7 combinations involving fibroblast growth factor 4 (FGF4) with CCND1. For CCND-2/3, there are 4 combinations each involving frizzled class receptor 5 (FZD5) with CCND-1/2. This depicts the possible synergy between CCND1-FGF4 and CCND-1/2-FZD5 that need to explored. Bao et al. [45] show that CCND1 co-localizes with FGF3, FGF4, and FGF19 at chromosome location 11q13. Brandt et al. [46] show that expression levels of previously described endothelial target genes of $\beta$-catenin were studied using qPCR, but no differences was observed in the expression of CCND1 after knockdown of FZD5. But this might not be the case with CCND-2/3. Gujral and MacBeath [1], CCND-1/2/3 was found to be up regulated ($+ive$ numbers) along with FGF4, while FZD5 was found to be down regulated ($-ive$ numbers). The search engine confirms the possible existence of the biological synergy between the above components by pointing out 3rd order combinations involving CCND1-FGF4 and CCND-2/3-FZD5 interactions.

7.4.9. Wnt Family Member (WNT)

One peculiarity that we can find in the table under the title of WNT1 is that the machine points to synergistic combinations of WNT1 with other families of WNT. This pattern emerged in all the top 10 ranked combinations. It might be of interest to investigate whether WNT1 works in tandem with other WNT family members, as if one observes in the other columns such behaviour is not observed.

8. Conclusion

This study demonstrates how biologists can use the machine learning based search engine to address the needle in a haystack problem of discovering meaningful combinations of higher order in a vast search forest, which on further wet lab test might assist in intervening the pathway at a combinatorial level, in time. The problem explodes combinatorially with even a small set of recorded genes in the above study, when one steps to explore 3rd order combinations. With the total number of $C_3^{71}$ (= 57155) combinations in this study, it becomes nearly impossible for any biologist to study the system wide dynamics of any pathway. The manuscript addresses these issues by enumerating and ranking a huge list of 3rd order combinations, demonstrating conserved machine learning rankings for wet lab established combinations across the different sensitivity methods used and presenting some of the patterns in the behaviour of some of the established combinations related to WNT3A response genes. In summary, the work presents a solution to the fundamental needle in a haystack problem of locating higher order gene combinations in a vast search forest, via use of powerful machine learning based search engine. Use of this engine is bound to assist many biologists/oncologists in search for meaningful higher order gene combinations that work in cell biology and make potential discoveries necessary for advancement in the study of cell/developmental biology as well as developement of therapeutics in diseased cells.