Data-driven spatio-temporal modelling of glioblastoma

Mathematical oncology provides unique and invaluable insights into tumour growth on both the microscopic and macroscopic levels. This review presents state-of-the-art modelling techniques and focuses on their role in understanding glioblastoma, a malignant form of brain cancer. For each approach, we summarize the scope, drawbacks and assets. We highlight the potential clinical applications of each modelling technique and discuss the connections between the mathematical models and the molecular and imaging data used to inform them. By doing so, we aim to prime cancer researchers with current and emerging computational tools for understanding tumour progression. By providing an in-depth picture of the different modelling techniques, we also aim to assist researchers who seek to build and develop their own models and the associated inference frameworks. Our article thus strikes a unique balance. On the one hand, we provide a comprehensive overview of the available modelling techniques and their applications, including key mathematical expressions. On the other hand, the content is accessible to mathematicians and biomedical scientists alike to accommodate the interdisciplinary nature of cancer research.


Introduction
Glioblastoma (GBM) is a malignant hierarchically organised brain cancer.It is both the most common and most aggressive type of primary brain cancer in adults [44].Not only does the diffusive invasion of glioma cancer cells into healthy tissue impede complete resection, but GBM harbours a subpopulation of highly therapy-resistant stem-like cells [113].Tumour recurrence is inevitable, resulting in a median survival time of 15 months for patients despite maximal treatment [191].The current gold standard of treatment is the Stupp protocol, which consists of maximal safe surgical resection, followed by radiotherapy and chemotherapy with temozolomide, an alkylating agent [174; 173].Mathematical cancer models have provided a deeper understanding of this immensely complex disease by unveiling the underlying mechanisms and offering quantitative insights [6; 3; 5; 37].Such models have entered all areas of GBM research, ranging from the classification and detection of brain tumours to therapy [207; 59].
This review provides the reader with an overview of existing mathematical and computational models that aim to simulate spatially resolved tumour growth.We discuss three main paradigms that have emerged for such in silico experiments.In Section 2, we introduce so-called continuum models that treat variables, such as the tumour cell density, as continuous macroscopic quantities based on conservation laws.Alternatively, one might represent each cell as an individual agent.Such discrete models are discussed in Section 3. Section 4 deals with hybrid multi-scale and multi-resolution models that merge and bridge different approaches.
Each of the three methods has its advantages and shortcomings.One must choose between them based on computational limitations and the level of detail required to answer the research questions of interest.With this review, we aim to assist researchers in choosing between the different methods by highlighting the drawbacks and assets of each approach and by showing how the different methods can complement each other.Moreover, we summarise the main concepts and the key mathematical expressions that lie at the core of each approach.We hereby aim to strike a balance between providing a brief overview and showing the mathematics involved.
Mathematical models of GBM draw on a wide variety of molecular and imaging data: histopathological data, computerised tomography (CT), positron emission tomography (PET), single-photon emission computerised tomography (SPECT), and magnetic resonance imaging (MRI), such as T1 weighted (+/-Gadolinium contrast), T2 weighted, T2-FLAIR, and diffusion-weighted imaging (DWI), and most recently spatial and single cell transcriptomics.The models have thus been employed to shed light on patient-specific data in vivo and ex vivo as well as on data from animal models and in vitro experiments.These analyses have provided invaluable insights and have deepened our understanding of glioma on a molecular and structural level.We address this issue in more detail in Section 5.The section also discusses the computational challenges posed by systematic inference of model parameters.Finally, Sections 6 and 7 provides a short overview of some clinical applications of these models and an overall summary, respectively.
While we discuss the different methods in light of GBM, we note that the same models are applied to other types of cancer.Indeed, the models build on concepts that are broadly used to study tissues.We do, therefore, not only cite sources that deal with GBM but occasionally refer the reader to illustrative papers from other areas of oncology and biology (see also [89; 99]).It is worth noting that the models are, in a broader sense, actually widely used across scientific disciplines.Throughout the paper, we thus present ideas and concepts that are also employed in other fields ranging from statistical mechanics to solid-state physics.For instance, the cellular Potts model presented in Section 3.1.3builds on the so-called Ising model used to describe ferromagnetism.We hence encourage the reader to think outside the box when exploring the literature, and, in this spirit, we provide a few citations to areas outside the realm of biology.
Before commencing, we would like to point the reader towards other reviews and papers for further details.
Lowengrub et al. [120] provide a detailed account of continuum models.For an elaborate discussion on discrete models, we refer the reader to Van Liedekerke et al. [116].Both Metzcar et al. [132] and Weeransinghe et al. [195] give a brief overview of the topic and include a list of recent references.For insights into hybrid multi-scale modelling, we recommend Deisboeck et al. [45] and Chamseddine & Rejniak [28].Falco et al. [58] present a concise overview highlighting their clinical implications.Ellis et al. [50] focus on mathematical models that intratumour heterogeneity and tumour recurrence based on next-generation sequencing techniques.Alfonso et al. [1] highlight the challenges that mathematical models face when dealing with glioma invasion.Finally, for an overview of the biology of the GBM from a clinical perspective, we refer the reader to the recent reviews by Finch et al. and McKinnon et al. [59;129].

Continuum models
When dealing with cancer treatment, we face questions related to tumour size, shape and composition.These questions all address cancer on a macroscopic scale.While macroscopic tumour dynamics emerge from interactions on a cellular level, it is possible to construct informative mathematical models without tracking individual cancer cells.Instead, the tumour and its environment can be represented as continuous variables that are governed by partial differential equations (PDE).Such continuum models capture many aspects of cancer in vitro, in vivo, and in patients.They can account for the impact of heterogeneous brain tissue on tumour growth, for different invasive tumour morphologies, and to some extent, even for potential tumour recurrence [e.g.62 ; 176].Moreover, they have successfully been applied in studies on the impact of chemotherapy and repeated immuno-suppression treatment [178; 63].Continuum models have thus been employed in diagnosis and treatment planning based on patient-specific data [e.g.177; 75; 35; 135; 86].However, PDEs smooth out small-scale fluctuations, which implies that continuum models do not apply to small cell populations, such as those found in the tumour margin.For such small cell numbers, stochastic events play a crucial role, and the applicability and predictive power of continuum models are limited.To properly understand cancer invasion, we must track individual cells.But to do so comes at a high or currently insurmountable computational cost (see Section 3).Thus the use of continuum models represents a trade-off that enables and supports scalability and mathematical insights.As a result, continuum models are widely used in the community [193; 120; 195].This section introduces the basic mathematical concepts of continuum models and their biological motivation.
Any such model in the literature builds on reaction-diffusion equations, describing variables such as the tumour cell density, the tumour volume fraction, the nutrient (oxygen, glucose) concentration, the neovasculature, the enzyme concentration, or other properties of the extracellular matrix (ECM) [e.g.211; 63; 42; 148; 66].Considering any such variable, ψ(x, t), which is a function of position x and time t, we have for its rate of change with time where J is the flux of the considered variable, while S is the sources and sinks for this variable.Thus, Equation (1) constitutes a conservation law.The exact expression for J and S, as well as the boundary conditions, will depend on the variable in question.For instance, when dealing with the quasi-steady diffusion of nutrients, Equation (1) generally takes the form [39; 168; 38; 124; 25] Here, D denotes a diffusion coefficient, while n(x) is the relevant nutrient concentration at the location x within the considered n-dimensional domain, which is oftentimes denoted by Ω.
Alternatively, let's consider the (normalized) cancer cell density, ρ(x, t), at a time t and location x.Many authors assume that the diffusion of cancer cells can be well-approximated by Fick's first law where D denotes a diffusion coefficient, which we discuss in detail in Section 2.1.As regards the sources and sinks of the cancer cell density, it is commonly assumed that cell proliferation, i.e. tumour growth, is well-described by a logistic growth term [e.g.175 ; 161].So, Equation (1) takes the form where λ is the growth rate of the tumour cell population.Other authors assume exponential tumour growth, substituting the second term on the right-hand side by λρ [e.g.177].Of course, more than one term might be necessary to summarise the relevant sources and sinks.Several more complex terms are, for instance, needed when considering differentiation between different interdependent tumour sub-populations (cf.Section 2.2).
It is worth stressing that Equation (1) cannot stand on its own.Other relations and constraints, including

