A stochastic hierarchical model for low grade glioma evolution

A stochastic hierarchical model for the evolution of low grade gliomas is proposed. Starting with the description of cell motion using a piecewise diffusion Markov process (PDifMP) at the cellular level, we derive an equation for the density of the transition probability of this Markov process based on the generalised Fokker–Planck equation. Then, a macroscopic model is derived via parabolic limit and Hilbert expansions in the moment equations. After setting up the model, we perform several numerical tests to study the role of the local characteristics and the extended generator of the PDifMP in the process of tumour progression. The main aim focuses on understanding how the variations of the jump rate function of this process at the microscopic scale and the diffusion coefficient at the macroscopic scale are related to the diffusive behaviour of the glioma cells and to the onset of malignancy, i.e., the transition from low-grade to high-grade gliomas.


Introduction
Gliomas are the most common type of primary brain tumours, accounting for 78% of all malignant brain neoplasia [43].They originate from mutations of the glial cells in the central nervous system and are classified by the World Health Organisation (WHO) into four grades according to the degree of malignancy (see [85] for a more detailed description).In this work, we mainly focus on the low grade gliomas (LGGs), which are a class of rarely curable diseases, often resulting in the premature death of the patient.Since in the last years some medical interventions have shown to improve the median survival time of the patients, the study of this class of tumour has become of great importance for the clinicians.
The development, growth, and invasion of gliomas in the brain is a very complex phenomenon, involving many interrelated processes over a wide range of spatial and temporal scales.As such, often the individual cell behaviours and the intracellular dynamics described at a microscopic scale are manifested by functional changes in the cellular and tissue level phenomena.Therefore, this multiscale nature of glioma evolution requires modelling techniques that are able to deal with different levels of description.
The first mathematical models for the study of brain tumours started to emerge in the early 1980s (see [26,28,27,80,81] for further details).Since then, the mathematical modelling of glioma evolution has evolved considerably and several different approaches have been proposed, going from discrete or hybrid microscopic models to macroscopic and multiscale frameworks.Discrete models at the microscopic scale, also called agent-based models, have been used to describe the dynamics of individual cells moving on a lattice (for some examples we refer the reader to [84,58,45], or, specifically, to [4] for cellular automata models and [39] for cellular Potts models).Further, stochastic discrete models for cell motion have also been proposed, e.g.describing 2D persistent random walk or 3D anomalous diffusion [29,57,5,72].In particular, recently in [72], the authors have presented the analysis of 3D cell tracking data, based on a persistent random walk model adapted into the context of glioma cell migration.At the macroscopic scale, several phenomenological models for glioma evolution stated in the form of reaction-diffusion-advection equations have been proposed and studied [78,44,80,77], also including patient-specific data (e.g. in the form of diffusion tensor imaging (DTI) information).This has allowed for a comparison between the real and the virtual tumour evolution [51,53,17,60].Concerning multiscale models, a broad and rich literature has been developed for the integration of microscopic and macroscopic dynamics (for some examples see [47,49,7,66,32,33,55,52,34]).In particular, in [32], a more detailed description of the migration process of individual cells, involving the dynamics of cell receptors and the interaction with the tumour microenvironment, is discussed.
A key aspect of modelling tumour evolution concerns cell movement, which is based on a combination of complex processes involving motility and migration: motility refers to the random movement from one location to another, while migration involves also the interactions between cells and the microenvironment [59].
The first description of particle movement, which uses a stochastic Markov process combining deterministic ordinary differential equations (ODEs) for the continuous movement with Poisson-like jumps for the random change of direction, was introduced in 1974 by Stroock [75] on the basis of the biological observations illustrated in [1].The concept of piecewise deterministic Markov processes (PDMPs) was introduced in 1984 in [24].An extension of [24] was then provided in [13,15], where the authors developed the extended generator and the differential formula for piecewise diffusion Markov processes (PDifMPs), showing that all the classes of proposed stochastic hybrid processes can be seen as a special case of their concept of a general stochastic hybrid system (GSHS).Further, in [10] a general class of continuous-time stochastic hybrid systems in which the continuous flow is the solution flow of a stochastic differential equation (SDE) was presented.These processes have been widely applied in different contexts, e.g. for interacting particle systems [11], air traffic management [14], or gene network [61]) and especially in biological modelling (for some examples, see [82,36,18,67,41,68,22,70]).However, it seems that the use of PDifMPs in the context of tumour growth, motility, and migration has not yet been investigated.In this article we extend the description of cell movement based on velocity jump processes with the use of PDifMPs in the context of glioma progression.In particular, we build a multiscale model, starting with a contact-mediated description of cell motion on the microscopic scale using PDifMPs.We use the extended generator for such processes to derive a generalised Fokker-Planck equation, including the description of the tumour-microenvironment interactions.The solution of this equation provides the joint density of the transition probabilities of this Markov process for all the involved variables.As the variables involved in these interactions are fast-acting compared to the macroscopic scale, we make use of a scale separation variable and the Hilbert expansion method to derive the corresponding macroscopic scale equation for the time and space variables (for a more general discussion of multiscale modelling and moment closure techniques, we refer the reader to [8,54,50]).
The paper is organised as follows.Section 2 contains a brief introduction to PDifMPs.In Section 3, we derive a stochastic multiscale model for glioma progression.Numerical simulations in a 2D scenario for the resulting macroscopic equation for the tumour cell density are presented in Section 4, including several studies on the effect of parameter variations.Finally, in Section 5, we review our results and discuss further directions of research.

