A Novel Prognostic Metric Resolves the MYCN enigma in Silico and Points to a Biosynthetic Regulatory Shift Driven by MYCN Amplification in Neuroblastoma

A novel prognostic metric resolves the MYCN enigma in silico and points to a biosynthetic regulatory shift driven by MYCN amplification in neuroblastoma | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results A novel prognostic metric resolves the MYCN enigma in silico and points to a biosynthetic regulatory shift driven by MYCN amplification in neuroblastoma Matteo Italia , Melody Parker , Rory Deignan , Fabio Dercole , Dawn Walker , View ORCID Profile Kenneth Y. Wertheim doi: https://doi.org/10.1101/2025.11.13.688307 Matteo Italia 1 Department of Electronic, Information, and Bioengineering, Politecnico di Milano , Milano, Italy . 2 Mathematical Oncology Laboratory (MOLAB), Instituto de Matemática Aplicada a la Ciencia y la Ingenieŕıa, University of Castilla-La Mancha , Spain . Find this author on Google Scholar Find this author on PubMed Search for this author on this site Melody Parker 3 Nuffield Department of Medicine, John Radcliffe Hospital, University of Oxford , Headley Way, Oxford, UK . 4 Insigneo Institute for in Silico Medicine, University of Sheffield , Sheffield, UK . Find this author on Google Scholar Find this author on PubMed Search for this author on this site Rory Deignan 4 Insigneo Institute for in Silico Medicine, University of Sheffield , Sheffield, UK . 5 Department of Physics, University of Sheffield , Sheffield, UK . Find this author on Google Scholar Find this author on PubMed Search for this author on this site Fabio Dercole 1 Department of Electronic, Information, and Bioengineering, Politecnico di Milano , Milano, Italy . Find this author on Google Scholar Find this author on PubMed Search for this author on this site Dawn Walker 4 Insigneo Institute for in Silico Medicine, University of Sheffield , Sheffield, UK . 6 Department of Computer Science, University of Sheffield , Sheffield, UK . Find this author on Google Scholar Find this author on PubMed Search for this author on this site Kenneth Y. Wertheim 4 Insigneo Institute for in Silico Medicine, University of Sheffield , Sheffield, UK . 6 Department of Computer Science, University of Sheffield , Sheffield, UK . 7 Centre of Excellence for Data Science, Artificial Intelligence, and Modelling, School of Digital and Physical Sciences, University of Hull , Kingston upon Hull, UK . 8 Computer Science, School of Digital and Physical Sciences, University of Hull , Kingston upon Hull, UK . Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Kenneth Y. Wertheim For correspondence: k.y.wertheim{at}hull.ac.uk Abstract Full Text Info/History Metrics Data/Code Preview PDF Abstract Neuroblastoma (NB) is the most common extra-cranial solid tumour in children. Although MYCN amplification is typically indicative of a poor prognosis, the expression level of MYCN has a non-monotonic relationship with the clinical outcome. This paper proposes an explanation for this phenomenon, which is called the MYCN enigma in the literature. The p53/MYCN metric, which is a measure of the p53 protein level relative to the MYCN protein level in a NB tumour, is the key concept. Our hypothesis is that the metric has a positive relationship with the outcome of a patient. It is presented as a binary classification system, where the prognosis is favourable when p53/MYCN > 3.5. The mathematical model presented in this paper describes the dynamics between MYCN, p53, ARF, and MDM2; their genes; and their mRNA transcripts. Simulations were carried out by solving the model numerically in a series of initial value problems. The results are aligned with a list of clinical and experimental observations. We extracted association rules with the Apriori algorithm and explored the parametric space stochastically. Assuming that MYCN enables biogenesis in NB tumours, our results support the prediction that MYCN amplification shifts the balance towards producing MYCN. When a wild-type tumour is deficient in MYCN, stress responses may restore biogenesis in general and preferentially produce MYCN in a negative feedback loop. Another prediction is that MYCN - amplified tumours, without treatment, require hard-to-attain biosynthetic rates and stress responses to achieve good outcomes. 1 Introduction Neuroblastoma (NB) is the most common extra-cranial solid tumour in children [ 1 ]: the median age of diagnosis is 573 days [ 2 ]. Impaired stem and progenitor cells in the neural crest (a transient structure in the human embryo) are responsible for its patho-genesis [ 3 – 5 ]. Due to mutations such as MYCN amplification and ALK activation [ 4 , 6 ], they cannot migrate and differentiate into more specific lineages, such as the sympathoadrenal lineage, which generates the sympathetic nervous system. Instead, they become NB cells and collectively, a malignant solid tumour [ 4 , 6 ]. A patient’s risk group is determined based on multiple criteria, including their stage defined by the International Neuroblastoma Risk Group (INRG) classification system, their age, their primary tumour’s histological type, and additional clinical and molecular parameters [ 7 ]. Patients from the low- and intermediate-risk groups are expected to have a five-year overall survival rate above 70 %, but high-risk patients are expected to succumb to the disease or relapse even after multi-modal therapy at a rate above 50 % [ 8 ], whereupon survival for more than three years is rare [ 9 , 10 ]. Emerging strategies such as personalised therapies [ 11 , 12 ], combination therapies [ 13 ], targeted inhibitors [ 14 , 15 ], and virotherapy [ 16 ] aim to improve survival outcomes. MYCN amplification in NB cells is considered a biomarker indicative of a poor prog-nosis. However, it is not well-defined. In some studies [ 17 , 18 ], a cell was considered MYCN -amplified (MA) when it had 10 or more copies of the gene (measured by FISH assay signals) and the signals surpassed at least three times the number of control signals (PAX3 probe, 2q35). In other studies [ 19 , 20 ], a four-fold increase was used as the threshold. According to one source [ 21 ], this mutation is generally present in around 25 % of all NB cases. It is also a defining characteristic of high-risk cases [ 22 ]. As an oncogene, MYCN can activate genes involved in metastasis, survival, proliferation, pluripotency, self-renewal, and angiogenesis [ 21 ]. It can also switch off genes responsible for differentiation, cell cycle arrest, and immune surveillance [ 21 ]. Importantly, [ 6 ] observed in their mouse xenografts that MYCN amplification and a specific activating mutation of ALK were sufficient to drive NB development from neural crest progenitor cells. Counter-intuitively, the abundance of MYCN mRNA in 91 primary NB specimens was found to vary non-monotonically with that of TrkA mRNA [ 23 ], a favourable prognostic indicator. It should be noted that being a favourable biomarker does not automatically mean TrkA mechanistically suppresses neuroblastoma. In other words, the correlation between the two quantities was found to be positive in some cases and negative in the rest. One possibility is that the relationship depends on another molecular species or even the broader cellular context. Out of the 91 specimens, the wild-type (WT) tumours without MYCN amplification displayed a much broader range of TrkA expression levels than their MA counterparts. Furthermore, the WT tumours with more MYCN mRNA also expressed more TrkA mRNA (favourable), but the MA tumours expressed less TrkA mRNA (unfavourable) than the WT tumours despite expressing even more MYCN mRNA. One interpretation is that an intermediate MYCN mRNA level leads to a favourable outcome for some unknown reasons. Motivated by these contradictory results, another study designed a gene signature called MYCN -157 [ 24 ]. It predicts poor outcomes independently of MYCN amplification and expression. In the same study, tumours with similar MYCN mRNA levels were found to have very different MYCN protein levels. The results reported in these two papers [ 23 , 24 ] form the basis of the MYCN enigma. This paper proposes a hypothesis to explain this enigma. Here, we define the MYCN enigma as the general observation that a NB patient’s clinical outcome does not depend monotonically on their tumour’s MYCN expression level or MYCN amplification status. We have compiled a list of specific clinical and experimental observations (section S1 in the supplementary file). A satisfactory explanation must be consistent with them. Our hypothesis is that the p53 protein level in a NB tumour, relative to the MYCN protein level therein, decides whether it is likely to progress or regress. Henceforth, we will simply call this metric p53/MYCN. This hypothesis is based on the known fact that MYCN (protein) activates p53 [ 25 ] and on the consensus that p53 is a tumour suppressor gene with diverse downstream effects, including DNA repair [ 26 ], cell cycle arrest [ 27 – 29 ], and apoptosis [ 30 , 31 ]. Although MYCN (gene) is typically considered an oncogene [ 21 ], it activates a tumour suppressor gene ( p53 ) too, suggesting a non-monotonic relationship mirroring the MYCN enigma. Intuitively, this tension could be the mechanism behind the phenomenon. If our hypothesis was true, a NB with low p53/MYCN would have a poor clinical outcome. For this hypothesis to be consistent with [ 23 ]’s dataset, p53/MYCN should follow a concave function of MYCN ’s mRNA level in a mixed collection of WT and MA NBs. Although we found numerous modelling studies involving MDM2 and p53 in the literature [ 32 – 38 ], and even some investigating their interactions with ARF and other interaction partners of p53 (such as ATM and WIP1) [ 39 – 43 ], our literature review led to the conclusion that a novel mathematical model coupling MYCN to the ARF/MDM2/p53 axis was needed to explore our hypothesis. We formulated eight non-linear ordinary differential equations (ODEs) to model the dynamics between MYCN, p53, ARF, and MDM2, as well as their corresponding genes and mRNA transcripts. The model describes a stream of genetic information flowing through a molecular net-work connecting MYCN to the ARF/MDM2/p53 axis, which is highly relevant to NB research [ 44 ]. It is also a suitable context for our task. TrkA and its associated products are not in the model because they do not interact with the network directly. There is also no evidence that MYCN activates TrkA like it activates p53 , so the p53/MYCN metric is mechanistically more relevant, making it a more direct proxy for a favourable prognosis than the mRNA level of TrkA . After parameterising the model, a series of computer simulations were carried out. Based on the results, a binary NB classification system was established: the prognosis is favourable when p53/MYCN > 3.5. It predicts whether a patient will have a good or poor outcome based on their p53/MYCN met-ric. The results are consistent with most of the clinical and experimental observations listed above too. The Apriori algorithm [ 45 ] and an exploratory data analysis helped us infer plausible mechanisms underlying the hypothesis. 2 Methodology We will begin this section by presenting a model comprising eight non-linear ODEs describing a gene regulatory network formed by MYCN , ARF , MDM2 , p53 , as well as their mRNA transcripts and proteins. After that and a summary of the computer simulations we performed, we will describe in detail the procedures we followed to establish the NB classification system based on our hypothesis and simulation results. We will explain how we extracted further insights from the simulation results with the Apriori algorithm and exploratory data analysis techniques too. 2.1 Model Description The mathematical model comprises four pairs of master equations. Each pair is based on the central dogma of molecular biology, which states that genetic information flows from a cell’s genome to its transcriptome to its proteome [ 46 ]. The first step of information transfer is called transcription, while the second is called translation. The master equations are presented and explained in subsection S2.1 in the supplementary file. In this paper, gene and mRNA names are italicised by default, but protein names are not. An exception is that a variable or parameter is always italicised. The first four equations below govern transcription and mRNA degradation, which can be derived from equation S1 . The four equations after them govern translation and protein degradation, which can be derived from equation S2 . Transcription, mRNA degradation, translation, and protein degradation are the four fundamental processes in this model. In the following subsections, we will provide the biological basis of the model. Briefly, the equations describe how environmental stressors influence the molecular network by activating or suppressing the protein complexation events illustrated in figure 1 , the protein-DNA interactions illustrated in figure 2 , and the protein-mRNA interactions illustrated in figure 3 , as well as the downstream effects on the four fundamental processes. For example, oncogenic stresses activate c-Myc to inhibit ARF degradation [ 47 – 49 ] (edge 1.i, figure 1 ). In this model, this stressor is denoted by s 6 . The detailed modelling procedures are explained in subsections S2.2, S2.3, S2.4, and S2.5 in the supplementary file. For an exposition of modelling gene regulatory networks, chapter 7 of [ 50 ] is recommended. Tables 1 and 2 summarise the model’s parameters, including their values, units, and biological descriptions. For a detailed account justifying the chosen values, please refer to subsection S2.6 in the supplementary file. Download figure Open in new tab Fig. 1: Protein complexation events. (a) MYCN protein’s ability to transcribe p53 and MDM2 is downregulated ( equations 2 and 3 ) when it forms complex A with ARF. (b) MDM2 protein is degraded by itself in Complex E and by ARF in complex B and complex F ( equation 7 ). The stressor s 5 promotes complex E formation. (c) p53 protein’s ability to transcribe p53 and MDM2 is downregulated ( equations 2 and 3 ) by MDM2 in complex D and complex F. It is also degraded by MDM2 in complex C ( equation 6 ). The stressors s 3 and s 4 prevent complex C and complex D from forming respectively. (d) The stressor s 6 reduces ARF protein’s degradation rate ( equation 8 ). Download figure Open in new tab Fig. 2: Interactions between proteins and DNA promoter regions. (a) The promoter region of MYCN is not regulated by these four proteins ( equation 1 ). (b) The promoter region of MDM2 is activated (more transcribed) by MYCN and p53 ( equation 3 ). (c) The promoter region of p53 is activated (more transcribed) by MYCN and p53 ( equation 2 ). (d) The promoter region of ARF is inhibited (less transcribed) when occupied by p53 and activated (more transcribed) when occupied by the stressor s 1 ( equation 4 ). Download figure Open in new tab Fig. 3: Interactions between proteins and mRNA transcripts. (a) MDM2 binds to and reduces the degradation rate of MYCN transcripts ( equation 1 ). (b) The transcripts of MDM2 are not regulated by these four proteins ( equation 3 ). (c) The rate at which p53 transcripts are translated is dependent on MDM2, which raises the rate when the stressor s 2 is on and lowers it when s 2 is off ( equation 6 ). (d) The transcripts of ARF are not regulated by these four proteins ( equation 4 ). View this table: View inline View popup Download powerpoint Table 1: Model parameters controlling the four fundamental processes in our model. DR denotes ‘degradation rate’. A denotes ‘assumed’. View this table: View inline View popup Download powerpoint Table 2: Equilibrium dissociation constants (DC) and Hill coefficient. P denotes ‘parameter’. S denotes ‘source’. A denotes ‘assumed’. 2.1.1 Environmental stressors A NB cell responds to a multitude of environmental stressors, which elicit intracellular molecular changes. For example, they may inhibit protein complexation events selectively. They can be thought of as the model’s inputs, but for the sake of simplicity, they are classified as parameters in this paper. There are six stressors, which are denoted by the parameters s 1 , s 2 , s 3 , s 4 , s 5 , and s 6 . Their effects on protein-protein, protein-DNA, and protein-mRNA interactions are illustrated in figures 1 , 2 , and 3 , respectively. They are dimensionless and lie on the interval [0, 1]. As a stressor increases from zero to one, the cell transitions from a stress-free environment to a maximum stress environment. Here are the intracellular changes induced by them and the molecular mechanisms that allow them to do so. s 1 is linked to an increased ARF transcription rate (edge 2.f, figure 2 ). This parameter represents oncogenic stresses that activate ARF transcription through E2F-1 [ 60 – 62 ] and Ras [ 63 , 64 ], as well as oxidative and heat shock stresses [ 65 ]. s 2 is linked to an increased p53 translation rate through MDM2 through two mechanisms (collectively represented by edge 3.c, figure 3 ). p53 translation is inhibited or induced by MDM2 depending on the cellular context. MDM2 naturally inhibits the translation activator L26 (scenario B, figure 4 or edge 3.d, figure 3 ), which would otherwise activate p53 translation (scenario A, figure 4 ). In the first mechanism (scenario D, figure 4 ), stresses in general prevent this negative action of MDM2 on p53 translation [ 66 – 68 ]. In the second mechanism (also scenario D, figure 4 ), specific stresses such as DNA damage induce MDM2 to bind to p53 mRNA transcripts. This process activates translation [ 69 , 70 ]. In other words, these specific stresses promote MDM2’s positive action on p53 translation. s 3 is linked to a reduced p53 degradation rate (edge 1.f, figure 1 ). In a stressed cell, both MDM2 and p53 are phosphorylated and acetylated to prevent complexation (complex C in sub section 2.1.2 ), which normally induces p53 degradation [ 71 , 72 ]. For example, ATM phosphorylates MDM2 in response to DNA damage induced by radiomimetic chemicals and ionising radiation, but not UV radiation [ 73 ]. ATM is a part of s 3 . s 4 prevents MDM2 from forming a complex with p53 (complex D and complex F in sub section 2.1.2 ) through phosphorylation and acetylation [ 72 ] (edge 1.g, figure 1 ). In the complexes, p53’s transactivation function is blocked, so s 4 enhances this function [ 72 ]. s 5 is linked to an increased MDM2 degradation rate (edge 1.h, figure 1 ). Cellular stresses signal through HAUSP to MDM2 [ 74 ], making MDM2 switch functionally from degrading p53 to dimerising with and degrading itself (complex E in sub section 2.1.2 ) [ 75 ]. s 6 is linked to a reduced ARF degradation rate (edge 1.i, figure 1 ). Oncogenic stresses activate c-Myc , which inhibits ARF’s ULF-dependent degradation pathway [ 47 – 49 ]. Download figure Open in new tab Fig. 4: p53 translation depends on MDM2 and whether the cell is stressed. (A) Naturally, L26 (a ribosomal protein) binds to p53 mRNA to promote translation. (B) In an unstressed cell, MDM2 binds to L26 to prevent it from promoting p53 translation. (C) In a stressed cell without MDM2, L26 acts as it does in scenario A to promote translation. (D) In a stressed cell with MDM2, MDM2 does not interact with L26, which can now promote p53 translation. MDM2 can also bind to p53 mRNA to promote translation directly. Note that this figure serves to illustrate the biological mechanisms, but in the mathematical model, the functions of L26 are treated in an abstract manner. In summary, s 1 targets the promoter region of ARF ; s 2 targets p53 mRNA transcripts; s 3 , s 4 , and s 5 target MDM2 and p53 to influence complexation; and s 6 targets ARF. 2.1.2 Protein complexation Complexation occurs when two or more proteins bind together to form a larger molecular structure called a complex. Complexation introduces non-linearity to the system because a binding event makes the concentrations of two or more proteins codependent and the resulting complex may exert further effects on the gene regulatory network by modulating the four fundamental processes. In the list below, we define six protein-protein complexes associated with the gene regulatory network. They are graphically depicted in figure 1 . A short description of the binding event(s) is provided for each complex. Note that it is possible for a protein to have multiple binding terminals. For example, complex C and complex D (edges 1.c and 1.b respectively, figure 1 ) share the same constituents (MDM2 and p53), but they bind at two different terminals in the two complexes. Unsurprisingly, they influence the gene regulatory network differently. Complex A: MYCN and ARF. In this complex, MYCN is bound to ARF with an equilibrium dissociation constant denoted by K A . ARF inhibits the transactivation function of MYCN [ 76 ], as represented by edge 1.d in figure 1 . Complex B: MDM2 and ARF. In this complex, ARF (either N-terminal or C-terminal) is bound to MDM2 (C-terminal) [ 77 , 78 ] with an equilibrium dissociation constant denoted by K B . ARF induces MDM2 degradation, as represented by edge 1.e in figure 1 . Complex C: MDM2 and p53. MDM2 (C-terminal RING finger domain) is bound to p53 (N-terminal) [ 75 , 79 , 80 ] with an equilibrium dissociation constant denoted by K C . MDM2 induces p53 degradation [ 81 ], as represented by edge 1.c in figure 1 . Complex D: MDM2 and p53. MDM2 (N-terminal) is bound to p53 (N-terminal transactivation domain) [ 72 , 82 – 84 ] with an equilibrium dissociation constant denoted by K D . This interaction switches off p53’s transactivation function, as represented by edge 1.b in figure 1 . Although p53 is thought to activate and inhibit transcription differently [ 85 , 86 ], it is assumed in our model that only unbound p53 can get close enough to inhibit transcription. Therefore, this complexation event prevents p53 from influencing transcription in all manners. Complex E: MDM2 and MDM2. Two copies of MDM2 form this complex with an equilibrium dissociation constant denoted by K E . In this dimer, MDM2 undergoes self-degradation [ 75 ], as represented by edge 1.a in figure 1 . Complex F: ARF, MDM2, and p53. ARF, MDM2 and p53 form a ternary complex. It is assumed that MDM2 is the bridge binding ARF with its C-terminal (as in complex B) and p53 with its N-terminal (as in complex D) [ 78 , 87 ]. Based on this assumption, the equilibrium dissociation constant is K F = K B K D . It is also assumed that, as in complex D, this complexation event prevents p53 from influencing transcription in all manners (transactivation and inhibition). How complex F is a combination of complexes B and D is illustrated in inset 1.j in figure 1 . The above complexation events are represented as quasi-equilibrium processes in the model, based on the assumption that protein-protein interactions occur on a much faster timescale than transcription, translation, and degradation. The mathematical details are provided in subsection S2.2 in the supplementary file. As explained in subsection S2.3 in the supplementary file, these interactions are modelled similarly to protein complexation, but the equilibrium dissociation constants ( table 2 ) are labelled numerically rather than alphabetically. 2.1.3 Interactions between proteins and DNA promoter regions In addition to protein complexation events, other binding events are non-linear phenomena too. Some of the four proteins can feedback to interact with the genome (DNA) by binding to certain promoter regions to control the relevant transcription rates. Figure 2 summarises how each protein activates or inhibits the promoters of the genes encoding the other three proteins. The promoter region of MYCN is assumed to be unbound since there is no evidence to the contrary. As a result, there is only one transcription rate for this gene. The p53 gene is transcribed by two transcriptional activators: MYCN [ 25 ] (edge 2.b, figure 2 ) and p53 itself [ 88 ] (edge 2.e, figure 2 ). They bind to non-overlapping sites in the gene’s promoter region [ 89 ]. As stated in sub section 2.1.2 , MYCN’s transcriptional activity is inhibited by ARF in complex A (edge 1.d, figure 1 ) and p53’s transcriptional activity is inhibited by MDM2 in complexes D (edge 1.b, figure 1 ) and F (inset 1.j, figure 1 ). In complex C (edge 1.c, figure 1 ), MDM2 induces p53 degradation through the process of ubiquitination. It is fair to assume that p53 cannot transcribe a gene whilst degrading, hence p53 can bind to promoter regions only when it is unbound. Taken together, the p53 gene is transcribed at a higher rate when its promoter region is occupied by unbound MYCN and/or unbound p53. The MDM2 gene is transcribed by two transcriptional activators: MYCN [ 90 ] and p53 [ 91 ] (edges 2.a and 2.c respectively, figure 2 ). Its promoter region has the same binding pattern as that of p53 , meaning the MDM2 gene is transcribed at a higher rate when its promoter region is occupied by unbound MYCN and/or unbound p53. However, the binding affinities differ between the two promoter regions. The ARF gene is down-regulated by p53 [ 49 , 92 , 93 ] and upregulated by a group of proteins activated by the environmental stressor s 1 , such as E2F-1, β -catenin, and Ras [ 60 , 64 , 65 ], as represented by edges 2.d and 2.f in figure 2 respectively. The gene is transcribed at a natural rate when its promoter region is unoccupied, at a higher rate when it is activated by s 1 , and the transcription rate drops to zero when p53 occupies the promoter region. 2.1.4 Interactions between proteins and mRNA transcripts Figure 3 graphically explains how MDM2 regulates the translation processes producing MYCN and p53. When MDM2 binds to MYCN mRNA transcripts, it reduces their degradation rate [ 94 ] (edge 3.a, figure 3 ). It is assumed that when MYCN mRNA is bound to MDM2, ribosomal subunits and initiation factors can still bind to it to initiate translation [ 95 , 96 ]. There is no evidence that the other three mRNA transcripts are protected from or encouraged towards degradation by proteins. On the other hand, MDM2 affects the rate at which p53 mRNA transcripts are translated (edge 3.b, figure 3 ). In an unstressed cell (edge 3.d, figure 3 ), p53 translation is indirectly inhibited because MDM2 is free to inhibit L26, which activates p53 translation [ 66 – 68 ]. In response to s 2 (associated with damaged DNA), L26 is free to promote p53 translation, which is also directly activated by MDM2 [ 69 , 70 ]; these positive effects are collectively represented by edge 3.c in figure 3 . The technical details regarding how these interactions are modelled mathematically, including how the effects of L26 are represented abstractly since its concentration is not a state variable, are provided in subsection S2.4 in the supplementary file. 2.1.5 Functional effects of complexation Our assumptions relating to the effects of complexation on transcription are as follows. MYCN is transcribed at one rate only. p53 is transcribed at two rates. The natural rate applies when its promoter region is unoccupied. The enhanced rate is achieved when uncomplexed MYCN (edge 2.b, figure 2 ) and/or uncomplexed p53 (edge 2.e, figure 2 ) occupies the region. MDM2 gene is transcribed at two rates. The natural rate applies when its promoter region is unoccupied. The enhanced rate is achieved when uncomplexed MYCN (edge 2.a, figure 2 ) and/or uncomplexed p53 (edge 2.c, figure 2 ) occupies the region. ARF is transcribed at three rates. When uncomplexed p53 occupies its promoter region, it is zero (edge 2.d, figure 2 ). When the promoter is unoccupied, transcription proceeds at the natural rate. When the gene is activated by s 1 , an enhanced transcription rate is achieved (edge 2.f, figure 2 ). Our assumptions relating to the effects of complexation on mRNA degradation are as follows. MYCN mRNA transcripts exist in two states, which degrade at a natural rate and a reduced rate respectively. In the uncomplexed form, they degrade naturally. When a transcript is bound by MDM2, it degrades at the reduced rate (edge 3.a, figure 3 ). The other three mRNA transcripts have a fixed degradation rate. Our assumptions relating to the effects of complexation on translation are as follows. Although MYCN mRNA exists in two forms (unbound or bound by MDM2), it is assumed that ribosomal subunits and initiation factors can bind to it in both cases to initiate translation [ 95 , 96 ]. Therefore, it has just one translation rate. p53 protein is produced at three rates: natural, enhanced, and reduced. When p53 mRNA is unbound, translation occurs at the natural rate. When MDM2 binds to it in response to s 2 (associated with damaged DNA), the enhanced rate applies (edge 3.c, figure 3 ). In an unstressed cell, MDM2’s action leads to a reduced rate of p53 production (edge 3.d, figure 3 ). MDM2 and ARF are produced at the same fixed translation rate. Finally, our assumptions relating to the effects of complexation on protein degradation are as follows. MYCN exists either in its unbound state or in complex A (bound to ARF), as represented by edge 1.d in figure 1 . The same natural degradation rate applies in both cases. p53 exists in four states: unbound, complex C (bound to MDM2), as represented by edge 1.c in figure 1 , complex D (bound to MDM2), as represented by edge 1.b in figure 1 , and complex F (bound to ARF and MDM2), as illustrated in inset 1.j in figure 1 . It degrades at an enhanced rate in complex C because MDM2 induces p53 degradation. In the other three states, it degrades naturally. MDM2 exists in six states: unbound, complex B (bound to ARF), as represented by edge 1.e in figure 1 , complex C (bound to p53), as represented by edge 1.c in figure 1 , complex D (bound to p53), as represented by edge 1.b in figure 1 , complex E (two copies of MDM2), as represented by edge 1.a in figure 1 , and complex F (bound to ARF and p53), as illustrated in inset 1.j in figure 1 . In complexes B and E, it degrades at an enhanced rate. In the other states, it degrades naturally. ARF exists in four states: unbound, complex A (bound to MYCN), as represented by edge 1.d in figure 1 , complex B (bound to MDM2), as represented by edge 1.e in figure 1 , and complex F (bound to MDM2 and p53), as illustrated in inset 1.j in figure 1 . In all four cases, ARF degrades naturally. How these functional effects are represented mathematically in the model is explained in subsection S2.5 in the supplementary file. 2.2 Numerical simulations We conducted three sets of in silico experiments. The configurations and computational details of these experiments can be found in subsection S2.7 in the supplementary file. First, we varied three groups of parameters with and without MYCN amplification in a grid search, meaning each parameter was set at predefined values within a range. The results were used to confirm the model’s stable long-term behaviour (equilibrium), which represents the clinical/treatment outcome; understand the outcome’s sensitivity to the three groups of parameters; and establish a NB classification system based on our hypothesis. Second, we focused on the two most sensitive groups of parameters, repeating the grid search without MYCN amplification and at five levels of MYCN amplification. Third, we managed to reach the conclusions of the first two experiments by exploring the parametric space randomly. 2.3 NB classification system based on our hypothesis In the first experiment, we tested 1259715 parametric combinations by using each to perform two sets of 10 simulations, once with g MYCN = 1 and once with g MYCN = 2. In each case, all 20 simulations ended in the same equilibrium state (eight-dimensional fixed point). Furthermore, the 2 · 1259715 = 2519430 fixed points were all found to be stable by perturbing the state vector from and checking that it returned to every fixed point. Then, we applied our hypothesis to the same dataset to create a NB classification system for the purpose of prognosing NB patients. 2.3.1 p53/MYCN: a new prognostic metric As stated in the introduction, our hypothesis is that the p53 protein level in a NB tumour, relative to the MYCN protein level therein, decides whether it is likely to progress (poor outcome) or regress (good outcome). This hypothesis rests on the fact that p53 and MYCN are, respectively, commonly considered to be a tumour suppressor gene and an oncogene with opposite effects. We applied the hypothesis to two simulations conducted in the first experiment. They were conducted with the default parameters and inactive stressors ( s i = 0), one with g MYCN = 1 and the other with g MYCN = 2. The arithmetic mean of the two equilibrium values of p53/MYCN was calculated and used to define a critical threshold. A patient or simulated equilibrium whose p53/MYCN is above the critical threshold is predicted to have a good outcome. Otherwise, a poor outcome is predicted. 2.3.2 MYCN mRNA and MYCN protein levels We wanted to check if the classification system and by extension our hypothesis is consistent with clinical datasets, especially [ 24 ]. Their dataset includes MYCN mRNA and MYCN protein levels, presented as a Boolean variable (high or low) and a discrete variable with four possible values respectively. In order to compare our simulation results with their dataset, we set up comparable quantities within our classification system. Working with the results of the first grid search, we noticed that the WT equilibria ( g MYCN = 1) are always associated with less MYCN mRNA than their corresponding MA equilibria ( g MYCN = 2). The mean value of the maximum and minimum MYCN mRNA levels in the combined set of WT and MA equilibria was used to define a critical threshold. An equilibrium above it was classified as having a high MYCN mRNA level and vice versa , in alignment with [ 24 ]. Working with the same set of simulation results, we merged the two sets of equilibrium MYCN protein levels ( g MYCN = 1 and g MYCN = 2) and set up three thresholds to divide them into four intervals. They were selected to keep the 1259715 WT equilibria ( g MYCN = 1) within the first three intervals and the 1259715 MA ( g MYCN = 2) ones within the last three intervals, again in alignment with [ 24 ]. This binary classification system predicts a good outcome if p53/MYCN≥3.5 and a poor outcome otherwise. It also describes a tumour along three dimensions: its MYCN amplification status, MYCN mRNA level, and MYCN protein level. Along these dimensions, there are two (WT or MA), two (high or low), and four possible discrete values respectively: 16 categories in total. Our classification system reduces three dimensions to one, but the original three dimensions enable comparison between different datasets. We created this classification system partly by aligning the results of the first grid search with [ 24 ]. Although the thresholds (except the one for p53/MYCN) were manually imposed to achieve alignment, the mere existence of suitable thresholds is a plus point for our hypothesis. 2.4 Data mining and visualisation To find evidence that our hypothesis is potentially correct in more general scenarios, we interpreted the simulation results of the three experiments within the classification system established on the basis of our hypothesis. Three-dimensional scatter plots were created to visualise the results in terms of the MYCN mRNA levels (pM), MYCN protein levels (pM), and p53/MYCN metric levels (dimensionless) at the simulated equilibria. The equilibria found in the second experiment were divided into four subsets along two dimensions. In one dimension, they were classified as WT or MA according to g MYCN . In the other dimension, they were classified as clinically good or poor according to the p53/MYCN metric. In each subset, we identified frequent itemsets of the tested parametric values and association rules between them with the Apriori algorithm [ 97 ]. The algorithm was implemented with the arules R package [ 98 ]. The technical details of this analysis are provided in subsection S2.8 in the supplementary file. The equilibria found in the third experiment were also divided into the same four subsets. Although we could have applied the Apriori algorithm again after dividing each parameter’s range into discrete intervals, we decided to visualise this dataset directly by generating augmented scatterplot matrices. They illustrate how the equilibria in each subset are continuously distributed in the parametric space. How we used the ggpairs function within the GGally R package for this purpose is explained in subsection S2.8 in the supplementary file. 3 Results This section begins by confirming the model’s stable long-term behaviour, presenting the results of a local sensitivity analysis (first experiment), and establishing the NB classification system. Then, the results of a global sensitivity analysis (second experiment) are presented within this classification system. After that, a more technical subsection delves deeper into the results of the second experiment by presenting the frequent itemsets and association rules extracted by the Apriori algorithm. The last subsection presents the results we produced in the third experiment. In all figures, MYCN m , MYCN , and p53/MYCN indicate the MYCN mRNA level (pM), MYCN protein level (pM), and the p53/MYCN metric (dimensionless), respectively. A blue data point represents a WT equilibrium, whilst a red data point represents a MA equilibrium. 3.1 Stability analysis, local sensitivity analysis, and binary NB classification system In the first experiment, we tested over 2.5 million parametric combinations, including two MYCN amplification statuses ( g MYCN = 1 and g MYCN = 2). In each case, the state vector equilibrated to the same eight-dimensional fixed point from 10 different initial conditions. Furthermore, all of them were found to be stable. The results suggest that periodic, quasi-periodic, and chaotic dynamics are not biologically relevant in this context even if they can in theory be reached from certain basins of attraction when given certain parametric combinations. Therefore, our model complies with the requirement that oncogenic mutations abolish p53 oscillations [ 32 ]. A subset of the same dataset is presented in figure 5 . Recall that there are three groups of parameters: the six stressors, h 5 , and the transcription and translation rates. Each panel presents the fixed points achieved by varying one group while keeping the other two at the default values, with or without MYCN amplification. The left panels indicate that the stressors principally affect the p53 protein level. The middle panels contain results that, relative to the results in the left and right panels, suggest that the model is not sensitive to h 5 . The fixed points in the right panels demonstrate large variations relative to the variations in the other four panels. Download figure Open in new tab Fig. 5: 3D visualisation of the fixed points reached in the local sensitivity analysis. They are a subset of the results of the first experiment. The left panels present the fixed points obtained by switching on/off the stressors. The middle panels present the fixed points obtained by varying h 5 with s 1 = 0.5, the other stressors switched off, and the rates set to the default values. Recall that s 1 was not set to zero because it is relevant to h 5 . The right panels present the fixed points obtained by varying the transcription and translation rates. In total, 64 + 3 + 6561 = 6628 fixed points are presented for each of g MYCN = 1 and g MYCN = 2. A blue data point represents a WT equilibrium, whilst a red data point represents a MA equilibrium. Both concentrations are in pM, but the p53/MYCN metric is dimensionless. Returning to the full dataset, we applied our hypothesis by setting thresholds in three dimensions. As a reminder, our hypothesis is that the p53 protein level in a NB tumour, relative to the MYCN protein level therein, decides whether it is likely to progress or regress. This metric is called p53/MYCN and a NB with low p53/MYCN would have a poor clinical outcome. A corollary of the hypothesis is that p53/MYCN should follow a concave function of MYCN ’s mRNA level in a mixed collection of WT and MA NBs. First, we focused on the two simulations with the default parameters and inactive stressors, one with g MYCN = 1 (WT NB) and one with g MYCN = 2 (MA NB). Using the p53 and MYCN protein levels at the two simulated equilibria, we calculated the p53/MYCN metric’s values there: approximately 4.526 for WT NB and approximately 2.483 for MA NB. Due to the fact that MYCN amplification is typically associated with poor outcomes, we calculated their arithmetic mean (approximately 3.5) and made this the critical threshold separating good outcomes (p53/MYCN≥ 3.5) from poor outcomes (p53/MYCN < 3.5). Second, a critical MYCN mRNA level set anywhere on [2.84, 2.98] would classify every MA fixed point in figure 6 as having a high MYCN mRNA level and every WT fixed point as having a low level, in agreement with [ 24 ]. The midpoint was chosen. Download figure Open in new tab Fig. 6: 2D visualisation of the fixed points found in the first experiment. The x-axis reports MYCN mRNA levels (pM). The y-axis reports MYCN protein levels (pM). In this experiment, the stressors, h 5 , and the transcription and translation rates were varied together to generate 64 · 3 · 6561 + 3 = 1259715 parametric combinations. Each combination was used to perform two simulations, once with g MYCN = 1 and once with g MYCN = 2. A blue data point represents a WT ( g MYCN = 1) equilibrium, whilst a red data point represents a MA ( g MYCN = 2) equilibrium. Third, the level of MYCN protein in [ 24 ] is a discrete variable with four possible values (zero, one, two, and three). Their WT levels are all smaller than or equal to level two, while their MA levels are all above level one. Like the level of MYCN mRNA, there are many solutions to this alignment problem, but we decided to respect the additional constraint of making the four intervals as evenly sized as possible. The following thresholds keep the simulated WT MYCN protein levels smaller than or equal to level two and the simulated MA levels greater than or equal to level one, with 88.89% of them being greater than or equal to level two. Level zero: from zero to 1500 Level one: from 1500 to 2500. Level two: from 2500 to 4000. Level three: above 4000. 3.2 Correct predictions based on our hypothesis Due to the results in the middle panel of figure 5 , we decided to set h 5 = 1 throughout the second experiment, but vary s i and the rates together. As a reminder, the range of p53/MYCN values in the middle panels is narrow relative to the ranges in the left and right panels in figure 5 . Since six values of g MYCN (one, two, three, four, six, and eight) were used to mimic MYCN amplification, we ran 2519424 (64 · 6561 · 6) simulations. Figure 7 presents the resulting fixed points within our classification system. The p53/MYCN metric predicts their clinical/treatment outcomes. Those above the critical threshold 3.5 are expected to have good clinical/treatment outcomes, whilst those below it are expected to have poor clinical/treatment outcomes. Download figure Open in new tab Fig. 7: 3D visualisation of the fixed points found in the second experiment, presented within our NB classification system. The stressors and rate parameters of the default configuration were varied together to produce 64 · 6561 = 419904 parametric combinations. Each combination was used to perform six simulations, with g MYCN set to one, two, three, four, six, and eight respectively. In the left panel, only the WT ( g MYCN = 1, blue) and MA ( g MYCN = 2, red) fixed points are presented. In the middle and right panels, the same fixed points are presented alongside the fixed points with even more copies of MYCN (three, four, six, and eight). Both concentrations are in pM, but the p53/MYCN metric is dimensionless. Based on the p53/MYCN metric, 75.66 % of the WT fixed points ( g MYCN = 1) belong to the category with good clinical/treatment outcomes. The remaining 24.34 % belong to the category with poor outcomes. The simulation results, when interpreted within our classification system and through the lens of our hypothesis, are in remarkably good agreement with [ 24 ]. Out of the 23 WT patients they studied, 18 achieved good outcomes (78.26 %) and five achieved poor outcomes (21.74 %). Within our classification system, 32.03 % of the MA fixed points with g MYCN = 2 belong to the category with good outcomes. The remaining 67.97 % belong to the category with poor outcomes. Our results are aligned with the general clinical observation that an elevated MYCN copy number correlates with a progressively heightened prevalence of adverse clinical and biological characteristics [ 99 ]. These simulation results are also not too far from the specific measurements taken by [ 23 , 24 ]. In Figure 2A in [ 23 ], although almost every Trk mRNA level reported for a MA specimen is close to zero (poor prognosis), three of the 15 reported levels are significantly higher, so 20% of the reported MA specimens achieved good outcomes. Figure 3C in [ 24 ] indicates that the six MA patients they studied all achieved poor outcomes. Figure 7 also shows the effects of increasing g MYCN beyond two on p53/MYCN. The poor outcome percentage continues to increase with the copy number. For instance, with g MYCN = 3, 83.64% of the tested parametric combinations led to fixed points classified as poor outcomes. This percentage increases to 93.32%, 97,61%, and 98.78% with copy numbers 4, 6, and 8, respectively. Along the other two dimensions ( MYCN mRNA and MYCN protein levels), the simulation results and [ 24 ] are aligned within our classification system too. The simulated WT fixed points all have low MYCN mRNA levels, similar to [ 24 ]. The simulated MA fixed points with poor outcomes (67.97 %) all have high MYCN mRNA levels, just like Figure 3 in [ 24 ]. Table 3 compares the same pair of datasets along the third dimension, the protein level of MYCN. The trends are similar. View this table: View inline View popup Download powerpoint Table 3: Comparison of the simulation results (Sim) of the second experiment and [ 24 ] (Exp) along one dimension of our classification system: the level of MYCN protein. The simulation results are all fixed points, whilst the experimental results describe tumours. In this table, both datasets are divided into three categories based on the MYCN amplification status and outcome of a case. WT denotes a WT category ( g MYCN = 1) and MA denotes a MA category ( g MYCN = 2). The outcome is/was either good or poor. For example, 14.57 % of the WT fixed points with good clinical/treatment outcomes have MYCN protein levels that belong to level 0. Since we set the thresholds for the MYCN mRNA and MYCN protein levels by comparing the results of our first experiment with [ 24 ], the fact that our model agrees with [ 24 ] is not a genuine emergent phenomenon. However, our classification system’s ability to prognose NB cases based on the p53/MYCN metric is an emergent property. The critical threshold of p53/MYCN was defined independently of [ 24 ] and at no point were the patients’ clinical outcomes in [ 24 ] used in our procedures. 3.3 MYCN amplification produces diminishing returns in terms of p53/MYCN In the second experiment, we increased g MYCN beyond two. The results are presented in the middle and right panels of figure 7 . There is a clear trend that a higher copy number of MYCN is associated with lower and less variable values of p53/MYCN (consistently poor outcomes). The mean and standard deviation (SD) of this metric both decrease as g MYCN increases in figure 7 . Numerically, the mean±SD (pM) is 6.04±3.63 with g MYCN = 1, 3.19±1.95 with g MYCN = 2, 2.24±1.40 with g MYCN = 3, 1.76±1.13 with g MYCN = 4, 1.27±0.85 with g MYCN = 6, and 1.02±0.71 with g MYCN = 8. This decreasing trend is consistent with the clinically observed trend that in general, as the copy number of MYCN goes up in a tumour, its clinical and biological characteristics get progressively more adverse [ 99 ]. More specifically, the WT tumours studied by [ 23 ] displayed a much broader range of TrkA expression levels than their MA tumours, most of whose Trk mRNA levels were found to be close to zero. As the level of TrkA mRNA is a favourable prognostic indicator, just like the p53/MYCN metric (but mechanistically less relevant), this agreement exemplifies our classification system’s predictive ability. By extension, it is favourable for our hypothesis. It is noteworthy that the p53/MYCN metric’s relationship with g MYCN is non-linear in the middle and right panels of figure 7 , demonstrating diminishing returns. As the copy number increases from one to two, three, four, six and eight, the marginal loss in p53/MYCN decreases progressively. This agrees with a corollary of our hypothesis: the mRNA level of MYCN should have a non-linear relationship with p53 and MYCN levels. Compared to the p53/MYCN metric, the levels of MYCN mRNA and MYCN protein vary with g MYCN in the opposite direction. The four means and standard deviations all increase with g MYCN . Consistently, [ 23 ] observed that their MA tumours almost uniformly expressed more MYCN mRNA than their WT tumours. However, in our simulations, the general trend (the mean) in the relationship between the p53/MYCN metric and the mRNA level of MYCN is a monotonically decreasing function. Also, a low copy number of MYCN is associated with a broad range of p53/MYCN values. Our results disagree with the concave function suggested by figure 2A in [ 23 ]. 3.4 Association rules extracted from the results of the second experiment We attempted to gain mechanistic insights into the results obtained in the second experiment. To do so, we divided the fixed points into four subsets corresponding to four genotype-outcome combinations: WT ( g MYCN = 1) or MA ( g MYCN = 2) and good (p53/MYCN above 3.5) or poor (p53/MYCN below 3.5) clinical/treatment outcomes. We applied the Apriori algorithm to the parametric combinations that gave us the fixed points in each subset. Altogether, we extracted 80 association rules, 20 rules for each genotype-outcome combination. For the sake of clarity and conciseness, table 4 reports four meta association rules, which summarise the association rules. The four sets of 20 association rules are presented in tables S1, S2, S3, and S4 in the supplementary file. It is important to bear in mind that the association rules are robust but not very common. For example, in table S1 , the top association rule is present in 1.1 % of the parametric combinations only (level of support). However, it is robust because of its high levels of confidence and lift (1 and 2.39), which clearly indicate that the antecedent (LHS) implies the consequent (RHS). View this table: View inline View popup Download powerpoint Table 4: Four meta association rules summarising the association rules extracted from four subsets of parametric combinations (second experiment). They describe wild-type and MYCN -amplified fixed points with good and poor outcomes. Each row presents the meta association rule summarising the dominant trait in 20 association rules. ‘/’ denotes OR. For example, H/High MYCN denotes H MYCN or High MYCN . Each rule follows the form: LHS ⇒ RHS, meaning that the left-hand side implies the right-hand side. L X , M X , H X denote the transcription rate of X is at 70%, 100%, 130% of the default value, respectively. Low X , Mid X , High X indicate the same for its translation rate. In the WT cases with good outcomes ( table S1 ), the meta association rule is that a low MYCN transcription or translation rate is required when the tumour already exhibits a high MYCN translation or transcription rate and a low p53 rate. Two rules indicate that when MYCN ’s transcription and translation rates are both high and one of p53 ’s biosynthetic rate is at the default value, the other biosynthetic rate should be high (level of confidence is 0.97). In the WT cases with poor outcomes ( table S2 ), the meta association rule is that a low p53 transcription or translation rate is required when the other p53 biosynthetic rate is already high and the stressor s 3 is active. Recall that the latter ( s 3 ) prevents MDM2 and p53 from forming complex C, where p53 degradation is enhanced. ATM is one protein that can perform this function in response to DNA damage (stress). The MA fixed points with good outcomes ( table S3 ) have very different association rules in their parametric combinations. 18 of the top 20 association rules pair a high MYCN biosynthetic rate and a default p53 biosynthetic rate with another high p53 biosynthetic rate. The MA fixed points with poor outcomes ( table S4 ) have less surprising association rules in their parametric combinations. 18 of the top 20 association rules pair a high p53 transcription or translation rate and a low MYCN biosynthetic rate with a low p53 translation or transcription rate. The antecedent also has one of the following: the other MYCN biosynthetic rate at the default value, an active s 3 , an active s 6 , an inactive s 4 , and a high ARF translation rate. 3.5 Random exploration of the parametric space To test the robustness of the conclusions drawn from the first two experiments, the parametric space was explored randomly in the third experiment. 10000 random parametric combinations were created and each was used to run six simulations with different values of g MYCN . The resulting fixed points are shown in figure 8 . Download figure Open in new tab Fig. 8: 3D visualisation of the fixed points reached in the simulations in the third experiment. The stressors and rate parameters of the default configuration were varied together and randomly, whilst keeping h 5 = 1, leading to 10000 parametric combinations. Each of these combinations was generated by sampling each stressor and rate randomly and uniformly between the minimum and maximum values tested in the second experiment. Each combination was used to perform six simulations, with g MYCN set to one, two, three, four, six, and eight respectively. In the left panel, only the WT ( g MYCN = 1, blue) and MA ( g MYCN = 2, red) fixed points are presented. In the middle and right panels, the same fixed points are presented alongside the fixed points with even more copies of MYCN (three, four, six, and eight). Both concentrations are in pM, but the p53/MYCN metric is dimensionless. In figure 8 , based on the p53/MYCN metric, 83.17% of the 10000 WT fixed points ( g MYCN = 1) have good outcomes and the remaining 16.83% have poor outcomes. On the other hand, only 24.17% of the 10000 MA fixed points ( g MYCN = 2) have good outcomes and the remaining 75.83% have poor outcomes. [ 24 ] studied 23 WT patients, 18 of whom had good outcomes (78.26%) and five had poor outcomes (21.74%). They only studied six MA patients and they all had poor outcomes (100%). Echoing the comments made in sub section 3.2 , our results agree with the general clinical observation that an elevated MYCN copy number correlates with a poor clinical outcome [ 99 ]. The WT fixed points all have low levels of MYCN mRNA, just like [ 24 ]. The MA fixed points all have high levels of MYCN mRNA, just like [ 24 ]. As can be seen in table 5 , the simulated MYCN protein levels mostly agree with [ 24 ] too. View this table: View inline View popup Download powerpoint Table 5: Comparison of the simulation results (Sim) of the third experiment and [ 24 ] (Exp) along one dimension of our classification system: the level of MYCN protein. The simulation results are all fixed points, whilst the experimental results describe tumours. In this table, both datasets are divided into three categories based on the MYCN amplification status and outcome of a case. WT denotes a WT category ( g MYCN = 1) and MA denotes a MA category ( g MYCN = 2). The outcome is/was either good or bad. For example, 10.67 % of the WT fixed points with good clinical/treatment outcomes have MYCN protein levels that belong to level 0. The results of increasing g MYCN beyond two are presented in the middle and right panels of figure 8 . There is a clear trend that a higher copy number of MYCN is associated with lower and less variable values of p53/MYCN (consistently bad outcomes). The mean and SD of this metric both decrease as g MYCN increases in figure 8 : 5.54±2.11 (mean±SD) with g MYCN = 1, 2.88±1.11 with g MYCN = 2, 1.98±0.75 with g MYCN = 3, 1.57±0.58 with g MYCN = 4, 1.14±0.43 with g MYCN = 6, and 0.90±0.34 with g MYCN = 8. Reflecting the decreasing trend of the p53/MYCN metric’s mean, 95.70% of the fixed points with g MYCN = 3 have poor outcomes; g MYCN = 4, 99.28%; g MYCN = 6, 99.94%; and g MYCN = 8, 100%. Compared to the p53/MYCN metric, the levels of MYCN mRNA and MYCN protein vary with g MYCN in the opposite direction. The two means and two standard deviations all increase with g MYCN . These results mirror those presented in figure 7 , so the comments made in sub section 3.3 are true for the results presented in figure 8 too. In agreement with [ 23 ], the WT fixed points display a much broader range of p53/MYCN than the MA fixed points. p53/MYCN displays a non-linear relationship with g MYCN . Aligned with [ 23 ], the levels of MYCN mRNA and MYCN protein become higher when g MYCN increases in figure 8 . Again in conflict with figure 2A in [ 23 ], p53/MYCN is a monotonically decreasing function of the mRNA level of MYCN and p53/MYCN is highly variable when g MYCN is small. 3.6 Visualising parametric distributions in the four scenarios of the third experiment Just as we did before learning association rules from the results of the second experiment, we divided the fixed points found in the third experiment into four subsets (genotype-outcome combinations): WT ( g MYCN = 1) or MA ( g MYCN = 2) and good (p53/MYCN metric above 3.5) or poor (p53/MYCN metric below 3.5) clinical/treatment outcomes. Figure 9 is a set of four ggplot2 generalised pairs plots. They visualise the four subsets of parametric combinations. Each plot was created using the ggpairs function, publicly available in the GGally R package. Download figure Open in new tab Fig. 9: Four ggplot2 generalised pairs plots comparing pairs of 14 parameters (six stressors and eight rates) in four subsets of the parametric combinations tested in the third experiment. The top two plots present the parametric combinations leading to WT fixed points. The bottom two plots present the parametric combinations leading to MA fixed points. The left two plots present the fixed points with good outcomes. The right two plots present the fixed points with poor outcomes. In each plot below the leading diagonal are 91 scatterplots comparing the unique pairs formed by the 14 parameters. On the diagonal are the probability density functions of the 14 parameters in the subset. Above the diagonal are the Pearson correlation coefficients of the 91 pairs of parameters. Each correlation’s statistical significance is shown too. In a matrix element, *** means the reported p-value is smaller than 0.001; ** means that it is smaller than 0.01; * means that it is smaller than 0.05; and a full stop (.) means it is smaller than 0.10. The top-left panel in figure 9 reveals interesting features shared by the parametric combinations giving us WT fixed points with good outcomes. The association rules in table S1 are robust but rare. They indicate the requirements (RHS) for such a fixed point given a set of conditions (LHS), but the conditions are rare. By contrast, this plot pertains to each parameter’s distribution in this sample, including its average value. The two distributions of MYCN rates are skewed towards their minima, whilst both distributions of p53 rates are skewed towards their maxima. The six stressors have relatively uniform distributions in this subset, implying that having specific combinations of stressors is unnecessary for a good outcome. The top-right panel pertains to the WT fixed points with poor outcomes. The distributions of the two MYCN rates are skewed towards their maxima, whilst the distributions of the two p53 rates are skewed towards their minima. Furthermore, the plot indicates that it is impossible for both of MYCN ’s rates to be low or both of p53 ’s rates to be high in a parametric combination. Interestingly, the distribution of s 3 is skewed towards its minimum. The bottom-left panel pertains to the MA fixed points with good outcomes. The distributions of both MYCN rates are skewed towards their minima, whilst the distributions of both p53 rates are skewed towards their maxima. Furthermore, it is impossible for both MYCN rates to be very high or for a low p53 rate to coexist with a high MYCN rate in a parametric combination. Compared to the top-left panel (WT, good outcomes), where the six stressors are almost uniformly distributed, s 3 and s 6 are skewed towards the maximum and minimum respectively. The bottom-right panel presents the parametric combinations giving us MA fixed points with poor outcomes. The distributions of both MYCN rates are skewed towards their maxima, but the distributions of both p53 rates are skewed towards their minima. It is impossible for both MYCN rates to be very low in a parametric combination. 4 Discussion The MYCN enigma is the general observation that a NB patient’s clinical outcome does not depend monotonically on their tumour’s MYCN expression level or MYCN amplification status. This paper presents a hypothesis inspired by MYCN’s ability to activate p53 despite the consensus that MYCN is an oncogene [ 21 ] and p53 is a tumour suppressor gene [ 25 ]. The key concept in this hypothesis is that the p53 protein level in a NB tumour, relative to the MYCN protein level therein, decides whether it progresses or regresses. We built a mathematical model comprising eight non-linear ODEs to model a molecular network connecting MYCN to the ARF/MDM2/p53 axis. The model itself is an addition to the literature on NB for the following reasons. The negative feedback loop formed by MDM2 and p53 has been modelled with ODEs in many studies, mostly with the aim of understanding how genetic damage triggers oscillations in p53’s level [ 32 – 38 ]. Some models even couple p53 dynamics with ARF and other species acting upstream (such as ATM) and downstream (such as WIP1) of p53 [ 39 – 43 ]. However, to the best of our knowledge, the presented mathematical model is the only one that couples MYCN to the ARF/MDM2/p53 axis. After parameterising the model, we tested the model to confirm that its long-term behaviour is stable for over 1.25 million parametric combinations with and without MYCN amplification. Only stable fixed points were found, so our model complies with [ 32 ]. Their theory states that oncogenic mutations abolish p53 oscillations. Then, we explored the parametric space systematically on a grid. We varied the pattern of active stressors (up to six), the Hill coefficient h 5 , and the eight transcription and translation rates. We varied the three groups independently and together at different levels of MYCN amplification, including the wild type without this mutation. Finally, we explored the parametric space randomly. 4.1 To what extent is the hypothesis correct? We applied our hypothesis to the simulation results of the first experiment to establish thresholds along three dimensions: MYCN mRNA’s level, MYCN protein’s level, and the p53/MYCN metric. This turned our hypothesis into a concrete binary classification system. It predicts whether a patient is likely to have a good or poor outcome based on the p53/MYCN metric level in their tumour. It also describes a tumour along three dimensions: its MYCN amplification status (Boolean variable), MYCN mRNA level (Boolean variable), and MYCN protein level (discrete variable with four possible values). These extra dimensions enable comparison between simulated fixed points and real tumours. Figure 5 presents the fixed points found in the first experiment. Compared to the six stressors and eight transcription/translation rates, the Hill coefficient h 5 is not a sensitive parameter. It controls the importance of cooperativity in s 1 ’s positive influence on ARF transcription. s 1 models oncogenic stresses that activate ARF transcription through E2F-1 [ 60 – 62 ] and Ras [ 63 , 64 ], as well as oxidative and heat shock stresses [ 65 ]. Unlike the other species, E2F-1 is a transcription factor, so it physically interacts with DNA and h 5 is most relevant to this aspect of s 1 . Therefore, our simulation results suggest that cooperativity plays a small role in how E2F-1 influences p53/MYCN. Interpreted within this classification system, the results of the second experiment ( figure 7 and table 3 ) are consistent with most of the clinical and experimental observations listed in section S1 in the supplementary file. Figure 7 clearly indicates a negative relationship between the copy number of MYCN and the p53/MYCN metric, a favourable biomarker in our classification system, suggesting a tumour’s clinical and biological characteristics get progressively more adverse with its degree of MYCN amplification. This is consistent with [ 99 ]. The non-linear relationship in figure 7 is a corollary that follows from our hypothesis too. When the MYCN copy number is higher, the marginal change to p53/MYCN associated with an additional copy of MYCN is lower. The p53/MYCN metric correctly predicts positive outcomes for around three quarters of the wild-type fixed points, just like the 23 WT patients [ 24 ] studied. The p53/MYCN metric indicates that most of the MA fixed points have poor outcomes. Three of the 15 MA patients studied by [ 23 ] and all six of the MA patients studied by [ 24 ] achieved poor outcomes. The discrepancies could be explained by two reasons. First, the two datasets [ 23 , 24 ] suffer from small sample sizes. Second, our results pertain to a MYCN copy number of two only and their MA patients might have acquired more copies of MYCN . Indeed, in figure 7 , the poor outcome percentage continues to increase when the copy number increases beyond two. Within our classification system, the simulated WT and MA fixed points all have low and high MYCN mRNA levels respectively, just like their clinical counterparts [ 24 ]. This trend is consistent with the observations of [ 23 ] too. Almost uniformly higher MYCN mRNA levels were measured in their MA specimens than their WT specimens. Table 3 indicates that the simulated protein levels match [ 24 ] very closely. In figure 7 , there is a clear trend that a higher copy number of MYCN is associated with not just lower but also less variable values of p53/MYCN (consistently bad outcomes). Interestingly, figure 7 does not resemble the graph of a concave function like Figure 2A in [ 23 ]. When the level of MYCN mRNA is very low, the p53/MYCN metric displays a broad range of values. To be aligned with [ 23 ], the WT fixed points with inadequate MYCN should all be under the critical threshold of 3.5. Conversely, the WT fixed points with enough MYCN should all be above 3.5. One possibility is that WT NB tumours have a mechanism (refer to sub section 4.3 ) triggered by a low level of MYCN to keep p53/MYCN under 3.5. It follows that the left-most blue fixed points above 3.5 in the middle panel of figure 7 are biologically irrelevant. Another explanation is that [ 23 ] investigated only seven WT cases (small sample size). 4.2 Expanding the hypothesis by going beyond the model Briefly, the results discussed in the last subsection support the conclusion that upon MYCN amplification, it becomes more difficult for a NB tumour to have p53/MYCN over 3.5. The conclusion favours our hypothesis as an explanation for the clinical and experimental results constituting the MYCN enigma. A patient’s clinical outcome depends neither on MYCN amplification nor MYCN expression alone. The level of p53/MYCN is the key. However, we remain ignorant of why it works and our hypothesis is clearly incomplete. It does not explain why the p53/MYCN metric should be kept under 3.5 in most WT tumours with inadequate MYCN and vice versa . This subsection aims to fill these gaps by examining the features shared by each group of parametric combinations that gave us a particular category of fixed points (genotype-outcome combination) in the second and third experiments. As a reminder, these features are the four sets of association rules in section S3 in the supplementary file and the parametric distributions in figure 9 . We will do so by speculating about what a real NB cell does to create the necessary conditions (parametric values). These speculations will take us beyond the scope of the mathematical model. In the WT cases with good outcomes ( table S1 ), the meta association rule implies that competitive dynamics are involved. When MYCN transcription is elevated, it takes resources away from MYCN translation and vice versa . That p53 ’s transcription or translation rate is low in 18 of the top 20 association rules suggests that the biosynthetic processes producing MYCN and p53 compete for the same resources. It is true that two rules suggest that MYCN ’s transcription and translation rates can both be high, but they are exceptions rather than the norm. For the parametric combinations giving us WT fixed points with good outcomes in the third experiment (top-left panel, figure 9 ), the two distributions of MYCN rates are skewed towards their minima, whilst both distributions of p53 rates are skewed towards their maxima. These features are consistent with the need to produce more p53 than MYCN to keep p53/MYCN high (good outcome). One possibility is that a transcription factor upregulates biogenesis in general, but it preferentially boosts p53 ’s processes. The relatively uniform distributions of the six stressors are consistent with a mechanism where these regulatory processes are independent of the six stressors, mirroring the predominantly good clinical outcomes observed in WT NB cases. We further hypothesise that MYCN itself is the regulator. It follows from our prediction that in a real WT tumour, the biosynthetic rates in these parametric combinations should not be possible when there is not enough MYCN. In figures 7 and 8, some of the blue points with p53/MYCN above 3.5 and MYCN at low levels are probably biologically irrelevant. In the WT cases with poor outcomes ( table S2 ), the meta association rule implies that the processes of p53 transcription and translation compete for the same resources. The distributions of the four MYCN and p53 rates (top-right panel, figure 9 ) are all consistent with the need to produce MYCN more efficiently than p53 to keep p53/MYCN low (poor outcome). s 3 is active and p53 has a high rate in all 20 association rules, whilst the distribution of s 3 in the whole sample is skewed towards its minimum. One possibility is that when MYCN is scarce, ATM signalling (or a functionally equivalent pathway) is triggered to restore biogenesis in general, but it preferentially boosts MYCN production. If ATM signalling gets too high, p53 production will increase too (association rules and [ 100 ]), but MYCN production will still be higher to maintain the WT-poor scenario. When there is enough MYCN, ATM signalling shuts down and MYCN becomes responsible for biogenesis again (back to the WT-good sce-nario). In addition to explaining why WT tumours have poor outcomes when they lack MYCN, our hypothesis also predicts that the biosynthetic rates giving us WT fixed points with poor outcomes should not be possible when there is enough MYCN (ATM signalling is off). In figures 7 and 8, some of the blue points with p53/MYCN below 3.5 and MYCN at relatively high levels may be biologically irrelevant. The distributions of MYCN ’s and p53 ’s rates in the subset of parametric combintions giving us MA fixed points with good outcomes (bottom-left panel, figure 9 ) are as expected. A MA tumour must produce p53 more efficiently than it produces MYCN to achieve a good outcome. However, the meta association rule extracted from this subset ( table S3 ) suggests that MA patients require extremely high biosynthetic rates to achieve good outcomes. We speculate that even MA tumours with abundant MYCN rarely achieve the required rates. The distributions of s 3 and s 6 are skewed towards the maximum and minimum respectively, whilst the corresponding distributions in the WT cases with good outcomes (top-left panel, figure 9 ) are almost uniform. The implication is that upon MYCN amplification, only certain combinations of active stressors can result in a good outcome. Recall our prediction made in the previous paragraph that when there is enough MYCN (as in a MA tumour), ATM signalling (a part of s 3 ) should become inactive. The distribution of s 3 in the bottom-left panel of figure 9 contradicts this prediction. We speculate that MA tumours with abundant MYCN rarely achieve the combinations of active stressors required for good clinical outcomes. Specifically, they cannot keep ATM signalling on. It follows that the red points (and others that are not blue) with p53/MYCN above 3.5 in figures 7 and 8 are biologically rare even if they are possible. The distributions in the bottom-right panel of figure 9 (MA, poor outcomes) are aligned with the intuition that a MA tumour must produce MYCN more efficiently than it produces p53 to achieve a poor outcome (p53/MYCN below 3.5). The almost uniform distributions of all six stressors are consistent with the consensus that MA tumours usually have poor outcomes, regardless of the combination of active stressors therein. We speculate that MYCN preferentially boosts p53 ’s rates in WT tumours and MYCN ’s rates in MA tumours. Maybe, MYCN amplification occurs with a regu-latory shift, shifting the preference from p53 production to MYCN production. MYCN amplification could be responsible for this shift, either directly by producing more MYCN or indirectly by triggering a stress response. Alternatively, a co-occurring mutation could be responsible. For example, ALK mutations frequently co-occur with MYCN amplification; activated ALK and MYCN collaborate in NB pathogenesis [ 101 , 102 ]. 4.3 Biological basis of the expanded hypothesis In sub section 4.2 , we went beyond the scope of our model to speculate on the mecha-nistic details of our hypothesis, including why our results differ from parts of [ 23 ] and [ 24 ]. The first component of the expanded hypothesis states that MYCN upregulates biogenesis in general in NB tumours. The second component predicts that ATM signalling is activated when MYCN is scarce to restore biogenesis in general and MYCN production in particular. When MYCN’s level is restored, ATM signalling becomes inactive again. The third component states that MA tumours require unrealistically high biosynthetic rates for good outcomes. These rates are hard to attain; so is the high level of ATM activity when MYCN is abundant. The final component is that there is a regulatory shift when MYCN amplification occurs, causing MYCN to shift from preferentially boosting p53 production to MYCN production. This expanded hypothesis is consistent with the NB literature, as discussed below. First component. Although MYCN is known for its role as a transcription factor, it has numerous interaction partners. It interacts with different partners in different environments [ 103 ]. As reviewed by [ 103 ], MYCN interacts with DNA topoisomerase II proteins to resolve torsional stress in supercoiled DNA, prevents R-loop formation, resolves R-loops, interacts with the RNA exosome to degrade aberrant RNA species, and binds to RNA species directly to enhance stress resilience. These functions allow transcription to proceed efficiently. MA NB cells express high levels of PARP enzymes, which are activated at sites of DNA damage and upon activation, recruit repair factors [ 103 ]. A signalling pathway linking MYCN through PARP enzymes to DNA damage response (DDR) genes such as BRCA1 was experimentally identified by [ 104 ]. In the experiment [ 104 ], suppressing MYCN, PARP enzymes, and DDR genes led to reduced cancer cell survival. Finally, MYCN regulates protein synthesis [ 105 ] and folding [ 106 ], as well as cellular metabolism [ 103 , 107 ]. These functions allow translation to proceed efficiently. Second and third components. Considering MYCN counteracts so many stresses, including DNA damage (previous paragraph), it is not unreasonable to speculate that MYCN deficiency activates compensatory cellular stress responses to restore biogenesis. For example, the HSF1-mediated proteotoxic stress response is a transcriptional programme activated by misfolded proteins; it produces translation machinery components [ 108 ]. The integrated stress response activates an adaptive gene expression programme [ 109 ]. The ATM (a part of s 3 ) and ATR DDR pathways are both promoted [ 110 ]. Interestingly, ATM signalling activates p53 [ 100 ], mirroring the 20 association rules in table S2 , which state that p53 has a high rate when s 3 is active. This idea is even more appealing if NB cells are addicted to MYCN. Many cancers are addicted to oncogenes [ 111 ]. Addiction is likely in this particular case because despite DDR defects, NB tumours can tolerate high levels of stress thanks to MYCN’s actions [ 15 ]. For example, MYCN downregulates ATM by miR-421 [ 112 ] and repairs damaged DNA (previous paragraph). In fact, the first action is exactly what the second component predicts—when MYCN’s level is restored, ATM signalling becomes inactive again. Its existence provides a biological basis for the third component too—MA tumours have abundant MYCN and are unlikely to experience a high level of ATM activity. Fourth component. In a MA tumour, which always has enough MYCN, MYCN upregulates biogenesis. Compared to a WT tumour, however, MYCN preferentially boosts MYCN production rather than p53 production in a MA tumour. This regulatory shift could be caused by MYCN amplification itself or by mutations that co-occur with MYCN amplification in neuroblastoma, perhaps in the ALK gene [ 101 , 102 ]. This explanation is plausible because [ 113 ] experimentally demonstrated the ability of both WT ALK and gain-of-function ALK to initiate MYCN transcription. 4.4 Validity In the first experiment, only stable fixed points were found as predicted by [ 32 ]. Using the model, we predicted two novel emergent phenomena, namely that E2F-1 does not transcribe ARF cooperatively and a non-linear relationship between the MYCN copy number of a patient and their outcome. Our model’s ability to reproduce and predict emergent phenomena makes it a contribution to the field in its own right. The fact that our simulation results are consistent with most of the clinical and experimental observations listed in section S1 in the supplementary file is encouraging. Also reassuring is that searching the parametric space systematically (experiment two) and randomly (experiment three) gave us similar results. Our hypothesis that the p53/MYCN metric decides a NB tumour’s clinical/treatment outcome is valid and deserves further investigation. Whilst our results differ from [ 23 ] and [ 24 ] in some respects, the expanded hypothesis proposed in sub section 4.2 reconciles the differences in theory. Furthermore, it is sup-ported by the publications mentioned in sub section 4.3 , so we believe that it deserves attention from experimental and clinical researchers. Nevertheless, it must be acknowledged that the model’s structure and the chosen parametric values form an incomplete representation of how things work within neuroblastoma cells. 4.5 Future work The four components of the expanded hypothesis, which we proposed in sub section 4.2 by adding mechanistic links not already in the mathematical model, are not testable even in silico . Therefore, it remains a tentative explanation of the MYCN enigma. For it to become a well-substantiated theory explaining how MYCN affects a heterogeneous population of NB cells by interacting, cooperating, and competing with p53 in different cellular contexts (combinations of mutations), we must build a more comprehensive mathematical model to guide an experimental study. The improved model would replace the eight static biosynthetic rates with functions dependent on the protein level of MYCN and the mutations present in a NB cell. It would include a negative feedback loop between MYCN and s 3 too. The static stressor would be replaced by a set of equations modelling a set of stress responses to MYCN deficiency. Based on sub section 4.3 , the HSF1-mediated proteotoxic stress response [ 108 ], the integrated stress response [ 109 ], and various DDR pathways (including ATM and ATR signalling) [ 100 , 110 ] are potential members of this set. Closing the loop, the new variables representing the selected stress responses would be connected to the biosynthetic rates. More fundamentally, the p53/MYCN metric is a population-level metric. When the hypothesis refers to its level in a cell, it is actually the average cell in a NB tumour. The hypothesis assumes that p53 is pro-apoptosis in a population of identical NB cells. However, as we proposed recently [ 114 ], a MA tumour may contain two distinct populations of cancer cells, where p53 is pro-survival and pro-apoptosis respectively. The basis is that p53 regulates hundreds of genes in a context-dependent manner, so it has diverse and contradictory downstream effects, including DNA repair, cell cycle arrest, and apoptosis [ 115 ]. The updated mathematical model should be incorporated into each NB cell agent in the multicellular model presented in our recent publications [ 114 , 116 ]. Currently, although MYCN and p53 are both multifunctional as required in the multicellular model and it naturally describes a heterogeneous population, the MYCN enigma is modelled simplistically as a conditional statement involving three Boolean variables. It states that if MYCN is amplified and MYCN is active in a NB cell agent, then p53 is inactive therein. The updated mathematical model, which would describe the molecular network more accurately, would represent the MYCN enigma better. One could use the modified multicellular model to test the population-level p53/MYCN metric’s ability to predict clinical/treatment outcomes in different scenarios. p53 has different functions, which could be modelled as a probability distribution over a heterogeneous population of NB cells. This change would generalise the p53/MYCN metric and make it applicable at the cellular level of resolution. 5 Conclusion The MYCN enigma—the observation that the outcome of neuroblastoma does not depend monotonically on MYCN expression or amplification—motivated our hypothesis that the amount of p53 protein relative to the amount of MYCN protein is a superior prognostic metric. Developing a novel mathematical model linking MYCN to the ARF/MDM2/p53 axis, we carried out numerical simulations to test our hypothesis. The results are aligned with numerous clinical and experimental observations. After extracting association rules and parametric correlations from the tested parametric combinations, we expanded the hypothesis by adding mechanisms outside the scope of the mathematical model. According to the expanded hypothesis, MYCN is responsible for biogenesis in NB cells and MYCN amplification and/or co-occurring mutations trigger a regulatory shift to reduce p53/MYCN (worse outcome). The expanded hypothesis also predicts a negative feedback loop between MYCN deficiency and ATM signalling (or equivalent), which restores biogenesis and favours MYCN production in a WT tumour deficient in MYCN. Finally, it predicts that a MA tumour requires very high biosynthetic rates and ATM signalling (or equivalent) to achieve a good outcome (p53/MYCN above 3.5). These conditions are very hard to attain biologically. The expanded hypothesis is supported by the literature, but biosynthesis and stress responses to MYCN deficiency are represented simplistically in the model. More mechanistic details about them must be added to the model. The multifunctional nature of p53 and tumour heterogeneity must be taken into account too. Such an improved model could enable more detailed, granular, and comprehensive simulations, which could in turn inform lab-based experiments and therapeutic design. 6 Data availability statement This manuscript has no associated data. Everything is included in the main manuscript and supplementary information section. 7 Author contributions and acknowledgments KYW conceived the project and secured funding. MP and KYW built the mathematical model with input from DW and RD. MI found the model parameters with input from KYW. KYW and MI proposed that p53/MYCN is a prognostic metric for neuroblastoma. MI designed and implemented the numerical simulations with input from KYW and FD. MI analysed the simulation results with input from KYW and FD. MP and MI created the visualisations with input from KYW and FD. KYW proposed the expanded hypothesis with input from MI, FD, and DW. KYW, MI, FD, and RD wrote the manuscript. KYW and DW supervised MP and RD. KYW and FD supervised MI. All six authors have read and approved the manuscript. KYW and RD are grateful for funding provided by the Insigneo Institute for in Silico Medicine. The authors want to thank Daniel Jordan for the preliminary results that inspired this project. Supporting Information S1 Clinical and experimental observations of relevance to the MYCN enigma A satisfactory explanation of the MYCN enigma must be consistent with the following clinical and experimental observations. In general, as the copy number of MYCN goes up in a NB tumour, its clinical and biological characteristics get progressively more adverse to survival [ 99 ]. Out of the 23 WT patients [ 24 ] studied, 18 achieved good outcomes (78.26 %) and five achieved poor outcomes (21.74 %). The six MA patients they studied all achieved poor outcomes. In another study [ 23 ], only three out of 15 MA patients achieved good outcomes (20 %) and the rest achieved poor outcomes (80 %). Together, the two datasets [ 23 , 24 ] suggest that an average MA patient has a 14.29% chance of achieving a good outcome. [ 24 ] measured the levels of MYCN mRNA and MYCN protein in their patients’ tumours. The results are reported in Fig. 3C in their paper. Tumours with similar MYCN mRNA levels were found to have very different MYCN protein levels. [ 23 ] measured the levels of MYCN mRNA and TrkA mRNA in 76 WT primary NB specimens and 15 MA ones. The results are reported in Figure 2A in their paper. They measured almost uniformly higher MYCN mRNA levels in the MA specimens than the WT specimens. It should be noted that whilst TrkA expression is indeed a favourable biomarker for neuroblastoma patients, it does not follow that TrkA mechanistically suppresses neuroblastoma. Again in Figure 2A in [ 23 ], almost every Trk mRNA level reported for a MA speci-men is close to zero (poor prognosis). The values reported for the WT specimens are much more diverse, showing a positive and linear relationship with the corresponding MYCN mRNA levels. Overall, the figure presents a concave function between the two sets of mRNA levels. NB cells are cancer cells. If the explanation involves p53 , it must agree with a recent theory [ 32 ]. It states that oncogenic mutations abolish p53 oscillations. S2 Methodological details S2.1 Master Equations Our gene regulatory network comprises four genes. Each pair of master equations describes the concentration dynamics of a gene’s mRNA transcript and protein. Equation S1 models transcription in terms of the transcript’s production and degradation rates. Equation S2 describes translation in terms of the protein’s production and degradation rates. In other words, each equation describes phenomena occurring at two timescales, namely synthesis and degradation. The two equations are linked by the mRNA transcript’s concentration level, on which the translation rate depends. The master equations for gene i , where m i ( t ) and p i ( t ) describe the concentration levels of its mRNA transcript and protein, are given below: In equations S1 and S2 , k i represents the transcription rate of gene i and l i represents its translation rate. In the same equations, δ mi and δ pi are the degradation rates of its mRNA transcript and protein respectively. After tailoring equations S1 and S2 to describe the four genes, we will use more readable variables than m i ( t ) and p i ( t ) to represent the eight concentrations (pM): [ MY CNm ] and [ MY CN ]; [ p 53 m ] and [ p 53]; [ MDM 2 m ] and [ MDM 2]; and [ ARFm ] and [ ARF ]. MYCN amplification can be simulated by increasing k MY CN . The key property of this pair of equations is that they are both linear, but the model pertains to a non-linear problem, so it should be non-linear. The model is non-linear because it also describes non-linear binding events. First, stressors outside the gene regulatory network can modulate the four fundamental processes. Second, the four proteins interact to form complexes to modulate these processes. Third, other complexation events contribute non-linearly too. S2.2 Modelling protein complexation Our model describes protein complexation in accordance with the procedures explained in Chapter 7 of [ 50 ]. The underlying assumption is that protein-protein interactions occur on a much faster timescale than the functions (transcription and translation) carried out by proteins and protein complexes, as well as their degradation. This leads to the simplifying assumption that protein-protein interactions are in equilibrium on the model’s slower timescale. The equilibrium dissociation constant of each complex describes the equilibrium of the relevant binding event. For example, applying this simplification to complex A, the fraction of unbound MYCN is as follows: where [ ARF ] is the protein concentration of ARF. Within our model, MYCN can only form one complex by binding to ARF (Complex A). Other proteins, especially MDM2, can exist in multiple complexes, sometimes with more than one complexation partner. The assumption underlying equation S3 can be generalised to this case, but the unbound fraction of such a protein is expressed in a more complex form. Again, the derivation procedures are explained in Chapter 7 of [ 50 ]. For example, the unbound fraction of MDM2 is as follows: Due to its ability to reside in multiple complexes, MDM2 cannot be fully understood in terms of a binary status (bound or unbound). The bound fraction of MDM2 is divided into further fractions. For example, the fraction of MDM2 in the dimer complex B is . Complex F is a trimer, so the associated fraction is more complex: . Again, the derivation procedures are explained in Chapter 7 of [ 50 ]. Equation S4 also exemplifies how five of the six environmental stressors (all except s 1 ) are represented in the model. As reviewed, s 3 prevents the formation of complex C, s 4 prevents the formation of complex D, and s 5 promotes the formation of complex E. Promotion and prevention are modelled by multiplying a term by the associated stressor variable (such as s 5 ) and its complement (such as 1 − s 3 ) respectively. Conceptually, these environmental stressors are like switches. They modulate some of the cellular processes described by the model in an abstract manner, but they do not represent specific molecular species. Derived similarly, the following equations are the fractions of unbound p53 and ARF: S2.3 Modelling protein-DNA interactions Mathematically, these binding events are handled similarly to protein complexation, but the equilibrium dissociation constants are labelled numerically rather than alphabetically. Cooperative binding is a mechanism that regulates protein-DNA interactions but not protein complexation. It is modelled by the Hill equation, which is parameterised by the Hill coefficient. When the Hill coefficient is bigger than one, it means that a copy of macromolecule A binds to a copy of macromolecule B more easily when the latter is already bound to other copies of A. Conversely, when the Hill coefficient is smaller than one, a copy of B loses some of its affinity for A with each binding event. When it is exactly one, cooperative behaviour is irrelevant. Most of the binding events described by our model are assumed to be non-cooperative, so they are represented by Michaelis–Menten terms. Examples are equations S3 and S4 . One major exception is the interaction between p53 and the genome. When p53 binds to DNA, it does so sequence-specifically as a tetramer and occupies two consecutive half-sites [ 117 , 118 ]. The Hill coefficient for this interaction is 1.8 [ 118 , 119 ]. We will illustrate how the Hill equation is incorporated into the terms associated with p53-DNA binding with an example. The promoter region of p53 is either free (unbound), bound to p53, bound to MYCN, or bound to both proteins. It is assumed in our model that when there are two transcriptional activators for a promoter, they bind to the promoter independently. Therefore, MYCN [ 25 ] and p53 [ 88 ] do not interact with DNA cooperatively in our model. The following equation models the fraction of p53 promoters that are unbound: Another exception is the influence of s 1 on the promoter region of ARF . In a stressed cell, ARF transcription is upregulated by proteins such as E2F-1, β -catenin, and Ras [ 60 , 64 , 65 ]. Earlier in subsection S2.2, we illustrated how five of the six environmental stressors are represented as simple switches in our model, that they do not represent specific molecular species. The stressor s 1 , which upregulates ARF transcription, is an exception. Within our model, it acts like an abstract transcription factor representing the set of proteins that upregulate ARF transcription. It can interact with the other molecular species in the model, but there is no evidence of cooperative binding in these interactions, so Michaelis–Menten kinetics is assumed. The fraction of ARF promoters that are unbound is given by the following equation: Derived similarly, the following equation is the fraction of MDM2 promoters that are unbound: S2.4 Modelling protein-mRNA interactions The fraction of unbound MYCN mRNA is given by the following equation: Although the positive interaction between MDM2 and p53 mRNA can be described by a Michaelis–Menten term as usual, the negative and indirect influence of MDM2 on p53 translation must be treated abstractly because L26 is not in the model. The solution is to assume that MDM2 either ‘occupies’ an activating site or an inhibiting site on a p53 mRNA transcript depending on whether the cell is stressed or not. As explained in subsection S2.2, s 2 acts like a simple switch modulating the two sites’ competition for MDM2. This representation is not mechanistically accurate, but it captures the functional consequences of MDM2’s dual role in a phenomenological sense (different inputs give different outputs). Based on these principles, the following equation gives us the fraction of unbound p53 mRNA: S2.5 Modelling functional effects of complexation Our model describes the functional effects of complexation using intuitive additive terms. When it is in a complex, a protein may degrade at a different rate or be incapable of transcribing another gene (transactivation). For example, MYCN degrades at two different rates, d MY CN (fast) and d MY CN−A (slow), in its unbound and bound states respectively. Therefore, the overall protein degradation rate of MYCN is [ θ MY CN d MY CN + (1 − θ MY CN ) d MY CN −A ][ MY CN ]. S2.6 Model parameters The model parameters can be divided into two major categories. The first category contains the rates of mRNA production, mRNA degradation, protein synthesis, and protein degradation in equations 1 to 8 . The second category contains the equilibrium dissociation constants parameterising the myriad of complexation events described in subsections 2.1.2, 2.1.3, and S2.4 in the main manuscript and subsections S2.2, S2.3, and S2.4 in the supplementary file. S2.6.1 Fundamental processes The natural mRNA and protein degradation rates, denoted by d m and d p , are 0.06 h -1 and 0.05 h -1 respectively [ 51 ]. It is assumed that MYCN mRNA transcripts degrade 10 times more slowly than the natural rate ( d MY CNm = 0.1 d m = 0.006 h -1 ). In complex C, p53 degrades faster than usual. In complexes B and E, MDM2 degrades faster than usual. It is assumed that these three enhanced protein degradation rates are 10 times higher than the natural rate ( d p 53 − C = d MDM 2 − BE = 10 d p = 0.5 h -1 ). By setting the derivatives in equations S1 and S2 to zero, it is possible to express the steady-state mRNA and protein abundance levels as and respectively. In order to estimate the transcription and translation rates of each species from the two expressions, we need the steady-state values of [ m i ] and [ p i ]. The Gene Expression Omnibus (GEO) repository contains mRNA sequencing [ 52 ] and ribosome profiling [ 53 ] datasets. Specifically, one dataset (accession number: GSE89413) records mRNA abundance levels in NB cells taken from different cell lines (including SH-SY5Y) as reads per kilobase million (RPKM), whilst a ribosome profiling dataset (accession number: GSE119615) records protein abundance levels in NB cells taken from the SH-SY5Y cell line as RPKM. First, we discuss the transcription rates. The steady-state mRNA abundance level ([ m i ]) of gene i is proportional to its RPKM value ( r i ) in the mRNA sequencing dataset divided by the sum of all the RPKM values (Σ r i ) therein. This relative transcription level can be made absolute by multiplication with the total copy number of mRNA in a cell ( N m = 2.3 × 10 4 copies per cell [ 51 ]). Therefore, [ m i ] = r i N m /Σr i . This quantity is multiplied by the natural mRNA degradation rate ( δ mi = d m = 0.06 h -1 ) to obtained k i = [ m i ] δ m . Considering the mRNA sequencing dataset describes the overall behaviour of NB cells, the model assumes that each estimate given by k i = [ m i ] δ mi corresponds to the highest transcription rate available to a gene when there are more than one rate. The estimates are denoted by α MYCN , β p53 , β MDM2 , and β ARF , which represent the natural transcription rate of MYCN and the enhanced transcription rates of the other three genes. In our model, MYCN is always transcribed naturally and the other three genes are transcribed faster when their activating mechanisms are active. In addition to having a fixed transcription rate, MYCN is a special case because the copy number of MYCN ( g MY CN ) is hidden in the estimate k MY CN , meaning k MY CN = g MY CN α MY CN . The SH-SY5Y cell line is not MA [ 120 ], so g MY CN = 1 and k MY CN = α MY CN for the purpose of parameterisation. The natural transcription rates of p53 , MDM2 , and ARF ( α p53 , α MDM2 , and α ARF ) are assumed to be 10 times slower than the corresponding enhanced rates. Moving on to the translation rates, the RPKM values ([ p i ]) in the ribosome profiling dataset can be used in a similar manner. First, [ p i ] /Σ [ p i ] estimates the relative translation level of gene i . Multiplying this estimate by the total number of proteins in a cell ( N p = 7.5 × 10 9 proteins per cell [ 51 ]) makes it absolute. Multiplying the absolute translation level by δ pi / [ m i ] gives the translation rate of gene i as l i = ( δ pi / [ m i ])( N p [ p i ] /Σ [ p i ]). The appearance of [ m i ] in the denominator means the rate is independent of the amount of mRNA. The natural protein degradation rate is d p = 0.05 h -1 . The steady-state abundance level of the corresponding mRNA [ m i ] is as described in the previous paragraph. Considering the dataset describes the overall behaviour of NB cells, the model assumes that each estimate l i = ( δ pi / [ m i ])( N p [ p i ] /Σ [ p i ]) corresponds to the highest translation rate available to a gene. Therefore, the four estimates are λ MYCN , µ p53 , λ MDM2 , and λ ARF . Apart from p53, each of them is produced at just one natural rate. p53 has a natural translation rate ( λ p53 ), a reduced translation rate ( κ p53 ), and an enhanced translation rate ( µ p53 ). The estimate is µ p53 and the natural rate λ p53 and reduced rate κ p53 are assumed to be 10 and 100 times slower respectively. S2.6.2 Equilibrium dissociation constants and Hill coefficient As reported in table 2 , literature values are available for many of these parameters. Two exceptions are K A and K B , so the model assumes a value of 100 nM based on the average order of magnitude of the other equilibrium dissociation constants. K 5 is the equilibrium dissociation constant for the interaction between the abstract transcription factor s 1 ; representing E2F-1, β -catenin, and Ras; and the ARF gene’s promoter region. Since s 1 is dimensionless and ranges from zero to one, K 5 = 1. S2.7 Computer simulations We conducted three sets of computer simulations. First, we varied three groups of parameters with and without MYCN amplification in a grid search. The results were used to confirm the model’s stable long-term behaviour (equilibrium), which represents the clinical/treatment outcome; understand the outcome’s sensitivity to the three groups of parameters; and establish a NB classification system based on our hypothesis. Second, we focused on the two most sensitive groups of parameters, repeating the grid search without MYCN amplification and at five levels of MYCN amplification. Third, we validated the conclusions of the first two experiments by exploring the parametric space randomly. S2.7.1 First grid search In the first experiment, we varied three groups of parameters together and for each parametric combination, we tested the model’s long-term stability. We tested 2 6 = 64 input combinations after binarising the six stressors so that s i = 0.9 or s i = 0. The unknown Hill coefficient h 5 was strategically varied from 0.2 (negative cooperativity) to one (no cooperativity) and then to four (positive cooperativity) to cover all three regimes. We tested three values for each of the four transcription rates ( α MYCN , β p53 , β MDM2 , and β ARF ) and four basal translation rates ( λ MYCN , µ p53 , λ MDM2 , and λ ARF ) based on experimental measurements. For each rate, the estimate in table 1 , 70 % of the estimate, and 130 % of the estimate were tested. For a process with two or three possible rates, all of them were varied to maintain the ratios between them. For example, when 70 % of the estimate of µ p53 was used, λ p53 was also reduced by 30 %. Within group three, we tested 3 8 = 6561 combinations of rates. The equilibrium dissociation constants were fixed under the assumption that they do not vary between two patients. Since s 1 is relevant to h 5 in the model, when h 5 was varied (0.2, one, and four) with the other stressors off (zero) and the rates at their default values, s 1 was set to an additional value: 0.5. It means three additional parametric combinations were tested. In total, 64 · 3 · 6561 + 3 = 1259715 parametric combinations were tested in the first experiment. Each parametric combination was used to perform two sets of 10 simulations within the MATLAB environment with the ode45 solver function, once with g MYCN = 1 and once with g MYCN = 2. For each combination and each MYCN amplification status, 10 simulations were conducted from 10 random initial conditions. An initial condition was created by sampling each of the eight variables randomly from the continuous uniform distribution in the range [0, 1]. A simulation was terminated at t = 10 6 or when an equilibrium was reached before t = 10 6 . Equilibria were detected by computing the ratio between the norm of the state variables’ derivatives and the norm of the state variables themselves at the end of each time step during the simulations. When the ratio fell below the cutoff threshold of 10 − 6 , the simulation was terminated. In mathematical terms, these equilibria are called stable fixed points, which correspond to clinically observable and relevant scenarios. The identified fixed points’ stability was tested numerically. Starting from each fixed point, we perturbed each state variable to a random value in the range between zero and 10 times the equilibrium value. The simulation was repeated from this perturbed condition and if the system (eight-dimensional vector comprising the state variables) returned to the fixed point, the fixed point was classified as stable. S2.7.2 Second grid search In the second experiment, we wanted to learn more about how the three groups of parameters interact to affect the system’s long-term behaviour. Informed by the results of the first experiment, we set h 5 = 1 and varied s i and the rates together. We did so at six levels of MYCN amplification, including the wild-type genotype. Within the MATLAB environment, 64 · 6561 = 419904 combinations of the parametric values tested in the first experiment were used to perform six sets of simulations with the ode45 solver function, one for each of g MYCN = 1, g MYCN = 2, g MYCN = 3, g MYCN = 4, g MYCN = 6, and g MYCN = 8. The same initialisation and termination procedures were adhered to. S2.7.3 Random exploration of the parametric space In the third experiment, in order to validate the conclusions of the two grid search experiments, we generalised them by sampling from the whole parametric space. We did so by varying s i and the rates together and randomly, whilst keeping h 5 = 1, leading to 10000 random parametric combinations. Each combination was generated by sampling each stressor and rate randomly and uniformly between the minimum and maximum values tested in the first two experiments. Then, like the second experiment, each parametric combination was used to perform six simulations with g MYCN = 1, g MYCN = 2, g MYCN = 3, g MYCN = 4, g MYCN = 6, and g MYCN = 8. The other procedures followed in the second experiment were adhered to too. S2.8 Data mining and visualisation Experiment two. The Apriori algorithm computes three metrics in order to identify frequent itemsets and association rules in each subset of parametric combinations: support, confidence, and lift. We will illustrate them using the parameters s 1 and s 2 in the subset of WT equilibria predicted to have good clinical/treatment outcomes. The support of s 1 being on (0.9) and s 2 being off (0) in this subset is the number of parametric combinations where s 1 = 0.9 and s 2 = 0 divided by the total number of combinations in this subset. If the support is above a user-defined threshold, the combination of s 1 being on (0.9) and s 2 being off (0) is a frequent itemset. The confidence of s 1 being on (0.9) with respect to s 2 being off (0) is the number of combinations where s 1 = 0.9 and s 2 = 0 divided by the number of combinations where s 2 = 0 in the subset. In other words, it is the support of s 1 being on (0.9) and s 2 being off (0) divided by the support of s 2 being off (0). This metric tells you how likely it is that s 1 = 0.9 in a parametric combination if s 2 = 0 therein in this subset. A major problem with this metric is that it is sensitive to the frequency of an item in a dataset. For example, if s 1 = 0.9 throughout the subset, the confidence of s 1 being on (0.9) with respect to s 2 being off (0) will of course be 100 %, but it does not buttress the claim that the likelihood of the former is boosted by the truth of the latter. The lift mitigates this problem of the confidence metric by taking the baseline occurrence of the item being boosted into account. For example, the lift of s 1 being on (0.9) with respect to s 2 being off (0) takes the baseline occurrence of s 1 being on (0.9) in the subset into account. It is the confidence of s 1 being on (0.9) with respect to s 2 being off (0) divided by the support of s 1 being on (0.9). If the lift is above one, it means that s 2 being off (0) increases the likelihood of s 1 being on (0.9) in a parametric combination in the subset. This association rule is denoted in this way: s 2 = 0 ⇒ s 1 = 0.9. Experiment three. Using the ggpairs function within the GGally R package, we created an augmented scatterplot matrix (also called a ggplot2 generalised pairs plot) for each subset of random parametric combinations. It is a 14×14 matrix because in the third experiment, we took random samples from a 14-dimensional parametric space (six stressors and eight rates). Below the leading diagonal are the 91 scatterplots comparing the unique pairs formed by the 14 parameters. On the diagonal are the probability density functions of the 14 parameters in this random sample. Above the diagonal, the Pearson correlation coefficients of the 91 pairs of parameters are reported. Each correlation’s statistical significance is shown too. In a matrix element, *** means the reported p-value is smaller than 0.001; ** means that it is smaller than 0.01; * means that it is smaller than 0.05; and a full stop (.) means it is smaller than 0.10. S3 Additional results Below are the four sets of 20 association rules (sub section 3.4 ) extracted from the results of the second experiment with the Apriori algorithm. We divided the fixed points found in the experiment into four subsets (genotype-outcome combinations): WT ( g MYCN = 1) or MA ( g MYCN = 2) and good (p53/MYCN above 3.5) or poor (p53/MYCN below 3.5) clinical/treatment outcomes. We applied the Apriori algorithm to the parametric combinations that gave us the fixed points in each subset. View this table: View inline View popup Download powerpoint Table S1: Top 20 association rules regarding the parametric combinations associated with WT fixed points with good outcomes in the second experiment. L X , M X , H X represent transcription rates of X that are 70%, 100%, 130% of the default (also called middle) value. Low X , Mid X , High X represent translation rates of X that are 70%, 100%, 130% of the default (also called middle) value. sion and sioff represent the active and inactive statuses of s i respectively. LHS denotes the antecedent (left-hand side) of a rule, whilst RHS denotes the consequent (right-hand side) of a rule. View this table: View inline View popup Download powerpoint Table S2: Top 20 association rules regarding the parametric combinations asso-ciated with WT fixed points with poor outcomes in the second experiment. L X , M X , H X represent transcription rates of X that are 70%, 100%, 130% of the default (also called middle) value. Low X , Mid X , High X represent translation rates of X that are 70%, 100%, 130% of the default (also called middle) value. sion and sioff represent the active and inactive statuses of s i respectively. LHS denotes the antecedent (left-hand side) of a rule, whilst RHS denotes the consequent (right-hand side) of a rule. View this table: View inline View popup Download powerpoint Table S3: Top 20 association rules regarding the parametric combinations associated with MA fixed points with good outcomes in the second experiment. L X , M X , H X represent transcription rates of X that are 70%, 100%, 130% of the default (also called middle) value. Low X , Mid X , High X represent translation rates of X that are 70%, 100%, 130% of the default (also called middle) value. sion and sioff represent the active and inactive statuses of s i respectively. LHS denotes the antecedent (left-hand side) of a rule, whilst RHS denotes the consequent (right-hand side) of a rule. View this table: View inline View popup Download powerpoint Table S4: Top 20 association rules regarding the parametric combinations associated with MA fixed points with poor outcomes in the second experiment. L X , M X , H X represent transcription rates of X that are 70%, 100%, 130% of the default (also called middle) value. Low X , Mid X , High X represent translation rates of X that are 70%, 100%, 130% of the default (also called middle) value. sion and sioff represent the active and inactive statuses of s i respectively. LHS denotes the antecedent (left-hand side) of a rule, whilst RHS denotes the consequent (right-hand side) of a rule. Funder Information Declared University of Sheffield , Insigneo Summer Research Programme References [1]. ↵ Matthay , K. , Maris , J. , Schleiermacher , G. , Nakagawara , A. , Mackall , C. , Diller , L. , Weiss , W. : Neuroblastoma. Nature reviews . Disease primers 2 , 16078 ( 2016 ) OpenUrl PubMed [2]. ↵ London , W. , Castleberry , R. , Matthay , K. , Look , A. , Seeger , R. , Shimada , H. , Thorner , P. , Brodeur , G. , Maris , J. , Reynolds , C. , et al : Evidence for an age cutoff greater than 365 days for neuroblastoma risk group stratification in the children’s oncology group . Journal of clinical oncology 23 ( 27 ), 6459 – 6465 ( 2005 ) OpenUrl Abstract / FREE Full Text [3]. ↵ Tomolonis , J.A. , Agarwal , S. , Shohet , J.M .: Neuroblastoma pathogenesis: dereg-ulation of embryonic neural crest development . Cell and tissue research 372 ( 2 ), 245 – 262 ( 2018 ) OpenUrl CrossRef PubMed [4]. ↵ Kameneva , P. , Artemov , A.V. , Kastriti , M.E. , Faure , L. , Olsen , T.K. , Otte , J. , Erickson , A. , Semsch , B. , Andersson , E.R. , Ratz , M. , et al : Single-cell transcriptomics of human embryos identifies multiple sympathoblast lineages with potential implications for neuroblastoma origin . Nature genetics 53 ( 5 ), 694 – 706 ( 2021 ) OpenUrl CrossRef PubMed [5]. ↵ Crane , J.F. , Trainor , P.A .: Neural crest stem and progenitor cells . Annu. Rev. Cell Dev. Biol . 22 , 267 – 286 ( 2006 ) OpenUrl CrossRef PubMed Web of Science [6]. ↵ Schulte , J. , Lindner , S. , Bohrer , A. , Maurer , J. , De Preter , K. , Lefever , S. , Heukamp , L. , Schulte , S. , Molenaar , J. , Versteeg , R. , et al : Mycn and alkf1174l are sufficient to drive neuroblastoma development from neural crest progenitor cells . Oncogene 32 ( 8 ), 1059 – 1065 ( 2013 ) OpenUrl CrossRef PubMed [7]. ↵ Monclair , T. , Brodeur , G.M. , Ambros , P.F. , Brisse , H.J. , Cecchetto , G. , Holmes , K. , Kaneko , M. , London , W.B. , Matthay , K.K. , Nuchtern , J.G. , et al : The international neuroblastoma risk group (inrg) staging system: an inrg task force report . Journal of clinical oncology 27 ( 2 ), 298 ( 2009 ) OpenUrl Abstract / FREE Full Text [8]. ↵ Esposito , M.R. , Aveic , S. , Seydel , A. , Tonini , G.P .: Neuroblastoma treatment in the post-genomic era . Journal of biomedical science 24 ( 1 ), 1 – 16 ( 2017 ) OpenUrl CrossRef PubMed [9]. ↵ Morgenstern , D.A. , Baruchel , S. , Irwin , M.S .: Current and future strategies for relapsed neuroblastoma: challenges on the road to precision therapy . Journal of pediatric hematology/oncology 35 ( 5 ), 337 – 347 ( 2013 ) OpenUrl [10]. ↵ Matthay , K.K. , George , R.E. , Yu , A.L .: Promising therapeutic targets in neuroblastoma . Clinical Cancer Research 18 ( 10 ), 2740 – 2753 ( 2012 ) doi: 10.1158/1078-0432.CCR-11-1939 OpenUrl Abstract / FREE Full Text [11]. ↵ Italia , M. , Wertheim , K.Y. , Taschner-Mandl , S. , Walker , D. , Dercole , F .: Mathematical model of clonal evolution proposes a personalised multi-modal therapy for high-risk neuroblastoma . Cancers 15 ( 7 ) ( 2023 ) doi: 10.3390/ cancers15071986 OpenUrl CrossRef PubMed [12]. ↵ Otero , J.G. , Bodzioch , M. , Belmonte-Beitia , J. : On the dynamics and optimal control of a mathematical model of neuroblastoma and its treatment: Insights from a mathematical model . Mathematical Models and Methsods in Applied Sciences 34 ( 07 ), 1235 – 1278 ( 2024 ) OpenUrl [13]. ↵ Tucker , E.R. , Jiménez , I. , Chen , L. , Bellini , A. , Gorrini , C. , Calton , E. , Gao , Q. , Che , H. , Poon , E. , Jamin , Y. , et al : Combination therapies targeting alk-aberrant neuroblastoma in preclinical models . Clinical Cancer Research , 1 – 15 ( 2023 ) [14]. ↵ Pucci , P. , Lee , L.C. , Han , M. , Matthews , J.D. , Jahangiri , L. , Schlederer , M. , Manners , E. , Sorby-Adams , A. , Kaggie , J. , Trigg , R.M. , et al : Targeting nras via mir-1304-5p or farnesyltransferase inhibition confers sensitivity to alk inhibitors in alk-mutant neuroblastoma . Nature Communications 15 ( 1 ), 3422 ( 2024 ) OpenUrl PubMed [15]. ↵ Southgate , H.E. , Chen , L. , Curtin , N.J. , Tweddle , D.A .: Targeting the dna dam-age response for the treatment of high risk neuroblastoma . Frontiers in oncology 10 , 371 ( 2020 ) OpenUrl PubMed [16]. ↵ Otero , J.G. , Belmonte-Beitia , J. , Jiménez-Śanchez , J. : Exploring neuroblas-toma’s cellular microenvironment: A novel approach using cellular automata to model celyvir treatment . Computers in Biology and Medicine 188 , 109782 ( 2025 ) OpenUrl PubMed [17]. ↵ Santiago , T. , Tarek , N. , Boulos , F. , Hayes , C. , Jeha , S. , Raimondi , S. , Rodriguez-Galindo , C .: Correlation between mycn gene status and mycn protein expression in neuroblastoma: a pilot study to propose the use of mycn immunohis-tochemistry in limited-resource areas . Journal of Global Oncology 5 , 1 – 7 ( 2019 ) OpenUrl CrossRef PubMed [18]. ↵ Wang , M. , Zhou , C. , Cai , R. , Li , Y. , Gong , L .: Copy number gain of mycn gene is a recurrent genetic aberration and favorable prognostic factor in chinese pediatric neuroblastoma patients . Diagnostic Pathology 8 , 1 – 6 ( 2013 ) OpenUrl CrossRef PubMed [19]. ↵ Campbell , K. , Naranjo , A. , Hibbitts , E. , Gastier-Foster , J.M. , Bagatell , R. , Irwin , M.S. , Shimada , H. , Hogarty , M. , Park , J.R. , DuBois , S.G. : Association of heterogeneous mycn amplification with clinical features, biological character-istics and outcomes in neuroblastoma: A report from the children’s oncology group . European Journal of Cancer 133 , 112 – 119 ( 2020 ) OpenUrl PubMed [20]. ↵ Yue , Z.-X. , Huang , C. , Gao , C. , Xing , T.-Y. , Liu , S.-G. , Li , X.-J. , Zhao , Q. , Wang , X.-S. , Zhao , W. , Jin , M. , et al : Mycn amplification predicts poor prognosis based on interphase fluorescence in situ hybridization analysis of bone marrow cells in bone marrow metastases of neuroblastoma . Cancer cell international 17 , 1 – 7 ( 2017 ) OpenUrl PubMed [21]. ↵ Huang , M. , Weiss , W.A. : Neuroblastoma and mycn . Cold Spring Harbor Perspectives in Medicine 3 ( 10 ) ( 2013 ) doi: 10.1101/cshperspect . a014415 OpenUrl CrossRef [22]. ↵ Mueller , S. , Matthay , K.K. : Neuroblastoma: biology and staging . Current oncology reports 11 , 431 – 438 ( 2009 ) OpenUrl CrossRef PubMed [23]. ↵ Tang , X.X. , Zhao , H. , Kung , B. , Kim , D.Y. , Hicks , S.L. , Cohn , S.L. , Cheung , N.-K. , Seeger , R.C. , Evans , A.E. , Ikegaki , N .: The mycn enigma: significance of mycn expression in neuroblastoma . Cancer research 66 ( 5 ), 2826 – 2833 ( 2006 ) OpenUrl Abstract / FREE Full Text [24]. ↵ Valentijn , L.J. , Koster , J. , Haneveld , F. , Aissa , R.A. , Van Sluis , P. , Broekmans , M.E.C. , Molenaar , J.J. , Van Nes , J. , Versteeg , R .: Functional mycn signature predicts outcome of neuroblastoma irrespective of mycn amplification . Proceedings of the National Academy of Sciences of the United States of America 109 ( 47 ), 19190 – 19195 ( 2012 ) doi: 10.1073/pnas.1208215109 OpenUrl Abstract / FREE Full Text [25]. ↵ Chen , L. , Iraci , N. , Gherardi , S. , Gamble , L.D. , Wood , K.M. , Perini , G. , Lunec , J. , Tweddle , D.A .: p53 is a direct transcriptional target of mycn in neuroblastoma . Cancer research 70 ( 4 ), 1377 – 1388 ( 2010 ) OpenUrl Abstract / FREE Full Text [26]. ↵ Smith , M.L. , Seo , Y.R .: p53 regulation of dna excision repair pathways . Mutagenesis 17 ( 2 ), 149 – 156 ( 2002 ) OpenUrl CrossRef PubMed Web of Science [27]. ↵ Hubbi , M.E. , Semenza , G.L .: Regulation of cell proliferation by hypoxia-inducible factors . American Journal of Physiology-Cell Physiology 309 ( 12 ), 775 – 782 ( 2015 ) OpenUrl [28]. El-Deiry , W.S. , Harper , J.W. , O’Connor , P.M. , Velculescu , V.E. , Canman , C.E. , Jackman , J. , Pietenpol , J.A. , Burrell , M. , Hill , D.E. , Wang , Y. , et al : Waf1/cip1 is induced in p53-mediated g1 arrest and apoptosis . Cancer research 54 ( 5 ), 1169 – 1174 ( 1994 ) OpenUrl Abstract / FREE Full Text [29]. ↵ Harper , J.V. , Brooks , G. : The mammalian cell cycle: an overview . Cell Cycle Control: Mechanisms and Protocols , 113 – 153 ( 2005 ) [30]. ↵ Greijer , A. , Wall , E .: The role of hypoxia inducible factor 1 (hif-1) in hypoxia induced apoptosis . Journal of clinical pathology 57 ( 10 ), 1009 – 1014 ( 2004 ) OpenUrl Abstract / FREE Full Text [31]. ↵ Yoon , M.-K. , Ha , J.-H. , Lee , M.-S. , Chi , S.-W .: Structure and apoptotic function of p73 . BMB reports 48 ( 2 ), 81 ( 2015 ) OpenUrl PubMed [32]. ↵ Xiong , L. , Garfinkel , A. : A common pathway to cancer: Oncogenic mutations abolish p53 oscillations . Progress in Biophysics and Molecular Biology 174 , 28 – 40 ( 2022 ) doi: 10.1016/j.pbiomolbio.2022.06.002 OpenUrl CrossRef PubMed [33]. Shin , Y.-J. , Hencey , B. , Lipkin , S.M. , Shen , X .: Frequency domain analysis reveals external periodic fluctuations can generate sustained p53 oscillation . PLoS One 6 ( 7 ), 22852 ( 2011 ) OpenUrl [34]. Hunziker , A. , Jensen , M.H. , Krishna , S .: Stress-specific response of the p53-mdm2 feedback loop . BMC systems biology 4 , 1 – 8 ( 2010 ) OpenUrl PubMed [35]. Ciliberto , A. , Novák, B., Tyson, J.J.: Steady states and oscillations in the p53/mdm2 network . Cell cycle 4 ( 3 ), 488 – 493 ( 2005 ) OpenUrl CrossRef PubMed Web of Science [36]. Lev Bar-Or , R. , Maya , R. , Segel , L.A. , Alon , U. , Levine , A.J. , Oren , M.: Generation of oscillations by the p53-mdm2 feedback loop: a theoretical and experimental study . Proceedings of the National Academy of Sciences 97 ( 21 ), 11250 – 11255 ( 2000 ) OpenUrl Abstract / FREE Full Text [37]. Xiong , L. , Garfinkel , A. : A common pathway to cancer: Oncogenic mutations abolish p53 oscillations . Progress in Biophysics and Molecular Biology 174 , 28 – 40 ( 2022 ) OpenUrl CrossRef PubMed [38]. ↵ Geva-Zatorsky , N. , Rosenfeld , N. , Itzkovitz , S. , Milo , R. , Sigal , A. , Dekel , E. , Yarnitzky , T. , Liron , Y. , Polak , P. , Lahav , G. , et al : Oscillations and variability in the p53 system . Molecular systems biology 2 ( 1 ), 2006 – 0033 ( 2006 ) OpenUrl CrossRef PubMed [39]. ↵ Luo , L. , Liu , H. , Yan , F .: Dynamic behavior of p53-mdm2-wip1 gene regulatory network under the influence of time delay and noise . Mathematical Biosciences and Engineering 20 ( 2 ), 2321 – 2347 ( 2023 ) OpenUrl [40]. Proctor , C.J. , Gray , D.A .: Explaining oscillations and variability in the p53-mdm2 system . BMC systems biology 2 , 1 – 20 ( 2008 ) OpenUrl PubMed [41]. Suarez , O.J. , Vega , C.J. , Sanchez , E.N. , Gonźalez-Santiago , A.E. , Rodŕıguez-Jorge , O. , Alanis , A.Y. , Chen , G. , Hernandez-Vargas , E.A. : Pinning control for the p53-mdm2 network dynamics regulated by p14arf . Frontiers in physiology 11 , 492386 ( 2020 ) OpenUrl [42]. Yu , H. , Liu , Q. , Bi , Y .: Ĺevy noise-induced phase transition in p53 gene regu-latory network near bifurcation points. Chaos , Solitons & Fractals 166 , 112885 ( 2023 ) OpenUrl [43]. ↵ Kim , D.-H. , Rho , K. , Kim , S .: A theoretical model for p53 dynamics: identifying optimal therapeutic strategy for its activation and stabilization . Cell Cycle 8 ( 22 ), 3707 – 3716 ( 2009 ) OpenUrl CrossRef PubMed [44]. ↵ Kim , E. , Shohet , J. : Targeted molecular therapy for neuroblastoma: the ARF/MDM2/p53 axis . Oxford University Press ( 2009 ) [45]. ↵ Agrawal , R.S. , Srikant , R. : R. fast algorithms for mining association rules . In: Proceedings of the 20th International Conference on Very Large Data Bases , VLDB , pp. 487 – 499 ( 1994 ) [46]. ↵ Crick , F .: Central dogma of molecular biology . Nature 227 ( 5258 ), 561 – 563 ( 1970 ) OpenUrl CrossRef PubMed Web of Science [47]. ↵ Chen , D. , Shan , J. , Zhu , W.-G. , Qin , J. , Gu , W .: Transcription-independent arf regulation in oncogenic stress-mediated p53 responses . Nature 464 ( 7288 ), 624 – 627 ( 2010 ) OpenUrl CrossRef PubMed Web of Science [48]. Hu , W. , Feng , Z. , Levine , A.J .: The regulation of multiple p53 stress responses is mediated through mdm2 . Genes & cancer 3 ( 3-4 ), 199 – 208 ( 2012 ) OpenUrl CrossRef PubMed [49]. ↵ Zindy , F. , Eischen , C.M. , Randle , D.H. , Kamijo , T. , Cleveland , J.L. , Sherr , C.J. , Roussel , M.F .: Myc signaling via the arf tumor suppressor regulates p53-dependent apoptosis and immortalization . Genes & development 12 ( 15 ), 2424 – 2433 ( 1998 ) OpenUrl Abstract / FREE Full Text [50]. ↵ Ingalls , B. : Mathematical Modeling in Systems Biology , ( 2013 ) [51]. ↵ Hausser , J. , Mayo , A. , Keren , L. , Alon , U .: Central dogma rates and the trade-off between precision and economy in gene expression . Nature Communications 10 ( 1 ) ( 2019 ) doi: 10.1038/s41467-018-07391-8 OpenUrl CrossRef PubMed [52]. ↵ Harenza , J.L. , Diamond , M. , Adams , R. , Song , M. , Davidson , H. , Hart , L. , Dent , M. , Fortina , P. , Reynolds , C. , Maris , J. : Corrigendum: Transcriptomic profiling of 39 commonly-used neuroblastoma cell lines . Scientific Data 4 , 170183 ( 2017 ) doi: 10.1038/sdata.2017.183 OpenUrl CrossRef PubMed [53]. ↵ Rodriguez , C.M. , Chun , S.Y. , Mills , R.E. , Todd , P.K .: Translation of upstream open reading frames in a model of neuronal differentiation . BMC Genomics 20 ( 1 ) ( 2019 ) doi: 10.1186/s12864-019-5775-1 OpenUrl CrossRef PubMed [54]. Lorenzin , F. , Benary , U. , Baluapuri , A. , Walz , S. , Jung , L.A. , Eyss , B. , Kisker , C. , Wolf , J. , Eilers , M. , Wolf , E .: Different promoter affinities account for specificity in myc-dependent gene regulation . eLife 5 ( JULY ) ( 2016 ) doi: 10.7554/eLife.15161 OpenUrl CrossRef PubMed [55]. Petty , T.J. , Emamzadah , S. , Costantino , L. , Petkova , I. , Stavridi , E.S. , Saven , J.G. , Vauthey , E. , Halazonetis , T.D .: An induced fit mechanism regulates p53 dna binding kinetics to confer sequence specificity . EMBO Journal 30 ( 11 ), 2167 – 2176 ( 2011 ) doi: 10.1038/emboj.2011.127 OpenUrl CrossRef PubMed Web of Science [56]. Zhou , S. , Gu , L. , He , J. , Zhang , H. , Zhou , M .: Mdm2 regulates vascular endothe-lial growth factor mrna stabilization in hypoxia . Molecular and Cellular Biology 31 ( 24 ), 4928 – 4937 ( 2011 ) doi: 10.1128/MCB.06085-11 OpenUrl Abstract / FREE Full Text [57]. Sang , P. , Shi , Y. , Lu , J. , Chen , L. , Yang , L. , Borcherds , W. , Abdulkadir , S. , Li , Q. , Daughdrill , G. , Chen , J. , Cai , J .: A-helix-mimicking sulfono-g-aapeptide inhibitors for p53-mdm2/mdmx protein–protein interactions . Journal of Medicinal Chemistry ( 2020 ) doi: 10.1021/acs.jmedchem.9b00993 OpenUrl CrossRef [58]. Shan , B. , Li , D.-W. , Brüschweiler-Li , L. , Brüschweiler , R. : Competitive binding between dynamic p53 transactivation subdomains to human mdm2 protein: Implications for regulating the p53mdm2 / mdmx interaction . Journal of Biologi-cal Chemistry 287 ( 36 ), 30376 – 30384 ( 2012 ) doi: 10.1074/jbc.M112.369793 OpenUrl Abstract / FREE Full Text [59]. Cai , X. , Yuan , Z.-M .: Stochastic modeling and simulation of the p53-mdm2/mdmx loop . Journal of computational biology : a journal of computational molecular cell biology 16 , 917 – 33 ( 2009 ) doi: 10.1089/cmb.2008.0231 OpenUrl CrossRef PubMed [60]. ↵ Bates , S. , Phillips , A.C. , Clark , P.A. , Stott , F. , Peters , G. , Ludwig , R.L. , Vous-den , K.H .: p14arf links the tumour suppressors rb and p53 . Nature 395 ( 6698 ), 124 – 125 ( 1998 ) OpenUrl CrossRef PubMed Web of Science [61]. Stanchina , E. , McCurrach , M.E. , Zindy , F. , Shieh , S.-Y. , Ferbeyre , G. , Samuelson , A.V. , Prives , C. , Roussel , M.F. , Sherr , C.J. , Lowe , S.W .: E1a signaling to p53 involves the p19arftumor suppressor . Genes & development 12 ( 15 ), 2434 – 2442 ( 1998 ) OpenUrl Abstract / FREE Full Text [62]. ↵ Aslanian , A. , Iaquinta , P.J. , Verona , R. , Lees , J.A .: Repression of the arf tumor suppressor by e2f3 is required for normal cell cycle kinetics . Genes & development 18 ( 12 ), 1413 – 1422 ( 2004 ) OpenUrl Abstract / FREE Full Text [63]. ↵ Palmero , I. , Pantoja , C. , Serrano , M .: p19arf links the tumour suppressor p53 to ras . Nature 395 ( 6698 ), 125 – 126 ( 1998 ) OpenUrl CrossRef PubMed Web of Science [64]. ↵ Berkovich , E. , Lamed , Y. , Ginsberg , D .: E2f and ras synergize in transcription-ally activating p14arf expression . Cell Cycle 2 ( 2 ), 127 – 134 ( 2003 ) OpenUrl CrossRef PubMed Web of Science [65]. ↵ Damalas , A. , Velimezi , G. , Kalaitzakis , A. , Liontos , M. , Papavassiliou , A.G. , Gorgoulis , V. , Angelidis , C .: Loss of p14arf confers resistance to heat shock-and oxidative stress-mediated cell death by upregulating β -catenin . International journal of cancer 128 ( 8 ), 1989 – 1995 ( 2011 ) OpenUrl PubMed [66]. ↵ Hu , W. , Feng , Z. , Levine , A.J .: The regulation of multiple p53 stress responses is mediated through mdm2 . Genes & cancer 3 ( 3-4 ), 199 – 208 ( 2012 ) OpenUrl CrossRef PubMed [67]. Ofir-Rosenfeld , Y. , Boggs , K. , Michael , D. , Kastan , M.B. , Oren , M .: Mdm2 regulates p53 mrna translation through inhibitory interactions with ribosomal protein l26 . Molecular cell 32 ( 2 ), 180 – 189 ( 2008 ) OpenUrl CrossRef PubMed Web of Science [68]. ↵ Takagi , M. , Absalon , M.J. , McLure , K.G. , Kastan , M.B .: Regulation of p53 translation and induction after dna damage by ribosomal protein l26 and nucleolin . Cell 123 ( 1 ), 49 – 63 ( 2005 ) OpenUrl CrossRef PubMed Web of Science [69]. ↵ Candeias , M.M. , Malbert-Colas , L. , Powell , D.J. , Daskalogianni , C. , Maslon , M.M. , Naski , N. , Bourougaa , K. , Calvo , F. , Fåhraeus , R. : P53 mrna controls p53 activity by managing mdm2 functions . Nature cell biology 10 ( 9 ), 1098 – 1105 ( 2008 ) OpenUrl CrossRef PubMed Web of Science [70]. ↵ Yin , Y. , Stephen , C.W. , Luciani , M.G. , Fåhraeus , R. : p53 stability and activity is regulated by mdm2-mediated induction of alternative p53 translation products . Nature cell biology 4 ( 6 ), 462 – 467 ( 2002 ) OpenUrl CrossRef PubMed Web of Science [71]. ↵ Meek , D.W .: The p53 response to dna damage . DNA repair 3 ( 8-9 ), 1049 – 1056 ( 2004 ) OpenUrl CrossRef PubMed [72]. ↵ Brooks , C.L. , Gu , W .: Ubiquitination, phosphorylation and acetylation: the molecular basis for p53 regulation . Current opinion in cell biology 15 ( 2 ), 164 – 171 ( 2003 ) OpenUrl CrossRef PubMed Web of Science [73]. ↵ Khosravi , R. , Maya , R. , Gottlieb , T. , Oren , M. , Shiloh , Y. , Shkedy , D .: Rapid atm-dependent phosphorylation of mdm2 precedes p53 accumulation in response to dna damage . Proceedings of the National Academy of Sciences 96 ( 26 ), 14973 – 14977 ( 1999 ) OpenUrl Abstract / FREE Full Text [74]. ↵ Brazhnik , P. , Kohn , K.W .: Hausp-regulated switch from auto-to p53 ubiquiti-nation by mdm2 (in silico discovery) . Mathematical Biosciences 210 ( 1 ), 60 – 77 ( 2007 ) OpenUrl PubMed [75]. ↵ Fang , S. , Jensen , J.P. , Ludwig , R.L. , Vousden , K.H. , Weissman , A.M .: Mdm2 is a ring finger-dependent ubiquitin protein ligase for itself and p53 . Journal of Biological Chemistry 275 ( 12 ), 8945 – 8951 ( 2000 ) OpenUrl Abstract / FREE Full Text [76]. ↵ Amente , S. , Gargano , B. , Diolaiti , D. , Della Valle , G. , Lania , L. , Majello , B .: p14arf interacts with n-myc and inhibits its transcriptional activity . FEBS letters 581 ( 5 ), 821 – 825 ( 2007 ) OpenUrl CrossRef PubMed [77]. ↵ Lohrum , M.A. , Ashcroft , M. , Kubbutat , M.H. , Vousden , K.H .: Contribution of two independent mdm2-binding domains in p14arf to p53 stabilization . Current Biology 10 ( 9 ), 539 – 542 ( 2000 ) OpenUrl CrossRef PubMed Web of Science [78]. ↵ Zhang , Y. , Xiong , Y. , Yarbrough , W.G. : Arf promotes mdm2 degradation and stabilizes p53: Arf-ink4a locus deletion impairs both the rb and p53 tumor suppression pathways . Cell 92 ( 6 ), 725 – 734 ( 1998 ) OpenUrl CrossRef PubMed Web of Science [79]. ↵ Honda , R. , Tanaka , H. , Yasuda , H .: Oncoprotein mdm2 is a ubiquitin ligase e3 for tumor suppressor p53 . FEBS letters 420 ( 1 ), 25 – 27 ( 1997 ) OpenUrl CrossRef PubMed Web of Science [80]. ↵ Shi , D. , Gu , W .: Dual roles of mdm2 in the regulation of p53: ubiquitination dependent and ubiquitination independent mechanisms of mdm2 repression of p53 activity . Genes & cancer 3 ( 3-4 ), 240 – 248 ( 2012 ) OpenUrl CrossRef PubMed [81]. ↵ Kamijo , T. , Weber , J.D. , Zambetti , G. , Zindy , F. , Roussel , M.F. , Sherr , C.J .: Functional and physical interactions of the arf tumor suppressor with p53 and mdm2 . Proceedings of the National Academy of Sciences 95 ( 14 ), 8292 – 8297 ( 1998 ) OpenUrl Abstract / FREE Full Text [82]. ↵ Kussie , P.H. , Gorina , S. , Marechal , V. , Elenbaas , B. , Moreau , J. , Levine , A.J. , Pavletich , N.P .: Structure of the mdm2 oncoprotein bound to the p53 tumor suppressor transactivation domain . Science 274 ( 5289 ), 948 – 953 ( 1996 ) OpenUrl Abstract / FREE Full Text [83]. Oliner , J.D. , Pietenpol , J.A. , Thiagalingam , S. , Gyuris , J. , Kinzler , K.W. , Vogelstein , B .: Oncoprotein mdm2 conceals the activation domain of tumour suppressor p53 . Nature 362 ( 6423 ), 857 – 860 ( 1993 ) OpenUrl CrossRef PubMed Web of Science [84]. ↵ Chen , J. , Marechal , V. , Levine , A.J .: Mapping of the p53 and mdm-2 interaction domains . Molecular and cellular biology ( 1993 ) [85]. ↵ Seto , E. , Usheva , A. , Zambetti , G.P. , Momand , J. , Horikoshi , N. , Weinmann , R. , Levine , A.J. , Shenk , T .: Wild-type p53 binds to the tata-binding protein and represses transcription . Proceedings of the National Academy of Sciences 89 ( 24 ), 12028 – 12032 ( 1992 ) OpenUrl Abstract / FREE Full Text [86]. ↵ Sullivan , K.D. , Galbraith , M.D. , Andrysik , Z. , Espinosa , J.M .: Mechanisms of transcriptional regulation by p53 . Cell Death & Differentiation 25 ( 1 ), 133 – 143 ( 2018 ) OpenUrl CrossRef PubMed [87]. ↵ Stott , F.J. , Bates , S. , James , M.C. , McConnell , B.B. , Starborg , M. , Brookes , S. , Palmero , I. , Ryan , K. , Hara , E. , Vousden , K.H. , et al : The alternative product from the human cdkn2a locus, p14arf, participates in a regulatory feedback loop with p53 and mdm2 . The EMBO journal 17 ( 17 ), 5001 – 5014 ( 1998 ) OpenUrl Abstract / FREE Full Text [88]. ↵ Wang , S. , El-Deiry , W.S .: p73 or p53 directly regulates human p53 transcription to maintain cell cycle checkpoints . Cancer research 66 ( 14 ), 6982 – 6989 ( 2006 ) OpenUrl Abstract / FREE Full Text [89]. ↵ Saldaña-Meyer , R. , Recillas-Targa , F. : Transcriptional and epigenetic regulation of the p53 tumor suppressor gene . Epigenetics 6 ( 9 ), 1068 – 1077 ( 2011 ) OpenUrl CrossRef PubMed Web of Science [90]. ↵ Slack , A. , Chen , Z. , Tonelli , R. , Pule , M. , Hunt , L. , Pession , A. , Shohet , J.M .: The p53 regulatory gene mdm2 is a direct transcriptional target of mycn in neuroblastoma . Proceedings of the National Academy of Sciences 102 ( 3 ), 731 – 736 ( 2005 ) OpenUrl Abstract / FREE Full Text [91]. ↵ Barak , Y. , Juven , T. , Haffner , R. , Oren , M .: mdm2 expression is induced by wild type p53 activity . The EMBO journal 12 ( 2 ), 461 – 468 ( 1993 ) OpenUrl CrossRef PubMed Web of Science [92]. ↵ Zeng , Y. , Kotake , Y. , Pei , X.-H. , Smith , M.D. , Xiong , Y .: p53 binds to and is required for the repression of arf tumor suppressor by hdac and polycomb . Cancer research 71 ( 7 ), 2781 – 2792 ( 2011 ) OpenUrl Abstract / FREE Full Text [93]. ↵ Robertson , K.D. , Jones , P.A .: The human arf cell cycle regulatory gene promoter is a cpg island which can be silenced by dna methylation and down-regulated by wild-type p53 . Molecular and cellular biology 18 ( 11 ), 6457 – 6473 ( 1998 ) OpenUrl Abstract / FREE Full Text [94]. ↵ Gu , L. , Zhang , H. , He , J. , Li , J. , Huang , M. , Zhou , M .: Mdm2 regulates mycn mrna stabilization and translation in human neuroblastoma cells . Oncogene 31 ( 11 ), 1342 – 1353 ( 2012 ) OpenUrl CrossRef PubMed [95]. ↵ Fredrick , K. , Ibba , M .: How the sequence of a gene can tune its translation . Cell 141 ( 2 ), 227 – 229 ( 2010 ) OpenUrl CrossRef PubMed Web of Science [96]. ↵ Wilkie , G.S. , Dickson , K.S. , Gray , N.K .: Regulation of mrna translation by 5 ′ -and 3 ′ -utr-binding factors . Trends in biochemical sciences 28 ( 4 ), 182 – 188 ( 2003 ) OpenUrl CrossRef PubMed Web of Science [97]. ↵ Han , J. , Pei , J. , Yin , Y .: Mining frequent patterns without candidate generation . SIGMOD Record (ACM Special Interest Group on Management of Data ) 29 ( 2 ), 1 – 12 ( 2000 ) doi: 10.1145/335191.335372 OpenUrl CrossRef [98]. ↵ Hahsler , M. , Gruen , B. , Hornik , K .: arules – A computational environment for mining association rules and frequent item sets . Journal of Statistical Software 14 ( 15 ), 1 – 25 ( 2005 ) doi: 10.18637/jss.v014.i15 OpenUrl CrossRef [99]. ↵ Campbell , K. , Gastier-Foster , J. , Mann , M. , Naranjo , A. , Ryn , C. , Bagatell , R. , Matthay , K. , London , W. , Irwin , M. , Shimada , H. , Granger , M. , Hogarty , M. , Park , Y. , DuBois , S. : Association of mycn copy number with clinical features, tumor biology, and outcomes in neuroblastoma: A report from the children’s oncology group: Mycn copy number in neuroblastoma . Cancer 123 ( 2017 ) doi: 10.1002/cncr.30873 OpenUrl CrossRef PubMed [100]. ↵ Karakostis , K. , Malbert-Colas , L. , Thermou , A. , Vojtesek , B. , Fåhraeus , R. : The dna damage sensor atm kinase interacts with the p53 mrna and guides the dna damage response pathway . Molecular cancer 23 ( 1 ), 21 ( 2024 ) OpenUrl PubMed [101]. ↵ De Brouwer , S. , De Preter , K. , Kumps , C. , Zabrocki , P. , Porcu , M. , Westerhout , E.M. , Lakeman , A. , Vandesompele , J. , Hoebeeck , J. , Van Maerken , T. , et al : Meta-analysis of neuroblastomas reveals a skewed alk mutation spectrum in tumors with mycn amplification . Clinical cancer research 16 ( 17 ), 4353 – 4362 ( 2010 ) OpenUrl Abstract / FREE Full Text [102]. ↵ Zhu , S. , Lee , J.-S. , Guo , F. , Shin , J. , Perez-Atayde , A.R. , Kutok , J.L. , Rodig , S.J. , Neuberg , D.S. , Helman , D. , Feng , H. , et al : Activated alk collaborates with mycn in neuroblastoma pathogenesis . Cancer cell 21 ( 3 ), 362 – 373 ( 2012 ) OpenUrl CrossRef PubMed Web of Science [103]. ↵ Müller , M. , Trunk , K. , Fleischhauer , D. , Büchel , G. : Mycn in neuroblastoma: the kings’ new clothes and drugs . EJC Paediatric Oncology , 100182 ( 2024 ) [104]. ↵ Zhang , W. , Liu , B. , Wu , W. , Li , L. , Broom , B.M. , Basourakos , S.P. , Korentzelos , D. , Luan , Y. , Wang , J. , Yang , G. , et al : Targeting the mycn–parp–dna damage response pathway in neuroendocrine prostate cancer . Clinical Cancer Research 24 ( 3 ), 696 – 707 ( 2018 ) OpenUrl Abstract / FREE Full Text [105]. ↵ Boon , K. , Caron , H.N. , Van Asperen , R. , Valentijn , L. , Hermus , M.-C. , Van Sluis , P. , Roobeek , I. , Weis , I. , Voute , P. , Schwab , M. , et al : N-myc enhances the expression of a large set of genes functioning in ribosome biogenesis and protein synthesis . The EMBO journal ( 2001 ) [106]. ↵ Kuzuoglu-Ozturk , D. , Aksoy , O. , Schmidt , C. , Lea , R. , Larson , J.D. , Phelps , R.R. , Nasholm , N. , Holt , M. , Contreras , A. , Huang , M. , et al : N-myc–mediated translation control is a therapeutic vulnerability in medulloblastoma . Cancer research 83 ( 1 ), 130 – 140 ( 2023 ) OpenUrl CrossRef PubMed [107]. ↵ Xia , Y. , Ye , B. , Ding , J. , Yu , Y. , Alptekin , A. , Thangaraju , M. , Prasad , P.D. , Ding , Z.-C. , Park , E.J. , Choi , J.-H. , et al : Metabolic reprogramming by mycn confers dependence on the serine-glycine-one-carbon biosynthetic pathway . Cancer research 79 ( 15 ), 3837 – 3850 ( 2019 ) OpenUrl CrossRef PubMed [108]. ↵ Dai , C. : The heat-shock, or hsf1-mediated proteotoxic stress, response in cancer: from proteomic stability to oncogenesis . Philosophical Transactions of the Royal Society B: Biological Sciences 373 ( 1738 ), 20160525 ( 2018 ) OpenUrl CrossRef PubMed [109]. ↵ Pakos-Zebrucka , K. , Koryga , I. , Mnich , K. , Ljujic , M. , Samali , A. , Gorman , A.M .: The integrated stress response . EMBO reports 17 ( 10 ), 1374 – 1395 ( 2016 ) OpenUrl Abstract / FREE Full Text [110]. ↵ Hamperl , S. , Bocek , M.J. , Saldivar , J.C. , Swigut , T. , Cimprich , K.A .: Transcription-replication conflict orientation modulates r-loop levels and acti-vates distinct dna damage responses . Cell 170 ( 4 ), 774 – 786 ( 2017 ) OpenUrl CrossRef PubMed [111]. ↵ Pagliarini , R. , Shao , W. , Sellers , W.R .: Oncogene addiction: pathways of thera-peutic response, resistance, and road maps toward a cure . EMBO reports 16 ( 3 ), 280 – 296 ( 2015 ) OpenUrl Abstract / FREE Full Text [112]. ↵ Hu , H. , Du , L. , Nagabayashi , G. , Seeger , R.C. , Gatti , R.A .: Atm is down-regulated by n-myc–regulated microrna-421 . Proceedings of the National Academy of Sciences 107 ( 4 ), 1506 – 1511 ( 2010 ) OpenUrl Abstract / FREE Full Text [113]. ↵ Schönherr , C. , Ruuth , K. , Kamaraj , S. , Wang , C.-L. , Yang , H.-L. , Combaret , V. , Djos , A. , Martinsson , T. , Christensen , J.G. , Palmer , R.H. , et al : Anaplas-tic lymphoma kinase (alk) regulates initiation of transcription of mycn in neuroblastoma cells . Oncogene 31 ( 50 ), 5193 – 5200 ( 2012 ) OpenUrl CrossRef PubMed Web of Science [114]. ↵ Wertheim , K.Y. , Chisholm , R. , Richmond , P. , Walker , D .: Multicellular model of neuroblastoma proposes unconventional therapy based on multiple roles of p53 . PLOS Computational Biology 20 ( 12 ), 1012648 ( 2024 ) OpenUrl [115]. ↵ Oren , M. , Prives , C .: p53: A tale of complexity and context . Cell 187 ( 7 ), 1569 – 1573 ( 2024 ) OpenUrl CrossRef PubMed [116]. ↵ Borau , C. , Wertheim , K. , Hervas-Raluy , S. , Sainz-DeMena , D. , Walker , D. , Chisholm , R. , Richmond , P. , Varella , V. , Viceconti , M. , Montero , A. , et al : A multiscale orchestrated computational framework to reveal emergent phenomena in neuroblastoma . Computer Methods and Programs in Biomedicine , 107742 ( 2023 ) [117]. ↵ McLure , K.G. , Lee , P.W .: How p53 binds dna as a tetramer . The EMBO journal ( 1998 ) [118]. ↵ Weinberg , R.L. , Veprintsev , D.B. , Fersht , A.R .: Cooperative binding of tetrameric p53 to dna . Journal of molecular biology 341 ( 5 ), 1145 – 1159 ( 2004 ) OpenUrl CrossRef PubMed Web of Science [119]. ↵ Brandt , T. , Petrovich , M. , Joerger , A.C. , Veprintsev , D.B .: Conservation of dna-binding specificity and oligomerisation properties within the p53 family . BMC genomics 10 , 1 – 14 ( 2009 ) OpenUrl CrossRef PubMed [120]. ↵ Huang , R. , Cheung , N.-K. , Vider , L. , Cheung , I. , Gerald , W. , Tickoo , S. , Holland , E. , Blasberg , R .: Mycn and myc regulate tumor proliferation and tumorigenesis directly through bmi1 in human neuroblastomas . The FASEB journal 25 , 4138 – 49 ( 2011 ) doi: 10.1096/fj.11-185033 OpenUrl CrossRef PubMed Web of Science View the discussion thread. Back to top Previous Next Posted November 14, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following A novel prognostic metric resolves the MYCN enigma in silico and points to a biosynthetic regulatory shift driven by MYCN amplification in neuroblastoma Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share A novel prognostic metric resolves the MYCN enigma in silico and points to a biosynthetic regulatory shift driven by MYCN amplification in neuroblastoma Matteo Italia , Melody Parker , Rory Deignan , Fabio Dercole , Dawn Walker , Kenneth Y. Wertheim bioRxiv 2025.11.13.688307; doi: https://doi.org/10.1101/2025.11.13.688307 Share This Article: Copy Citation Tools A novel prognostic metric resolves the MYCN enigma in silico and points to a biosynthetic regulatory shift driven by MYCN amplification in neuroblastoma Matteo Italia , Melody Parker , Rory Deignan , Fabio Dercole , Dawn Walker , Kenneth Y. Wertheim bioRxiv 2025.11.13.688307; doi: https://doi.org/10.1101/2025.11.13.688307 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Systems Biology Subject Areas All Articles Animal Behavior and Cognition (7617) Biochemistry (17633) Bioengineering (13856) Bioinformatics (41841) Biophysics (21399) Cancer Biology (18529) Cell Biology (25422) Clinical Trials (138) Developmental Biology (13352) Ecology (19860) Epidemiology (2067) Evolutionary Biology (24281) Genetics (15582) Genomics (22461) Immunology (17700) Microbiology (40293) Molecular Biology (17140) Neuroscience (88413) Paleontology (666) Pathology (2823) Pharmacology and Toxicology (4813) Physiology (7632) Plant Biology (15107) Scientific Communication and Education (2042) Synthetic Biology (4284) Systems Biology (9808) Zoology (2267)