The Expanded p53 Interactome as a Predictive Model for Cancer Therapy

The tumour suppressor gene TP53 is implicated in the majority of all human cancers, thus pivotal to genomic integrity. Even though over 72,000 PubMed publications are linked with the keyword p53 and this number is continuously increasing, due to the complexity of its interactions we are still far from fully elucidating p53’s role in tumorigenesis. Computational methodologies are novel tools to depict and dissect complex disease networks. The Boolean PKT206 p53–DNA damage model has previously demonstrated good predictive capability for p53 wild-type and null tumours in various in silico knockouts. Here, we have expanded PKT206 to generate a more clinically robust representation of p53 dynamics. The new PMH260 model incorporates 260 nodes representing genes, with 980 interactions between them representing inhibitions and activations. Additional biological outputs, including angiogenesis, cell cycle arrest and DNA repair were also amalgamated into the model. Three in silico knockouts of highly connected nodes (p53, MDM2 and FGF2) were generated and logical steady state analysis and dependency relationships determined. 71 % of predictions were considered true from superimposition of human osteosarcoma and HCT116 microarray profiles. In silico knockout analysis revealed 98 potential novel predictions, of which 13 were validated by literature; 83 % of them were overlapping with PKT206. Thus the expanded Boolean PMH260 model offers a promising platform for clinical potential in targeted cancer therapeutics.


INTRODUCTION
The important functions of the tumour suppressor gene TP53, involved in pathological states ranging from neurological to cardiovascular pathologies, underscores the requirement for elucidation of p53 complex regulatory networks.In particular this is important for cancer, where p53 mutations are implicated in over 50 % of human malignancies, and defects within its pathways are linked to tumorigenesis [1].Many of these perturbations are associated with aberrant levels of important anti-proliferative processes such as apoptosis [2] and angiogenesis [3] to which p53 has a well-defined role.Cancer at the cellular level is dynamic, involving both epigenetic and genetic alterations resulting in genomic instability and inducing malignant change [4].Treating cancer patients is not effective mostly because different molecular perturbations and tumour signatures exist dependent on the individual patient.
Even though over 72,000 PubMed publications are linked with the keyword p53 and this number is increasing, we are still far from understanding the molecular mechanisms of the p53 system, which contribute to the malignant phenotype.This myriad of information along with p53 network intricacy is challenging to amalgamate.Relevant information about p53 actions at the molecular and physiological level must be integrated with the cell environment that is pivotal to development of optimal individualized targeted therapies.In order to gain increased knowledge of how p53 transcriptional regulatory networks govern genomic aberrations and malignant transformation, a systematic approach is fundamental to assimilate the diverse information into a coherent and holistic framework.
Cancer systems biology is a relatively novel, yet promising field that can utilize genome wide data to analyse intracellular network perturbations [5].These integrative approaches can systematize and interpret experimental data of multiple dimensions to elucidate and predict molecular principles and deviations of the underlying cancer phenotype.These models may be manipulated in silico to mimic external signals and mutations observed in vivo.Predictive models can be of both clinical and pharmaceutical relevance for the development of novel effective therapies.Indeed, several computational studies of different mathematical complexities have already been developed in an effort to address this [6].
Concerning mathematical models, Boolean networks are a promising predictive framework and have been applied successfully to model various biological phenomena [7, 8 and 9].Boolean logic states that a node is assigned discrete expression values which qualitatively predict the temporal evolution of the system, where at each time point a node state is determined by the state of other upstream nodes by transfer of a Boolean function.This contrasts with continuous ranges of kinetic models, which incorporate more parameters that can be difficult to measure [10].Boolean states of '1'/'ON' or '0'/'OFF', referring to for example activation of transcription or not, can effectively recapitulate gene expression in silico.Furthermore, feedback loops which are frequent and important signalling processes may also be considered.
The PKT206 Boolean model, comprising 206 nodes with 738 interactions [8], has previously demonstrated predictive potential of p53 network dynamics, with 51-75 % correct predictions derived when compared to various microarray profiles.Furthermore, application of the novel signal transduction score flow algorithm to PKT206 resulted in greater predictive ratios (82 %) when compared to identical 'omics' data [11].Both these qualitative (STSFA) [11], and quantitative (Boolean), algorithms [8] have demonstrated that PKT206 successfully predicts overall system properties.However, the p53 network is extensive with a plethora of p53 responsive genes that govern multivariate pathways.Consideration of this is fundamental for a better molecular depiction of the p53cancer network.Therefore in this study we have expanded PKT206 to PMH260 (p53 Michelle Hussain, 260 nodes), which includes 260 nodes with 938 interactions between them of inhibition or activation.For clinical relevance additional biological outputs were also introduced into the model to include: apoptosis, angiogenesis, cell cycle arrest, cellular senescence and DNA repair.In order to test the predictive capability of the extended p53 model we superimposed 'omics' data to PMH260.Furthermore, three highly connected nodes of the network were subjected to in silico knockout tests of different p53 and DNA damage statuses to predict network perturbations.

