An in silico hybrid continuum-/agent-based procedure to modelling cancer development: Interrogating the interplay amongst glioma invasion, vascularity and necrosis

This paper develops a three-dimensional in silico hybrid model of cancer, which describes the multi-variate phenotypic behaviour of tumour and host cells. The model encompasses the role of cell migration and adhesion, the in ﬂ uence of the extracellular matrix, the e ﬀ ects of oxygen and nutrient availability, and the signalling triggered by chemical cues and growth factors. The proposed in silico hybrid modelling framework combines successfully the advantages of continuum-based and discrete methods, namely the ﬁ nite element and agent-based method respectively. The framework is thus used to realistically model cancer mechano-biology in a multiscale fashion while maintaining the resolution power of each method in a computationally cost-e ﬀ ective manner. The model is tailored to simulate glioma progression, and is subsequently used to interrogate the balance between the host cells and small sized gliomas, while the go-or-grow phenotype characteristic in glioblastomas is also investigated. Also, cell – cell and cell – matrix interactions are examined with respect to their e ﬀ ect in (macroscopic) tumour growth, brain tissue perfusion and tumour necrosis. Finally, we use the in silico framework to assess di ﬀ erences between low-grade and high-grade glioma growth, demonstrating signi ﬁ cant di ﬀ erences in the distribution of cancer as well as host cells, in accordance with reported experimental ﬁ ndings.


Introduction
Gliomas are a common type of primary central nervous system tumours [1,2]. The understanding of their origin has evolved from the idea that they could be initiated by glia [3] to the current conjecture of expansion from restricted progenitors, or neural stem cells [4,5], whence tumour characteristics and behaviours are somewhat bound by hysteresis to its origins. The median overall survival of patients affected by gliomas varies from 6.5 to 8 years for low-grade gliomas, to a measly 1.25 years for glioblastoma multiforme [6,7]. The most important hallmark of this kind of tumours is its phenotypic plasticity, underpinning both its ability for malignant progression from low-to highgrade types, and its invasive behaviour regularly resulting in a widespread infiltration of otherwise healthy brain tissue by glioma cells. Especially the latter is largely responsible for the lack of effective therapeutic approaches, and ultimately the aforementioned poor prognosis.
Despite the accumulated clinical wisdom and biological data, the invasion dynamics and its underlying mechanisms remain unclear as a result of a lack of ethical means to investigate the connections between clinical observables and in vitro data, for as instance the geometry of tumour growth and shape can affect the biological data we can sample in vivo [8]. The question about how to account for spatial characteristics when analysing cancer growth dynamics is novel and largely unexplored by current standards. Mathematical and computational modelsalso collectively referred to as in silico modelscan play an important role in extrapolating findings from in vitro experiments to in vivo. Thus, in silico modelling as a means to synthesise, to test the compatibility of different data speaking to the same theory, and to discover known unknowns [9] has a long history in this field, as recently reviewed by Alfonso et al. [10] and Magi et al. [11].
Herein we contribute a new tool to serve the community trying to clear this complicated puzzle and open up a path towards better prognosis: a novel method that combines the advantages of agent-based and continuum models. In particular, the agent-based approach allows to model detailed cellular dynamics (encoded either as logics or as simulated genetic expression patterns), cell-cell and cell-matrix interactions. Importantly, agent-based modelling allows to take into account highly heterogeneous system dynamics. This is problematic to reconcile with the continuum approach. However, continuum models allow to cover a large spatial volume in 3D.
Specifically, we develop an interface between two software frameworks, namely the agent-based modelling software platform BioDynaMo [12] and the finite element (FE) analysis software for multiscale modelling FEB3 [13][14][15]. The resulting hybrid toolkit offers the advantages of both types of simulations, enabling highly efficient models of large-scale 3D biological systems. Notably, it can take into account macroscopic dynamics, which can be measured with medical imaging for example, as well as cell-cell interactions, allowing the experimental assessment of model predictions on multiple scales [16].
Models that combine discrete cellular dynamics with continuousscale resolution models (such as growth-consumption dynamics) have already emerged in the past decade [17] usually referred as hybrid models. In general, hybrid models are built by building interfaces across different modelling solutions, each representative of the best practices accumulated by the research community most concerned with the specific problem and scale (tissue growth, electrical activity, mechanical action, metabolism, etc.) [18][19][20]. Our approach builds on our predecessors' intuition by offering ways to model intracellular signalling and the transport of growth factors and oxygen by differential equations and the finite element method (Section 2.1), and combining this with discrete cell dynamics simulated by an agent-based approach (Section 2.2).
The agent-based simulation engine BioDynaMo is a member of a powerful breed of computational biology solutions [21], together with for instance PhysiCell [22], concerned with but not limited to (i) simplifying access to simulation by providing toolkits and instructions; (ii) translating the best practices in parallelization and code optimisation from enterprise software to open-source software for research; (iii) promoting seamless scalability of work from single researcher's exploration of a conjecture, to large scale simulations and collaborations across the cloud. The ability to scale from laptop to heterogeneous multicore distributed architectures within the same framework would promote better testing and reproducibility all along the pipeline of work of the community, while maintenance and upgrading of the environment would ensure that the community could shoot for "HPC on the cloud" performances while greening and securing their practices.
Importantly, what makes our proposed in silico hybrid procedure stand out is that it solves the multiscale modelling problem by coupling the tissue scale with the cell scale by volume averaging (Section 2.3). To the best of our knowledge, this approach is the first one to take into account detailed genetic rules, cellular behaviours and cell-cell interactions in 3D, while simulating substance diffusion using the highly efficient volume averaging technique. All the above said, we believe the proposed hybrid solution to hold a great potential interest for the mathematical and computational oncology communities, and following we will detail how the hybrid simulation environment works.

