Integrating Intracellular Dynamics Using CompuCell3D and Bionetsolver: Applications to Multiscale Modelling of Cancer Cell Growth and Invasion

In this paper we present a multiscale, individual-based simulation environment that integrates CompuCell3D for lattice-based modelling on the cellular level and Bionetsolver for intracellular modelling. CompuCell3D or CC3D provides an implementation of the lattice-based Cellular Potts Model or CPM (also known as the Glazier-Graner-Hogeweg or GGH model) and a Monte Carlo method based on the metropolis algorithm for system evolution. The integration of CC3D for cellular systems with Bionetsolver for subcellular systems enables us to develop a multiscale mathematical model and to study the evolution of cell behaviour due to the dynamics inside of the cells, capturing aspects of cell behaviour and interaction that is not possible using continuum approaches. We then apply this multiscale modelling technique to a model of cancer growth and invasion, based on a previously published model of Ramis-Conde et al. (2008) where individual cell behaviour is driven by a molecular network describing the dynamics of E-cadherin and -catenin. In this model, which we refer to as the centre-based model, an alternative individual-based modelling technique was used, namely, a lattice-free approach. In many respects, the GGH or CPM methodology and the approach of the centre-based model have the same overall goal, that is to mimic behaviours and interactions of biological cells. Although the mathematical foundations and computational implementations of the two approaches are very different, the results of the presented simulations are compatible with each other, suggesting that by using individual-based approaches we can formulate a natural way of describing complex multi-cell, multiscale models. The ability to easily reproduce results of one modelling approach using an alternative approach is also essential from a model cross-validation standpoint and also helps to identify any modelling artefacts specific to a given computational approach.


About Multiscale Modelling
Computational models of complex biomedical phenomena, such as tumour development, are becoming an integral part of building our understanding of underlying cancer biology. Mathematical models which are generated from biological data and experiments, e.g., in vivo or in vitro, through phenomenological observations in real patients help in explaining the mechanisms of this complex phenomenon. Quantitative, predictive models have the potential to significantly improve biomedical research by allowing virtual, in silico modelling.
Experimentalists and theoreticians have agreed that cancer progression involves processes that interact with one another and occur at multiple temporal and spatial scales. The time scales involved vary from nanoseconds to years: signalling events in the cell typically occur over fractions of a second to a few seconds, transcriptional events may take hours, cell division and growth and tissue remodelling require days, tumour doubling times are on the order of months, and tumour growth occurs over years, etc. Typical spatial scales range from nanometres for protein-DNA interactions to centimetres for a the development of a solid tumour mass, tumour-induced angiogenesis, tissue invasion, etc. These scales are strongly linked with each other. A phenomenon cannot be completely considered using a single scale, fully isolated without taking into account what happens at other smaller or larger scales.
In general, when incorporating different temporal and spatial scales into mathematical models, there are three commonly used viewpoints: the subcellular level, the cellular level, and the tissue level. Or, from a modelling point of view these levels can also be referred to as the microscopic scale, the mesoscopic scale, and the macroscopic scale, respectively. Cancer usually starts at the subcellular level marked by events that occur within the cell, such as genetic mutations, transduction of chemical signals between proteins, and a large number of intracellular components that regulates outward activities at the cellular level such as uncontrolled cell division, and cell detachment that leads to epithelial-mesenchymal transition (EMT), etc. The main activities of cell populations, such as interactions between tumour cells and host cells, intravasation and extravasation processes, proliferation, apoptosis, aggregation and disaggregation properties, are all viewed from a larger scale, that is the mesoscopic scale. The macroscopic scale concerns activities that occur at the tissue level such as cell migration, convection and diffusion of chemical factors, all of which are typical for continuum processes [1].
During the last decade or so many approaches to multi-cell, multiscale modelling of cancer growth and treatment therapy have been developed. For example, see articles by [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] for modelling details and [21][22][23] for reviews on multiscale modelling. The goal of each approach is, in the first instance, to be able to replicate observed experimental results and data. Since the biology of cancer is very complex, models have to focus on ''first order'' effects and introduce certain simplifications to make them computationally feasible. These simplifications often introduce modelling artefacts i.e., observed model behaviours or side effects which are due to the particular choice of the mathematical/ computational method. Isolating the source of modelling artefacts is very difficult and quantifying the impact of such modelling artefacts on model predictions is a daunting task. Therefore, in order to identify deficiencies and limitations of modelling methods currently in use, we have to be able to routinely conduct rigorous model cross-validation to ensure that predictions of different modelling approaches for a single biological system are in agreement, at least qualitatively, with each other and with experimental data. Since in many situations experimental data is hard to find or simply unavailable model, the issue of model crossvalidation is even a more important issue.
For mathematical models of biomedical systems to be credible and usable on a larger scale by a variety of biomedical researchers, they have to be: a) easy to set up, b) easily reproducible, c) transparent and open to peer review and challenge, d) publicly accessible and able to run on multiple operating systems without the need to recompile, and e) interactive and easily modifiable.
In this paper we present a case study on model cross-validation. We reproduce a cancer invasion model, originally described in [14], using a CC3D-based implementation and compare our simulation results to those of the original paper (in which a centrebased implementation was used). We document the details of model building based on the published article, highlight obstacles in reproducing published results and suggest a streamlined, systematic approach to cell-based model cross-validation.