Definition and notation
In this section, we provide a brief introduction to PDifMPs and the construction of their paths.We refer the reader to [15] and [61] for a general description of stochastic hybrid systems.
For the couple of non-exploding processes (St, Vt), we assume that the first stochastic component (St) t∈[0,T ] possesses continuous paths in E1 and the second component (Vt) t∈[0,T ] is a jump process with right continuous paths and piecewise constant values in V.The times (Ti) i∈N at which the second component jumps form a sequence of randomly distributed grid points in [0, T ].The motion of the PDifMP (Ut) t∈[0,T ] on (E, B(E)) is defined by its characteristic triple (φ, λ, Q) as follows: At the end point Ti+1 of each interval, si+1 is set to the current value of φ( • , ui) to ensure the continuity of the path.Further, a new value vi+1 is chosen as fixed parameter for the next interval according to the jump mechanism described below.We define also the function b with values in R d 1 , which represents a family of drift coefficients, and the d1 × m matrix σ with real coefficients.
Assumption 2.1.We assume that b : E → R d 1 and σ : E → R d 1 ×m are linearly bounded and globally Lipschitz continuous for all s ∈ E1.
For any vi ∈ V, this assumption ensures the existence and uniqueness of the solution to (1) (see Theorem 5.2.1 in [62]).Moreover, the stochastic flow satisfies the semi-group property, i.e., • λ : E → R+ is the jump rate, i.e, it determines the frequency at which the second component of (Ut) t∈[0,T ] jumps.
is the transition kernel that determines the new values of the second component after a jump occurs.For all u ∈ E, it satisfies Q(u, {u}) = 0, meaning that the process cannot have a no-move jump.
Moreover, for all t ∈ [Ti, T ], i ≥ 0, we define the survival function of the inter-jump times as This function states that there is no jump in the time interval [Ti, t) conditional on the process being in the initial state ui.Let U be a uniformly distributed random variable on [0, 1], thus Assumption 2.2.Let λ : E → R+ be a measurable function such that ∀ui ∈ E and T > 0 Moreover, there exists a measurable function ψ : [0, 1] × E → E such that for ui ∈ E and A ∈ B(E) ψ represents the generalised inverse function of Q.For a fixed t, ψ(U(ω), U (ω)) is a random variable describing the post-jump locations of the second component of U .
Assumption 2.3.For all A ∈ B(E), Q( • , A) is measurable, while for all u ∈ Ē the function Q(u, • ) is a probability measure.
Summarising, the first component of the triple (φ, λ, Q) describes the continuous evolution of the trajectories of the process (Ut) t∈[0,T ] between jumps in time intervals defined by the survival function S, while the couple (λ, Q) yield the jump mechanism.All three components of (φ, λ, Q) are coupled.

Construction
From the local characteristics (φ, λ, Q), it is possible to iteratively construct the sample path Ut as follows.Let (Un) n≥1 be a sequence of iid random variables with uniform distribution on [0, 1] and u0 = (s0, v0) ∈ E the initial value of (1) at T0 = 0, such that u0 can be either an F0-measurable random variable (independent from the Wiener process) or a deterministic constant, for some ω ∈ Ω.We apply the survival function S(t, u0) defined in (2) and use its generalised inverse ζ with the first element U1 to determine T1 = ζ(U1, u0), i.e., the first jump time of the second component of Ut.We then define the sample path Ut up to the first jump time as The trajectory of Ut follows the stochastic flow φ given in (1) starting from U0 = u0 until a first jump occurs at the random time t = T1.The post-jump state UT 1 is determined through the measurable function ψ.For all A ∈ B(E), the distribution of ψ(U2, u0) is given by where τ1 is the waiting time until the first jump occurs, i.e. τ1 = T1.
Restarting the process from the post-jump location UT 1 , we define the next waiting time before a jump occurs from the survival function (2).In this way, we find the next jump time T2 = T1 + τ2.
Consequently, the state of the process in the interval [T1, T2) is given by We proceed recursively to obtain a sequence of jump times (Ti) i≥1 , such that the generic sample path of Ut, for t ∈ [Ti, Ti+1), is defined accordingly by The number of jump times that occur between 0 and t is denoted by Assumption 2.4.For all t > 0 and for every starting point ui ∈ E, E[Nt|u = ui] < ∞.
This assumption ensures the non-explosion of the process Ut.Under the Assumptions 2.1-2.4 the piecewise diffusion process can be constructed as a strong càdlàg Markov process (see [15] for further details), called then a Piecewise Diffusion Markov Process (PDifMP).

Extended generator of the PDifMP
The notion of infinitesimal generator is an extremely important tool for the study of Markov processes [9,24].In the following, we adopt the definition in [61,15], and, for the reader's convenience, we recall the theorem that fully characterised the extended generator (see [9] and references therein for further details about the difference between extended and classic generators).
Theorem 2.1.Let Ut be a PDifMP with characteristics (φ, λ, Q).The domain D(A) of the extended generator A consists of all bounded, measurable functions f on E ∪ ∂E satisfying: Then, for f ∈ D(A), u = (s, v) ∈ E, the extended generator Af is given by where for s = (s1, . . ., s d 1 ).Here, ∇sf (s, v) • b(s, v) is the inner product in R d 1 , σ T is the transpose matrix of σ, ∇ T s is the transpose operator of ∇s.We refer to [15] for the definition of L loc 1 (p) and the proof of this theorem.