Methods
The methodological presentation of the in silico hybrid modelling platform is outlined in the following paragraphs. The first and second subsections will overview the macroscopic and microscopic in silico methodologies respectively, while the third will describe the coupling of the former approaches into the hybrid platform.

Macroscopic domain: continuum-based finite element solver
Let us denote the volume of the tissue domain of analysisreferred hereafter as the macroscopic domainby Ω, which is the total volume of the biological tissue involved: both host and non-host (cancerous) tissue, and is bounded by surface Γ. The boundary Γ is set sufficiently distant from the focal region of interest (e.g., the tumour) in order to avoid the imposed boundary conditions impacting the solution. Fig. 1A shows a schematic representation of the three-dimensional domain of analysis, with the total volume Ω being equal to (20 mm) 3 .
At the macroscopic (length scale) level, the continuum-based model accounts for three compartments, each consisting of various species: the biochemical components compartment, the cells compartment and the extracellular matrix (ECM). The biochemics compartment encompasses the balance of oxygen and nutrients, the concentration of which is denoted by ξ , growth factors and enzymes, each type denoted by g j ( = … j N 1, , g ) and μ respectively. The cells compartment encompasses the cells involved in the biophysical problem under consideration, each cell species being denoted by c i ( = … i N 1, , c ), while the ECM compartment collectively accounts for the structural aspects of the stroma (again both host and non-host tissue), denoted by ∊. N c and N g above denote the number of individual cell phenotypes and growth factors considered in the model, which can take any value. As described below, each component is described at the continuum level through a boundary value problem: the balance of each species is mathematically Fig. 1. A. Cut-through illustration of the 3D unstructured finite element mesh (outline surface mesh is shown slightly transparent) B. From the same perspective, the representative averaging volumes (RAVs) are illustrated in the 3D domain (scaled up for visualisation purposes, while FEs are hidden) with each agent represented as a small spheroid, coloured with respect to their phenotype (purple for the neurons, green for the healthy glial cells, whereas orange for the normoxic cancerous cells and red for the hypoxic ones). C. From bottom to top, consecutive zoom-in pictures (2.5×, 9×) depict the cells, shown as spherical agents and coloured with respect to their phenotype. The yellow arrows point out the scattered cancerous cells, located at adjacent RAVs at the centre of the 3D domain of analysis. Light blue thick line depicts a one millimetre scale, while the cells size in the RAVs is scaled up by a factor of five for illustration purposes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) modelled via coupled differential equations, accompanied with proper boundary and initial conditions.
Assuming the presence of a microvascular network that supports tissue with a uniform source of oxygen and other nutrients, oxygen diffuses into the interstitial space and is consumed by the cells; thus, the balance equation of the oxygen saturation level ξ , expressed in a Lagrangian frame of reference, reads where D ξ is the isotropic diffusion coefficient for oxygen (given in m 2 day −1 ), α ξ is the constant oxygen production rate owing to the supply from the microvascular network (expressed in day −1 ), and δ ξ i , is the oxygen consumption rate by each of the i-th cell phenotype species (in ) controls the supply of oxygen from the blood vessels to the adjacent tissue and, consequently, the cells. Also, "d dt . / " denotes the material derivative (= ∂ ∂ + ∂ ∂ t v xv . / · . / ; the interstitial biofluid flow velocity vector (expressed in spatial coordinates at a material point in the extracellular matrix). The present macroscopic model accounts for chemical cues that are responsible for cell activation, growth, mitosis, migration etc. The balance of these growth factors (each species denoted by g j in dimensions: mg/mm 3 ) is described by a set of reaction-diffusion equations [23,24] that account for random spatial diffusion, growth factor secretion and uptake by the various cell species, c i , and the natural decay. The corresponding partial differential equations (PDEs) are expressed as where D j is the isotropic diffusion coefficient for each of the chemical agents involved, g j (in m 2 day −1 ), δ j represents the decay rate of the agent, and α j i , is the gross rate of g j due to both secretion and uptake by the cells (given in day −1 ). In addition to this, the proposed model accounts for the dynamics of the concentration of matrix metalloproteinases (MMPs), μ. MMPs are enzymes that degrade the extracellular matrixalso referred to as matrix-degrading enzymes. The balance of the MMPs concentration, μ, in the extracellular space of the tissue obeys an identical PDE to that of Eq. (2) [16], with corresponding model parameters: D α , μ μi , and δ μ . This model treats the ECM as a porous solid medium of porosity, e p , with the volume fraction occupied by stromal components (collagen, fibronectin, laminin) is: ∊ = − e 1 s p . Therefore, the sum of the volume fraction for all phasesthis includes the ECM, ∊ s , the cells, ∊ c , and the interstitial fluid, ∊ fgives unity: ∊ + ∊ + ∊ = 1 s c f . The present model assumes fixed volume fraction for the interstitial fluid, whereas the structural integrity and composition of the ECM, and hence ∊ s , is assumed to change with time. We describe structural changes at the stroma of the tissue (at the macroscopic level) using a first-order ordinary differential equation for the ECM volume fraction that accounts for the remodelling of the ECM, and the degradation of the matrix due to the presence of chemical cues (here metalloproteinases are modelled to cleave ECM fibres) at the interstitium [19] where ∊ δ is the ECM degradation rate due to the presence of matrixdegrading biochemicals μ (given in days −1 ), ∊ α is the ECM self-regeneration rate and ∊ α i , is the cell-stimulated ECM remodelling (e.g. deposition of collagen or/and fibronectin as a result of cell activity). Also, Finally, the population of each cell species is controlled by the following PDE: where α i is the cell production (if negative-valued, the cell decay) term, while the first term of the equation in the right hand side describes the potential for migration or invasion of the i-th cell phenotype. Both terms are calculated by the solver at the microscopic domain and, subsequently, projected to the macroscopic domain solverin this case the finite element solver as explained below. In all differential Eqs. (1)-(4) a zero-flux boundary condition for all species is considered, which is applied at the external boundary of the macroscopic domain of analysis, Γ, see Fig. 1. Also, proper initial conditions are selected based on the simulated problem as explained in the Results section. The weak form of all coupled differential Eqs.
(1)-(4) can be subsequently derived using calculus of variations, and the (continuous Galerkin) Finite Element (FE) method is employed as a discrete method [25]. As in conventional three-dimensional FE problems, the domain of analysis is discretised here with 3D elements (e.g., tetrahedrons, hexahedrons) that employ Lagrange nodal basis functions. Time discretization and numerical solution of the above equations was carried out through an explicit Euler scheme, where the time increment was properly selected such that explicit time-integration produces stable numerical solutions after the Courant-Friedrichs-Lewy (CFL) condition is satisfied. More precisely, assuming that the diffusion coefficient for the oxygen, D ξ , takes the highest value compared to the corresponding coefficient of the growth factors, and provided the smallest edge length of all FEs in the mesh, Δℓ, then, the time step must comply with the CFL limit: The model parameters adopted for Eqs. (1)-(3) are listed in Table 1 of the supplementary material.

Microscopic domain: agent-based modelling solver
At the microscale, cell behaviour and their mechanics is described using an agent-based modelling (ABM) approach. Opposed to the use of (partial or ordinary) differential equations in continuum-based modelling, in ABM cells are treated as discrete objects namely individual agents, which obey certain rules and mechanisms with respect to their phenotype, microenvironment and other external (biochemical or biomechanical) factors. For an overview and comparison between continuum and discrete models, the reader is referred to the review of Schaller and Meyer-Hermann [26]. In the present ABM, cells are described as spherical entities residing in a control volume, the representative averaging volume (see Section 2.3). For simplicity and as it will be clarified in the following subsection, the control volume of a microscopic domain is considered cubic with fixed dimensions: (100 μm) 3 , which contains various types of cells, C i , host and cancerous ones. The present brain tumour agent-based model assumes the presence of neurons (Ne) and glial cells (G), as well as cancerous glial cells and cancer-associated fibroblasts (CAF). Also, the model distinguishes between normoxic-state cancerous glial cells (CG) and hypoxic-state cancerous glial cells (hCGs). For Ne and G species, the cells densities considered are 5 × 10 4 per mm 3 and 10 5 per mm 3 respectively to match observed density in human cortex [27]. Moreover, assuming avascular tumour growth of a primary brain cancer, the initial population of the CG species has a density of 10 4 cells per mm 3 , which for simplicity are evenly distributed at the centre of the domain of analysis (within radius~0.5 mm). Fig. 2A illustrates in a flowchart the phenotypic behaviour of host (Ne, G), cancerous (CG, hCG, CAF) and necrotic cells (NC) with respect to the level of oxygen saturation, chemical cues and the extracellular matrix (ECM). The present model assumes a passive role for the host cells: balance between cell apoptosis and duplication is preserved in the simulation unless oxygen drops below a certain level, then they simply die and corresponding cell agents are removed from the simulation.
However, host fibroblasts trans-differentiation is implicitly described by explicitly increasing the CAF population in a logistic fashion and in proportion to the CGs population. CGs are prone to unregulated cell growth and division, provided the oxygen/nutrients concentration in the matrix doesn't decrease below a certain threshold (oxygen saturation level must be ⩾0.9) The growth rate value of CGs was set by matching simulation results of our in silico hybrid model to in vitro data of cancerous cells growth rate [28]; however, cancer cell mitosis is modelled as a stochastic event governed by a probability density function having uniform distribution.
The increase of the CG population size leads to a decline of oxygen concentration in the medium, resulting in a hypoxic microenvironment, since blood vessels oxygen support is insufficient. This model assumes host cells (Ne, G) incapable to survive under hypoxia, thus, as explained above, they diein other words they are simply removed from the ABM solver. The present model, however, can assume that only cancerous cells (normoxic or hypoxic) can transform permanently into necrotic cells (NC). However, parts of the tissue with sustained hypoxia may lead to anoxia where cancerous cells can survive and they become necrotic, whereas adjacent blood vessels seize functioning and they stop providing any oxygen/ nutrients (i.e., production rate parameter in Eq. (1) is set: . Also, CGs under hypoxic conditions (hCG) can express a migratory phenotype dictated by local gradients of the oxygen saturation level (chemotaxis), while ECM microstructural properties become important to the cells' preference to migrate (durotaxis).
In this model, tumour cells degrade the ECM locally (as outlined briefly in Section 2.1) and hence they prefer moving towards a denser matrix. Through transdifferentiation of native fibroblasts to CAFs, CGs secrete matrix metalloproteinases (MMP) that in turn degrade the ECM by breaking the structural fibres of the stroma (collagen, fibronectin)resulting in a heterogeneous ECM. Relevant observations have been reported in gliomas in vitro and in vivo [29][30][31]. In addition, CGs through CAFs also discharge platelet-derived growth factors (PDGF) for vasculogenesis, which subsequently fuel the unregulated proliferation of CGs [32,33]. Transition between cell states (normox-ia↔hypoxia→necrosis) is modelled as a stochastic event depending on a given cell's internally encoded rules in interaction with the microenvironment. Due to the motile phenotype of some cells, cells are bound to interact with one another mechanically which, at each time point of the simulation, is computed numerically via the repulsion force between any two adjacent (neighbour) cells. The resultant force acting on a given cell is evaluated by taking the vector sum of the concentrated forces originating from all neighbouring cells. Subsequently, the relative displacement of cells is then computed according to Newtonian mechanical dynamics, where inertial force effects are assumed negligible compared to viscous forces [34].
The agent-based model has been implemented using the BioDynaMo platform [12]. BioDynaMo allows to simulate cell mechanics and their interactions in 3D, while enabling the user to specify various model properties such as the genetic "rules", diffusable substance properties and mechanical forces.

Hybrid macro-to-micro model
The in silico hybrid multiscale cancer model couples the macroscopic scale (i.e., the FE mesh representing the tissue and the corresponding macroscopic state-variables; see Section 2.1) with the microscopic scale (of the agent-based discrete model representing the cells and the extracellular space) using the representative averaging volume (RAV) technique. The RAV technique has been successfully applied to date studying various problems in porous media, including cancer modelling [14]; however, a thorough treatise of volume averaging methods can be found in manuscripts [35,36]. Briefly, each quadrature point of a (macroscopic) finite element 1 is associated with a unique representative volume, the RAV, which contains a cloud of randomly distributed points: the cells. For simplicity, the RAV takes a cubic shape, L Δ 3 , having dimensions as described in Section 2.2 and oriented to align with the global frame of reference of the macroscopic domain. Initial estimates for the amount of each cell population, their physiological condition and state (i.e. oxygen levels, chemical cues, etc.) at each representative volume are determined by linear interpolation of the nodal (macroscopic) data to every Gauss point within the finite element (i.e., ∊ ξ g μ c , , , , j s i ), see Fig. 2B. In turn, these nodal data are determined by the macroscopic initial conditions imposed at the beginning of the simulation, or by the applied boundary conditions (where applicable), or obtained from the corresponding nodal data evaluated from the previous time increment of the FE solver, FEB3. Subsequently, the cells' dynamics within each RAV is simulated by the agent-based modelling solveropposed to employing a constitutive equation or any other sort of a deterministic modelat each quadrature point. After several internal iterations, of fixed time step T Δ , of the microscopic/cell-level simulator BioDynaMo, cell populations C i are re-evaluated for the various species (see Fig. 2A; the neurons, Ne, the glial cells, G, the normoxic and hypoxic cancerous glials, CG and hCG, and the necrotic cells, NC), and for each representative volume. Subsequently, the average cell rate at each RAV and the average direction vector of the cells escaping the bounds of the representative (cubic) volume is calculatedthe latter computations being applicable to cells with a motile phenotype only. These data are evaluated for each species C i and then up-scaled back to the macroscopic (tissue-level) FE simulator, FEB3, through quantities α i and ∂ ∂ c X / i respectively. More specifically, at each RAV and for each cell species, the average population increase (or decrease), C Δ i , is evaluated between successive increments of the macroscopic FE solver, separated by a fixed time step t Δ , and the rate is . Also, between successive increments of the macroscopic FE solver, the net flux of inwards or/and outwards cell migration is computed for each face of the RAV, C Δ i , thus, giving the local gradient (that describes the cells migration potential) from C L Δ /Δ i . In summary, as illustrated in Fig. 2B, the simulation procedure of the proposed in silico hybrid procedure is performed in a partitioned manner, with FEM and ABM working alternatively and passing necessary data interchangeably by projection and averaging until all simulation steps have been completed.
The proposed in silico cancer model distinguishes two layers of implementation that are identified from the two spatial scales involved. The agent-based modelling (ABM) method has been implemented in a multi-threaded, object-oriented, C++ code and incorporated in the open-source framework BioDynaMo. BioDynaMo employs several opensource libraries: ROOT, ParaView, and OpenMP. Moreover, the multiscale finite element procedure has been implemented in a scalable C+ + code, and incorporated into the existing numerical analysis framework FEB3. FEB3 (pronounced Phoebe) is founded on several highperformance, open-source numerical libraries: PETSc, libMesh, GSL, and blitz++. FEB3 has been designed in an object-oriented manner and facilitates parallel computations using the message passing interface technology of the MPICH library. The tissue domain FE three-dimensional meshes have been generated using Gmsh, and decomposed and distributed across multiple processors using the ParMETIS library. The coupling and communication between the (tissue-level) FE solver and the (cell-level) ABM solver, BioDynaMo, has been fully facilitated within FEB3. All simulations presented in this work were carried out on a desktop machine having an Intel i9-7940X CPU (@3.1 GHz ×14) and 128 GB main memory, operating Linux (Ubuntu 18.04, kernel version: 4.15.0-54-generic). Simulations of a 30 days tumour development (44,000 simulation steps) required on average 4 h to be completed while the main memory usage was less than 3 GB. However, it is worth highlighting here that the code of our solvers was not fully optimisedfuture code development could yield further improvements in efficacy. Furthermore, some numerical libraries used by FEB3 and BioDynaMo were built without optimisation flags.

Convergence and sensitivity analysis
A convergence and sensitivity analysis of the multiscale methodology is presented here. This involved two separate sets of simulation tests: (i) In the first set, the finite element (FE) mesh density 2 varied while the size of the representative averaging volume (RAV) of the microscopic solver was fixed. Specifically, the coarsest FE meshes consisted of 1626 elements and 425 nodes, and 4322 elements and 768 nodes respectively, the medium density FE mesh consisted of 3628 elements and 875 nodes, whereas the densest meshes consisted of 5526 elements and 1276 nodes, and 7594 elements and 1387 nodes respectively.
The purpose of this parametric analysis was to investigate the stability and convergence of the present multiscale framework with respect to the finite element discretisation.(ii) In the second set, the size of the RAV varied: (80 μm) 3 , (100 μm) 3 and (120 μm) 3 , while the size of the FE mesh was kept fixed. Thus, the purpose of this parametric analysis was to investigate the sensitivity of the hybrid FEM/ABM procedure to the scaling effects of averaging volume.
It is important to note here that in all simulation experiments, the initial cell density for each population (G: 5 × 10 4 cells/mm 3 and Ne: 10 6 cells/mm 3 everywhere in the domain of analysis; CG: 10 4 cells/ mm 3 were distributed spherically within a millimetre radius; zero for all other cell phenotypes) and macroscopic state variable (O 2 saturation level set to unity; ECM volume fraction set to vary 0.299 ± 0.1) was considered the same, while the model parameters associated with the macroscopic solver (oxygen uptake and vascular supply rate, growth factors rates, etc.) were identical and the solution and space/time discretisation settings for the FE solver were the same throughout the analysis. Also, for the purposes of this analysis, the model parameters have been selected to simulate a slow tumour growth condition in order to make cell-scale phenomena (i.e. low duplication rate, suppressed cancer cell migration) easier to monitor within the adopted time frame.
Cut-through FE meshes of the (2 cm) 3 domain of analysis is depicted at the middle row in Fig. 3; whereas at the bottom row of the same figure, the representative averaging volumes is shown as it increases from 5.12 × 10 5 μm 3 left to 17.28 × 10 5 μm 3 in right. All three figures of the RAVs are produced by taking a cut-through of the same FE mesh, while the finite elements are made here transparent for illustration purposes.
For the first simulation test, Fig. 3A shows the tumour volume growth (in mm 3 ) over a period of thirty days while Fig. 3B shows the evolution of the total tumour cell density (the population is calculated by taking the sum of CG, hCG and NC cell densities). Evidently from the former plot, having fixed all model parameters, the numerically predicted tumour volume appears to converge towards 750 mm 3 approximately as the mesh density increases. Despite of this, no clear conclusion can be made from the time plots in Fig. 3B where a relative difference between consecutive simulation end-point tumour cell populations is observed at about 5,000 cells/mm 3 . Videos 1-3 from the supplementary material animate (in clockwise order) using colour contours the cell density for CG, hCG and NC as well as the tumour cell density (in ×10 3 cells/mm 3 ) and the distribution of the oxygen level for the FE (macroscopic) solver. Also, in the same videos, the ABM (microscopic) solver simulation results of a selection of representative averaging volumes is shown where different colours depict the cell phenotypes: dark blue for the necrotic cells, light blue for neurons, green for glial cells, and orange and red for normoxic and hypoxic cancerous glial cells.
Furthermore, the computational time of the multiscale methodology has been assessed with respect to the spatial resolution (FE mesh density), while keeping the time resolution fixed (i.e., the time integration step, T Δ ) and the ABM parameters also fixed. The computational times recorded for each FE mesh discretisation are illustrated in a plot in Section 1 of the Supplementary Material. As expected, the computational cost increases with the density of the mesh (of the macroscopic domain) which in turn translates in an increase to the amount of ABM simulations required at the cell level. However, more importantly, we observe that the ABM solver took up to approximately 70%-90% of the total computational time required by the hybrid (multiscale FE-ABM) simulator.
For the second simulation test, Fig. 3C shows the tumour volume growth time plots for the three RAV sizes considered while Fig. 3D shows the evolution of the total tumour cell density for the corresponding averaging volumes. From Fig. 3C it can be seen that with a 25% increase of the RAV size (from 80 μm to 100 μm RAV side length) the final tumour size prediction increases by 66%, while a further 20% RAV increase (from 100 μm to 120 μm) gives an approximate 13% increase in the final tumour prediction. Similarly in Fig. 3D, we observe that the predicted tumour cells population prediction becomes less sensitive to the size of the microscopic domain of analysis, in other words the predicted values converge with larger RAVs.
Finally, Videos 4-6 from the supplementary material show animations of the simulation results for different averaging volume sizes, which are evident by comparing amongst visualisations 4A, 5A and 6A. Overall, the sensitivity analysis demonstrates convergence of simulation results in these exemplary cases. Therefore, in summary of the convergence and sensitivity analysis results above, a RAV size of (100 μm) 3 with a moderate FE mesh density of (3628 elements and 875 nodes) has be selected for the brain glioma growth simulations. This combination provides the best compromise between spatial resolution and computational cost that permitted us to conduct in-depth investigations of tumour growth.

Low-grade tumour growth
Using our in silico hybrid model, we simulated the development of a low-grade tumour over the period of a month. These low-grade tumours are well-differentiated and characterised by a slower growth rate compared to high-grade tumours. Thus, in this simulation, cancerous cells express a slow growth rate and no migratory behaviour. In the resulting simulation, tumour expansion is observed as shown in Fig. 4B, where the tumour diameter increases from about 37 mm at day 10 to about 64 mm at day 30. Similarly, the total tumour (TC = {CG + hCG + NC}) volume increases from 180 mm 3 ( ± 7.1 mm 3 ) at day 10 of the simulation to 784 mm 3 ( ± 28 mm 3 ) at day 30, as shown in Fig. 4C. This tumour growth is supported by cell divisions of the cells in the microscopic domains, and the cumulative number of cells composing the tumour mass in all microscopic domains rises up to 54,800 ( ± 1120 cells/mm 3 ) at the end of the simulation. This link between macroscopic domain and individual cell elements is represented in Fig. 4A. As CGs proliferate, and as we assume a uniform microvascular network providing a uniform source of oxygen, an hypoxic region emerge at the centre of the tumour (Fig. 4D), leading to CGs expressing hCGs phenotype. Later on, as oxygen becomes scarcer, hCGs survival becomes impossible thus creating a necrotic centre, as shown in Fig. 4G. This oxygen deprivation also impacts healthy hosts (Ne and G), whose densities decrease as they get closer to the tumour centre, even reaching a density of 0 at the centre. These characteristics can be observed at both macroscopic and microscopic scale in Fig. 4A, where the tumour centrecomposed of NC (coloured in blue) and hCGs (coloured in red)is surrounded by CGs (cells coloured in orange). The tumour obtained at the end of this simulation is shaped as a spheroid and clearly defined, with no isolated CGs outside of the tumour mass as shown in Fig. 4D.

High-grade tumour growth
High-grade brain tumours are more malignant tumours and thus have worse prognosis than low-grade tumours. They are known to have a high regrowth behaviour, even after complete surgical removal of the tumour mass, in some cases due to their more diffused organisation. In order to simulate such more aggressive tumour development, three modifications have been made to our model, the first one being the secretion of MMPs by CAFs. These enzymes, by degrading the ECM, facilitate cell migration and so give rise to more invasive behaviour of cancerous cells. Secretion of a growth factor has also been added, recapitulating the beneficial effect of PDGF on tumour growth. Finally, in this simulation cancerous cells adopt a "go-or-grow" phenotype [37], where hCGs express a high migratory behaviour without the possibility to grow or divide, while CGs grow and divide, without exhibiting any migratory behaviour. In other words, when cancerous cells are in an hypoxic environment, they express a migratory behaviour. Oxygen supply is assumed to be uniform across all the tissue. Simulation results for this case are illustrated in Fig. 5 and in Video 7 of the supplementary material. During tumour development, we can observe that the early growth behaviour, and the resulting shapes, are still similar between low-grade and high-grade tumours, as they have a well defined spherical shape (Fig. 4B day 10 and Fig. 5B day 10 respectively). Likewise, their volumes are still in the same range, being 151 mm 3 and 212 mm 3 , respectively. This means that under the specified conditions, high-grade tumours exhibit a volume 1.4 times bigger than low-grade ones.
However, differences emerge as the development continues, with a bigger volume difference between low and high-grade at day 20 (2.1 times bigger for the high-grade) and day 30 (2.4 times bigger for the high-grade), tumours reaching a final volume of respectively 784 mm 3 and 1890 mm 3 after 30 days of development.
Moreover, differences in shape can be observed at day 30, highgrade tumours exhibiting a more diffused phenotype with outgrowth of lower density at the surrounding of the main tumour mass, while lowgrade tumours remain well defined and more spherical (Fig. 4B and Fig. 5B).
In addition, healthy cells (glials and neurons) are significantly affected by the tumour, as cancerous cells spread covering a larger area when compared to the low-grade case (see Figs. 5B,D and 4B,D respectively). As evident in Fig. 5C, due to the pronounced hypoxia and anoxia, healthy cells within the tumour region are susceptible to abrupt apoptosis. Fig. 5G, depict the absence of cells from both healthy populations respectively. Finally, secretion of MMPs by the proliferative CGs at the rim of the tumour (Fig. 5D) lead to degradation of the ECM both at the tumour localisation and surrounding (not shown here). As a consequence, these changes of the ECM density facilitate migration of cancerous cells, hence, the expansion and spread out of the glioma. All these factors give rise to more aggressive tumours, exhibiting high growth rate and migratory behaviour, with a more diffused phenotype characterised by outgrowth of low CGs density at the surrounding of the main tumour mass.

PDGF secretion promotes tumour development and invasive behaviour
As experimental studies showed the promoting effect of PDGF on tumour growth [32,33], we investigated in our model the role of PDGF on tumour development. To achieve this, we modified the cellular response to growth factor. Two simulation conditions are considered, with or without beneficial effect of growth factor on cancerous cells (CG and hCG), all other parameters remaining identical. More concretely, this growth factor enhances the survival, growth and division behaviour of cancerous cells.
Surprisingly, although the obtained tumours exhibit a similar volume at day 10 to the default development condition, important differences emerge later on. Indeed, we can notice a smaller volume at day 30 for the tumour obtained without growth factor beneficial effect than the tumour with growth factor beneficial effect (1300 mm 3 and 1890 mm 3 respectively). These results are summarised in Fig. 6A. In addition, we can notice that the absence of the growth factor effect, the tumours are more homogeneous. In particular, their shapes are more spherical, with less outgrowth of low CG density around the tumour mass as shown by Fig. 6B. Thereby, the promoting effect of PDGF on cancerous cells gives rise to a bigger and more aggressive tumour, exhibiting more heterogeneity and a quicker development, and so worsens patients prognostic outcomes.
In addition, we demonstrate here that the shape and characteristics of the tumours are strongly influenced by the behaviour of individual cells. Results of the simulation without the promoting impact of the growth factor on cancerous cells are presented in Video 8 in the supplementary material. Moreover, we have also investigated the sensitivity of the tumour growth to changes of hCG migration direction, dependent upon the ECM and oxygen gradients (Supplementary Material: Section 2 and Video 9). We find that tumour shape and volume are not strongly dependent on such changes.

Discussion
This paper presents a novel in silico hybrid modelling procedure that couples tissue biomechanics using the finite element method and cell mechanics via the agent-based method. To achieve this, we have integrated two open-source numerical analysis platforms: FEB3 and BioDynaMo respectively. Although the hybrid continuum-/agent-based modelling concept is not new [19], to the best of our knowledge, this is the first that the two modelling methods are coupled using the volume averaging technique.
The proposed in silico hybrid platform demonstrates great potential to simulate the development of brain gliomas in a detailed but computationally cost-effective manner. For the working cases presented in this paper, each simulation took about four hours to finish on a multicore desktop machine. In particular, the ABM solver required the simulation of~5.5 × 10 5 -7 × 10 5 agents coupled via the RAV technique with a FE mesh (of a 8 cm 3 tissue volume) consisting of~4 × 10 3 tetrahedral elements in total. Had we employed ABM alone to simulate glioma growth in a domain of the same size, it would have required for the solver to simulate approximately~1.3 billion agentsthe scale of data involved in such an analysis is numerically prohibitive even for a moderate sized computer cluster. Furthermore, the proposed in silico hybrid platform was capable to recapitulate the dynamics of a growing tumour in the brain that may exhibit either a benign or an aggressive behaviour. In addition to this, the platform recapitulated the main points of brain tumour necrosis and oxygen deprivation, as well as the dynamics of a heterogeneous developing tumour matrix.
In spite of the successes demonstrated, there is significant scope for improvement and alleviating model simplifications involved in the proposed in silico hybrid platform. For instance, we plan to circumvent the existing coarse-grained description of the extracellular matrix heterogeneity, and incorporate explicitly the ECM solid and fluid mechanics directly at the microscale. To achieve this, the present model will be coupled with our biphasic multiscale model that encompasses the matrix microstructure and interstitial fluid (micro) flow [14], and therefore directly model the biomechanics of proteoglycans (i.e. hyaluronan and collagen) and the parenchyma.
Moreover, a natural extension of this hybrid model is to include tumour-induced angiogenesis, which will be accomplished by plugging the proposed procedure in the coupled angiogenesis model we have recently developed [13]. Moreover, our hybrid model can be further used as a platform for the study of the efficacy and delivery of drugs to tumour spheroids in in vitro experiments. This can be achieved by adding extra equations in the macroscopic solver FEB3 and, hence, model the balance of the therapeutic agents at the interstitiumit can be seamlessly accomplished by plugging into our modular reaction-diffusion solver the extra PDEs required [15]. In a similar fashion to adding more biochemical factors and drugs into the  simulator, one could also incorporate extra agents in BioDynaMo to introduce more cell types in the in silico experiments. Last but not least, our hybrid model offers great promise being used for the design and development of organ-on-chip devices in cancer, neurodegenerative diseases, tissue regeneration, etc.
To conclude, our simulations yield differences in the spatial distribution of CGs relative to the core of the tumour. These results are explanatory powerful and yield experimentally verifiable predictions. Histological analysis of the density of cancerous and host cells in the adjacent environment of the tumour could provide further support for these simulation results. The possibility for hypothesis generation further supports the case of mechanistic computer models to complement experimental cancer research. Ultimately, such models could generate generate predictions of tumour progression to inform on surgical and/ or radiotherapy planning [38,39]. The practicality of the proposed hybrid method is an important feature for such a clinical setting.
We anticipate our modelling approach can be further employed for a wide range of biomedical problems, other than cancer. In particular, models of healthy biological development, biofilm growth, immune system dynamics, or synthetic tissue growth are based on relevant dynamics.

Data availability
All data and computer codes used in the simulations of this paper are freely available through the figshare repository: In silico hybrid FEM-ABM model (https://figshare.com/projects/In_silico_hybrid_FEM-ABM_ model/67343).