Anisotropic Diffusion
By using a scalar for the diffusion coefficient in Equation ( 4), we assume isotropic tumour growth.However, GBM spreads anisotropically, primarily expanding along pre-existing structures, such as blood vessels and white matter tracks [164; 70; 153; 110; 74; 40; 22].To take the heterogeneous structure of the brain into account, Swanson et al.
[177] hence proposed to adopt different values for the diffusion coefficients in grey and white matter.Concretely, based on CT scans by Tracqui et al. [183], Swanson et al. [177] found the diffusion coefficient in white matter to be more than five times larger than in grey matter.
Taking the idea of heterogeneous diffusion further by including anisotropy, other authors [e.g.90; 84; 110; 143; 9] substitute the diffusion coefficient with an n-dimensional diffusion tensor, D(x, t) ∈ R n×n .To evaluate D(x, t), Painter & Hillen [143] deploy diffusion tensor imaging (DTI) data.DTI is a magnetic resonance imaging (MRI) technique that measures the anisotropic diffusion of water molecules and hereby maps highly structured tissue.
This technique provides the diffusion tensor for water molecules, DΠ(x, t), throughout the brain.Of course, due to the size difference, the movement of cancer cells is more restricted than that of water molecules, which means that DΠ(x, t) does not adequately describe glioma growth, i.e.D(x, t) = DΠ(x, t).However, based on a transport equation for individual cell movement, Painter & Hillen [143] establish a relation between D(x, t) and DΠ(x, t) expressed in terms of the Fractional Anisotropy (FA) that is commonly used to quantify DTI data [16].
where the diffusion tensor D ∈ R n×n is symmetric and positive-definite, as it is related to the variance-covariance matrix of the probability distribution function that describes the velocity changes of individual cells [see 82 ; 143].
In other words, when dealing with three-dimension data, D(x, t) is a symmetric and positive-definite 3 × 3 matrix that incorporates the impact of the local environment on cell migration.The model by Painter & Hillen [143] has been extended and applied by other authors [53; 54; 176].
Models of the form of Equation ( 6) are referred to as Fokker-Planck models, while Equation ( 4) is an example of a Fickian model.The additional advection term of the Fokker-Planck model has a demonstrable impact on the solution [see also 15].We also note that Equation 6 would correspond to the Fisher's equation if the first term on the right-hand side were substituted by D∇ 2 ρ.Indeed, some authors employ a diffusion term of this kind to describe the cancer cell density [114].

Mechanical interactions, cell types, lineage, and feedback
Tumours are hierarchically organised.GBMs harbour stem-like cells (GSC) as well as proliferating (GCP) and differentiating (GTP) subpopulations [cf.169; 188; 113; 22].Since these three cell types exhibit very different behaviours, it is insightful to differentiate between these subpopulations, even when dealing with continuum models.
Models that distinguish between viable and necrotic tumour tissue are the first step in this direction.Examples of such models can be found in the papers by Wise et al. [199;198] and Frieboes et al. [26; 60].Their work is based on reaction-diffusion equations of the form where the index i runs over all subpopulations, ui denotes the velocity of the considered cell species, and Jmec,i is the flux that arises from mechanical interactions Note that while Equation (7) appears to differ from the other reaction-diffusion equations listed above, this is merely a matter of notation.It's a Fickian model that can be derived by inserting Equations ( 3) and (8) into Equation (1).The advantage of phrasing the problem in this manner is that the mechanical flux reflects the mechanical interaction energy that can be obtained from an understanding of the underlying cell biology [see also 101; 200; 102; 62].By introducing Jmec,i, it is thus possible to inform the model about the properties of the tumour and host without the necessity of constructing a suitable diffusion tensor.
In papers that employ the Equation ( 8), ui is computed by imposing relations similar to Darcy's law that links the velocity to a gradient in pressure [e.g.62].For glioma, the relevant expression often takes the form ui = −∇p + Fi, where p is the solid pressure arising from the tumour proliferation, whilst Fi reflects mechanical interactions.
We exemplify the source functions that enter Equation ( 7) by listing the relevant terms for the necrotic tissue according to Wise et al. [199], for which Here, the indices 'd' and 'v' refer to the dead and viable cancer cells, respectively, while λA, λN, and λC denote the rates of apoptosis, necrosis, and the clearance of dead cells.Moreover, H is a Heaviside step function, and nN is a viability limit for the nutrient concentration below which cells die.For comparison, the source function for the viable tissue takes a similar form: where λM denotes the rate of mitosis, and n∞ is the far-field nutrient level.
By distinguishing between GSCs, GCPs, and GTPs, Kunche et al. [112] and Yan et al. [204; 205; 203] have incorporated the GBM lineage and hereby taken the discussed models one step further.Kunche et al. [112] use this to investigate feedback regulation of cell lineage progression.Furthermore, while previous papers only consider adhesion when computing the mechanical interactions, Chen et al. [33; 31; 32] include the impact of elastic membranes and the implications of the calcification of dead tumour cells [144].