Generalised Fokker-Planck equation
The adjoint of the generator is used to derive the generalised Fokker-Planck equation, describing the time evolution of the probability distribution g(t, s, v) of the process.The equation is given by where the adjoint operator of A dif reads We refer to [6,40] for further details on the derivation of Fokker-Planck equations for general Markov processes.

Application to tumour modelling
Gliomas can be considered as dynamical ecosystems where cells undergo constant changes due to many cellular processes, e.g.migration, proliferation, death, or creation of new blood vessels [79,3].We focus on the process of cell movement, which is responsible for the global diffusive features that characterise glioma evolution.Cell movement can be divided into motility and migration.Motility refers to the random or spontaneous motion of cells from one location to another, while cell migration involves many interconnected biological aspects, such as environmental cues driving it.Thus, methods that take into consideration the stochastic nature of this phenomenon (i.e., motility) while accounting for environmental cues influencing it (i.e., migration) are important for providing a more complete understanding of the entire process.Following [52,50], we model the process of cell movement under the influence of subcellular scale interactions, considering the effects of the amount of bound receptors located on the cell membrane.Specifically, we consider the role of integrins in this dynamics [23,30].Referring to cell migration, we take into account the alignment of the tissue as a cue enhancing the efficiency of cell invasion [83,25], as cells tend to attach to the fiber and crawl along them, a phenomenon referred to as contact guidance.However, since the direction that cells decide to follow remains random, there is a need to consider a stochastic description for the motility component.
Inspired by particle movement models [75,63,64], we propose piecewise diffusion Markov processes for the modelling of cell movement.In the context of persistent random motion, the continuous stochastic component of the PDifMP describes the contact guidance phenomenon, while its second component describes the random motility dependent on the velocity jump process.This approach makes it possible to describe the cellular migratory response to environmental signals while keeping the random aspect of cell motility.Moreover, it also allows us to show how several well-established methods proposed in the literature (e.g.see [63,64,50]) can be cast into a rigorous PDMP framework.

Interactions between cells and microenvironment
In order to migrate through the complex brain structure, glioma cells must adapt quickly to the physical characteristics of the environment.Their interactions with the extracellular matrix (ECM) [37] are mediated by the binding between the integrins and the ECM fibrillar proteins.These bindings allow them to exert the forces necessary for them to migrate [69,30].As these processes happen at a sub-cellular level, we describe the mechanism behind cell motion modelling the dynamics of the receptors on the tumour cell membrane.
Let y(t) ∈ (0, 1) be the concentration of bound integrins and let us assume that the binding between integrins and tissue occurs in areas of highly aligned fibers [32].The binding process can be described with the following general reaction where R0 defines the total number of cell surface receptors, Q(x) the macroscopic volume fraction of tissue (including ECM and brain fibers), depending on the position x ∈ X ⊂ R 3 , and k + and k − the rates of attachment and detachment between cell and tissue [34,32].Within this framework, denoting by x = x0 + vt, we look at the path of a single cell moving from an initial position x0 with velocity v ∈ V ⊂ R 3 .V = αS 2 is the closed set for cell velocities, where S 2 denotes the unit sphere on R 3 and α the mean speed of a tumour cell, which is assumed to be constant.Since we are interested in the interactions between cell surface receptors and the ECM, and this binding process takes place for fixed position x, we ignore any type of randomness resulting from the velocity change.The mass action kinetics for the concentration y(t) is governed by the following ODE: Since the integrin dynamics are much faster than the macroscopic time scale phenomena, we assume that they equilibrate rapidly [33,50,19].Thus, after rescaling y/R0 → y, we consider the unique steady state y * of (10), given by and we define a new internal variable z := (y − y * ) ∈ Z = (y * − 1, y * ) ⊂ R, which measures the deviation of y from its steady state [32,50].Considering the piecewise location of a single cell x = x0 + vt through the density field Q(x), z satisfies where The internal variable z is bounded as long as ∇xQ(x) is bounded and its sign depends on the current orientation of the cell w.r.t the gradient of Q(x).

