Multicellular model of neuroblastoma proposes unconventional therapy based on multiple roles of p53

Abstract

Neuroblastoma is the most common extra-cranial solid tumour in children. Over half of all high-risk cases are expected to succumb to the disease even after chemotherapy, surgery, and immunotherapy. Although the importance of MYCN amplification in this disease is indisputable, the mechanistic details remain enigmatic. Here, we present a multicellular model of neuroblastoma comprising a continuous automaton, discrete cell agents, and a centre-based mechanical model, as well as the simulation results we obtained with it. The continuous automaton represents the tumour microenvironment as a grid-like structure, where each voxel is associated with continuous variables such as the oxygen level therein. Each discrete cell agent is defined by several attributes, including its cell cycle position, mutations, gene expression pattern, and more with behaviours such as cell cycling and cell death being stochastically dependent on these attributes. The centre-based mechanical model represents the properties of these agents as physical objects, describing how they repel each other as soft spheres. By implementing a stochastic simulation algorithm on modern GPUs, we simulated the dynamics of over one million neuroblastoma cells over a period of months. Specifically, we set up 1200 heterogeneous tumours and tracked the MYCN-amplified clone's dynamics in each, revealed the conditions that favour its growth, and tested its responses to 5000 drug combinations. Our results are in agreement with those reported in the literature and add new insights into how the MYCN-amplified clone's reproductive advantage in a tumour, its gene expression profile, the tumour's other clones (with different mutations), and the tumour's microenvironment are inter-related. Based on the results, we formulated a hypothesis, which argues that there are two distinct populations of neuroblastoma cells in the tumour; the p53 protein is pro-survival in one and pro-apoptosis in the other. It follows that alternating between inhibiting MDM2 to restore p53 activity and inhibiting ARF to attenuate p53 activity is a promising, if unorthodox, therapeutic strategy. The multicellular model has the advantages of modularity, high resolution, and scalability, making it a potential foundation for creating digital twins of neuroblastoma patients.

Fig 1. Model structure.

(A) There are two agent populations: neuroblastoma and Schwann cell agents. They communicate using intercellular juxtacrine (contact-dependent) and paracrine (diffusive) signals, such as apoptotic signals. (B) The tumour’s microenvironment is represented by a continuous automaton, which comprises voxels populated by the cell agents. (C) The intracellular mechanisms mediating cell cycling and cell death in response to the microenvironment and intercellular interactions are represented by Bernoulli trials with conditional probabilities (full details provided in S1 Text). (D) Cell-cell overlap is resolved by a centre-based mechanical model (also known as an off-lattice soft sphere model), wherein cells repel each other like soft spheres. (E) A neuroblastoma cell agent belongs to one of four clones: wild-type (green), MYCN-amplified (magenta), TERT-rearranged (red), and ATRX-inactivated (blue). Each clone comprises six smaller clones (or subclones) with different combinations of the following: ALK activation/amplification, other mutations activating the MAPK/RAS signalling pathway, and p53 inactivation. At the highest resolution, a neuroblastoma cell agent’s clonal identity is defined in terms of the subclone to which it belongs. For the sake of simplicity, it has a clone ID.

Fig 2. Sequence of events in one realisation of the stochastic process.

At the start of each time step (one hour), any neuroblastoma cell agents are sequentially evaluated. Each agent senses and responds to its microenvironment by updating its attributes. Then, it attempts to progress in the cell cycle. If it is apoptotic or necrotic, it may get removed from the system. Otherwise, it may divide into two daughter agents. Schwann cell agents are then updated similarly. All agents are then considered collectively to minimise the total overlap in the system, thus updating their spatial coordinates. This involves solving a centre-based mechanical model by applying Euler’s method iteratively. Each step finishes with an update of the continuous automaton. This accounting step does not change the agents’ attributes; it simply records their latest spatial distributions, any changes in the supply rate of oxygen and its abundance, and matrix production.

Fig 3. Calibration pipeline.

We generated 3000 near-random combinations of 20 unconstrained parameters by Latin hypercube sampling and systematically eliminated unrealistic sets of parameters in a series of calibration studies. In the first study, they were assessed by their ability to reproduce neuroblastoma’s in vitro growth kinetics [39]. In the second study, their ability to mimic its hypoxic response was evaluated [40]. The third study was designed based on the regulatory dynamics between neuroblastoma and Schwann cells, observed in vitro[41]. The fourth study selected the parametric combinations that reproduced the clinical outcomes associated with different histological categories [42]. The last two studies were designed to reproduce the clinical outcomes associated with different mutations [20].

Fig 4. Violin plots illustrating the initial macroscopic properties of 1200 heterogeneous virtual tumours, namely each tumour’s dimensionless oxygen level (O2) and the fraction of its total cell population assigned to be Schwann cells (SC fraction).

(A) Each violin plot contains 1200 data points (all cases), one for each virtual tumour. (B) Each violin plot contains 45 data points (regressed cases only), one for each regressed tumour. (C) Each violin plot contains 1155 data points (progressing cases only), one for each progressing tumour.

Fig 5. Violin plots illustrating the initial fractional compositions of the wild-type (WT, green) and MYCN-magnified (MA, magenta) clones in 1200 heterogeneous virtual tumours.

The fractional composition of a clone or subclone is simply the number of living neuroblastoma cell agents inside that clone or subclone divided by the total number of living neuroblastoma cell agents in the virtual tumour containing the clone or subclone. (A) and (D) Each violin plot contains 1200 data points (all 1200 cases), one for each virtual tumour. (B) and (E) Each violin plot contains 45 data points (regressed cases only), one for each regressed tumour. (C) and (F) Each violin plot contains 1155 data points (progressing cases only), one for each progressing tumour.