Modelling the macroscopic environment
Understanding the microenvironment is essential for understanding GBM since the brain region and other properties, such as the patient's age, have been shown to play a key role in tumour development and heterogeneity [e.g.22; 154].
The concentrations of different chemicals, including (but not limited to) nutrients, drugs, matrix-degrading enzymes, and ECM macromolecules, are commonly modelled using reaction-diffusion equations.The associated diffusion is often in the form D∇ 2 ψ, but other second-order spatial derivatives can be found in the literature [7; 92; 147].The source terms reflect the processes at play.For instance, the rate of oxygen consumption is often assumed to be proportional to the local oxygen concentration and might be proportional to the local cell density.
Overall, the use of continuum models to model, say, nutrient flows is justified through the different scales that define cancer.Continuum models reliably capture the spatial gradients and temporal variation of nutrients in a way that allows us to assess the behaviour of tumour cells.For some purposes, it might even be adequate to assume that the nutrient levels stay constant over time since cell proliferation takes place on a much longer time scale than nutrient diffusion (cf.Equation 2).
Initially, tumours do not possess their own vasculature but rather rely on the diffusion of nutrients and waste products.During this so-called avascular phase, the tumour may lack some features of malignancy and appear more benign.But beyond a certain size (1-3 mm in radius), the nutrient inflow can no longer sustain the growing cell population [126].Hypoxia, i.e. oxygen shortage, sets in.Without stimulating angiogenesis, i.e. recruiting vasculature, tumour growth will stagnate.However, in response to the hypoxic conditions, the cancer cells release tumour angiogenic factors (TAFs), which stimulates the migration and proliferation of endothelial cells (ECs).New blood vessels sprout towards the tumour, and tumour growth resumes.Understanding the transition from avascular to vascular tumour growth is essential since it is a critical step toward malignancy.To study angiogenesis, many authors model both the diffusion of TAFS and ECs using reaction-diffusion equations [e.g.
Table 1: Overview of the reviewed continuum models.The references are not exhaustive.

Model type Features References Continuum model (CM), tumour
Concept: The tumour cell density can be modelled as a fluid using PDEs based on conservation laws.Assets: CMs capture macroscopic tumour growth at low computational cost.Drawbacks: They cannot capture heterogeneity or the behaviour at low cell densities (e.g. at the tumour front).
The concentration of and sensitivity to the ith components are denoted as ci and χi, respectively [see 126].
Like continuum models for tumour cells, continuum models for ECs cannot capture the behaviour of individual cells.Some authors, therefore, rely on agent-based models to study the ECs during angiogenesis [166; 43; 202; 29].
We refer to Section 3 for a detailed discussion of agent-based models [see also 4].Table 1 summarises the properties of all continuum models discussed in this section.

Discrete (agent-based) models
Continuum models are based on the assumption that the behaviour of tumour cells is well-approximated by macroscopic averages.However, this assumption breaks down for small cell populations since stochastic events dominate these.Therefore, continuum models cannot reliably capture the onset of tumour growth, and they do not give a complete picture of the invasive tumour front or margin, i.e. the tumour-host interface, where the cancer cell density is low -although it is possible to address this drawback by linking different scales, as shown by Trucu et al. [184] (see also Section Discrete models that describe spatial tumour growth or other cell populations, such as ECs (cf.Section 2.3), fall into two categories: lattice-based and off-lattice methods.While the former operates on a fixed lattice, the latter does not impose the same restrictions.Lattice-based models can be divided into three sub-categories [116; 132].First, there exists a group of so-called cellular automata (CA) models that treat all processes, including cell movement, as stochastic processes (cf.Section 3.1.1).For these models, each lattice site (e.g.pixel or voxel) can either at most be occupied by a single cell or a population of a limited size [97; 96; 125; 165].Each cell or cell cluster is thus characterised by its position on the grid.Secondly, there exist so-called lattice gas cellular automata (LGCA).While these models consider single cells and are conceptually very similar, they differ from other CA models by attributing a velocity to each cell (cf.Section 3.1.2).After accounting for stochastic events, these models thus include a deterministic evaluation of cell propagation at every time step.Thirdly, there are socalled cellular Potts models (CPM) that phrase cell interactions within the framework of statistical mechanics (cf. Section 3.1.3).CPMs attribute several lattice points to each cell, whereby these models capture cell morphologies.
Off-lattice models, also known as lattice-free models, come in two flavours.Centre-based models (CBM) describe cells as simple particles whose interactions can be expressed as physical forces.Deformable cell models (DCM) or vertex models extend this picture by resolving the cells using several nodes.Like CPMs, DCMs can thereby account for cell morphologies.
One drawback of discrete models is their relatively high computational expense due to the large number of cells that typically need to be simulated.After all, tumours reach 10 5 − 10 6 cells per mm 3 .However, one benefit of discrete models is their ability to capture emergent properties.The simple rules that govern the (nonlinear) interactions between individual agents will naturally give rise to the dynamics of macroscopic tumour growth.As discussed in Section 3.1.2,reaction-diffusion equations might thus naturally follow from up-scaling agent-based models [80; 79; 87].Due to the ability of discrete models to capture emergent properties, agent-based models can be calibrated based on patient-specific data and be used to predict macroscopic parameters [122; 85].
Concerning the implementation of discrete models, it is worth mentioning that some software tools are publicly available.There is no need for a research group to start from scratch.For a recent overview of open-source toolkits, we kindly refer the reader to Table 1 in the paper by Metzcar et al. [132].In this connection, we note that one must always select and adapt the numerical approach based on the problem at hand.For instance, not all problems require the same geometry.While colon cancer cells and cells in vitro move in two dimensions, glioblastoma tumour growth in vivo is intrinsically three-dimensional, which affects the growth pattern [67].In the following two sections, we will discuss lattice-based models in more depth (Section 3.1) and then turn to a detailed discussion of off-lattice models (Section 3.2).

Lattice-based models
In lattice-based models, we follow a population of cells on a rigid grid.If a regular (e.g.Cartesian) lattice is chosen, artefacts in the established tumour growth pattern might arise, reflecting the lattice symmetries [116].
While unstructured meshes are harder to implement and combine with existing PDE solvers, they do not suffer from this shortcoming [18].

Cellular automata with one or several cells per lattice site
In this section, we discuss two types of cellular automata (CA) models: those with at most one and those with multiple cells per lattice site.Both approaches have their advantages and drawbacks.Below, we will start by discussing the single-cell CA models and then continue to frameworks that include many cells per lattice site.
When considering single tumour cells, we follow each of these through a sequence of time steps.During every time-step (t → t + ∆t), a cell might migrate to a neighbouring lattice site, it might die, it might divide, and it might grow, whereby it temporally occupies two sites before dividing.Each of these events can be modelled as stochastic processes or be imposed through deterministic conditional statements that encompass our knowledge about cancer -additional parameters, such as age, mutation, or phenotype, might be included in these rules [146].For instance, while cells might die spontaneously, external factors can make cell death increasingly likely.
Thus, tumours develop a necrotic core due to nutrient deprivation.One might, therefore, set the probability for individual cell death to depend on the distance to the tumour front [165], on the local nutrient concentration, or the local level of toxic metabolites [125].
Tumour cell migration is also affected by the location of the cell within the tumour.If we consider a cell that is deeply buried within the tumour, all surrounding lattice sites are already occupied, and there is nowhere to migrate.This notion also plays a role in cell proliferation.To divide, a cell deep within the tumour must push its neighbours away to make space for the daughter cell, which becomes increasingly difficult as the distance to the tumour edge increases.For simplicity, many authors account for this notion by introducing a sharp cut-off: If there are no free sites within a certain distance, the cell will be unable to divide [47].Other models assume a smooth decline in proliferation with depth [87; 88].On a macroscopic level, this suppression of cell division leads to a constant proliferative rim, below which we encounter a quiescent cell population surrounding a hypoxic region with a necrotic core at its centre.In accordance with data, this model predicts that a tumour spheroid enters linear growth after initially growing exponentially [83].As regards cell division, we also note that experiments suggest that the cell cycle time follows a Γ-like distribution, e.g. an Erlang distribution [152].The rules imposed on the individual cells can account for a broad range of cell properties.For instance, if the CA model accounts for the hierarchical nature of tumour cells (GSCs, GCPs and GTPs), the probability of cell proliferation will also depend on the attributed cell-type [146].
When solving the system numerically, one must choose a suitable time-step, ∆t.For this purpose, it is helpful to consider the cell population as a whole.Every time a cell dies, moves, divides, or grows, the system as a whole will transition into a new state.The time-step should be chosen accordingly such that only one event is likely to occur within ∆t, and thus ∆t is itself a function of time.To achieve this, one might distil the temporal dynamics of the system into a master equation and employ kinetic Monte Carlo algorithms [71; 20; 87].Moreover, the order in which the individual cells are updated should be random to avoid grid artefacts [120].
Rather than considering individual cells, one might track cell clusters of a limited size, keeping track of the number of cells at each lattice site.The underlying concepts stay the same.The advantage of reducing the resolution in this manner is the reduced computational cost, which allows for modelling tumours of centimetre size [152].Meanwhile, a coarser resolution might lead to artefacts affecting the predicted tumour growth.Another drawback of both single-and multiple-cell CA models is the fact that they do not properly account for cell morphologies.For models that address this issue, see Section 3.1.3.

Lattice gas cellular automata models
Like the models in Section 3.1.1,lattice gas cellular automata (LGCA) models follow a set of stochastic rules to determine whether individual cells die, grow, or divide.But in addition to specifying the location of the cells on the lattice, LGCA models attribute a discrete set of velocity channels to each spatial location.Hence, the cells roam a discretised phase space.
Analogously to migration in the CA models presented above, cells might stochastically migrate in phase-space: They can move to neighbouring velocity channels following probabilistic rules.At most, one cell can occupy any given velocity channel at a given location.This restriction is not to say that these channels are unique, i.e. at each location, there might be one or more channels with the same velocity, including channels at which the cell is at rest.

When using
LGCA models, the time evolution of the system is split into two distinct steps.First, the model allows for a number of stochastic events that correspond to the number of cells.Afterwards, the model updates the position of every cell based on their individual velocities complying with momentum and mass conservation.
In this manner, the model alternates between accounting for probabilistic interactions and deterministic cell movement.These basic concepts of LGCA models are summarised in Fig. 1.
Within the mean-field approximation, Hatzikirou & Deutsch [80] show that the Boltzmann equation governs the macroscopic dynamics of LGCA models (The Boltzmann equation can be used to describe the evolution of thermodynamic systems that are out of equilibrium.Applications range from quantum mechanics over cosmology to magnetic resonance imaging [e.g.46]).Furthermore, through Taylor expansions and scaling, they arrive at a reaction-diffusion equation with a logistic growth term, giving a microscopic justification for the success of the models presented in Section 2.
LGCA models offer an intuitive implementation of cell migration.By comparing model predictions to experimental data, the group of Deutsch and Hatzikirou has thus used such models to gain insights into cell-cell interactions and mechanisms that underlie tumour invasion [34; 182; 79].In particular, they focus on the socalled "Go or Grow" (GoG) hypothesis.This hypothesis addresses the fact that glioma cells exhibit an inverse correlation between cell migration and proliferation [69].It states that migration and proliferation are mutually exclusive events, i.e. that moving cells cannot divide.This dichotomy stems from shared signalling pathways [see also 72].Since the transition from benign to malignant tumour growth is coupled with a transition from a highly proliferative to a highly migrative phenotype, it is paramount to understand the link between these two phenomena.Specifically, the group of Deutsch and Hatzikirou investigate the role of hypoxia [see also 133].For simplicity, let's consider a Hamiltonian of the form [185; 166]

Cellular Potts models
The first term on the right-hand side of Equation ( 12) summarises the contact energy at the cell interfaces.The sum runs over all pairs of adjacent lattice sites (x,x').Each cell (and the ECM) has a unique identifier σ, i.e.
lattice sites with the same value of σ are associated with the same cell (or the ECM).Moreover, each cell has a type τ [e.g.GSC, GCP, or GTP, see 65].E(τ (σ(x)), τ (σ(x'))) denotes the contact energy per unit surface area (in 3D) between a cell of type τ (σ(x)) and a cell of type τ (σ(x')).The Kronecker delta, δ (σ(x), σ(x')), is needed since we only encounter an interface if σ(x) = σ(x'), while there will be an contact energy associated with cells where the Boltzmann temperature, β, describes the motility of the cells.
One can include additional terms in Equation ( 12) to add more information about cell biology.For instance, in analogy to the volume constraint, Ouchi et al. [142] suggest including constraints on the surface area of the cell (in 3D) [see also 19].Other authors include chemotaxis, motility, haptotaxis, and haptokinesis in this manner [163; 43; 179; 115].One immediate drawback of CPMs is that they are limited to phenomena and cell properties that can be formulated as terms in a Hamiltonian.
Meanwhile, when accounted for, cell division and death are treated as stochastic events or occur when certain conditions are met.For instance, Gao et al. [65]

Off-lattice models
One of the main drawbacks of lattice-based models is that the lattice introduces a lower length scale.Off-lattice models overcome this problem at the expense of increased computational cost.We give a schematic overview of the two different types of off-lattice models in Fig. 3.

Center-based models
Centre-based models (CBMs) treat cancer cells as physical particles whose trajectories can be deduced from their equations of motion [see 181].The cells themselves are usually represented as spheres (in 3D) or viscoelastic ellipsoids that deform when subjected to external forces [41].During mitosis, the mother cells are oftentimes modelled as dumbbells, consisting of two overlapping spherical cells [48].Constraints on areas and volumes might lead to additional forces (F i ), while the movement of the node (v i ) gives rise to viscous forces.The cell has a cytoskeleton (CSK, dashed lines).
The equations of motion can be derived within an energy-based picture or by modelling the relevant interactions as physical forces that act on the particles.In the latter approach, the models draw on Newton's second law and commonly assume that inertial forces and any average acceleration can be neglected (cf.Brownian dynamics).
Put another way, we are dealing with friction-dominated overdamped motion, which means that friction forces balance out all the other forces [47; 83; 116].For cell i on a substrate, we can then write [cf.116] The first term on the left-hand side expresses cell-substrate friction forces, while the second term deals with friction between cells.Here, Γ cs ji and Γ cc ji are tensors that account for anisotropies and inhomogeneities.For instance, for spherical cells in a homogeneous and isotropic environment, Γ cs ji = γI, where γ is a damping coefficient, and I is the identity matrix [47].Thus, the drag force that acts on a cell in an isotropic viscous environment is proportional to the cells' velocity [127].
On the right-hand side of Equation ( 14), F sub i denotes forces that emerge from adhesion and repulsion between cell i and the substrate.F mig i describes the migration forces.For instance, chemotaxis might introduce a preferred direction for migration [47].One might account for the memory of cells, i.e. persistence of a cell to move in a given direction.Some authors do so using inertial forces.In addition, some authors include cell polarity giving the cell a preferred direction of motion.F mig i becomes a noise term that describes a random movement in the absence of such forces and assumptions.The last term on the right-hand side summarises the adhesion and repulsion between cell i and all other cells.These forces will typically depend upon the distance between the cell centres or the amount of overlap, δji.To exemplify the use of Equation ( 14), let's consider spherical cells, whose movement is dominated by isotropic friction on the substrate and cell-cell interactions.Meineke et al. [131] and Drasdo [47] use this model to describe cellular mono-layers in the intestinal crypt and in vitro.For mono-layers, we might assume that all cells have roughly the same velocity, which leaves us with the first term on the left-hand side and the last term on the right-hand side of Equation ( 14): which amounts to a set of coupled ODEs that can be solved using explicit or implicit schemes, such as a forward Euler scheme [131; 127].In general, Equation ( 14) can be written as a set of linear equations and be solved using matrix manipulation, for which standard software packages are available.
The computational challenge of CBMs is two-fold: First, all interacting cell pairs must be identified.Secondly, the relevant interactions take place on chemical and dynamic time scales that are much shorter than the time scale that governs tumour growths.Besides their high computational cost, it should be noted that CBMs might suffer from artefacts when cells are densely packed [116].Despite these shortcomings, CBMs have widely and successfully been applied to study tumour-related issues, such as drug response [e.g.64; 145].

Deformable cell (vertex) models
Deformable cell models (DCMs) are conceptually similar to CBMs: Cell dynamics are distilled into equations of motion.However, rather than representing each cell by simple geometrical objects, DCMs depict cells as a number of connected nodes (cf.Fig. 3).In contrast to CBMs, DCMs hereby capture cell morphology and intracellular mechanics.The nodes within a cell are connected by viscoelastic elements that are either part of the cell membrane or the cytoskeleton.It is essential to include a cytoskeleton since the modelled cells might buckle without an internal scaffolding structure [24].
Assuming that inertial forces can be neglected, the equation of motion for any given node takes the same form as that of a single cell in CBMs (cf.Equation 14).Thus, frictional forces balance out all other forces acting on the node [see also 187].Node-node interactions can be modelled by letting the viscoelastic elements act as damped linear springs.To better depict the behaviour of biopolymers, one might model the elements using nonlinear forces [117; 138].Besides node-node interactions, we must introduce forces that express area and volume constraints on the cell as a whole or on subdivisions spanned by several nodes.By considering subdivisions of the cell membrane, we might furthermore introduce forces that express the resistance of the cortex to bend [138].
Simulating realistic cell dynamics requires a large number of nodes per cell, and the models can become rather complex.Due to the resulting high computational cost, it is often only feasible to simulate single cells or small cell populations on short time scales.
For completeness, it should be mentioned that there is a subgroup of DCMs called vertex models (VMs).In such models, the cells are represented as polygons spanned by a grid of nodes, i.e. adjacent cells share vertices and edges.There is no space between cells.VMs are useful when dealing with densely packed cells and have thus found applications in tissue mechanics but have not been applied to GBM [141; 2].
Table 2 gives an overview of all discrete methods discussed above.

Hybrid multi-scale models
In the literature, the term hybrid or multi-scale model might denote any combination of the different methods discussed above.A large variety of approaches thus fall into this category.Adopting the nomenclature from Deisboeck et al. [45], we distinguish between three types: composite hybrid models, adaptive hybrid models, and calibrated models, such as continuum models with functional parameters.Composite hybrid models track and connect phenomena across various scales while keeping the scale associated with each individual aspect fixed (cf. Section 4.1).Thus, all tumour cells are treated on equal footing.
In contrast, adaptive hybrid models dynamically adjust the local resolution based on the required level of detail (cf.Section 4.2).That is to say that not all tumour cells are represented in the same manner.Finally, the third type of hybrid model denotes, for instance, continuum models, whose parameters have been calibrated based on agent-based models or biophysical considerations of microscopic or mesoscopic phenomena (cf.Section 4.3).

Model type Features References Cellular automata (CA) models
Concept: All CA-type models are discrete and lattice-based.Cell migration, death, movement, division, and growth are modelled as stochastic processes.Particularity: Multiple cells or at most one are allowed per lattice site.Assets: Centimetre-scale simulations are feasible when tracking cell clusters.Individual cells can be represented.Drawbacks: They do not represent cell morphology, have a limited spatial resolution and are computationally heavy.For multiple cells per site, the individual cell positions are not tracked.

Lattice gas cellular automata (LGCA) models
Concept: They are CA-type models.There are multiple velocity channels (and hence cells) per lattice site.A deterministic step evaluating movement is included (see Fig. 1).Assets: Individual cells are represented.Convergence to CM has been shown.Drawbacks: They do not represent cell morphology, have a limited spatial resolution and are computationally heavy.

Cellular Potts models (CPM)
Concept: They are CA-type models.There are multiple sites per cell.A Hamiltonian function governs stochastic events.Assets: CPMs draw on the physics behind the cell dynamics and represent cell morphology.Drawbacks: They have a limited spatial resolution.Centimetre-scale simulations are not feasible.
Authors have included a broad range of phenomena through a variety of terms in the Hamiltonian [142; 166; 65; 43; 179].

Centre-based models (CBM)
Concept: CBMs model cells as physical particles governed by their equations of motion.No lattice is used.Cells are simple geometrical objects.Assets: The physics behind the cell dynamics is accounted for.Cells and their movement are resolved.Drawbacks: CBMs are computationally heavy.Centimetre-scale simulations are not feasible.Cell morphology is not modelled.

Deformable cell models (DCM)
Concept: Cells are depicted as a collection of connected nodes described by their equations of motion.Assets: DCMs model cell morphology and stresses on a subcellular level.There is no lower length scale since no lattice is used.Drawbacks: Only simulations of small cell populations or single cells are feasible.

Composite hybrid models
When using a discrete model to represent tumour cells, the extracellular biochemical players, such as oxygen, are commonly modelled as continuous fluids [150; 151; 149; 103].Moreover, intracellular changes, i.e. processes on a sub-cellular scale, might be tracked for each individual cell using ODEs as discussed in Section 4.1.1.These different microscopic, mesoscopic, and macroscopic mechanisms are all coupled.For instance, the nutrient levels both depend on the cell density and impact the probabilities that are used to predict cell actions.When this interdependence is taken into account, the discrete tumour cell model becomes one of many components in a larger framework that spans and links multiple temporal and spatial scales.The framework as a whole then constitutes a composite hybrid model [cf.10; 7; 11; 194; 108; 106; 104; 27; 180; 9; 28; 64].Such hybrid frameworks have been used to study a wide variety of phenomena, including drug resistance [e.g.52 ; 145].In each case, the focus of the study prescribes the resolution of different features.For instance, in some studies, discrete models of the vasculature are combined with continuum models of the tumour to better understand angiogenesis or drug delivery [e.g.29].As another example, May et al. [128] couple discrete tumour models with biomechanical calculations of stresses and strains to arrive at a more accurate prediction of the tumour shape.Finally, we note that composite hybrid models have been employed in Bayesian inferences to establish posterior probability distributions for relevant model parameters [88] (see also Section 5).

Ordinary differential equations
When dealing with agent-based hybrid models, some authors include sub-cellular dynamics, modelling molecular networks, pathways, and reactions [10; 11; 210; 105].Such processes are governed by a system of coupled nonlinear Ordinary differential equations (ODEs) that represent mass-balance equations and rely on the stoichiometry of the underlying biochemical reactions.
ODEs thus play a vital role in some composite hybrid models.In general, ODEs are a common tool in mathematical biology and oncology.For instance, they can capture global tumour growth patterns [162].While such models do not spatially resolve the tumour, they have successfully been used to study the response to drugs and radiation treatments, recovering experimental constraints on the overall growth pattern [e.g.167 ; 111].The ODE for tumour growth take the generic form where N (t) is the total cell number or volume, and S(N ) denotes an appropriate source function.

Adaptive hybrid models
Detailed modelling of individual cells is required to capture tumour growth at low cell numbers.Agent-based models thus give unique insights into glioma invasion.However, at high cell densities, a lower resolution suffices.
Continuum models realistically capture all relevant aspects of large chunks of bulk tumours.Adaptive hybrid models take advantage of this notion by describing different parts of the tumour tissue using distinct modelling approaches.Thereby, it becomes possible to draw on the strengths of discrete models when depicting centimetresized tumours without the insurmountable computational cost that would otherwise be involved.Bearer et al. [14] and Frieboes et al. [61] couple an off-lattice approach with the continuum models by Wise et al. [199] discussed in Section 2.2.To study glioma invasion, they include dynamic transitions between the continuum and cell-based representations based on the microenvironment.In their model, both mutations and hypoxia might induce the development of a migratory phenotype [see also 159].
The adaptive hybrid models discussed above allow for mass transfer between the different components and are carefully constructed to obey laws, such as mass conservation.However, it is worth noting that the continuum models are pre-assumed without guaranteeing that the imposed functional form emerges from up-scaling the cell-based representation (cf Section 4.3).Moreover, due to the assumptions that enter the continuum and agentbased models, single-cell measurements might lead to different parameter values than those required for continuum models to fit patient-specific data [109].
Rather than coupling agent-based and continuum models, other authors [208; 209; 116] combine agent-based models with different resolutions.Such approaches are known as multi-scale agent-based modelling.For instance, one might couple CA models that consider single cells with models that deal with cell clusters at a coarser lattice resolution [see also 160, for a multi-resolution study of polymers].

Calibrated models
Various authors encode information from cell-based models or biophysical considerations into heterogeneous and time-varying parameters of continuum models.This goal can be achieved in different ways.As briefly discussed in Each method mentioned above relies on simplifying approximations to derive an analytical expression, which might limit their predictive power.To circumvent the need to explicitly state expressions for the macroscopic variables altogether, Kavousanakis et al. [100] use a form of equation-free modelling called coarse projective integration [see also 12].This method invokes the so-called coarse time-stepper that consists of four phases: lifting, simulation, restriction, and extrapolation.These four phases are repeated for every macroscopic time step.In the lifting step, an ensemble of cell distributions is drawn from the macroscopic tumour cell density.
Each cell distribution is then evolved over a short time span using an agent-based model -Kavousanakis et al.We note that the general idea of calibrating coarse models based on more resolved representations is not only limited to continuum models.Recently, Van Liedekerke et al. [186] have thus calibrated the mechanical interaction forces of CBMs based on DCMs.

Data-driven modelling
Researchers draw on a wide range of experimental data to inform their mathematical models of GBM and do so following several different approaches.One can use data to select the underlying mechanism or inform the model parameters.For example, as discussed above, some researchers use DTI MRI data to determine the diffusion properties of cancer cells.Moreover, data allow for qualitative comparisons and constraints on the simulations.
For instance, most models are set to recover generic properties of cancer, such as spherical avascular tumour growth.Finally, one can use data to infer model parameter values using a statistical framework, either maximum likelihood estimation or a Bayesian approach.Systematic incorporation of imaging and molecular data into the mathematical models of cancer via statistical inference is a promising and timely area of future development.This area of development is bolstered by the latest progress in imaging and omics approaches which are producing quantitative data on an unprecedented scale.Furthermore, the latest progress in computational power, simulationbased inference methods, and machine learning make this task computationally feasible.Several modelling studies aim to fit specific data from patients, animal models, or in vitro experiments.Here, we provide a list of illustrative examples and summarise our perspective on important areas of future research.We also refer the reader to some recent reviews on the topic [57; 81].
Two decades ago, Stamatakos and collaborators [171; 170] used agent-based models to investigate the impact of chemotherapy in vivo based on patient-specific PET, SPECT, T1-weighted MRI, histopathologic and genetic data.They reduced the computational cost by clustering the cells into macroscopic (1 mm 2 ) regions rather than simulating individual cells.More recently, Gallaher et al. [64] published a study in which they fit MRI and ex vivo cell tracking data from rats to off-lattice agent-based models, simulating individual cells.They accomplish this fit by using random sampling.However, due to the high computational cost, agent-based models are most commonly found in studies focusing on small cell populations.For instance, Oraiopoulou et al. [139] use agent-based models to recreate the invasive morphologies of GBM observed in vitro.Meanwhile, many papers on agent-based models do not include a fit to real-world data but rather perform in silico experiments to make qualitative and quantitative predictions.For instance, based on such analyses, Perez-Velazquez et al. [145]  The parameters used in state-of-the-art continuum models [109] to capture the migration and proliferation of GBM are derived from MRI and CT images, and the application of continuum models to patient-specific data in vivo is widespread in the literature.Thus, Swanson et al. [177; 178] already used such data to inform their models twenty years ago [see also 77].Later studies that build on their work have included further imaging data, such as PET [190; 75; 135].Isotropic diffusion reaction equations yield a detailed picture based on T1Gd and T2 Flair MRI images [35].Other studies have employed histopathology [62].The anisotropic continuum models discussed in Section 2.1 likewise built on MRI data, drawing on DTI images to compute the diffusion matrix [90; 143; 53; 8; 176].
Due to the relatively low computational cost of continuum models, they can be employed in Bayesian sampling schemes as illustrated by Lipkova et al. [118], who infer model parameters based on MRI and PET images using a Markov Chain Monte Carlo algorithm.Lipkova et al. demonstrate that their Bayesian approach can be used to improve personalised radiotherapy.Ezhov et al. [56] address the same inverse problem using machine learning.Like imaging data, molecular data offers unique insights into GBM [21; 130].In particular, different sequencing techniques, including bulk RNA-seq [23], scRNA-seq [136; 17; 36; 197; 201; 158] (recently reviewed in [95] and [98]) and spatial RNA-seq [155], have helped researchers to uncover different aspects of brain tumours.Sequencing data have thus contributed to our understanding of important cell types, the nature of the invasive tumour front, the spatial cellular organisation of the tumours, and the molecular basis of the heterogeneity of GBM.For instance, Ravi et al. [155] recently identified spatially distinct clusters of genes via spatial transcriptomics, and Neftel et al. [136] identified four main cellular states among malignant tumour cells by analysing data at the single-cell level.Furthermore, Neftel et al. demonstrated that the relative frequency of these four states varies between GBM samples.On the one hand, such information on the heterogeneity of GBM is valuable in its own right and can help to inform the underlying assumptions used in the mathematical models.As exemplified in Section 2.2, mathematical models thus present a unique tool to investigate the implications of such information on the tumour phenotype.On the other hand, the vast number of spatial GBM sequencing data that is becoming available can directly be used in parameter inference, model selection or validation of model predictions.
6 Clinical application of spatio-temporal modelling in Glioblastoma Glioblastoma is a biologically complex and dynamic tumour that exists within the intricate and responsive brain environment.Parcelling and modelling discreet aspects of these biological processes and systems holds the potential for numerous direct clinical care applications and has the potential to address several areas of unmet need.
Mathematical diffusion-tensor modelling is currently widely used in the context of MRI diffusion-tensor-based tractography for the modelling of white-matter tracts.Understanding patient-specific anatomy, and anatomicopathological relationships, is critical for effective surgical planning and avoidance of iatrogenic injury [196].Future refinements in data modelling can be expected to improve spatial resolution and reliability, which would directly translate to improved validity of neuro-imaging with direct clinical applications (cf.Section 5).
A key concern in clinical practice is the prediction of tumour growth patterns, sites of progression, and the location of tumour recurrence.This triad is important both for therapeutic planning and prognostication [91].
Mathematical models have the potential to refine these facets of care and hence to offer highly-tailored conformal radiation dosing and targeted surgical supra-marginal resection [137; 189].These therapeutic adaptations might maximise the prognostic benefits and lead to the destruction of the tumour while sparing normal neural tissue.
Similarly, when combined with the modelling of growth rates, this information could be assimilated to generate more accurate prognostic maps with the prediction of anatomico-pathological functional deficits and overall prognosis.
Recent years have been characterised by increasing recognition of glioblastoma's molecular intra-and intertumoural heterogeneity and its functional consequences.This is exemplified in the evolution of the WHO classification of tumours from a purely histopathological to a fully-integrated molecular classification [119].As we have summarised, mathematical models of glioblastoma hold the potential to illustrate novel prognostic aspects of tumour biology and the micro-environment not currently captured.These might include structural metrics, such as tissue density and cellular linage, or functional ones, including oxygenation, nutrient supply and metabolism.
The potential for applying spatio-temporal mathematical modelling of glioblastoma in clinical care is substantial but thus far remains largely unrealised.As summarised in this review, the foundation of knowledge in this area is now considerable.The stage is set for those who wish to advance the field from theoretical models to practical applications.

Summary and discussion
This paper presents an overview of mathematical models that provide means to scrutinise and predict the spatial and temporal evolution of tumours.While we focus on models for glioblastoma and highlight characteristics of this aggressive cancer type, we also attempt to give a rounded picture of state-of-the-art approaches across the field of mathematical oncology.
Spatially resolved tumour models fall into two main categories: They either treat the tumour cell density as a continuous variable (cf.Section 2) or deal with autonomous agents representing individual cells or cell clusters (cf.Section 3).
The computational cost associated with continuous models is relatively low.This property allows for the modelling of the entire tumour volume and the application to a wide variety of (patient-specific) data.Furthermore, continuous models can capture many of the characteristics of GBM.We summarise how continuous models achieve this goal in Table 3. GBM spreads anisotropically following existing structures in the brain, including white matter tracks and blood vessels.On a macroscopic level, this behaviour is well-described by reaction-diffusion equations when including an anisotropic diffusion term (cf.Section 2.1).
Moreover, by including source terms that couple a set of such reaction-diffusion equations, continuous models can link different cell sub-populations and account for interactions with the environment (cf.Section 2.3).
Continuous models are hereby able to encapsulate the hierarchical nature that is a defining feature of GBM (cf. Section 2.2).The models also show great versatility: By including additional biologically motivated terms into By allowing for feedback between the cells and environment, one can introduce preferential migration, leading to anisotropic growth along structures.

Cell hierarchy and lineage
A set of coupled PDEs described subpopulations with different properties.
Individual cells underlie rules that are affected by cell type.Mitosis and apoptosis E.g. logarithmic or exponential source terms are included in the PDE.These might link different sub-populations and be associated with properties of the environment.
Both phenomena underlie stochastic or deterministic rules that are invoked for each agent.
the reaction-diffusion equations, continuous models can be painlessly extended to include phenomena ranging from immune responses to the impact of radiotherapy.
Tumour cells must always be understood in relation to their surroundings.In this connection, aspects, such as nutrition flows, the extracellular matrix, or the vasculature, are often (if not always) represented using continuous models.Even when using agent-based models, continuous models thus play a vital role.However, continuous models have their drawbacks.They do not account for the stochastic events that dominate at low tumour cell densities.Consequently, they are not ideal for addressing some questions regarding the invasive tumour front or tumour recurrence that is a prevalent obstacle in treating GBM.
At low cell densities, agent-based models have an edge.Such models come in many different flavours.Some agent-based models borrow from the description of other phenomena.They thus lean on theoretical models, e.g. from solid-state physics, posing biological events in terms of physical forces and potentials (cf.Sections 3.1.3,3.2.1, and 3.2.2).Other models rely on purely stochastic rules to encode biological information (cf.Sections 3.1.1and 3.1.2).Moreover, agent-based models differ in the details that they encompass.While some models address aspects of cell morphology, others do not.While some models define a lower threshold for their spatial resolution by employing a lattice on which the agents move, others do not.
Like continuous models, agent-based models are able to capture the essential properties of GBM (see Table 3).
Anisotropic diffusion along existing structures can be included by altering the rules that govern cell migration based on the environment.Moreover, by altering these rules based on the cell type, the hierarchical organisation of GBM can be incorporated into the model.
The high computational expense of agent-based models constitutes their main drawback.It can be a significant stumbling block on the way to unleashing their full potential.However, this shortcoming can be partly circumvented.One might thus combine different approaches into a hybrid model bridging across the scales involved (cf.Sections 4.1 and 4.2).Alternatively, one might study small cell populations to derive emergent macroscopic properties that can be used to calibrate models on a coarser scale (cf.Section 4.3).Indeed, one of the main assets of agent-based models is their ability to link an understanding of individual cells to macroscopic events.
To further lower the computational cost when confronting the models with data, one might use sophisticated statistical techniques or machine learning methods that demonstrably reduce the number of models that need to be constructed [see also 121; 94], as discussed in Section 5.
Each of the available mathematical models has its strength and weaknesses.Which model to employ will hence be dictated by the research question and the available data.First, some approaches are better suited for depicting certain phenomena than other methods are.One might drastically simplify the work associated with the practical implementation by choosing the method wisely.Secondly, not all models are equally informative.
The question is whether the additional information can be exploited and for what purpose [see also 51].On the one hand, the chosen model sets limits on the insights that one might gain.This notion favours more complex models.On the other hand, one should avoid cracking nuts with a sledgehammer.To exemplify this, it is thus worth noting that non-spatial models make full use of certain types of (patient-specific) data whilst being unable to address any questions on tumour morphology.They are undoubtedly very valuable tools in their own right.
The same holds for each spatio-temporal cancer model presented above.With this in mind, we hope that our review provides a useful guide and helps the reader select a suitable mathematical model for their research.
Mathematical models are only useful if they describe and relate to experimental data.Our review thus ends with a summary and discussion of the application of mathematical models to experimental data and its clinical relevance.A wealth of imaging and molecular data is becoming available, particularly due to advances in single-cell and spatial transcriptomics.At the same time, researchers have begun to implement spatial models into statistical frameworks to perform systematic inference based on these data.While few such inference analyses have yet been published, advances in machine learning and computing power promise to further enable the systematic incorporation of imaging and molecular data into the spatial models of GBM.We thus envision a greater impact of mathematical modelling in unravelling the basic biology of GBM and its applications in clinical treatment in the near future.
30; 140; 126; 120; 75].As regards the ECs, we again encounter logistic or exponential source terms.Meanwhile, the gradient of the TAF concentration, as well as that of adhesive molecules (e.g.fibronectin), enter through a flux term encapsulating chemotaxis and haptotaxis: 4).Due to the aggressive invasion of healthy tissue by GBM, this shortcoming is especially problematic for understanding brain tumours.Moreover, continuum models cannot adequately deal with heterogeneity.To address these issues, one needs to model cancer on a mesoscopic scale following individual tumour cells or cell clusters [see also76].Mathematical models that do so are referred to by various names: discrete models, agent-based models, individual-based models, or cell-based models.In all cases, cells or cell clusters are modelled as autonomous agents that follow a set of (stochastic) rules prescribing cell movement, death, division, and growth.By doing so, discrete models capture the variability among cells that might, for instance, arise as a response to the microenvironment [e.g.7].Note that while discrete models treat cells as autonomous quantities, continuum models are commonly used to depict the microenvironment (cf.Section 4.1).PDEs hence still come into play when dealing with, e.g.nutrient and drug concentrations [e.g.120; 123].

Figure 1 :
Figure1: Schematic overview of the concepts behind lattice gas cellular automata (LGCA) models.The left panel summarises the cell dynamics that are handled as stochastic processes: cell migration, death (and lysis), growth, and cell division.Filled (orange) circles denote occupied states, while there are no cells in the empty circles.On its own, the left panel corresponds to the cellular automata models with at most one cell per lattice site that are discussed in Section 3.1.1.However, as illustrated in the right panel, several velocity channels are associated with each location, and multiple cells can occupy the same position.Filled orange circles denote occupied channels, while empty channels are indicated with empty circles.Here, there are five velocity channels per location: four channels that lead to migration and one channel at which the cell is at rest.After determining the stochastic processes, the model computes the movement of the cells based on the velocity channels that they occupy.
Potts models (CPM) define a Hamiltonian function, H, to incorporate cell movement, growth, interactions between adjacent cells, and interactions between cells and their microenvironment [see 181, for a detailed discussion of Hamiltonian mechanics ].Moreover, CPMs attribute multiple lattice sites to each cell, which allows the models to account for cell morphology.Again, CPMs do so via the Hamiltonian function.The basic concepts are illustrated in Fig. 2

Figure 2 :
Figure 2: Schematic overview of cellular Potts models.Each cell (and the ECM) has a unique identifier.Here, we consider cells 1, 2, and 3 and the ECM denoted by 0. Each cell occupies several lattice sites.Randomly selected sites at the cell borders (*) might swap affiliation depending on the energy change associated with this swap allowing the cells to grow and move.

Figure 3 :
Figure 3: Schematic overview of off-lattice models.The left panel illustrates spherical cells in a centrebased model (CBM), highlighting forces on cell 1.Due to the overlap (δ 21 ) between cell 1 and 2, cell 1 is repulsed by cell 2.Overlapping cells forming a dumbbell also occur during cell division.F 31 exemplifies an adhesive force occurring without direct contact.In an isotropic environment, the drag force will be proportional to the velocity, v 1 .The right panel illustrates a cell made of multiple nodes in a deformable cell model (DCM).Node i at the cell membrane is subject to forces from the surrounding nodes (F ji ).Constraints on areas and volumes might lead to additional forces (F i ), while the movement of the node (v i ) gives rise to viscous forces.The cell has a cytoskeleton (CSK, dashed lines).
Contact mechanics offers different approaches to model these cell-cell interactions, including Hertzian formulations and the Johnson-Kendall-Roberts (JKR) model of elastic contact [cf.48; 80].Some authors simply assume that all forces involved in cell-cell interactions are well-approximated by linear springs, generalised linear springs, or polynomial expressions [e.g.131; 25; 127, and references therein].
The hybrid model presented by the group ofStolarska and Yangjin [107;172;105] treats necrotic tumour zones, quiescent tumour tissue, and the surrounding microenvironment as continuous fluids.The relevant PDEs are solved on a regular grid.The cells in the proliferative tumour rim (100-200 µm), on the other hand, are modelled using a CBM.They are represented as autonomous ellipsoids.At the boundary between the cell-based and continuous components, interpolation is used to establish the boundary conditions for the reaction-diffusion equations and the forces exerted on the individual agents [see also 55, for a related example from molecular dynamics].

Section 3 . 1 . 2 ,
one might up-scale agent-based models using the mean-field approximation [cf.80].Other authors draw on the emergent macroscopic properties of agent-based models to derive phenomenological relationships [cf.120; 122].
use a CA model.In the restriction step, the macroscopic density is updated by computing the ensemble average of the agent-based simulations.Under the assumption that the macroscopic cell density evolves more slowly than the microscopic variables, the observed change in the cell density yields an estimate of its time derivative.Using the forward Euler method, this information can be employed to project the cell density further into the future without considering individual cells.
Ezhov et al. construct a library of 100,000 simulations with different parameter values and use their atlas to train a neural network.The trained neural network can infer the patient-specific model parameters based on MRI images within minutes, allowing for effortless model personalisation.Today, continuum models are thus at the point where they have the potential to enter clinical settings and contribute to personalised treatments.While the same does not yet hold for agent-based models due to their high computational cost, machine learning might likewise help to overcome this obstacle in the future, as discussed by Joergensen et al.[94].Machine learning algorithms would play the same role as they do in continuum models by circumventing the need to produce further simulations once a surrogate model or inference algorithm has been trained.Similar strategies are applied when dealing with agent-based models outside of cancer research[156].
cf. Sec.2.2.Different properties of the environment, such as nutrient flow or vasculature, are modelled on a macroscopic scale.Assets: They are computed at a low computational cost.Spatial gradients and time variations, e.g. in nutrient levels, are sufficiently well resolved to study tumour cells.Drawbacks: CMs of ECs do not capture the behaviour of individual cells.
addressed the development of drug resistance, Kim and collaborators [106; 104] investigated the key mechanisms behind molecular switches, and Schmitz et al. studied growth patterns.

Table 3 :
Overview of how different model types can deal with properties that define glioblastoma growth.