PDifMP description for glioma cell movement
To model cell movement under the influence of external signals, we assume that the sample path of an individual cell starting in position x0 and moving in a certain direction due to contact guidance for a random period of time is given by Here, the second term in the r.h.s represents the stochastic variability in the velocity, with σ ∈ R being the diffusion coefficient and Wt the standard Wiener process.Due, for instance, to collisions with other cells in their surrounding [56,65], during the movement a cell stops for a negligible duration and reorients its path [56].This causes the cell to adopt a new velocity to continue migrating in the new direction until another obstacle is encountered.To describe this process, we rely on the introduced PDifMP framework.We set E1 = X × Z ⊂ R 4 and we denote by St := (Xt, Zt) the continuous component describing cell motion.Their evolution is characterised through the SDE (12) for cell motility and the ODE (11) for the interactions with the microenvironment.Both processes are affected by spontaneous velocity changes induced by the jump process Vt.Then, we denote by E = E1×V the state space of the piecewise process Ut = (Xt, Zt, Vt) for cell motility and migration and by φ : [0, T ] × E → E1, the solution to the coupled system (11)- (12).
As the duration of reorientation is negligible, we describe the direction of a cell at a given instant.Moreover, under the additional assumption that the motion is Markovian in the state space, we state that cell direction is described with an inhomogeneous Poisson-like process [38], whose intensity depends on time, position on the scaled sphere V, and internal state.Thus, the cell reorientation rate referring to the jump rate function λ : [0, T ] × E → R+ of the stochastic process ut depends on the integrin state z.This means that the binding process is seen as the onset of reorientation.In particular, following [76], we assume that, if many integrins are bounded, cells tend to change direction frequently in order to escape the densely packed areas, resulting in an increased rate λ.Thus, following [33,34], we set λ(ut) := (λ0 − λ1zt) ≥ 0 with λ0 and λ1 positive constants.In particular, λ0 refers to the basal turning frequency of an individual cell [74] accounting for the "spontaneous" cell motility, while the term λ1z represents the variation of the turning rate in response to environmental signals.
Following the construction described in Section 2.2 with initial state u0 = (x0, z0, v0), we use the jump rate function λ defined in (2) to determine the duration of movement before any reorientation of direction occurs.Moreover, considering that the velocity jump process vt is of Markovian type, we have that cells retain no memory of their velocities before the reorientation.Thus, we define the Markov transition kernel Q, determining the post-velocity jump state of the process ut, using K(x, v, v ), which describes the distribution of newly chosen velocities, having that K(x, v, v ) = K(x, v).Definition 3.1.Let ν be the standard Lebesgue measure on (V, V) and K : X × V → [0, ∞] be a measurable function with respect to the σ-algebra X ⊗ V such that Then, the mapping defines a Markov transition kernel over V, where ν(dv) = dv.
Denoting by q(x, v) the fiber distribution function over V, with v = v v ∈ S 2 , and by a scaling constant [49,48], we assume that the dominant directional cue leading cell migration is given by the fiber network.Thus, the transition probability kernel is given by For the fiber distribution function q(x, v), different expressions can be found in the literature, such as the Von Mises-Fisher Distribution, the Peanut Distribution Function, or the Orientation Distribution Function (ODF) [66,2].A comparison among these distributions have been proposed in [20], in both 1D and 2D scenarios.We rely on this analysis and we choose the ODF for describing q(x, v), i.e., we set Here, v stands for the fiber direction, x for spatial position within the brain, while D is the diffusion tensor taking into account information about the water diffusivity in the brain [20].We also assume that fibers are not polarised, i.e., q(x, v) = q(x, −v) for all v ∈ S 2 .It is straightforward to verify that q is a probability distribution on S 2 [33,32,34].From (2), it is possible to construct the sequence of jump times (Tn) n≥1 , with Tn = τ1 + • • • + τn for all n ≥ 1 (and T0 = 0 by convention), such that the process Ut describing cellular movement is piecewise constructed on each interval [Ti, Ti+1), i = 1, . . ., n, via the characteristics (φ, λ, Q) given by Here, zt is the solution of (11) and vt is a piecewise constant over each interval of random length Ti+1 − Ti.As proven in [15], this construction leads to a càdlàg strong Markov process, describing cell motion in an anisotropic environment.In summary, the overall system describing a contact-mediated movement of glioma cells at the microscopic scale reads The solution of ( 17) is a triple and hereafter we will refer to (Xt, Zt, Vt) as (xt, zt, vt) as we are talking about the sample path of Ut.