Fig 6. Violin plots illustrating the initial fractional compositions of the TERT-rearranged (TR, red) and ATRX-inactivated (AI, blue) clones in 1200 heterogeneous virtual tumours.

The fractional composition of a clone or subclone is simply the number of living neuroblastoma cell agents inside that clone or subclone divided by the total number of living neuroblastoma cell agents in the virtual tumour containing the clone or subclone. (A) and (D) Each violin plot contains 1200 data points (all 1200 cases), one for each virtual tumour. (B) and (E) Each violin plot contains 45 data points (regressed cases only), one for each regressed tumour. (C) and (F) Each violin plot contains 1155 data points (progressing cases only), one for each progressing tumour.

Fig 7. Two sets of 95% confidence intervals for 45 regressed and 1155 progressing virtual tumours’ medium initial clonal fractional compositions.

The six columns on the left record the 45 regressed virtual tumours’ mediums, the lower bounds of their confidence intervals (CI LB), and the upper bounds (CI UB). The six columns on the right record the same statistics for the 1155 progressing virtual tumours. Subclones 1, 2, 3, 4, 5, and 6 constitute the wild-type (green) clone. Subclones 7, 8, 9, 10, 11, and 12 constitute the MYCN-amplified (magenta) clone. Subclones 13, 14, 15, 16, 17, and 18 constitute the TERT-rearranged (red) clone. Subclones 19, 20, 21, 22, 23, and 24 constitute the ATRX-inactivated (blue) clone. At the highest resolution, a neuroblastoma cell agent’s clonal identity is defined in terms of the subclone to which it belongs. For the sake of simplicity, it has a clone ID. For each subclone, we resampled from its 1200 initial fractional compositions to build two confidence intervals. First, we resampled the sample 12000 times, taking 45 values with replacement and recording their medium on each occasion. Using the 12000 mediums, we built a 95% confidence interval for the subclone’s medium initial fractional composition. Second, we repeated the procedures with a resample size of 1155 values to build another 95% confidence interval.

Fig 8. Outcomes of clonal competition in 1155 progressing heterogeneous virtual tumours.

(A)–(D) The four scatter plots pertain to four groups of clones: wild-type (WT, green), MYCN-amplified (MA, magenta), TERT-rearranged (TR, red), and ATRX-inactivated (AI, blue). Each point on a scatter plot relates a clone’s initial fractional composition to its enrichment at the end of the corresponding simulation: its final fractional composition divided by its initial fractional composition. The four clones’ fractional compositions in each virtual tumour, both initial and final, always add up to one. (E)–(H) The violin plots present the final fractional compositions of the six subclones in each clone. The 24 fractional compositions of every virtual tumour always add up to one.

Fig 9. Sensitivity analysis evaluating 1000 MYCN-amplified (MA) clones with different gene expression profiles.

(A) The violin plots present how the expression levels of each gene are distributed over the profiles that allowed their corresponding MYCN-amplified (MA) clones to expand in the simulations. The genes are MYCN (M), the genes encoding the MAPK/RAS signalling pathway (MR1 and MR0), p53, p73, and HIF. MR1 denotes enhanced MAPK/RAS signalling in a MYCN-amplified neuroblastoma cell agent. MR0 denotes enhanced MAPK/RAS signalling in a neuroblastoma cell agent whose MYCN is not amplified. (B) Each bar quantifies the amount of variance in these gene expression profiles that can be explained by a particular principal component. (C) The lollipop chart presents the gene expression profile that parameterised the simulation that matched what is known [20] most closely. In this simulation, the three mutated clones dominated the wild-type clone and in each clone, the subclones with mutations in p53 and the genes encoding the MAPK/RAS signalling pathway dominated their peers. Furthermore, the number of living neuroblastoma cells in each clone increased during the simulation: the clones expanded. (D) The four time series confirm that, in the simulation parameterised by the gene expression profile presented in (C), all four clones (wild-type or WT, green; MYCN-amplified or MA, magenta; TERT-rearranged or TR, red; and ATRX-inactivated or AI, blue) expanded.

Fig 10. Results of an in silico drug trial for a virtual tumour containing a singular MYCN-amplified (MA) clone.

The violin plots from (A) to (E) illustrate how the inhibitory activity levels with respect to each drug target are distributed over the 26 effective drug combinations that led to regression in the corresponding simulations. MR denotes the MAPK/RAS signalling pathway, Bcl denotes Bcl-2 and Bcl-xl collectively, BB denotes Bak and Bax collectively, and telo denotes telomerase. (F) Six time series illustrating how the six subclones within the MA clone responded to an effective drug combination in the corresponding simulation.

Fig 11. Results of an in silico drug trial for a virtual tumour containing a singular MYCN-amplified (MA) clone.

(A) The outcome of a principal component analysis of a dataset relating to 305 effective and 297 ineffective drug combinations against the virtual tumour; each bar quantifies the amount of variance in the dataset that can be explained by a particular principal component. (B) to (D) The weights of the first principal component of the dataset. Each bar quantifies the importance of inhibiting a drug target, including MYCN; the MAPK/RAS signalling pathway (MR); JAB1; CHK1; ID2; IAP2; HIF; BNIP3; VEGF; p53; p73; p21; p27; Bcl-2 and Bcl-xl (Bcl); Bak and Bax (BB); CAS; CDS1; CDC25C; ALT; and telomerase (telo). (E) The scatter plot projects the entire dataset on its first two principal components and colours two clusters predicted by hierarchical clustering. (F) The scatter plot presents the final outcomes of the simulations associated with the 602 drug combinations in the dataset and colours the same two clusters.