Modelling genotypes in their physical microenvironment to predict single- and multi-cellular behaviour

A cell’s phenotype is the set of observable characteristics resulting from the interaction of the genotype with the surrounding environment, determining cell behaviour. Deciphering genotype-phenotype relationships has been crucial to understand normal and disease biology. Analysis of molecular pathways has provided an invaluable tool to such understanding; however, it has typically lacked a component describing the physical context, which is a key determinant of phenotype. In this study, we present a novel modelling framework that enables to study the link between genotype, signalling networks and cell behaviour in a 3D physical environment. To achieve this we bring together Agent Based Modelling, a powerful computational modelling technique, and gene networks. This combination allows biological hypotheses to be tested in a controlled stepwise fashion, and it lends itself naturally to model a heterogeneous population of cells acting and evolving in a dynamic microenvironment, which is needed to predict the evolution of complex multi-cellular dynamics. Importantly, this enables modelling co-occurring intrinsic perturbations, such as mutations, and extrinsic perturbations, such as nutrients availability, and their interactions. Using cancer as a model system, we illustrate the how this framework delivers a unique opportunity to identify determinants of single-cell behaviour, while uncovering emerging properties of multi-cellular growth. Contact Francesca M. Buffa, University of Oxford, francesca.buffa@imm.ox.ac.uk


Introduction
A comprehensive understanding of living organisms, including their development and the occurrence and progression of disease, requires systematic efforts into deciphering the link between the multitude of different genotypes and phenotypes, that is the set of observable characteristics resulting from the interaction of the genotype with the surrounding environment, determining cell morphology and behaviour.
Efforts to characterise such relationships have proliferated in recent years, thanks in part to the increased capability to efficiently collect the necessary data in an ever higher number of organisms and individuals.
Such efforts have spanned from cataloguing genetic variation in thousands of individuals [see e.g. (1)] and searching for genotypes-phenotypes associations, to silencing or inactivating thousands of genes in laboratory high-throughput screens to study their function [see e.g. (2)]. However, efficient instruments to achieve a comprehensive understanding of the causal nexus between a given genotype and the observed phenotype are still lacking, and this is particularly true when such a nexus is complex and multifactorial (3)(4)(5). One promising means to achieve such an understanding has been the characterisation and modelling of biological pathways (5).
Several of the key biological pathways regulating cellular function are increasingly understood, along with their dysregulation in multiple diseases (6)(7)(8)(9). However, much less is known about how these pathways interact and determine the behaviour of individual cells and multi-cellular systems. This has fuelled methodological development of efficient representations of such interactions in order to facilitate the study of the underlying potential mechanisms. A number of different approaches have been proposed, including modelling by differential equations, and network modelling methods, such as petri nets and logical networks (10)(11)(12). Such methods have produced encouraging results (13)(14)(15)(16), but it is becoming increasingly evident that modelling molecular pathways and signalling, or gene, networks in isolation, dissociated from the cellular context, does not reflect the crucial impact of the microenvironment in determining the phenotype (17)(18)(19). Additionally, as single-cell sequencing and imaging technologies are providing new in-depth information about the genotype and phenotype of single, or small groups of, cells (20), modelling approaches that enable consideration of cells both as independent entities and as a population, are becoming increasingly attractive.
To address the above needs, we have developed a novel computational framework which combines Agent-Based Modelling (ABM) and gene network modelling. A detailed description of this method is provided in the Methods section. Briefly, microC addresses the fundamental challenge of integrating genotype with phenotype data, while accounting for the effect of the physical microenvironment, a key determinant of the phenotype. The most innovative aspect of this framework is that it enables to build models of the genotype-phenotype relationship in a 3D spatially-aware microenvironment, including aspects such as signalling to and from the microenvironment, and signalling between cells. Importantly, the collective behaviour and evolution of cellular populations emerges from the properties and behaviour of individual cells, which in turn are governed by the underling dynamics of specific signalling networks, and the interactions with the surrounding cells and microenvironment. This permits prediction not only of the behaviour of individual cells, and the entire population of cells, but also study of the possible causative mechanisms of such behaviour.
In this paper, we introduce the microC framework and we illustrate the range of potential applications that it enables. Specifically, we perform perturbation experiments of increasing complexity, in which we monitor over time the three-dimensional growth and evolution of populations of pre-malignant and cancer cells. We illustrate how microC can be used to gradually increase the genetic and microenvironmental heterogeneity to monitor clonal competition, signalling to and from the microenvironment, and cell-cell interactions.