CC3D-Bionetsolver framework for multiscale simulation
Modelling methodologies that explicitly represent individual cells are particularly appropriate for modelling and simulation of cancer invasion. There are important events and physical phenomena associated with cancer invasion on the single-cell level that can only be suitably captured in computational simulations by accounting for individual cell properties and important aspects of cell-cell interactions, such as changes in cell-cell contact area.
In modelling the various stages of cancer progression, certain computational and mathematical methodologies are more suitable than others. For example, in the case of solid avascular tumour growth, continuum models are well-suited since they capture bulk properties of tissues. Instead of explicitly treating individual cells, collective properties of the whole tumour tissue are modelled, such as cell density and oxygen concentration. An advantage of such an approach is that systems with a large number of cells, such as on the order of 10 6 or higher, can be handled. On the other hand, explicit representation of individual cells and their properties (e.g., locations, radii, morphology, surface area, volume, etc.) can become computationally burdensome when trying to model on the order of 10 4 to 10 6 cells. Nevertheless, such individual cell-based modelling approaches are capable of capturing phenomena and behaviour in multicellular systems that continuum strategies cannot capture.
Systematic development of biomedical models may be divided into the following distinct stages: a) creating a conceptual biomedical model, b) developing a formal description of the model based on an established modelling language such as the Systems Biology Markup Language or SBML, c) translating the formal language into a set of mathematical representations, for example, SBML is translated into a set of ordinary differential equations or ODEs, and d) developing a computational implementation of c).
''Traditional'' biomedical model building usually skips intermediate stages and jumps from a conceptual model description directly into low-level code. This is often convenient from the perspective of a modeller but it greatly impedes model crossvalidation, reuse or sharing. Problem solving environments, such as CC3D, Mason, or Flame, greatly reduce the amount of effort necessary to build models which rigorously follow stages a)-d) and at the same time offer the same level of flexibility in model construction as low-level programming languages. To build and run our models we used CC3D -an open source simulation environment based on the Glazier-Graner-Hogeweg (GGH) model which allows simulating cell behaviours on an individual cell basis, where individual cells can interact with each other or with the underlying medium. Several models of tumour growth and angiogenesis have already been simulated using CC3D environment. See, for example, articles by [24][25][26][27][28].
Multiscale models in CC3D-Bionetsolver are described using a combination of the CompuCell3D Markup Language (CC3DML) and Python scripting. Such a combined approach allows one to build complex biomedical models and does not require recompilation when running them. In a typical CC3D simulation ''static'' aspects of the model, such as lattice size, simulation runtime, list of cell types, initial conditions or cadherin affinities, are usually described using CC3DML. We can replace CC3DML with equivalent Python syntax. The ''dynamic'' part of the CC3D model is described using Python scripting. Since Python is a fullfeatured programming language, modellers are able to express complex cell type differentiation rules, couple cell properties to concentrations of diffusive chemicals or to cell-cell signalling or parameterise cell adhesive properties in terms of underlying molecular or gene regulatory networks.

Comparison of center-model and GGH-model for multicellular simulation
Here we briefly discuss the main differences and some similarities between the centre-based model of [14] and our model based on the GGH model. As indicated, CC3D is a software application that implements the GGH model, allowing lattice-based simulation of multicellular systems. Each biological cell is represented as a set of contiguous sites on a lattice and the system evolves in time through an energy minimization procedure. On the other hand, the centre-based model represents each biological cell in terms of the location of its centre of mass and its radius. This fundamental distinction between the two methodologies is illustrated in Fig. 1.
The cells of the centre-based model behave as elastic spheres and equations describing their behaviour and interactions are derived on the basis of classical mechanical concepts. The centrebased model approximates cell-cell contact areas using the radii of neighbouring cells and the distance between their centres. In contrast, the concept of cell neighbour has an explicit representation in the GGH model since two cells share one or more lattice edges (for 2D simulations) or faces (3D simulations). Because of these differences, each modelling approach has relative strengths and weaknesses with respect to capturing different biophysical processes and phenomena. On the other hand, the GGH and centre-based models also have some important similarities. Both methodologies use continuum, reaction-diffusion equations to model extracellular chemical fields and they both incorporate cellcell adhesion and mechanical constraints on cell shape. In each case, extracellular chemical fields can both modify and be modified by cell behaviours or properties such as cell growth rates, secretion, absorption and chemotaxis.

An application to multiscale modelling of cancer growth and invasion
The multiscale model of epithelial-mesenchymal transition (EMT) developed by [14] incorporates important aspects of Ecadherin-b-catenin signaling and its coupling to cell-level properties of intercellular contact and adhesion. This model requires explicit representation (on a cell-to-cell basis) of localised and spatially heterogeneous changes in cell-cell adhesion strength and contact areas. It is at this level of granularity that invasive cancer cells sense and respond to their environment. In terms of biological processes, the model of [14] captures cell-contact-dependent recruitment of Ecadherin and b-catenin to the cell membrane and reincorporation of both back into the cytoplasm. Computationally, the simulations incorporated (1) time-varying changes in cell-cell adhesion as a function of a system of ordinary differential equations (ODEs) for intracellular reaction kinetics of E-cadherin-b-catenin signalling and (2) changes in rate parameter values in the reaction kinetic model as a function of changing contact areas between neighbouring cells.

Results
We ran three sets of 3D simulations to model: (1) detachment waves of b-catenin in a thin layer of epithelial cells, described in subsection 0.5, (2) tumour growth and detachment of cells from a layer of epithelial cells, and (3) tumour growth and detachment of cells from a multicellular tumour spheroid, both described in subsection 0.6.
Initially, all cells were individually created in the shape of a cube of size 7|7|7 pixels, with gaps of 1 pixel length between them. From 0 MCS to 20 MCS we allow the cells to grow, during which time the volumes and surface areas of the cells increase and the cells become more spherical. During this period of the simulations, cell-cell contact areas undergo an equilibrating transient that does not reflect natural phenomena. Thus, we did not start the numerical integration of the differential equations (corresponding to the subcellular biochemical networks) until 20 MCS. Keeping in mind that the subcellular model is sensitive to changes in intercellular contact areas, if numerical integration occurred during the initial cell shape changes, unrealistic subcellular dynamics could occur as an artefact of these changes. Starting the integration at 20 MCS helped avoid this. All parameter values used in the computational simulations are listed in Table 1, unless stated otherwise.