Derivation of the mesoscopic equation and its macroscopic limit
We rely on the definition of the extended generator of (Ut) t∈[0,T ] given in Section 2.3 to obtain a mesoscopic equation describing the evolution of the joint probability density function of all microscopic variables.In the specific, for all test functions f ∈ D(A), the extended generator A of the above defined process Ut reads where λ and Q are given in (16).Notice that the integral term in ( 18) is defined over V as the transition kernel Q has a density defined on V.
Let g(t, x, z, v) be the joint pdf of the microscopic variables at time t ∈ [0, T ], position x ∈ X, internal state z ∈ Z, and velocity v ∈ V.In this context, we refer to as glioma density function.The adjoint operator A * g is given by Thus, following the analysis of Section 2.4, the generalised Fokker-Planck equation for the evolution of g(t, x, z, v) reads where, from ( 14), the turning operator reads [63,46,66,33] Remark 3.1.Note that for σ = 0, (20) coincides with the kinetic transport equation derived in [33,50,32].This means that the PDMP resulting from setting σ = 0 in ( 20) is the formally defined mathematical model underlying the description in [33,50,32].
We introduce the notations for the mean fiber orientation and the variance-covariance matrix of the fiber orientation distribution, respectively.Notice that the symmetry on the fiber distribution implies Eq = 0.
Following [33,50], we model proliferation as an effect of cell-tissue interactions via integrin binding Here, M (t, x) denotes the macroscopic cell density, that is, the marginal distribution of g(t, x, z, v) over all possible velocities and internal states, i.e., Moreover, µ(M ) is the growth function and the kernel X (x, z, z ) is a probability density in the second variable z characterising the transition from z to z during the proliferation process at position x.For X we only assume that the operator P(g) is uniformly bounded in the L 2 -norm, a reasonable biological condition related to the space-imposed limits on cell division.Thus, for the evolution of g(t, x, z, v) we obtain the following equation Due to the high dimensionality of ( 23), numerical simulations of this equation would be too expensive.Moreover, clinicians are more interested in the macroscopic dynamics of the tumour rather than in the lower scale interactions.Thus, we derive the macroscopic equation for the evolution of the tumour density, based on the definition of the moments of g with respect to v and z: Notice that we do not consider higher order moments of g with respect to z as the subcellular dynamics are much faster than the events taking place on the other scales, so that the deviation z is close to zero.Dropping the (t, x) notation for simplicity, the moment equations reads and Following [32,33], we consider a parabolic scaling of the moment equations setting x → x and t → 2 t for space and time variables, respectively.In particular, we scale the growth rate function µ(M ) with 2 as it accounts for faster dynamics.Thus, we obtain and We consider the Hilbert expansion methods [31,50] expanding the moments of g as By equating the same powers of in ( 26) and ( 27), we derive the equation for the leading order coefficient M0 of the Hilbert expansion of M .Thus we obtain With classical scaling arguments (see [32] for more details), we obtain M z 0 = m z 0 = 0 and m0 = q(x,v) w M0.On account of that, using the symmetry assumption, i.e., Eq = 0, from (31), we obtain M z 1 = 0, and Moreover, considering (30) and following the analysis in [63,32], we get M1 = 0, and Replacing it into (32) and integrating over V, we get: where Therefore, the evolution equation for M0 reads where denotes the function that carries the information about the influence of the subcellular dynamics, while refers to the macroscopic tumour diffusion tensor.In addition, the tumour drift velocity is given by In view of the results obtained in [32], the -correction terms for M can be left out and, after ignoring the higher order terms and discarding subscripts, we obtain the following evolution equation characterising the macroscopic glioma density: Using the theory of monotone operators for nonlinear parabolic equations and following the approach in [73,71], it is possible to prove the existence, uniqueness and non-negativity of the solution of the following parabolic problem with homogeneous Neumann boundary conditions. where We refer the reader to Appendix A.1 for more details about the necessary assumptions on the operators and for an outline of the proof of the well-posedness of the macroscopic problem (38).