Rationale for the framework
Previous studies have modelled cells as computational agents and started to consider replacing ABM conditional statements, that drive cellular behaviour, with gene networks, that can represent more realistically the dynamical features of the intra-cellular system. In particular, small scale predefined logical networks have been used to model the cell-cycle arrest in a model of avascular spheroid growth (21), and to study cell differentiation in a hyper-sensitivity reaction model (22). However, these networks were considered as static entities, providing simple rules and not fully embedded in the ABM simulation. In another example, cells have been equipped with relatively more complex decision making models, represented as differential equation (23), but this approach was also not fully embedded in the ABM modelling, and is further limited by the requirement of prior knowledge for the many kinetic parameters to represent pathways accurately. microC fully exploits ABM technology. Specifically, the cellular behaviour is determined by a gene network, encapsulated within the cell, and by the interactions of the gene network with the surrounding microenvironment ( Figure 1). Both the cells and the inner networks are simulated using ABM. As a result, cellular behaviour may be studied both at the single-cell level, decoding the link between a cell genotype and its phenotypic realization in a given contextual environment, but also at the more global level, allowing the search for emerging properties of the multi-cellular system.  The clonal types in a simulation are defined by the user via mutation profiles. A mutation profile defines a (sub)population of cells where specific gene mutations are present. This mutations are introduced as constraints on the gene status, such as a constantly active/inactive status (activating/inactivating mutations) or a change in function (e.g. amplified or conditional behaviour). There is no limitation with respect to the number of mutation profiles or the number of mutated genes for each profile. In the extreme case, each cell can have a distinct mutation profile, which regulates the mutation level in each gene of the encapsulated network.

Subcellular gene networks
Networks are encapsulated within cell agents, and drive their decisions. Although any type of network model may be used, our current implementation exploits logical Boolean networks. The latter have been shown to preserve key dynamical characteristics of the network (24) and they can be designed quickly, without the need for accurate estimates for a large number of parameters. Briefly, the nodes of such logical networks can have only two states (active or inactive). Nodes may be connected with other nodes via links. All nodes are assigned logical rules that determine the current and future state of the node.
We use asynchronous network update, first described by Thomas (25), and widely adopted since in Boolean gene networks. We distinguish between four types of nodes: genes, receptors (input), output nodes, and fate-decision nodes. Genes have both incoming and outgoing links, receptors have only outgoing links, and output and fate-decision nodes have only incoming links. Cell-fate decision nodes have a crucial role in the simulations as they trigger the actions which determine cell behaviour.

Cell status and cell-fate decision rules
An activated cell-fate decision node is associated with a specific action, our current implementation includes the following actions:  Apoptosis. The cell dies and is removed from the simulation after a specified time.
 Necrosis. The cells dies; it remains in the simulation and occupies space, but it does not otherwise actively participate.
 Growth Arrest. The subcellular network simulation for a cell in growth arrest is "paused" for a given period of time. Furthermore, the cell interacts with the microenvironment at a reduced rate, for example the consumption rate drops to half of the normal rate. Cells are modelled as sinks that consume oxygen at a rate proportional to the local oxygen concentration. In particular, oxygen consumption is modelled through the equation: where R0 is the initial consumption rate, C is the concentration of oxygen in the specific grid cell, CN,f is a threshold value that determines the lowest possible oxygen concentration (currently fixed at 80% of the oxygen activation threshold), and finally CO2,opt is an optimal oxygen concentration, currently set to 0.28mM. The latter two parameters have predetermined values in microC, whereas the initial consumption rate (R0) and the oxygen activation threshold may be set by the user. The latter is a precondition that triggers the necrotic cell-fate decisions. We set the necrosis threshold at 0.02mM (21). Cells in growth arrest, consume oxygen at half the initial rate, whereas necrotic cells do not consume oxygen. Parameters such as the diffusion coefficient and the initial/boundary conditions of oxygen concentration can be adjusted as required by the specific application.