Detachment Waves of Epithelial Layer Simulations
To simulate detachment waves of b-catenin in a thin layer of epithelial cells, we performed the simulation on a domain or a lattice of 264|224|60 pixels in x, y, and z directions, respectively, with the z-axis being perpendicular to the page. In the lattice, we place a sheet of cells with 30 cells along the x-axis (horizontal), 25 cells along the y-axis (vertical), and 1 cell along the z-axis. As mentioned, initially each cell occupies a cube 7|7|7 pixels and we insert a gap of 1 pixel between each cell, as can be seen in the top left figure of Fig. 2. The aim of giving a 1-pixel gap for this simulation is to give space for the cells to grow where cell volume increases followed by increasing cell surface area until the cells become spherical and tightly attached to each other, as can be seen from the top second left figure (MCS 30) of Fig. 2. The initial target volume for cells is set to 1:2 times the cell volume, making the average volume of each cell about 410 pixels. We set 1 pixel equal to 2mm. Therefore one tumour cell has a volume of about 3280mm 3 . The sheet represents a thin layer of tissue with a volume of 0:48|0:408|0:014mm 3 .
In the intracellular model, summarised in Eqs. (15)-(18), disruption of cell-cell adhesion occurs when there is an increase in the concentration of free b-catenin in the cytoplasm, that is when b-catenin concentration exceeds a specified threshold value as a result of disassociation of E-cadherin-b-catenin complex at the cell membrane. The threshold we specified for our simulations is 50:0. As explained previously, for cell detachment to occur, nuclear b-catenin must exceed this threshold value.
Depending upon whether b-catenin is above or below the critical EMT-MET threshold, the terms A 1 and A 2 in Eqs. (11) and (12) are changed appropriately. In each case, if b-catenin is below the threshold, the first terms in Eqs. (11) and (12) are used. On the other hand, if b-catenin is above the threshold, the second terms in each of the equations are used. However, in the SBML  implementation of our subcellular model, we do not actually implement two separate sets of equations for attached (belowthreshold) and detached (above-threshold) cells. Instead, we include both terms from Eq. (11) and both terms from Eq. (12) in the same SBML file. We effectively ''include'' or ''omit'' one term or the other (depending on whether cells are below-threshold or above-threshold) by either (1) setting a equal to 0 and n equal to a non-zero value (see n in Table 1) for the case of a belowthreshold cell or (2) setting n equal to 0 and a equal to a non-zero value (see a in Table 1) for the case of an above-threshold cell. In our CC3D-Bionetsolver implementation (i.e., our Python script), the increase of b-catenin concentration above threshold is deliberately initiated by decreasing the value of k z at a specified time (70 MCS) from k z~1 :5 to k z~1 :0. This parameter influences the association rate of b-catenin with the proteasome. When k z~1 :5, b-catenin-proteasome complex formation is sufficiently rapid to keep the b-catenin concentration of all cells well below the threshold of c T~5 0:0. However, when k z is decreased to a value of 1:0, b-catenin accumulates in the cytoplasm as a result of decreased proteasomal degradation.
We check the b-catenin concentration for every cell at each MCS. If the b-catenin concentration for a cell of type ''Low-BetaCat'' increases above a threshold value of 50:0, the cell type is changed to ''HighBetaCat'', n is set to 0:0 instead of 100:0 and a is set to 2:0 instead of 0:0. Similarly, when the b-catenin concentration of a ''HighBetaCat'' cell decreases below the threshold, the value of n for that cell is set to 100:0 and a is set to 0:0 and the cell type is switched to ''LowBetaCat''.
In the case of an EMT event (i.e., a cell type change from ''LowBetaCat'' to ''HighBetaCat''), changing the values of n and a as described is equivalent to swapping the expressions in A 1 and A 2 , between the below-threshold (½bvc T ) expressions and the above-threshold (½bwc T ) expressions. Physically, this corresponds to (1) a cessation of E-cadherin-b-catenin complex formation in the membrane (n~0:0) and (2) an accelerated dissociation of Ecadherin-b-catenin complex (i.e., the dissociation rate parameter d i (t), is increased by a~2:0) to form cytoplasmic (free) E-cadherin and free b-catenin. Together, the effects of these two phenomena are (1) an increased concentration of b-catenin in the cytoplasm and (2) a significantly reduced adhesion strength between the transformed cell and its neighbouring cells due to the loss of Ecadherin-b-catenin complex in the membrane.
In Fig. 2 an increase of b-catenin above threshold occurs in several cells, randomly. When one cell is induced with a high bcatenin concentration above the threshold c T , the cell becomes vulnerable to a loss of cell-cell attachment resulting in EMT. The event propagates outward from this localised event, affecting neighbouring cells. When a given cell detaches, the neighbouring cells in turn become vulnerable to EMT because of increased free b-catenin concentration inside the neighbouring cells. These cells detach from surrounding cells and the effects propagate throughout the layer of cells. At 130 MCS, we observe a small group of cells that start to detach. By 200 MCS, detachment waves have spread outward to adjacent cells. As time evolves, some cells at other positions also show detachment waves independently. Eventually around 500 MCS all cells in the layer have been affected and have detached from each other. This is the hallmark of EMT events. Due to the stochastic nature of the GGH model, regular waves of cell detachment which originate from one cell and then spread radially and regularly outward as seen in [14] can not be produced using CompuCell3D. Nevertheless, the results are qualitatively the same between our implementation and that of [14]. Also, our aim in this paper is to illustrate differences between the two approaches. They are different in different ways and we may biologically conjecture that neither is superior to the other.
In order to see how the concentrations of proteins inside individual cells vary over time, we wrote functions in our CC3D-Bionetsolver code to record values to output files of all concentrations each MCS. In Fig. 3, we plot concentrations of b-catenin, E-cadherin-b-catenin complex, and b-catenin-proteasome complex for a typical cell undergoing EMT and MET (see the simulation results shown in Fig. 2). Because of the stochastic nature of the GGH model, the concentrations fluctuate in response to fluctuations in contact area between cells. When the concentration of b-catenin increases significantly (due to loss of contact area between cells and complete detachment) the transition curve (during the period of detachment) becomes smooth (i.e., fluctuations cease). After the cell regains contact with other cells, the curve is observed to fluctuate again. The top figure of Fig. 3 shows the concentrations of b-catenin, E-cadherin-bcatenin complex, and b-catenin-proteasome complex when running the simulation up to 5000 MCS. Here we see three cycles of detachment and attachment, as shown from the repeated cycles of high and low concentrations of b-catenin (yellowish-green line).
Once the cells have detached from each other, they are free to randomly migrate from their original positions. However, in this simulation we did not apply any source of attractants that should cause the cells to migrate away from the layer; the cells detach but stay in their position or slightly move due to the stochastic nature of CC3D. Allowing the simulation to proceed all the way to 5000 MCS, we observe that b-catenin concentrations in detached cells gradually decrease back toward the threshold. When the concentration reaches the threshold, a in the internal model is again set to zero and n is set to a non-zero value. This alters the internal kinetics such that b-catenin is no longer rapidly degraded. Instead, it accumulates inside the cells and is reincorporated into E-cadherin-b-catenin complex. This process occurs rapidly so that near-zero concentrations of free b-catenin are observed in some cells as seen in the bottom left figure of Fig. 3 (bottom figures are plots of the concentrations to 1000 MCS or for one cycle of detachment).
The increase of E-cadherin-b-catenin complex increases the adhesiveness of cells and they undergo mesenchymal-epithelial transition (MET) resulting in the reattachment of neighbouring cells. Thus, cells again exhibit an epithelial phenotype, but this time with an irregular configuration of the cell layer, or loss of epithelial configuration. This is because of the random migration of cells away from their original positions that occurred when they were detached. The results we report here resemble those in Fig. 7 of [14]. We also observe from the simulation results that, after the first stage of EMT events, a few cells migrate so far that they cannot reattach to other cells. These cells remain as mesenchymal cells.
As for cells that cannot reattach after the first detachment (because they have migrated too far from other cells and thus remain mesenchymal), the concentrations of the subcellular proteins immediately reach their own steady states, as shown by plots of data in Fig. 4.