Numerical simulations
We perform 2D simulations of the macroscopic equation for the tumour cells (37) to study the impact of both the subcellular dynamics and the stochastic parameter σ on the overall tumour evolution.With this aim, we firstly specify parameters and coefficient functions involved in the equation.Concerning the tumour diffusion tensor D T (x) in (35), we numerically compute it using the orientation distribution function given in (15), where D(x) represents the water diffusion tensor obtained from processing (patient-specific) DTI data.Taking advantage of this DTI information, for the macroscopic tissue density Q(x) we assume the following expression where F A refers to the fractional anisotropy of the tissue.We refer to [32] for its definition.This choice is motivated by the fact that the fractional anisotropy represents a measure of the fiber alignment and, since in this setting fiber alignment is guiding cell migration, it is reasonable to assume that the function Q(x) expresses higher values where the tissue is more anisotropic.Following several previous works (e.g.see [33,20]), for the growth rate µ(M ) we employ a logistic growth term defined as with µ 0 the constant growth coefficient and K M the tumour carrying capacity.Finally, we report in Table 1 the range for the constant parameter values involved in the macroscopic setting (37).The values for the stochastic parameter σ are proposed based on the ranges of the other parameters.We present 2D numerical simulations performed with a self-developed code in Matlab (MathWorks Inc., Natick, MA).The computational domain is a horizontal brain slice reconstructed from MRI scans.The DTI dataset used to compute the D T (x) was acquired at the Hospital Galdakao-Usansolo (Galdakao, Spain), and approved by its Ethics Committee: all the methods employed were in accordance to approved guidelines.A Galerkin finite element scheme for the spatial discretisation is considered, together with an implicit Euler scheme for the time discretisation.For the initial condition, we consider a Gaussian-like aggregate of tumour cells centered at (x 0 , y 0 ) = (−35, −41), situated in the left-bottom part of the brain slice.To be specific, M 0 = e (x−x 0 ) 2 +(y−y 0 ) 2 8 .  Moreover, Figure 2 shows the initial tissue density estimated with (39).In particular, yellow areas refer to regions where the fibers are highly aligned and, thus, the value of F A(D(x)) is closer to one, while black-red areas refer to more isotropic regions, where the fibers are randomly distributed [32].
We present different sets of simulations to obtain insight into several features characterising the proposed approach.In detail, (A) we consider the model for σ = 0 and we evaluate the effects of the variation of λ 1 and λ 0 on tumour evolution; (B) we fix the value of λ 0 and λ 1 and we assess the effects of the variation of σ on tumour evolution, i.e., the role of the stochastic parameter in the overall dynamics; (C) we consider different combinations of λ 1 and σ and we show how their respectively effects merge; (D) following the approach proposed in [12], we discuss the effects of λ 0 , λ 1 , and σ on the estimation of the onsets of malignant transformation from low grade to high grade gliomas.
Starting from the numerical test (A), we analyse the effects of varying λ 1 (referring to it as experiment A.1) and λ 0 (referring to it as experiment A.2).These experiments are motivated by the fact that obtaining a clear biological estimation for λ 0 and, especially, for λ 1 is quite difficult.Thus, understanding the impact of their variation becomes a fundamental point to address.As described in Section 3.1.2,λ 0 refers to the basal turning frequency of an individual cell, while λ 1 takes into account the role of the receptor dynamics in the evolution.Recalling the expression of the turning rate λ(z), we could describe the constant parameters λ 0 and λ 1 as the weights of the receptors-independent and receptors-dependent cell turning, respectively.Starting from the analysis on the parameter λ 1 and in line with some studies concerning the effects of its variability [50] on tumour evolution, we consider the range λ 1 ∈ [−5, 5] (s −1 ) and we assess the effects of changes in both its sign and modulus.Considering that the turning rate λ(z) = λ 0 − λ 1 z has to be non-negative, we should ensure that λ 0 ≥ λ 1 z, meaning that • if λ 1 ≥ 0, the non-negativity is ensured for λ 0 ≥ λ 1 /2; • if λ 1 ≤ 0, the non-negativity is ensured for λ 0 ≥ λ 1 .
Thus, to obtain reasonable values of the turning rate, we should assume λ 1 ≤ λ 0 .Although we are aware that negative values of these parameters are not sustained by biological observations, we also include them in our analysis because we want to assess the sensitivity of our results to these parameter changes.In Figure 3, we firstly show the evolution of the tumour density over time in the limit case in which λ 0 = λ 1 = 0.8 (s −1 ).37) with the parameters listed in Table 1 and for λ 0 = λ 1 = 0.8 (s −1 ).The tumour evolution is shown after 200, 400, and 600 days.
We notice how cells spreading is highly influenced by the underlying fiber structure.Cells clearly tend to move along preferential directions, determined by the fiber bundles, and this gives rise to a heterogeneous tumour mass with an irregular shape, which is a common characteristic for this kind of brain tumours.
Referring to the tumour situation at the last time step, i.e., after 600 days, we compare the tumour evolution for different values of the parameter λ 1 ∈ [−5, 5] (s −1 ), as described in experiment A.1.Results are shown in Figure 4.The main effect of varying λ 1 consists in obtaining a greater or lower level of heterogeneity in the distribution of the tumour cells inside the tumour mass.The external border of the neoplasia, in fact, does not seem to be particularly affected, while the internal dissemination of the cells shows evident changes when λ 1 varies from largenegative values to large-positive values.In particular, clear differences with respect to the case λ 1 = 0 can be observed for quite large values of the parameter (|λ 1 | > 1), while the evolution is qualitatively similar in the cases |λ 1 | < 1.Such differences can be better observed in Figure 5, where the differences between the solution of system (37) for λ 1 = 0 (s −1 ) and the solution of the same system for the different values of λ 1 used in Figure 4 are shown.The impact of λ 1 variation can be immediately grasped.There is a clear difference in the spreading inside the tumour mass and in the cell response to the anisotropy of the brain tissue.The impact becomes stronger when λ 1 increases in modulus, and especially for |λ 1 | > 1.In this case, in fact, the haptotactic component of the dynamics is stronger (in an attractant or repellent way, depending on the sign of λ 1 ) and, thus, the heterogeneity of the underlying brain tissue have a larger impact on the dynamics.The mechanism that drives cell migration along the tissue structure can be visualised in details in Figure 6, where the leading eigenvector of the tensor D T (x) (related to the fiber direction) is plotted together with the differences in the tumour density at 600 days for λ 1 = 5 (s −1 ) and λ 1 = −5 (s −1 ).Recalling the expression given for the tissue density (39), from the left plot of Figure 6 we notice that, where the fibers are strongly aligned (e.g.along the central vertical bound), we obtain negative values of the difference M λ1=−5 − M λ1=0 .Here, in fact, the gradient of tissue Q driving the haptotactic movement is bigger and, due to the negative value of λ 1 , cells tend to avoid this area, moving away from it.Conversely, looking at the right plot of Figure 6, we obtain exactly the reverse behaviour.In fact, the positive value of λ 1 leads to a much stronger haptotactic movement towards these fiber bundles.Thus, the difference shows positive values in the same regions described above.
We then test the effect of varying the parameter λ 0 , as described in experiment A.2. Results of this test are shown in Figure 7, where the difference between the solution of Figure 4: Experiment (A.1).Numerical simulations of equation (37) with the parameters listed in Table 1, λ 0 = 0.8 (s −1 ) and for different values of λ 1 (s −1 ).The tumour evolution is shown after 600 days.
(37) for λ 0 = 0.8 (s −1 ) and the one obtained for λ 0 varying in the interval [0.25, 5](s −1 ) are illustrated.We observe two different trends for λ 0 ≥ 0.8 or λ 0 ≤ 0.8.Smaller values of the parameter lead to a larger spreading of the tumour cells with respect to the case λ 0 = 0.8, while larger values of it lead to a reduced invasion of the tumour mass.In fact, smaller values of λ 0 mean a reduced random turning of the cells, thus a greater persistence in their migration, which macroscopically translates into a large spread.Instead, larger values of λ 0 imply a larger frequency of cell turning and, thus, a macroscopic lower degree of persistence and spread in the tissue.In particular, the main difference is in the region of the outer rim of the neoplasia.Concerning the numerical test (B), we fix λ 0 = λ 1 = 0.8 (s −1 ) and we vary the value of the parameter σ relating to the variability of the cell velocity in the microscopic model (12) and, thus, leading the additional diffusion term appearing in the macroscopic model (37).Results of the simulations for σ ∈ [0.01 − 0.2] (mm 2 •s −1 ) are shown in Figure 8.As expected from equation ( 37), the effect of the parameter σ consists of a larger spread of the tumour cells inside the brain tissue.In particular, the larger the value  .1).Differences between the solution of system (37) with λ 1 = 0 and the solution obtained for λ 1 varying in the interval [−5, 5] (s −1 ).Results are shown after 600 days.Here λ 0 = 0.8 (s −1 ), while the remaining parameters are taken from Table 1.
of σ is, the stronger the diffusion phenomenon characterising glioma cells appears.
For large values of σ, we observe more regular tumour borders and a more isotropic cell migration because the additional diffusion term does not depend on the diffusion tensor (35).These features can be better appreciated in Figure 9, where the differences between the solution of equation ( 37) for σ = 0 and for σ = 0 are shown.This figure clearly depicts an extensive and more homogeneous diffusion of the tumour mass for large values of σ.We obtain, in fact, negative values of the differences only in areas inside the tumour core (due to the balance between a faster spread and the same cell proliferation rate), while positive differences in the areas around the tumour border.In particular, comparing the first rows of Figures 9 and 7, we notice that the increase of σ values has an effect similar to the decrease of λ 0 values, i.e., a larger tumour spread in the area of tumour outer rim.It is interesting to observe how the same macroscopic cell behaviour is obtained from two different microscopic processes.
In fact, increasing σ allows for a stronger effect of the stochastic component related to the variation of cell velocity, while decreasing λ 0 reduces the random turning of the cells and determines a greater persistence in their direction of migration.
Referring to test (C), we analyse the interplay between the effects of the parameters  .1).Differences between the solution of system (37) with λ 1 = 0 and the one obtained for λ 1 = −5 (s −1 ) (left plot) and λ 1 = 5 (s −1 ) (right plot) after 600 days.The differences are plotted against the fiber direction.).This is still an effect of the higher value of λ 1 , here set at λ 1 = 5 (s −1 ).Finally, the combination of low values for both σ and λ 1 used for scenario (C.3) determines a smoothness of the internal distribution of tumour cells as well as a reduced cell spread in the healthy tissue.
For the last test (D), we discuss the onsets of malignant transformation from low grade glioma (LGG) to high grade glioma (HGG) in relation to the possible variations of the parameters λ 0 , λ 1 and σ.
LGGs are usually slowly-growing, infiltrative tumour with a very unpredictable clinical course.Most LGG patients face transformation of their tumour into higher grade one, with a worse prognosis.This process is known as malignant transformation and it is usually defined on the basis of contrast enhancement on MRI scans or histopathological evidences.In line with the approach proposed in [12], we estimate the time instant τ OSM of the onset to the malignant transformation of cells into a more aggressive high grade tumour.The main aim of the proposed experiment (D) consists in showing how our approach is able to replicate the same qualitative behaviours of [12] (where a comparison with patient data is proposed), but with a more detailed and precise description of the microscopic processes related to cell migration.Specifically, τ OSM is defined as the first time instant at which the LGG cell density becomes greater than a certain threshold M crit , which we set to 0.6K M [12].We run several numerical tests varying one parameter at the time and estimating the resulting time of onset of the malignancy.Table 2 collect the results of these experiments.We observe that the parameter λ 1 seems to not have such an evident impact on the time of onset of malignancy.In fact, τ OSM varies only of ± 28 days.Instead, both λ 0 and σ strongly affect the estimation of τ OSM .In Figure 11 the estimated values of τ OSM with respect to λ 0 and σ are plotted together with the corresponding interpolant curves, showing the trends of τ OSM = τ OSM (λ 0 ) (left plot of Figure 11) and τ OSM = τ OSM (σ) (right plot of Figure 11).Increasing the value of λ 0 leads to a reduction of the time τ OSM at which LGG turns into HGG, while increasing σ has the reverse effect, i.e., it leads to an increase of τ OSM .The parameter λ 0 is, in fact, related to the tumour responsiveness to the tissue structure, and large values of this parameter refer to a loss of responsiveness, which is a common characteristic in HGG.Moreover, observing that the overall diffusion coefficient of tumour cells in equation ( 37) is proportional to 1 λ0 + σ 2 2 , increasing σ (or equivalently decreasing λ 0 ) corresponds to an increase of this diffusion coefficient.Thus, comparing these results with the ones shown in [12] (e.g.see Figure 7 in there), we notice a good qualitative agreement between them and a similar behaviour for the evolution τ OSM .We would like to remark that this is only a first possible approximation for the estimation of τ OSM and we are aware that there are several other factors involved in the definition of the transformation from LGG to HGG, apart from the increase in the tumour density.Surely, the tumour density values have an evident impact on the definition of τ OSM , however, from a mathematical point of view, it is difficult to provide a formal definition for it.Thus, as a first attempt, we decide to rely on the definition given in [12] for τ OSM , leaving its possible extensions for future works.