Additionally, we are assuming a conditional consumption and production rate of diffusible substances that determine cell-cell interaction. In case, the corresponding receptor or output node is activated, we assume a consumption or production of constant rate, that has been defined by the user (see

Geometric measures
We have used sphericity to assess the roundness of spheroids: where ψ is sphericity, V is the volume of the object and A its surface area.
The radius of a spheroid was determined at each stage of growth as the average distance between the coordinates of the initial centre point of the simulation, and the outermost cells of the growing spheroids.

Cloud Implementation
The models presented in this study are accessible via a web interface (https://www.microc.org). This interface also enables modification of the models and input parameters to conduct experiments other than those discussed in this paper. We have prepared a detailed protocol ( is a widely-used XML-based standard exchange format for sharing data; it is a flexible data model that can be used for object-relational data and a wide variety of graphs (26). GraphML is another XML-based, widely-used, data sharing format for graphs (27). GINML is an extension of GXL and can be produced for example by the logical model editor GINsim (14). The web server converts any of the above formats to the

Results
In-silico growth of a heterogeneous population of cells communicating within a threedimensional microenvironment To illustrate microC, we report six different ways in which it can be applied ( Figure 1). In each case, we illustrate the modelling; we report the model predictions and their agreement or disagreement with results from wet-lab studies; and we conduct a sensitivity analysis assessing the impact of parameter choice on those predictions. We chose the example of cancer, a complex disease where methods to study the link between genotype and phenotype are particularly and urgently needed (4). For the simulations, we focused on small cancer spheroids, and on pathways underlying the main hallmarks of cancer, namely sustained proliferative signals, resistance to cell death and evasion of growth suppressors (9).
As noted above, the proposed framewor enables simulations of individual cells or group of cells, where each cell contains an inner gene network that simulates the cascade of signalling events occurring in each cell upon stimulation, and determines the cell behaviour ( Figure 1). First, we describe how we evaluated the ability of microC to import the required gene networks. We then confirm the faithfulness of our ABM encoding of these networks, and show that we can reproduce the expected network behaviour. We then scale up the simulation by increasing the genetic and microenvironmental heterogeneity in a stepwise manner. Specifically, we consider a population of homogenous pre-cancerous cells growing in a 3D environment; then, we gradually introduce mutations in onco-and tumour suppressor genes ( Figure 1A) to study how they affect the multi-cellular growth ( Figure 1B). Subsequently, we allow the resulting clones (carrying different single or multiple mutations) to grow in competition ( Figure 1C), and study the parameters affecting clonal evolution. Finally, we investigate the additional effects of introducing extrinsic perturbations such as lack of oxygen and presence of growth factors, and enabling cell-microenvironment ( Figure 1D) and cell-cell interactions ( Figure 1E).

Evaluation of the dynamics properties of the inner-cell gene networks
Whilst there is a large number of methods, from traditional statistics to emerging deep-learning approaches, which permit accurate prediction of the behaviour of a population of cells given some initial inputs; such methods often constitute a black-box when it comes to interpreting mechanistically the results.
They emphasize the predictive ability of the model, with respect to the study of the possible causative mechanisms of the predicted behaviour. In this first example, we show how for each single cell microC enables monitoring of the paths which have followed to arrive at a given prediction. In particular, we show how microC enables understanding of the dynamical properties of the gene networks within each cell.
We consider a previously developed large mitogen-activated protein kinase (MAPK) network (29). This is a Boolean network that has been assembled in a knowledge driven, bottom-up, manner by considering gene-gene relations (e.g. activation/repression) as reported in the published literature. As such, this network it is not specific to a given context, cell line or genotype; however, it is likely to be somewhat biased towards signalling and interactions observed in cancer cell lines, as the pathways modelled are key to the hallmarks of cancer. First, we verify that microC can correctly import the gene network, reformat it, recodify it, and execute it using our ABM framework. We then ask how intrinsic perturbations, namely mutations commonly occurring in cancer cells, affect the functioning of the network.

Effect of mutations on gene network dynamics and multi-cellular growth
The range of abilities, or hallmarks, that a cell needs to acquire in order to progress to become a cancer cell have been extensively described, and the role of somatic mutations in such processes is well documented (9). However, the chain of events from a single mutation occurring in a normal cell to the acquisition of such hallmarks, and then to the occurrence and progression of cancer, is extremely complex and not fully understood. Here, we used microC to begin to dissect such questions by asking how the occurrence of single and multiple mutations might impact on the behaviour of the inner signalling networks, and how this results in different patterns of single-cell and multi-cellular growth.
We consider the growth of single-clone populations with loss-of-function mutations in the well-known tumour suppressor genes p53 and PTEN, and the activation of the known cancer driver EGFR. We observe considerable differences in the 3D growth patterns of cell types carrying different mutations, including WT, single mutations carrying clones (EGFR activating mutation, EGFR+; p53-or PTEN-loss-offunction mutations), and clones with co-occurring mutations in two or three of these genes. Importantly, cooccurring mutations show a significantly more aggressive phenotype with faster growth ( Figure 2C) with respect to cases in which single mutations occur ( Figure 2D, and Figure S3A). This is in agreement with experimental data showing that co-existing alterations in these specific genes in pre-malignant breast cancer cell lines significantly increases growth with respect to the single alterations in EGFR, P53 and PTEN genes (30,31). Furthermore, we confirm that the presence of active growth factor signalling is determinant for rapid growth. Specifically, all EGFR+ clones (EGFR+, EGFR+PTEN-, EGFR+p53-, EGFR+p53-PTEN-) exhibit initial exponential growth while clones with no activated EGFR signalling grew at a much slower rate ( Figure 2C).
We then compare the geometrical properties of spheroids growth in a 3D environment. This shows that spheroids consisting of the most aggressive clones (higher proliferations rates), namely EGFR+p53and EGFR+p53-PTEN-, grow in a much more symmetrical manner and have higher sphericity values ( Figure 2E, Figure S8), a measure of how near an object is to a perfect sphere (see Methods). Although there is not much published evidence to confirm this prediction, as spheroids are typically grown for just one cell line (one genotype) at any given time, and morphology is either not assessed or not reported, this is an extremely intriguing result. It highlights that simple physical space-constraint rules can dictate profound morphological differences between the growth of multi-cellular spheroids, hence tumours; and it reveals how this is linked with the cell genotype. Importantly, it also highlights how to make the correct assumptions about physical constraints when modelling cellular growth is paramount in order to predict the correct behaviour. In this simulation, we assume that proliferating cells can randomly occupy any of the free neighbouring slots, and, in the case that no slots are available, proliferation is postponed until a free slot becomes available. Whilst this is a reasonable approximation in the case of small spheroids grown in suspension, it needs to be evaluated further in, and it may not be suitable for, cases of more complex organoids cultures or simulations of tumours in-vivo.

Figure 3. Paths of clonal evolution: the order in which mutations occur affects spheroids' growth.
A

Scrutiny of clonal evolution paths: significance of the order at which mutations occurs
In the previous examples, we considered the simultaneous occurrence of multiple mutations. This reflects reasonably well experimental conditions where mutations are artificially introduced, but it does not represent clonal evolution in real tumours. Scrutiny of clonal evolution paths is an emerging field of research, however this has typically been done without taking into account the signalling and context in which mutations occur. In this example, we ask how the order in which mutations occur affects both individual cells and the overall multi-cellular growth.
Focusing on the aggressive EGFR+p53-PTEN-clone ( Figure 3A), we examine the six possible evolution paths for this clone ( Figure 3A-B). This revealed that the order in which mutations occur has a significant effect on a spheroid's growth ( Figure 3C). Specifically, the time of occurrence of the EGFR+ mutation was crucial: clones that acquire this mutation before the loss of a tumour suppressor resulted in larger spheroids, followed by those clones that acquire the EGFR+ mutation, and finally those that acquire the EGFR+ mutation ( Figure 3C).
These results confirm first that cancer can occur from multiple evolutionary paths, but they also suggest that a proliferation stimulus, namely EGFR activation, followed by the loss of a tumour suppressor, p53 loss-of-function, and then PTEN loss-of-function, results in the most rapid evolution. This appears to support the preponderance of experimental data indicating that p53 mutations is a relatively late, rather than cancer initiating, event in a number of cancers (see e.g. (32)). However, the evidence on this point is contrasting, and it is well known that p53 is affected by multiple mutations, with different functional implications (for a review see e.g. (33)), so further studies are needed to include such differences.
Importantly, and differently from other approaches to study clonal evolution, microC enables exploration of how differences in growth rates between cells carrying different mutations might result from the underlying characteristics of the network (Figure 1, 2 and S4). For example, EGFR activation directly affects the status of a large group of genes including ELK1, CREB, MYC and RAS that promote proliferation and block apoptosis ( Figure S5), so to acquire this mutation early would be very advantageous.
Finally, growth rates depend also on dynamical network parameters, such as the speed of the specific intercellular processes. We performed a thorough sensitivity analysis by changing the value of the parameters, starting from values suggested in the literature and expanding the range far from this initial choice (see Methods and Supplementary Material, and Figure S6). This identified the temporal decision window as one of the critical parameters, which is tightly linked with the temporal ratio between intracellular to intercellular processes. This analysis shows that changes in this parameter can significantly affect the resulting growth curves, and thus growth rate, and the effect is different when different mutations are considered. In our case, a value of 100 was the optimal choice in order to preserve dynamical behaviour of the model, while minimizing computational intensity of the simulations ( Figure S6).

Emerging competition patterns impact the growth of multi-clonal cell populations
Next, we asked how competition between different clones, grown together in a 3D multi-cellular spheroid, affects the 3D growth dynamics. To this end, we introduce multiple clones in the same environment ( Figure   4A and 4B), and compare the growth curves of the resulting spheroids ( Figure 4C), and their final size ( Figure 4D), with those observed when the same clones were grown in isolation, that is in single-clone spheroids.
In the multi-clonal simulations, we reveal that clones with aggressive phenotypes systematically take over the free surface area of the spheroids, thus restricting the rest of the clones in the central parts of the spheroid. This is particularly evident in the 2D geometry ( Figure 4A). This is a striking finding and it implies that the population of the aggressive phenotypes not only increases rapidly against the rest of the clones due to faster growth, but it also creates a physical barrier to for any further proliferation of the inner, slower growing, population.
This effect significantly affects both the growth curves ( Figure 4C) and the final population fractions Interestingly, these findings can be examined in light of the proliferation rates and the dynamical properties of the different clones ( Figure 2B). In particular, EGFR+p53-and EGFR+p53-PTEN-clones, are the only ones characterized by a single possible outcome for the cell fate decision, which is proliferation. This is a strong advantage when competing with clones that are characterized by multiple decision outcomes. For example, EGFR+ clones have a higher probability of undergoing apoptosis or growth arrest, which constitutes a strong disadvantage under direct competition since apoptotic events free space that are likely to be rapidly occupied by more aggressive clones.

Interaction between the genotypes and the physical microenvironment affects multicellular growth and clonal selection
Although the mechanisms are not yet fully elucidated, increasing evidence points to a key role of the microenvironment in clonal selection, hence multi-cellular growth of heterogeneous populations. To illustrate how microC can address this, we focused on hypoxia simulations as this is one of the major microenvironmental differences between cancer and normal tissue (34)(35)(36)(37). Specifically, we aim to study the formation of necrotic cores in larger spheroids ( Figure 5A-B), and their growth under artificially uniform well-oxygenated conditions (our Control configuration) with respect to growth when oxygen consumption and diffusion are enabled (we refer to this as the Hypoxia configuration) (see also Methods section).
We introduce the following into our spheroid model: a) we define environmental agents that can simulate the diffusion of oxygen in the microenvironment, b) we upload a new extended network which includes a hypoxia responsive module ( Figure S5), and c) a new cell fate decision, necrosis, to be triggered as response to extremely low oxygen concentrations (0.02mM O2,) (see below and Methods). We then simulate the modified model in our control and hypoxia configurations. Interestingly, we observe that EGFR+ clones are the only ones with the potential to outgrow a critical spheroid size that triggers necrosis.
For a configuration setup with initial consumption rate 0.005mM.s -1 and diffusion coefficient 10 -9 m 2 .s, the necrotic core that was formed ranged between 20 -360 cells ( Figure 5C) in agreement with previous reports (38). Accordingly, spheroids with necrotic cores had on average 4.8% -23.6% less living cells than spheroids without necrotic cores.
We perform a sensitivity analysis on the parameters involved in the diffusion-reaction equation (Equation 1), namely we change the oxygen consumption rate and the diffusion coefficient to one of a broad range of values ( Figure S10). We found that the consumption rate was the only parameter which, when changed, significantly affects oxygen concentration ( Figure S10A), which may impact spheroid growth indirectly by triggering necrosis. This prediction is in agreement with previous evidence showing oxygen consumption to be the most critical kinetic parameter (38).
We then study the combined effect of hypoxia (as defined above), and clonal competition (see section 3). Namely, we introduce all eight clones in the simulation environment and compared growth in isolation and competition, under either hypoxic and control conditions ( Figure S9B). Overall, our simulations showed that spheroids where hypoxia was enabled had reduced clonal diversity ( Figure 5D). In particular, we noticed that the EGFR+p53-and EGFR+p53-PTEN-clones increased their presence in the total population. This was additional to the initial enrichment of this population due to competition (white bar series in Figure 5D). On the contrary, the rest of the clones decreased their presence in the spheroids by an average of 4 -5%. This selection pressure, and consequent reduction of clonal diversity due to hypoxia can be explained by the spatial distribution of clones ( Figure 4A). Less aggressive clones are more likely to be segregated to the central part of the spheroids, that eventually will become necrotic under hypoxia. Of note, in this example we have not introduced possible mutations which might occur as a result of hypoxia, which could impact on the observed clonal diversity. 6. Impact of cell-cell signalling on multi-cellular growth in a heterogeneous microenvironment Our final example asks how accounting for cell-cell interaction affects the growth of spheroids. This experiment highlights one of the most innovative aspects of microC. As cell signalling and its relevance for normal development and disease are increasingly understood, a framework which enables such modelling, and study of the consequences of such signalling on cell behaviour, is correspondingly advantageous.