Tumour Growth and Invasion
For simulations involving tumour growth, the GGH target volume is incremented each MCS during growth phases at a constant rate of 0:02 times the current cell volume and GGH target surface area is also incremented at a constant rate of 0:02 times the current cell surface area. This results in a doubling of cell number approximately every 40 MCS. Cell division was set to occur when the volume of a cell exceeded 2 times its initial volume. This rate of growth was not necessarily intended to reflect in vivo rates of tumour cell growth. Rather, the purpose in our simulations is simply to let the tumour grow to a specified size so that we can then initiate EMT events and observe the subsequent dynamics of cell detachment and migration.   cells start out cube-shaped with size 7|7|7 pixels and a 1 pixel gap between each of them. In this simulation, we apply a linear concentration gradient of chemoattractant in the x-axis direction to generate cell migration.
From a single thin layer, the tumour grows and becomes a bulky layer as a result of rapid cell division. In the implementation of CC3D-Bionetsolver, it is possible to let the tumour grow indefinitely, but in this simulation we limit the cell division. Cells are permitted to grow and divide until the total number of cells in the tumour mass exceeds 500 cells. After 200 MCS we initiate an increase of free b-catenin concentration, as previously described, by reducing the value of k z from 1:5 to 1:0. Beginning around 500 MCS, some cells in the outer layer show high concentrations of free b-catenin. These cells eventually break away from the primary tumour mass and migrate in the direction of increasing chemoattractant concentration (away from the tumour mass).
As EMT events propagate over the tumour surface and more cells begin to detach from the outer layer, cells underneath the surface are exposed to the medium. The reduced amount of cellcell contact area that these underlying cells experience destabilises them and makes them vulnerable to EMT. The b-catenin concentrations in these cells increase above threshold and eventually the cells undergo EMT and detach from the tumour. In this way, the effects of early EMT events propagate into the tumour surface as the tumour mass grows and a continual series of detachment events are observed to occur.
To show the distribution of free b-catenin inside the cells that remain within or attached to the primary tumour mass, we provide a cross sectional view of the tumour mass along the yz plane in Fig. 6. Cells that are bound to other cells inside the tumour are roughly blue in colour. This indicates a free b-catenin concentration lower than the threshold c T . On the other hand, cells in the outer layer are exposed to medium and have less cell-cell contact area. The colour of these cells and those immediately underneath them range from yellowish green to dark orange. This indicates higher concentrations of free b-catenin near to or greater than c T . These results are in good agreement with the simulation results of [14] and experimental data of [29]. While our mathematical model does not explicitly model (or make a distinction between) the two types of b-catenin, we assume that the concentration of free b-catenin inside the cytoplasm (which we explicitly model) provides some indication of the concentration of nuclear bcatenin.
To study the sensitivity of multiscale dynamics to the b-catenin degradation rate parameter k 2 , [14] performed simulations with different values of k 2 (corresponding to different degrees of tumour cell invasiveness in tumour invasion assays). The invasion assay has been used in vitro as a measure of invasive potential of tumour cells. In our simulations, if k 2 is small this results in high concentrations of free b-catenin. If concentrations exceed threshold, then cells are susceptible to cell-cell detachment and may become invasive by breaking away from the primary tumour mass. In other words, sufficiently small k 2 can be thought of as a marker for malignant or invasive tumour cells. [14] used b-catenin degradation rate values of k 2~1 0 min {1 (fast degradation rate), k 2~1 min {1 (medium degradation rate), and k 2~0 min {1 (no degradation).
Our CC3D-Bionetsolver implementation is, for some reason, very sensitive to small changes in k 2 . In other words, tumour cell invasiveness in our simulations varies significantly with only small variations in k 2 values (much smaller than those used in [14]). Because of this sensitivity, we only varied k 2 within a very small range using a value of k 2~0 :95 for the low degradation rate, k 2~1 :0 for the medium degradation rate and k 2~1 :05 for the fast degradation rate. The resulting data that we collected from our simulations are summarised in Fig. 7, where, qualitatively, the results are the same as those in [14]. This illustrates the differences between the two approaches, which is the main aim of this paper. We have plotted the number of cells that reached a fixed distance over time. In the implementation, we remove cells that reach a certain distance from the main tumour mass. For this, we chose a distance 70 pixels. The maximum number of cells in the simulations (and therefore the maximum number of cells that can be removed) is 500 for all simulations. The curve obtained from the slow degradation rate simulation (k 2~0 :95) increases exponentially over a short period of time (purple line), while that obtained using k 2~1 :0 shows a more gradual increase in the number of removed cells (blue line). Finally, the curve corresponding to k 2~1 :05 (the fast degradation rate invasion assay) increases very slowly, indicating that only a small number of cells detached and were removed beyond the distance threshold of 70 pixels (green line). 0.6.2 Multicellular Spheroid Tumour (MTS). It was also of interest to see how our CC3D-Bionetsolver implementation could mimic the growth and invasion of multicellular tumour spheroids or MTS. These are spherical aggregations of (malignant) cells that can be grown in vitro. MTS are particularly used in cancer research for studying multicellular resistance or chemo-or radiotherapy assays [30]. They can be used to study cell-cell and cell-matrix adhesion in vitro as well as the influence of the environment on many cellular functions including differentiation, cell death, apoptosis, gene expression and regulation of proliferation. MTS exhibit the characteristics of threedimensional solid tumours.
For MTS simulations, we use a cubic lattice with size 240|240|240 pixels in x, y, and z directions. The simulations begin with one cube-shaped cell (size 7|7|7 pixels) placed at the centre of the cubic lattice. To maintain tumour compactness as cells divide and to prevent undesirable effects before we trigger detachment, we set the threshold value of b-catenin (c T ) to a relatively high value (c T~7 0). This ensures that no cells undergo EMT during the growth phase (in which the tumour is permitted to grow and become spherical in shape). An image of the tumour during this stage in the simulation is shown in the top right figure in Fig. 8. At 400 MCS the value of k z is decreased from 1:5 to 1:0  and at 500 MCS cells at the surface of the tumour spheroid can be seen with a high concentration of free b-catenin. In these simulations, we apply a radial chemoattractant gradient increasing outwardly in all directions from a minimum value at the center of the cubic lattice. After losing cell-cell adhesion with neighbouring cells (due to EMT resulting from above-threshold concentrations of free b-catenin), detached cells migrate radially outward in the direction of increasing chemoattractant concentration.
An interesting feature of the data collected from MTS simulations is that it gives an indication of the size of the tumour. In Fig. 9 we show cell positions for a single cell over time in simulations with different values of k 2 . Cell positions with respect to an initial position (where the cell was created as a result of mitosis) were written to an output file for selected cells. All data in Fig. 9 were taken from cell ID 1, which actually was not created by mitosis, but instead was present in the initial lattice configuration at 0 MCS. The data indicated by the red line were generated using a value of k 2~0 :95, the blue line represents data using k 2~1 :0, and the black line resulted from a simulation using k 2~1 :02. All data initially show an identical change in cell position from the centre of the lattice toward the same fixed position at 40 pixels. The cell resides here for an extended period of time before migrating quickly toward the edge of the lattice. This position of 40 pixels can be assumed to be the radius of the MTS. All three simulations (using different values of k 2 ) indicate the same value for tumour radius. On the other hand, for each k 2 value, the cell detaches from the primary tumour mass at a different time. This can be seen in the latter portions of each of the curves. In each case, there is a portion of the time-course that increases linearly (indicating the cell has detached from the main tumour). This linearly increasing portion occurs at a different point in time for each of the three simulations.
In a study by [31] of the growth and invasion of glioblastoma multiforme (GBM) in 3-dimensional collagen I matrices, invasive distance is defined as the radius of the entire GBM system minus the radius of the MTS. Thus, in our simulations, invasive distance corresponds to the distance that cells move radially outward from the MTS after detaching from the primary tumour mass as seen in Fig. 9.
Plots of invasion distance obtained from our MTS simulations show patterns similar to the data obtained from simulations using a layer of cells. This can be seen in Fig. 10. In the case of low bcatenin degradation rate (k 2~0 :95), invasion assay data, indicated by the purple line, show an exponential increase in the number of cells that have reached a distance of 70 pixels (i.e., have been removed from the simulations). For simulations using k 2~1 :0 and k 2~1 :02, cell removal rates are slower than for the simulation using k 2~0 :95, thus suggesting less invasive tumours.
Our simulation results have verified in vitro and in vivo experiments, that the level of invasiveness of tumour cells can be assessed from the extent of the loss of cell-cell adhesion. We can see in our simulations that high level of invasiveness is achieved by down-regulation of cell-cell adhesion, that is by decreasing the values of k 2 . We use k 2~0 :95 to simulate more invasive scenario and k 2~1 :0 for less invasive scenario as shown by the bottom right and bottom left figures in Fig. 11, respectively. This then must be followed by up-regulation of cell-matrix adhesion, another component that is required for successful invasion. This ''discrete analogy'' can be related to the inverse relation between cell-cell and cell-matrix adhesion, that is in order to invade and migrate through the surrounding tissue, cell-cell adhesion should be sufficiently low and cell-matrix adhesion should be sufficiently high.
In the study by [31] of glioblastoma multiforme (GBM) growth and invasion, it was shown the effects of increasing collagen concentration on the level of invasiveness of GBM cells, which is similar to increasing cell-matrix adhesion. GBM implanted in a high collagen concentration at early times shows growth patterns typical of malignant tumours where invasive cells gradually accumulate from the centre of MTS, invading outwardly in all directions, as shown in the top right figure of Fig. 11 of the experimental data. On the other hand, GBM that has been implanted in a low collagen concentration shows relatively few invasive cells that tend to invade along distinct branches, as shown in the top left figure. We can relate this to invasion assay simulations that we have performed using a slow b-catenin degradation rate k 2 (where a low k 2 value implies more invasive tumour cells). Here, we compare the results of our simulations in Fig. 11 with experimental data from [31]. Using k 2~0 :95 to simulate the invasive scenario and k 2~1 :0 to simulate the less invasive scenario, our simulations show different growth patterns of MTS that are strikingly noticeable between MTS with more invasive cells (bottom right figure) and MTS that is less invasive (bottom left figure). Qualitatively, our simulations are comparable to the experimental data.
Although we cannot directly compare our simulation results (bottom right and left figures) with the experimental results of [31] (top right and left figures), the patterns of invasion from decreasing cell-cell adhesion (our simulation results) show similarities with the patterns of invasion from increasing cell-matrix adhesion (experimental results). We note that GBM is a sarcoma and likely not use E-cadherin/b-catenin signalling as it is not originated from epithelial tissues. Instead, sarcomas along with other types of brain tumours, express N-cadherin that also mediate calciumdependent intercellular adhesion. Nevertheless, another paper by [18] developed another multiscale model of transendothelial migration (TEM) involving N-cadherin in which the pathway that they developed is not far different than the pathway using Ecadherin, based on their literature study. Hence, there may be possibility that the kinetics of intracellular proteins of GBM similar to the kinetics we have described here.

Discussion
In this paper, we have developed a multiscale individual cellbased model to study the roles of intracellular E-cadherin and bcatenin dynamics in cell-cell adhesion within tumours and tumour cell invasion. To model the intracellular biology, we used a mathematical model developed by [14]. We used CC3D, a latticebased simulation environment, for modelling the cellular level and Bionetsolver, a programming library, for modelling the subcellular (or intracellular) level. The integration of CC3D and Bionetsolver modelling tools enables us to study cell behaviours that are driven by the dynamics inside cells. It allows us to tune the level of detail at the intracellular level, without switching the simulation framework, and examine the effects of changing details at the cellular level.
In the model presented here, we examined invasive behaviours of cancer cells by modifying key parameters that are responsible for cell adhesion. Studies have suggested that nuclear b-catenin upregulation may characterise invasive cell populations in many types of cancer [29,[32][33][34]. It is possible to tune parameters that regulate the concentration of free b-catenin (including nuclear bcatenin) to study cancer invasiveness in silico. Two parameters considered by [14] and that we considered in our simulations are k z and k 2 . The parameter k z influences the association rate of bcatenin with the proteasome. A sufficiently high k z helps maintain appropriate cell-cell adhesion and tumour compactness because it keeps the b-catenin concentration of all cells well below a threshold value. However, when k z is decreased to a sufficiently low value, free b-catenin accumulates in the cytoplasm as a result of decreased b-catenin-proteasome complex. This leads to EMT events, in which cells lose cell-cell adhesion, break off from the primary tumour body, and migrate through and invade surrounding tissue. Varying the parameter k 2 (which influences the rate of b-catenin degradation) affects the invasive potential of cells as demonstrated by our simulated invasion assays. Sufficiently low k 2 results in cells that are more invasive than cells with a comparatively high k 2 . Our simulation results obtained by varying k 2 are qualitatively comparable to experimental data obtained in a study of multicellular tumour spheroids.
While we were able to qualitatively reproduce results from [14], there were noticeable discrepancies that are likely due to fundamental differences in the two simulation methodologies. In contrast to the centre-based implementation of [14], where it is possible to manipulate a single cell and thereby initiate detachment waves, our CC3D-BionetSolver framework does not easily permit a similar level of control. In other words, it was difficult to control cell properties in such a way to cause detachment waves to appear from a single cell in an epithelial layer and propagate radially outward in a regular manner (as shown in Figures 5 and 6 in the paper by [14]). Instead, by reducing k z from 1:5 to 1:0 in our GGH model, detachment waves randomly arose from localised groups of cells within epithelial cell layers and propagated outward irregularly. The discrepancies indicate that our approach and the centre-based approach are quantitatively different which was one of the aims of this paper. Nevertheless the results are qualitatively the same. See Fig. 12 for a comparison.
Another difference, is that the stochastic nature of the GGH model results in fluctuations of intracellular variables (concentrations) because of the fluctuating contact areas between cells. This can be seen from the plots of concentration data shown in Figs. 3 and 4. It should be noted that the question of what fundamental differences exist between these two simulation methodologies is distinct from the question of how well the simulation results collectively (of either methodology) reflect or correspond to actual experimental observations. This latter issue, while centrally important in the field of biological modelling, does not fall within the scope of the current study. Primary contributions of our study include the following: (1) It brings to light important differences that exist between two major individual cell-based modelling methodologies (the centre-based model and the GGH model) within the context of cancer biology and (2) it provides an introduction to CC3D-Bionetsolver, a recently developed multiscale framework for multicellular simulation.

Glazier-Graner-Hogeweg or GGH Model
The GGH model contains description of objects (e.g., cells, ECM, diffusible fields), interactions (e.g., cell-cell adhesion, morphogen-dependent cell growth), initial conditions (e.g., initial configuration of cells based on a time-lapse microscopy image), and the time evolution of cell properties (e.g., b-catenin concentration dynamics driving adhesive cell properties or rulebased cell type differentiation).
In the GGH model cells are represented as spatially extended domains on a fixed lattice, usually 3D Cartesian lattice or 3D hexagonal lattice. Each cell is simply a collection of lattice pixels having the same index (also referred to as cell id) s i ð Þ where i denotes lattice pixel, see Fig. 13. The GGH also allows compartmentalised cells where domains represented cells are further subdivided into subcompartments representing distinct parts of a biological cells (e.g., membrane, organelles, etc) [35].
The dynamics of cells in the GGH model is described by effective energy formalism and implemented as a Monte Carlo algorithm. At each step we randomly select a pixel i as a target pixel and randomly select one of its neighbouring pixels i 0 (in this paper we use consider pixels up to fourth nearest neighbour) as a source pixel. Then we attempt to change its index from s i ð Þ to the index of s'~s i 0 ð Þ. For each pixel copy attempt, we calculate the change in the overall system effective energy DE and accept the attempted pixel reassignment with probability P DE ð Þ: for DEƒ0 where T m is a parameter representing the effective cell motility. If i and i 0 belong to the same cell i.e. when s i ð Þ~i 0 ð Þ we do not copy the index.
The net result of this algorithm is that the cellular pattern in the GGH model evolves to minimise effective energy. We use this property of the GGH model to construct energy terms in such a way that their minimisation mimics actual cellular behaviour.
The simulation is subdivided (temporally) into so-called Monte Carlo Steps (MCS) which correspond to a unit of physical time. By convention, each MCS consists of number of pixel copy attempts equal to the total number of lattice sites. The conversion between pixel and physical distance (or MCS and physical time) depends on model parameters. In a simple case for example, in Bionetsolver we set timestepBionetwork to 0.03 and if Bionetsolver gets called every MCS then 1 MCS corresponds to 0:03 hours. In this paper we do not specifically set a relationship between MCS and the physical time because in the computational simulations we also incorporate cell mitosis or cell division which in the process itself also requires another time convention. In the mitosis process we do not apply any intracellular pathway, but instead we use a builtin mitosis function provided by CC3D.
The physical distance is recovered by converting pixels into units of length. This conversion is more straightforward than the correspondence between MCS and physical time and in our simulations we set 1 pixel to correspond to 2mm.  The effective energy, also called the Hamiltonian and denoted by either H or E, is the core of the GGH model. The Hamiltonian is typically expressed as a sum of terms, each term representing different cellular behaviours, interactions, mechanics, etc. The effective energy mixes true energies such as cell-cell adhesion with terms that mimic energies (e.g., the response of a cell to a chemotactic gradient of a chemical field). In our simulations we have used Hamiltonian containing adhesion energy term, two terms implementing constraints on cellular shapes (volume and surface) and one term implementing chemotactic force: The first term of Hamiltonian in Eq. 2 represents variations in energy due to adhesion between cells of different types, a boundary energy, that depends on J(t(s(i)),t(s(j))) between two cells (s(i),s(j)) of given cell types (t(s(i)),t(s(j))) at a link (the interface between two neighboring pixels). The sum is over all neighbouring pairs of lattice sites i and j (note that the neighbour range may be greater than one), the boundary energy coefficients are symmetric, where s(s) is the surface area of cell s, S t (s) is cell's target surface area, and l surface (t(s)) is cell's inverse membrane compressibility. The chemotaxis term E chemo represents interaction of cells with external chemical gradient c. It is easier, and somewhat less confusing to actually write expression for DE chemo than for E chemo . For a given pixel copy attempt we define change in the chemotaxis energy term as: where j denotes target pixel and i denotes source pixel for an attempt to copy s i ð Þ to pixel j. Strictly speaking E chemo is quasi energy term which is used to produce biased cell motion up (or down) the concentration gradient depending on sign of l chemo . The l chemo determines how strong a given cell responds chemotactically to the external gradient c.
In Supporting Information S1 we provide a brief CC3D tutorial which covers the basic usage of the GGH model. More detailed information about the model can be found in [28,36,37] and tutorial documentation on http://www.compucel3d.org/ Manualhttp://www.compucel3d.org/Manual.

Bionetsolver programming library
Bionetsolver is a Czz library with a high-level Python API that permits easy definition of sophisticated models coupling reactionkinetic models described in the SBML with GGH objects for execution in CompuCell3D. Bionetsolver makes use of the SBML ODE Solver Library (SOSlib) to implement reaction-kinetic network dynamics which can regulate the cell dynamics generated by the GGH core. For further information on SOSlib, the reader may refer to the paper in [38]. SOSlib provides functionality both for reading SBML models and solving them as a system of ODEs. In addition to this functionality, there are three classes -BionetworkSBML, BionetworkTemplateLibrary and Bionetwork -that provide some additional convenience in storing and manipulating SBML models as well as creating ODE integrators and time-stepping the integrators.
The Python API of Bionetsolver provides a set of 7 core functions that can be called from within a CC3D Steppable. These 7 functions are used for initialisation and manipulation of Bionetsolver objects from within the steppable. In this way, the entire specification of a multiscale (cell-subcell-level) simulation can be written in Python and executed in the CC3D player.
After the Bionetsolver API is imported and initialised, SBML models are loaded with a loadSBMLModel function and each SBML model can be added to one or more template libraries using the function addSBMLModelToTemplateLibrary. When loadSBMLModel is called, a string argument is required that signifies a name for the SBML model. Similarly, when ad-dSBMLModelToTemplateLibrary is called, the user provides the SBML model name (specified when loadSBMLModel was called) as well as a string argument that signifies the name of the template library. A single SBML model may be added to several template libraries and each template library may contain one or more SBML models. The code example below shows the use of these functions. Both of them are called from within the start function of a CC3D steppable. Notice that when loading an SBML model, both a model name and a model key are provided as arguments. The model key is used to reference specific SBML models in certain function arguments. In the example below, a single SBML model is loaded and is added to two bionetwork templates, ''LowBetaCat'' and ''HighBetaCat''. To associate bionetworks with CC3D cells, the template name must be the name of a CC3D cell type.
# Create a bionetwork SBML model named SimpleModel sbmlModelName = ''SimpleModel'' sbmlModelKey = ''SM'' sbmlModelPath = os.getcwd()+''/MultiScaleModels/sbmlModels/SimpleExample.sbml'' bionetAPI.loadSBMLModel(sbmlModelName, sbmlModelPath, sbmlModelKey) # Add ''SimpleModel'' to templates called ''LowBetaCat'' and ''High-BetaCat'' bionetAPI.addSBMLModelToTemplateLibrary(''SimpleModel'', ''LowBetaCat'') bionetAPI.addSBMLModelToTemplateLibrary(''SimpleModel'', ''HighBetaCat'') In addition, a setBionetworkInitialCondition function can be used to specify initial conditions for parameters and state variables in any SBML model within a template library. As arguments to this function, the user provides a template library name, a variable or parameter name and the corresponding initial numerical value of the variable or property. As indicated above, a set of SBML models stored within a template library can be associated with a CC3D cell type by providing the cell type name as the name of the template library. When the function initializeBionetworks is called, a separate bionetwork object is created for each cell of the given cell type and the previously specified initial conditions (specified using setBio-networkInitialCondition) are set for each of the bionetworks. Any parameters or state variables for which setBionetworkInitialCondition was not called are simply initialised, by default, to values specified in the original SBML models. A code excerpt below shows the use of these functions. Note that the variable name provided as the second argument to setBionetworkInitialCondition is prefixed with the SBML model key (in this case SM for SimpleModel). In case the same variable or parameter name, k1, appears in another SBML model in the same template library, this prefix uniquely identifies the parameter k1 found in the model indicated by the key SM. The initializeBionetworks function accepts a single numerical argument which corresponds to the timestep length used in the simulation to update the SBML models.
# Set initial conditions for templates ''LowBetaCat'' and ''HighBetaCat'' bionetAPI.setBionetworkInitialCondition(''LowBetaCat'', ''SM k1'', 0.9 bionetAPI.setBionetworkInitialCondition(''HighBetaCat'', ''SM k1'',0.1 # Create bionetwork instances from templates and initialise states initializeBionetworks(0.05) All of the functions mentioned above are initialisation functions and are called within the start function of the CC3D steppable. In addition, there are three more functions that are called within the step function of the CC3D steppable. These are (1) time-stepBionetworks, for time-stepping the ODE integrators, (2) getBionetworkValue, for retrieving SBML property and state variable values from the integrators, and (3) setBionetworkValue, for setting SBML property values of the integrators. The code excerpt below shows an example of how getBionetworkValue and setBionetworkValue can be called inside a Ô forÕ loop that iterates over CC3D cells. The cell ID (referenced by cell.id) is used to retrieve and set variable and parameter values for bionetworks (i.e., SBML models) associated with each cell.
for cell in self.cellList: S1 = bionetAPI.getBionetworkValue(''SM S1'', cell.id) bionetAPI.setBionetworkValue(''SM S1'', 0.1, cell.id) CC3D cell-level properties can be retrieved using procedures described in the CC3D documentation and the CC3D demo simulations. Finally, SBML property values can be set as a function of CC3D cell properties and, likewise, CC3D cell properties can be set as a function of SBML state variable values. This is how mechanistic coupling can be established between SBML (subcellular) and CC3D (cell-level) properties and dynamics.
Finally, Bionetsolver has a function copyBionetworkFromParent that can be used inside the updateAttributes function of a CC3D mitosis steppable to copy a parent cell bionetwork to a child cell that has just undergone mitosis. Copies of the parent SBML model integrators are created and the states of the original integrators are copied to the new integrators. The use of this function is illustrated in the code excerpt below.

E-cadherin and b-catenin kinetics
In [14], the kinetics of E-cadherin and b-catenin in a cell are conceptually modelled as follows. After being synthesised, Ecadherin is released to the cytoplasm as free E-cadherin (the concentration is denoted by ½E c ). Free b-catenin (concentration ½b) is assumed to be distributed in the cytoplasm and near the cell membrane. When there is signalling for cell-cell contact, free Ecadherin in the cytoplasm (½E c ) is transported to the cell membrane (concentration ½E m ) where its cytoplasmic domain binds to free b-catenin to form E-cadherin-b-catenin complex (concentration ½E=b) and the extracellular domain binds to Ecadherin-b-catenin complex of adjacent cells. If cell detachment occurs, E-cadherin-b-catenin complex dissociates, releasing free bcatenin and sequestering E-cadherin into the cytoplasm by endocytosis. The free b-catenin is then degraded and downregulated after binding with proteasome. These intracellular interactions are summarised in a schematic diagram shown in Fig. 14.
The mathematical model we use to describe the dynamics of these key chemical species concentrations (including the component influencing cell-cell adhesion i.e., E-cadherin-b-catenin complex) in each individual cell i formulated precisely as in [14] i.e.
The parameter k m is the rate of production of b-catenin, c i (t) is a time-dependent function describing the amount of cadherin stimulated to form bonds in response to increased cell-cell contact, and d i (t) is a function that describes the amount of cadherin released as a result of broken bonds during cell detachment. These functions (c i (t) and d i (t)) depend on the rate of change in contact area between adjacent cells. They are defined as follows: and d i (t)~X newdetachments a d,j (t)r d : where r c and r d are rate constants associated with E-cadherin translocation between the cell membrane and the cytosol. The rate constant r c reflects how quickly E-cadherin is transported from the cytoplasm to the cell membrane in response to cell-cell contact signalling and r d reflects the rate of the reverse action (from the membrane to the cytoplasm) when cell detachment occurs. In [14], a c,j (t) and a d,j (t) are time-dependent functions quantifying the rate of increase in contact area (when E-cadherin is transported from the cytosol to the membrane) and the rate of loss of contact area (when Ecadherin is reincorporated from the membrane to the cytosol), respectively. These functions are defined as follows [14]: whereâ a(t) j is the approximated contact area between cells i and j at time t divided by the surface area of cell i.
The condition for attachment is assumed valid as long as the concentration of free b-catenin ½b is below a threshold value, c T . For detachment to occur, the amount of free b-catenin in the cytoplasm must be sufficiently high with an additional contribution from dissociated E-cadherin-b-catenin complex. In other words, free bcatenin must be higher than the threshold value (½bwc T ). This claim is based on several studies that have found upregulation of b-catenin in the cytoplasm (or termed nuclear b-catenin) at the invasive front of colorectal carcinomas [29,33], invasive breast cancers [34], fibrosarcoma, clear cell sarcoma and carcinosarcoma [32]. Nuclear accumulation of b-catenin initiates the loss of epithelial differentiation and gain of mesenchyme-like capabilities of the tumour cells at the invasive front, while in the central areas of the primary tumours, nuclear b-catenin was found to be localised to the cell membrane. Nuclear accumulation of b-catenin has been the most powerful predictor of liver metastasis in colorectal cancer. This may be an important marker for adjuvant therapy or other treatment modalities.
In order to obtain a nondimensional system of equations, we nondimensionalise Eqs. (7a)-(7d) in the usual way by setting the following dimensionless variables: The dimensionless parameter values used in the simulations of this chapter can be found in Table 1. These values are based on the parameters used in [14] where we nondimensionalise the parameters by assuming T~1 min and E~1 nM, unless stated otherwise.

Supporting Information
Supporting Information S1 Information about how to configure, run, and modify two simple CC3D GGH-based simulations.