PhysiBoSS: a multi-scale agent-based modelling framework integrating physical dimension and cell signalling

Abstract Motivation Due to the complexity and heterogeneity of multicellular biological systems, mathematical models that take into account cell signalling, cell population behaviour and the extracellular environment are particularly helpful. We present PhysiBoSS, an open source software which combines intracellular signalling using Boolean modelling (MaBoSS) and multicellular behaviour using agent-based modelling (PhysiCell). Results PhysiBoSS provides a flexible and computationally efficient framework to explore the effect of environmental and genetic alterations of individual cells at the population level, bridging the critical gap from single-cell genotype to single-cell phenotype and emergent multicellular behaviour. PhysiBoSS thus becomes very useful when studying heterogeneous population response to treatment, mutation effects, different modes of invasion or isomorphic morphogenesis events. To concretely illustrate a potential use of PhysiBoSS, we studied heterogeneous cell fate decisions in response to TNF treatment. We explored the effect of different treatments and the behaviour of several resistant mutants. We highlighted the importance of spatial information on the population dynamics by considering the effect of competition for resources like oxygen. Availability and implementation PhysiBoSS is freely available on GitHub (https://github.com/sysbio-curie/PhysiBoSS), with a Docker image (https://hub.docker.com/r/gletort/physiboss/). It is distributed as open source under the BSD 3-clause license. Supplementary information Supplementary data are available at Bioinformatics online.

4 Figure S4: Effects of simulation parameters in spheroids.
A: Growth without TNF. Snapshots of a simulation of a spheroid growth at t = 0, 8, 16 and 24 h. B: Growth under a 0.5 ng/mL continuous dose of TNF. Snapshots of a simulation at t = 0, 8, 16 and 24 h. C: Growth under a pulsed injection of TNF, 0.5 ng/mL for 5 min every 300 min. D: Variation of secretion rate value. Time evolution of the cell population in each fate (three leftmost panels) for a secretion rate of: 0 (1st panel), 0.1 (2nd) and 10 (3rd). TNF was injected at the beginning of the simulation for 5 min at a dose of 1 ng/mL. Note that when secretion rate is 10, TNF accumulates in the environment (grey shading).

Note on the different rates of transition
Regarding the different rates of transition used in our model, two sets of nodes were established: first, most of the nodes present in the model represent protein-to-protein interactions and signal transduction such as Apoptosis pathway activation, TNF activation of Survival pathway and ROS activation. The value for these nodes was 1 MaBoSS rate that corresponds roughly to 10 minutes of PhysiBoSS time. Additionally, a subset of nodes, the ones that represent the activation of BCL2 and cFLIP as well as translation of different mRNAs (mROS, mcIAP and mXIAP) (all of them activated by NFkB) and the translation of cIAP from mcIAP, were considered to be activated in much more time (1/24 MaBoSS rate that corresponds roughly to 4 hours of PhysiBoSS time). We considered that protein-to-protein interactions and signal transduction happened 24 times faster than the transcriptional and translation steps. This was an educated guess based in previous works: Hirata et al. (2002) and Shimojo et al. (2008) found that transcription is 12 times faster than translation in mouse cells (and Lewis (2003) considered this as well for zebrafish). In our case, most interactions are protein activations, and those are faster than translation, therefore considering a factor of 24 between protein activations and translation.

Note on spheroid size comparison to experiments
In Figure 4A, we present simulations without TNF injections: spheroids with an initial radius of 100 µm with ≈1000 cells reached a size of approximatively 150 µm radius with ≈4500 cells in 24h, a growth of 4.5x factor. In this case, these results consider that oxygen is freely available to cells, even in the centre of the spheroid. Additionally, in Figure S6A in Supplementary Information we present results without TNF injection and with oxygen and nutrients limitations, a closer set-up to experimental data. In this case, an initial 9000-cells spheroid has a 200 µm radius and after 24 hours grows to a final ≈16500 cells spheroids with a 245 µm radius, a growth of 1.8x factor.
Experimentally, comparable population growths have been observed in Freyer et al. (1984) work, which measured a fivefold increase in two spheroids (from 1000 to 5000 cells in one and from 600 to 3000 cells in another). Nonetheless, the number of cells in the spheroids are quite different, which is likely due to the different cell strains we are comparing: V-79 Chinese hamster fibroblasts in Freyer et al. (1984) and 3T3 Swiss albino mouse embryo fibroblasts simulated in present work with data taken from Tay et al. (2010) and Kellogg et al. (2015). This density discrepancy is also a probable cause for the differences we see in terms of volume changes. Freyer et al. (1984) 1000-cells spheroid's radius goes from 210 µm to a final 230 µm for the 5000-cells spheroid (9.5% radius increase) and 600-cells spheroid's radius goes from 150 µm to a final 220 µm for the 3000-cells spheroid (32% radius increase), while our simulated spheroid has a 50% radius increase without oxygen limitation. However, in simulations with oxygen and nutrient limitations, we observed that these are impeding factors for growth once the spheroid goes beyond a threshold volume. Indeed, considering these limitations, we measured a radius increase of approximatively 25%, which is in range of the experimental values. These considerations highlight the importance of multi-scale frameworks that take into account both environmental and genetic information.
5 Figure S5: Genetically heterogeneous population under TNF treatment. Colour code: green, Proliferative cells; red, cells committed to Apoptosis; black, cells committed to NonACD. Grey level background in graphs indicates presence of TNF in continuous injection at 0.5 ng/mL. All simulations are in oxygen-limited regime, initial spheroid radius is 200 µm, which accounts for roughly 9000 cells, + stands for over-expression and -stands for knock-out.

Simulations of heterogeneous population composed of 75% of WT cells (orange) and 25% of mROS+ and cIAP-mutated cells (blue) or with CASP3+ and Cytc+ mutated cells (light blue). A: Snapshots of a simulation for each case at initial and final time (24 h), with cells coloured by
7 Table S1: Different multi-scale modelling frameworks Examples of different multi-scale modelling frameworks that could be used to address in different manners problems similar to the ones presented in this work. Intracellular framework can be differential equations, ordinary (ODE) or partial (PDE), or logical modelling (Boolean or multilevel). Extracellular framework can be based on a lattice (cellular automata (CA) or Cellular Potts (CP)) or independent from a lattice description (agent-based (AB) or vertex-based (VB)). Dimensions (Dim.) can be 1: 1D (lines), 2: 2D (areas), 3: 3D (volumes) or 3m: 3D objects in a monolayer. OS compatibility can be Linux (L), Mac OS (M) or Windows (W).