Discussion
To the best of our knowledge, this is the first hierarchical stochastic model in which piecewise diffusion Markov processes are used to describe glioma cell motion within a multiscale framework.We start with the description of glioma cell movement at the microscopic scale using a PDifMP, which combines a stochastic model for cell motility and a deterministic one for cell migration.The latter looks at the response of glioma cells to external and environmental cues.The extended generator of the formulated PDifMP takes the form of an integro-differential equation in all the involved variables.Its solution yields the density of the transition probability of the Markov process.Using scaling arguments, we then obtain the equation describing the evolution of the tumor density at the macroscopic level.In this way, our approach allows us to take into account the macroscopic level properties as well as the features characterising the microscopic processes.Using numerical simulations of the macroscopic setting we analyse the role and influence of both the parameters involved in the jump rate function λ of the PDifMP and the parameter σ related to the stochastic variability in the cell velocity.In particular, we observe how the parameter λ 0 at the microscopic scale promotes a major spreading of the tumour mass inside the brain tissue, regardless of the specific brain structure, while λ 1 relates to cell responsiveness to the guided movement along the brain fibers.The fully detailed formulation of glioma cell motion with the PDifMP allows us to observe that the jump rate function determines the distribution of the waiting times of the process being in a particular state.Thus, for a constant jump rate (λ = λ 0 ) there is no influence of the microenvironment on the motion and a larger frequency of cell turning determines at the macroscopic scale a reduced migration along the fibers.Instead, including the term λ 1 z results in an increase in reorientations in response to the brain structure and, thus a visible heterogeneity inside the tumour bulk.A particularly interesting result is obtained by comparing the numerical experiments A.2 and B. In fact, we show how a similar macroscopic behaviour -large cell spreading around the outer rim of the tumour -can result from two different sources at the microscopic level: either from increasing the value of σ and thus the diffusion of cells, or from reducing the value of λ 0 and thus the random cell rotations, resulting in higher cell persistence .
With respect to well-known multiscale models of this type [32,34,33], in the present note we also include a further novel aspect concerning the transition to malignancy of the tumour mass.In particular, by accepting the hypothesis that the loss of responsiveness of glioma cells to the tissue structure can be seen as a sign of the transition from LGG to HGG, we numerically show that the time at which this transition happens can be estimated with our approach and it is highly influenced by the parameters σ and λ 0 .The obtained results are perfectly in line with the ones presented in [12], confirming the reliability of the proposed approach.
With our work, we aim at emphasising how the use of PDMP or PDifMP for the description of the phenomena leading cell movement is of paramount importance for rigorously modelling the cellular scale processes.An interesting point would concern a numerical comparison of the cell behaviours at the different scales (microscopic and macroscopic) with either the deterministic and the stochastic formulation.Moreover, in the present notes, glioma cell motion is described in relation to the binding with the tissue, but the proposed approach can be extended in order to incorporate other biologically relevant aspects of tumour progression.For instance, following [21], the influence of microenvironmental acidosis on glioma cell migration and the consequent pH-repellent chemotactic process can be considered.This could be done assuming different expressions for the jump rate function of the PDifMP, e.g.allowing its dependence on different interactions between cells and microenvironment or relating it to the tumour response to treatments.Another interesting direction for future development concerns the modification of the jump process using stochastic differential equations to model not only jumps in the cell velocity, but also jumps in the position, trying to recover the typical feature of tumour recurrence in different (and quite far from the original tumour location) regions of the brain.Finally, here we propose a first possible way to analyse the transition to malignancy.However, as stated in the above section, this process is much more complex and we are working towards the development of an interdisciplinary study in which an extension of our approach could be used to shed light on the intricate biological processes underlying this transition.
is the stochastic flow of the continuous first component of (Ut) t∈[0,T ] .Starting at T0 = 0 with initial value u0 = (s0, v0) ∈ E, the process φ(t, u) represents the solution of a sequence of SDEs over the consecutive intervals [Ti, Ti+1) of random length.At each random point Ti ∈ [0, T ], i ≥ 1, there are newly updated initial values ui = (si, vi) ∈ E, where si serves as the initial value and vi as a parameter in the following SDE defined on the interval [Ti, Ti+1):

Figure 8 :
Figure 8: Experiment (B).Numerical simulations of equation (37) with parameters listed in Table 1 and for different values of σ.The tumour evolution is shown after 600 days.Values of σ are expressed in mm 2 • s −1 .The figures referring to the cases σ = 0.15 and σ = 0.2 are shown on a less zoomed region to better assess the tumour invasion in the tissue.

Figure 9 :
Figure 9: Details of experiment (B).Differences between the solution of system (37) for σ = 0 and the one obtained for σ ∈ [0.01, 0.2] (mm 2 • s −1 ).Results are shown after 600 days.The remaining parameters are taken from Table 1.The figures referring to the cases σ = 0.15 and σ = 0.2 are shown on a less zoomed region to better observe the tumour invasion in the tissue.

Figure 11 :
Figure 11: Experiment (D).Estimation of the time of onset of malignant transformation τ OSM for different values of the turning rate λ 0 and the free parameter σ.The remaining parameters are taken as in Figure 3.

Table 2 :
Estimations of the onsets of malignant transformation τ OSM for different values of λ 1 , σ, and λ 0 .