Generation of PMH260
PKT206 generation was described in detail elsewhere [8].Briefly, PKT206 considers 203 nodes representing genes / proteins, connected with 738 interactions of inhibition or activation derived from the STRING database.These binary interactions were connected to two downstream outputs: apoptosis and cellular senescence, and one upstream DNA damage input.All nodes were connected to each other through direct interactions.The PKT206 model includes all direct interactions with p53 and all direct interactions between genes/proteins that interact with p53.The model was analysed using CellNetAnalyzer as described below [13].
For construction of PMH260, interactions and nodes were derived from the STRING database (9.1), available at (http://string-db.org) [12].Additional nodes (n = 54), and interactions between new and previous nodes (n = 242) from PKT206 were amalgamated into the new model.Interactions were filtered by a confidence score over 0.7, which is a threshold for high confidence assigned by STRING.Manual curation of these interactions was also undertaken by extensive scientific literature search in order to ensure model accuracy.

Addition of biological outputs
Three additional biological outputs of angiogenesis, cell cycle arrest and DNA repair were included into the PMH260 model.Dependent on the biological processes they regulate, nodes (excluding input and outputs) (n = 254) were connected to these outputs by edge functions of inhibition, activation, or ambivalent factor.These were initially selected using Gene Ontology (GO) terms and were additionally validated by literature to confirm the interactions.As the PKT206 model included only one input node (DNA damage) and two output nodes (apoptosis and cellular senescence), the additional biological outputs of angiogenesis, cell cycle arrest and DNA repair were added into PMH260.Then, relevant nodes downstream of p53 within the PKT206 model (n = 203) were also linked to the above mentioned three outputs.

In silico analysis of PMH260
CellNetAnalyzer (CNA) is an independent software platform designed for MATLAB enabling the user to run various algorithms for structural and functional analysis of biochemical, metabolic, and signalling networks [13].The model is represented by an interaction hypergraph where nodes are connected by interactions of binary effects (activation or inhibition).Interactions can be combined by operators of: 'AND' (where a hyperarc connects several input nodes), 'OR' (several arcs connect to the same output node) and 'NOT' (opposite state of the source node activates the target node).For in silico analysis the PMH260 network was generated in CNA and two techniques were utilized: Dependency Matrix calculations and Logical Steady State Analysis (LSSA).
Dependency relationships between species may be analysed with the construction of a dependency matrix (M) that displays all pair-wise node dependencies.
These relationships are calculated by the shortest distance between node pairs, represented by six different values indicating relationships between nodes as defined by: i.If no negative or positive path exists between node i and node j, i has no effect on j.
ii.If both a negative and positive path exist from node i to node j, i is an ambivalent factor of j.
iii.If only negative paths exist from node i to node j with negative feedback loops present, i is a weak inhibitor of j.
iv.If only positive paths exist from node i to node j with negative feedback loops present, i is a weak activator of j.
v. If only negative paths exist from node i to node j with negative feedback loops absent, i is a strong inhibitor of j.
vi.If only positive paths exist from node i to node j with negative feedback loops absent, i is a strong activator of j.

Knock out analysis of highly connected nodes
To investigate node dependency relationships in response to various in silico knockouts (KOs), the connectivity of all nodes was determined using CellNetAnalyzer (3.1) [13], and visualised in Cytoscape (2.8) [14], and the three most highly connected nodes in the network chosen.Selected nodes and corresponding edges were deleted from the model to represent in vivo KO or mutation.Depending on the effect of a particular node deletion, the relationships between remaining nodes within the network can be defined by the six above-mentioned values.

Application of Logical Steady State Analysis
Logical Steady State (LSS) analysis is a method for predictive inputoutput relationships in signalling networks [14].A steady state is reached when the state of nodes is fully consistent with its Boolean function, and as such once a network has reached a LSS it will remain there over time.The state of some nodes may remain undetermined if multiple LSS solutions are possible for the same input conditions.
Files listing proteins and interactions were generated for PMH260 in accordance with [14].Four different scenarios were constructed which represent in vivo processes such as loss of p53 function due to mutation (mimicked in silico by removal of the p53 node) in the presence and absence of DNA damage (Table 1).
Deletion of p53 from the network was undertaken by removing the node and all associated edges to generate a p53 KO scenario.The DNA damage input was switched to '0'/'OFF' or '1'/'ON' depending on input required.

Genome wide validation of PMH260
Predictions that were derived from LSSA of the various in silico scenarios were compared to microarray data of human osteosarcoma U2OS (p53 + / + ) and SaOS2 (p53 -/ -) cell lines exposed to etoposide-induced DNA damage and untreated conditions from [8].To further test the model's predictive power, 'omics' profiles of HCT116 p53 wild type and null human cell lines available from the Gene Expression Omnibus (Accession number GSE10795) were also investigated.This analysis was conducted in accordance with [8,11] and is described below.The predicted state of a node, (i) for p53 wild-type and p53 mutant were defined as S(i)wt and S(i)mu respectively.Both could take values of 0, 1 or NaN.A variable Emod represented the predicted change of gene state from p53 wild type to mutant where: The expression level of a gene (i) derived from experimental validation was defined by the parameter Eexp where: The fold change (FC) between experimental and simulated data was calculated for each gene, defined by the equation: Where M1(i) is the median of expression values in the target condition and M2(i) is the median of expression values in the source condition.Threshold values (θ) were applied to normalise expression profile distributions, using (x) and (σ) of the log10 FC scores, of (θmax and θmin), defined as: In Gene activity obtained from microarray data was further determined by the equation: log10 FC > x + σ: Up-regulated log10 FC < x -σ: Down-regulated x -σ <= log10 FC <= x + σ : No change Where x represents the mean and σ the standard deviation.
For comparison of both datasets, we analysed changes in each gene response between experimental and in silico conditions.The predicted change of gene activity may be defined by a variable Emod of three states.Differences between model simulations (Emod) and experimental observations (Eexp) were defined as (Emod -Eexp) and could take on one of three values (0, 1 or 2): a correct prediction corresponds to the value of 0, which means that both experimental and simulated outcomes were the same; a small error prediction corresponds to the value of 1, which means that one outcome was 'no change' and the other 'up or down-regulated'; a large error prediction corresponds to the value of 2, which means that one outcome was 'up-regulated' and the other was 'down-regulated'

PMH260 network topology
The PMH260 network comprises 260 nodes in total with 980 interactions of inhibition or activation.The DNA damage input to the model serves as environmental stress of which 42 nodes are linked to.Of these, five are inhibited by DNA damage and the remaining activated.34 feedback loops were identified to which p53 participated in over 50 % (n = 18).117 nodes were located upstream of p53, with 164 downstream.Of the total 254 internal nodes, 27 nodes were considered as ambivalent and these were located both up and downstream of p53.Among the outputs, 28 nodes were linked to DNA repair, 62 to cell cycle arrest, 68 to cellular senescence, 42 to angiogenesis and 118 to apoptosis (Figure 1).

PMH260 network analysis
We investigated the LSS of the PMH260 model by considering four in silico comparative scenarios of different DNA damage and p53 status to identify network perturbations resulting from loss of p53.These scenarios are summarized in Table 1.
204 node states out of 260 (78.5 %) were determined in p53 wild-type background with DNA damage OFF, meaning that 149 nodes were upregulated and 56 nodes were down-regulated (Figure 2A, lanes 4 and 5).Thus, 55 nodes had undetermined states ('NaN'), (Fig 2A , Lane 6).These undetermined states are obtained where several stable states may exist for the same input conditions.Similarly, 205 states were determined out of 260 nodes (78.8%) in p53 wild-type background with DNA damage ON (Figure 2A, lanes 1 and 2).In comparison, reduced system stability was observed in p53 KO models with over half of all nodes undetermined from the total: 147 (57 %) and 148 (57.1 %) for DNA damage ON and OFF respectively, when compared to p53 wild type (Figure 2A; lanes 9 and 12).All output signals remained active across all scenarios.
In the presence of p53, a greater number of upregulated nodes were observed compared to p53 KO backgrounds: 164 (63 %) when DNA damage was ON in the presence of p53, compared to 101 (39 %) when p53 was deleted with DNA damage ON (Figure 2A; compare lanes 1 and 7).
Nearly a 4-fold increase was observed in the number of down-regulated nodes in the presence of DNA damage in p53 wild type background, 41 (15.7 %), when comparing to p53 KO, 11 (4 %), (Figure 2A, compare lanes 2 and 8).Only a 2-fold increase in the number of down-regulated nodes was seen when DNA damage was OFF in the presence of p53, 56 (21.5 %), compared to when DNA damage was OFF in p53 KO background, 30 (11.6 %), (Figure 2A; compare lanes 5 and 11).

LSSA of angiogenic and apoptotic node states
As both angiogenesis and apoptosis are pivotal to tumour formation and progression, we focused on investigating the effect of their node state changes in response to the four in silico comparative scenarios presented in Table 1.The Gene Ontology database was used to annotate biological processes involved in angiogenesis and apoptosis for each node.Of the apoptotic nodes, 45 were anti-apoptotic and 73 pro-apoptotic nodes, constituting 28 % and 17 % of the total network respectively (Figure 2C).Globally, increased up-regulation of pro-apoptotic nodes was observed in p53 wild type, compared to p53 KO backgrounds (Figure 2C; compare lanes 1 and 7 to lanes 13 and 19).Across all scenarios a greater shift of nodes to undetermined states was observed when p53 was removed from the network (Figure 2C; compare 6 and 12 to lanes 18 and 24).The numbers and percentages of expressed angiogenic and apoptotic nodes under these comparative conditions are summarized in Supplementary Table 2.

Distribution of pro-angiogenic node states
Increased pro-angiogenic node up-regulation (n = 21) was observed in p53 wild type backgrounds compared to when p53 was deleted (n = 17) (Figure 2B; compare lanes 1 and 13).This was regardless of DNA damage presence as it was also noted in scenario II (Table 1 for p53 wild type with DNA damage OFF, with p53 knockout with DNA damage OFF having 21 and 16 nodes up-regulated respectively (Figure 2B; compare lanes 7 and 19).
When DNA damage was ON, only a single node (STMN1) was down-regulated in the p53 KO scenario.This was compared to 4 (14.8 %).when p53 was present.(Figure 2B; compare lanes 14 and 2).The growth factor PDGFRB was up-regulated when p53 was deleted from the network when DNA damage was ON.
Some nodes were considered DNA damage dependent.For example, the oncoprotein STMN1 was exclusively up-regulated in the absence of DNA damage and down-regulated in its presence.In fact, STMN1 down-regulation was observed across all scenarios regardless of p53 status.Similarly, the proto oncogene JUN was repressed only when DNA damage was OFF and enhanced in its presence.
Four pro-angiogenic nodes (14.8 %) were down regulated in p53 wild-type backgrounds regardless of DNA damage (BDKRB1, EPHB4, IGF1R and PDGFRB).Up-regulation of nodes ELAVL1 and PDGFRB along with JUN was observed in response to DNA damage, but they were repressed in its absence in p53 KO backgrounds.The ambivalent factor FGF2 was down-regulated in p53 KO background in response to DNA damage (see Supplementary Document 1).

Distribution of anti-angiogenic node states
There were no major differences in antiangiogenic node state changes across all comparative scenarios.No down-regulation of nodes was observed in the presence of p53 (Figure 2B; see lanes 5 and 11).Two anti-angiogenic nodes were down-regulated (tumour suppressors PML and ELAVL1) when DNA damage was OFF in the absence of p53 (Figure 2B, lane 23).The number (n=11) (73.3 %) and particular nodes up-regulated remained identical in p53 wild type backgrounds, regardless of DNA damage input, (Figure 2B; lanes 4 and 10).Anti-angiogenic node upregulation was slightly increased in the p53 KO background (n = 6) (40 %) when DNA damage was ON to 4 (27 %) when DNA damage was OFF (Figure 2B; compare lanes 16 and 22).The only differentially up-regulated node across all scenarios was the tumour suppressor ING5 which was enhanced only in p53 KO backgrounds in response to DNA damage

LSSA of apoptotic node states caused by changes in p53 and DNA damage
The majority of pro-apoptotic nodes were upregulated when p53 was present compared to p53 deletion backgrounds.In fact, nearly a 3-fold increase of up-regulated pro-apoptotic nodes (n = 53) (72.6 %), was seen in the presence of p53, compared to when p53 was deleted from the network (n = 20) (27.4 %), (Figure 2C; compare lanes 1 and 7 to lanes 13 and 19).This was observed in both scenarios I and II (Table 1, compare columns 1 and 2).Thus, this was regardless of DNA damage input.Some nodes were exclusively repressed when p53 was deleted and only enhanced in its presence, such as the DNA damage response gene GADD45A.
Investigating the effects of DNA damage OFF and ON in the presence of p53 (scenario II of Table 1), the majority of pro-apoptotic nodes exhibited no change; the number of up-regulated nodes for DNA damage ON was 56 (77 %), and 53 (72.6 %) for DNA damage OFF (Figure 2C, see lanes 1 and 7).The majority of up-regulated nodes were identical in both conditions (Supplementary Document 1).

Distribution of anti-apoptotic node states
Interestingly, a greater number of anti-apoptotic nodes were up-regulated in the presence of p53 than in p53 KO backgrounds.This was regardless of DNA damage.However, over 50 % of these anti-apoptotic up-regulated nodes are also pro-apoptotic including: DAXX, DDIT4, DUSP2, DUSP4, EGFR, ESR1, FHL2, and PRKCA.In addition, several of these nodes also have roles in angiogenesis, and include: EGFR, PRKCA and PTGS2 (Supplementary Document 1).
Comparing up-regulation of anti-apoptotic to proapoptotic nodes, nearly a 2-fold increase among proapoptotic nodes is observed when p53 is present (Figure 2C; compare lanes 1 and 7 to lanes 4 and 10).Of the 31 anti-apoptotic nodes up-regulated in p53 wild type scenario, the majority were also considered pro-apoptotic, with additional roles in angiogenesis.
Down-regulation of anti-apoptotic nodes was greater in p53 wild-type backgrounds compared to p53 deletion backgrounds (Figure 2C; compare lanes 5 and 11 to lanes 17 and 23): 7 (15.5 %) in p53 wild type with DNA OFF, and 5 (11 %) in p53 WT with DNA damage ON, compared to 2 (4.4 %), and 4 (8.8 %), for DNA damage ON and DNA damage OFF respectively in p53 deletion backgrounds.Of these in p53 WT backgrounds, the majority of nodes were identical apart from the proto oncogene BCL3 and the growth factor regulated gene EPHB4, which were downregulated in the presence of p53 (Supplementary Document 1).Some nodes were up-regulated across all conditions such as the anti-apoptotic factor BCL2.
Comparing the effect of DNA damage ON and OFF in p53 wild type backgrounds (scenario III of Table 1), the number of up-regulated nodes was similar across both conditions: 32 (71 %) in response to DNA damage and 31 (68 %) in its absence (Figure 2C; see lanes 4 and 10).Of these, the majority of nodes were identical (Supplementary Document 1).

Dependency changes in response to in silico removal of highly connected nodes
The most highly connected three nodes within the model, being p53, MDM2 and FGF2, were chosen for in silico KO analysis to investigate network perturbations in their absence, mimicking in vivo loss of function mutations.The effect of a particular node deletion can result in relationship changes within the network, as defined by the six values given in the methods section.These changes are presented in Table 2.
The total number of dependencies (n = 67,600) (260×260) represents elements in the dependency matrix of the p53 wild-type model.Of these, 40653 correspond to interactions with no effect, 20988 are ambivalent factors, 2761 are weak inhibitors, 3048 are weak activators, 48 are strong inhibitors and 102 strong activators.For any knockout (KO) scenario the total number (n = 67,081) represents the elements within the network after a particular node deletion (259x259).We focused on changes in strong activators or strong inhibitors as these have the greatest effect on the cell.p53 KO had the strongest effect when compared to the wild-type interactome.When p53 was deleted from the network, the number of strong inhibitors nearly doubled from 48 in p53 wild type to 71 in p53 knockout, (Table 2).The p53s negative regulator MDM2, which is contained within an intricate feedback loop with p53, induced only two changes when deleted.
The majority of the dependency changes in p53 KO scenario were altered to strong activators (n = 65) and strong inhibitors (n = 23), from ambivalent factors in PMH260 wild type.As noted in Supplementary Table 1, two changes from weak activators in wild type to strong activator in p53 KO were noted, being DYRK2 and MAPK1WA.Interestingly, several nodes changed from having no influence on angiogenesis in the wild type scenario to having a strong inhibitory (COL181A, PPARG and PML), or activating (MAPK1WA, MMP2, MMP13, NTN1, S100B, SGK, PRAK and STMN1) angiogenic effect in the absence of p53.
Similarly, TFDP1, GTSP1 and EZH2 had no influence on cell cycle arrest in the presence of p53, however when p53 was excluded from the network all three became positive regulators of cell cycle arrest.Only one weak inhibiting interaction to strong activator was observed in p53 wild type to KO, namely for PPMID targeting CHK1.For MDM2 KO, one dependency change was observed within the whole network, namely for ATM to DYRK2.These relationship changes are summarized in Supplementary Table 1.

Genome wide validation of PMH260
The PMH260 model has derived potential novel predictions during in silico analyses of LSSA and dependency changes.Some of these predictions were confirmed either by previous literature or verified in the laboratory.Even so, only a minority of predictions could be experimentally verified.To investigate the predictive efficiency of PMH260, microarray profiles of etoposide treated and untreated osteosarcoma SaOS2 (p53 -/ -) and U20S p53 (p53 + / + ) were compared against LSSA data.For further investigation, microarray profiles of HCT116 colon cancer cell lines, (p53 null and wild type) were also compared.
The fold change and expression profiles of both microarray and in silico datasets were calculated according to the formula described in [8, 10 and 25].Genes classified as 'small error' or 'large error' (see methods) were ranked based upon the number of times their response was incorrectly predicted.A 'small error prediction', was assigned a score of 1, a 'large error prediction' a score of 2. Incorrect gene scores for all predictions across all simulations were added to produce a total score.For example, if a gene was given a large error prediction in three conditions, and a small error prediction in one condition, its total score would be 3*2 + 1 = 7.A score of 4 was defined as the threshold to consider that the response of a gene was consistently incorrectly predicted by the network.To achieve this score, a gene must receive a small error prediction in over half of the simulations, or a large error prediction more than once, or a large error prediction and two small error predictions.
The true prediction percentage ranged from 55 to 71 % with an average of 61.2 % across all simulations when compared to microarray data.Large errors were in the minority, corresponding to less than 5 % of the total predictions, with a mean of 2.9 %.Small errors comprised 35 % for all simulations (Table 3).

DISCUSSION
The tumour suppressor p53 maintains genomic integrity [27].Dissecting the role of p53 and its regulatory networks is pivotal to a greater understanding and subsequent implementation of successful therapies.Systems biology teamed with traditional reductionist approaches are novel, yet promising tools to describe the complexity of diseased systems such as cancer [5].
PKT206 has demonstrated good predictive capability for p53 -DNA damage pathways.Even so, the p53 network is extensive, with over 1000 p53 responsive genes described [28].To capture these, a larger model is necessitated.Extensive models allow for a global system overview and thus better representation of biological phenomena, which can effectively capture sub networks and mutation drivers.For correct representation, models must incorporate and depict the heterogeneity of diseased states and the processes that govern them.Angiogenesis and apoptosis for example are pivotal properties to the tumour environment.These processes are involved in cell signalling and thus consideration of these is crucial to capture tumour dynamics [5].Here we have incorporated an additional three outputs allowing for the many processes involved in p53 governed tumorigenesis to be modelled and analysed systematically.Thus, the expanded PMH260 model should have better clinical relevance.
To mimic in vivo mutations and elucidate the role of deleted nodes on PMH260 network dynamics, three in silico KO tests of highly connected nodes were performed.The greatest effect was observed when p53 was removed from the network deriving 98 predictions.In comparison, deletion of the p53 negative regulator MDM2, resulted in two changes into strong activator dependency changes, whereas no change in strong dependencies occurred for FG2 removal; numerous changes in other dependencies occurred across all knockouts.This is due to the fact that p53 is the most highly connected 'hub' node of the model and participates in the majority of feedback loops.The majority of dependency changes are potential novel predictions.13 of these 98 predictions were confirmed by literature search or experimentally verified [8, 15-19 and 21-26].Of the 24 dependency changes described in PKT206 for p53 deletion, 22 (83 %) in PMH260 were in agreement.This highlights the accuracy and reproducibility of the p53 interactome.Most of the prediction differences derived from dependency simulations were from new nodes introduced into PMH260.Several nodes were altered from no angiogenic effect to having a positive influence on angiogenesis, similarly two nodes changed from no influence in the wild type to positive cell cycle regulators in the p53 mutant.This indicates the role of p53 in maintaining the stability of the whole network, highlighting its well defined role as a tumour suppressor, and offers therapeutic potential for p53 null and wild type tumours.
As apoptosis and angiogenesis are clinically significant, we investigated the effects of both pro and anti-apoptotic and angiogenic changes in response to different p53 and DNA damage statuses using LSSA.
Increased up-regulation of pro-apoptotic nodes were observed compared to anti-apoptotic in p53 wild type backgrounds.This suggests that in the presence of p53 there would be a favourable apoptotic increase.Many of these up-regulated anti-apoptotic nodes are ambivalent factors towards apoptosis: for example, over 50 % of up-regulated anti-apoptotic nodes were also pro-apoptotic.PMH260 at present is a global cancer model and thus differential expression for a particular gene may be cell / tissue specific.Indeed, p53 is known to transcriptionally regulate target genes in a context specific manner [29].This however provides a window of therapeutic opportunity by inhibiting anti-apoptotic nodes of specific tumour phenotypes.Those anti-apoptotic plus pro-angiogenic nodes are ideal candidates for targeted inhibition.For example, SGK, TGFA, TPT1, and VEGFA display such activity.Indeed several anti VEGFA drugs have been described with clinical trials undertaken [30].
As PMH260 revealed both important established and novel predictions that could be of therapeutic interest, investigation of its predictive accuracy was needed.Thus we superimposed and compared osteosarcoma and HCT116 expression profiles to the p53 interactome.Good prediction ratios were derived with 58-71 % of correct predictions across all simulations, which significantly exceeds an expected probability of 33 % for a random model which has three possible prediction outcomes with equal probability.Importantly, combining true predictions and small errors resulted in over 90 % of successful predictions.
Predictive ratios here are similar to those obtained from PKT206.This is promising, as the model increases so does network complexity, which is often difficult to simulate and analyse.Comparing differential expression profiles of apoptotic and angiogenic factors in silico to 'omics' p53 wild type and null profiles has derived several important predictions of therapeutic relevance.For example, pro apoptotic GADD45A was down-regulated only in p53 KO (SaOS2) backgrounds.Whilst anti-apoptotic JUN and ambivalent factor DAXX were DNA damage responsive and up-regulated in both p53 wild type and KO backgrounds.In parallel with [8], growth factors were also uncovered as contributing factors to osteosarcoma; the proangiogenic and apoptotic ambivalent factor FGF2 was down-regulated in SaOS2 cells under DNA damage, whilst anti-apoptotic PDFGRB and IGF1R, along with growth factor regulated genes EPHB4 and PTGS2, were up-regulated in SaOS2 cells under DNA damage and HCT116 cells, suggesting cell specific regulation.Furthermore, several novel predictions were obtained which may be of relevance contributing to the p53 null tumour phenotype, For example, TFDP1, GSTP1 and EZH2 altered to positive cell cycle regulators upon p53 deletion.Similarly, genes changed to pro-angiogenic factors in p53 KO scenario.
In summary, the PMH260 model provides increased coverage of potential expression changes in cancer.Validation of its predictive power through comparison to microarray profiles showed good results.This highlights the predictive performance of large-scale Boolean models for the investigation of DNA damage inducible pathways.Several important prospective predictions have been reported here, which offer putative anti-cancer therapies for the inhibition or enhancement of target genes or perturbed pathways in both p53 null and wild type tumours in osteosarcoma and colon cancers.This offers a promising platform for personalized targeted therapies where individual 'omics' patient's profiles may be superimposed for single or multiple mutations.

Figure 1 .
Figure 1.The PMH260 p53-DNA damage model.Five layers are observed in accordance with node function; the input signal of DNA damage (green), upstream of p53 (yellow), the network hub, p53 contained within the crucial feedback loop with MDM2 (red), downstream of p53 (turquoise) and five downstream outputs (orange).Red lines signify inhibition, blue activation.

Figure 2 .
Figure 2. Distribution of all PMH260 nodes under LSSA.A) All nodes in each condition.B) Nodes regulating angiogenesis under LSSA.C) Nodes regulating apoptosis under LSSA.DD: DNA damage, WT: Wild type.KO: Knockout.

Table 2 . Distribution of changes in dependency matrix elements in response to three in silico KO tests
. A greater effect on the network was observed in the absence of p53 than other KOs.Of the 6 defined relationships, the majority of changes were from ambivalent to strong activators.

Table 3 . Percentage of all predictions from comparison of in silico data to microarray profiles.
Correct predictions were in the majority across all simulations with large errors occupying a small percentage of results when comparing osteosarcoma and colon cancer 'omics' profiles to in silico data under LSSA.The expected correct prediction rate for a random model is 33 % since there are three possible prediction outcomes (up, down or unchanged) with equal probability in a random model.