A-B. Cells with different mutation profiles (colours corresponds throughout
We study EGF signalling ( Figure 5E) We observed a significant increase in the population of cells under the Hypoxiasignalling condition.
These effects were consistent with the aggressiveness of the clones that we have seen in the previous examples; namely, p53-and p53-PTEN-clonal populations increased significantly more than the PTENand WT clones under hypoxia when EGF signalling was accounted for. Interestingly, enabling EGF signalling between cells further increased the tendency of hypoxia to reduce clonal diversity ( Figure 5D and 5G). This reveals that the reduction in clonal diversity attributable to central necrosis ( Figure 5D), is further increased due to the different proliferation rates of clones which are sensing EGF released under hypoxia.
A sensitivity analysis of the parameters determining the strength and the length of the EGF response, namely the activation threshold of the stimulus receptor and consumption rate of the growth factor, showed that low activation threshold values are more likely to activate the EGF receptor for any given EGF concentration above the threshold value, whereas lower consumption rates of EGF are more likely to activate the EGF receptor for longer ( Figure S10B). Both events increase the probability of proliferation, that in turn affects the size of the spheroid ( Figure S10C).

Discussion
To respond to the need to deepen understanding of the intricate relationship between cell genotype and phenotype, we have developed microC, which combines a powerful stochastic method, ABM, and gene networks modelling, into a novel framework for in-silico experimentation. Specifically, microC addresses the challenge of modelling and simulating the dynamics and evolution of a heterogeneous population of cells within a changing microenvironment. ABM is increasingly used to simulate the dynamics and evolution of complex systems in applications ranging from engineering to ecology (39)(40)(41)(42); this approach makes microC a naturally multiscale framework, enabling assessment of both multi-cellular systems and individual cells, and their physical and molecular properties.
One of the most innovative aspects of microC is that it substantially extends the capacity of ABM simulations of living cells, considering each individual cell as a meta-agent. Each cell is considered as a community of computational agents, the genes and molecules, acting and interacting within each cell. This enables the user to modify and/or replace gene networks in the cells, to define new constrains representing different types of mutations in different cells or in the same cell, and to customize cell-cell and cellmicroenvironment interaction parameters. As a consequence, the proposed framework is naturally suited to study how the behaviour of a specific perturbation, such as a mutation, occurring in individual cells within a three dimensional, dynamic microenvironment, affects the collective behaviour of other, similar or different, cells and the whole system. This reveals in some cases unexpected patterns of collective behaviour, not a-priori defined in the model, and which could not be predicted by observing the individual elements. Instead, there are emerging properties of the system as a whole and affect the evolution of the cell population and the behaviour of the single elements in return.
microC comes with example networks that can be used, modified or replaced with new ones by users.
Using one of these networks, a MAPK signalling network, we illustrated microC features through the study of growth patterns of 3D in-silico spheroids, affected by perturbation of intrinsic (mutations) and extrinsic (oxygen and growth factor availability) factors. We demonstrated that clones carrying multiple mutations have a growth advantage with respect to clones carrying single mutations, and that the order in which these mutations are acquired has a significant impact on the spheroid growth rate and final size. This was true both when the clones were grown in isolation and when they were allowed to compete with other clones. We also showed that this effect increased under more extreme microenvironmental conditions, namely lack of oxygen.
We also generated a hypothesis linking the geometric properties of spheroids with their genotype. The morphological characteristics of spheroids of the same clone vary considerably and may be a source for high data variability in in-vitro experiments (43). Here, we demonstrated how the dynamics of a gene network which is specified by the genotype and microenvironmental characteristics, affects the proliferation rate, which in turn has a significant impact on the overall spheroids' shape, irrespective of other elements such as adhesion forces which are long known to affect morphology (44).
Finally, we observed expected emerging properties, such as the formation of necrotic cores in the larger spheroids, and new unpredicted emerging properties, such as the effect of hypoxia and EGF signalling on clonal diversity of the spheroids. These properties were not encoded at the cellular level, but rose from the cell agents competition for nutrients, space and cellular interaction.
In summary, in this paper we assessed the capabilities of a new modelling framework that links genotype with phenotype, via gene networks and signalling pathways. We have provided a number of examples of using the framework, which not only enables a broad range and new types of modelling studies, but it also delivers a microenvironment for in-silico experimentation built using well-recognised formats and shared standards, thus enabling widely used model representations. This provides an environment for prediction, experimentation and reasoning using existing gene networks, as well as a starting point for gene network development.