Treatment-induced shrinking of tumour aggregates: A nonlinear volume-filling chemotactic approach
Luis Almeida, Gissell Estrada-Rodriguez, Lisa Oliver, Diane Peurichard, Alexandre Poulain, Francois Vallette
TTreatment-induced shrinking of tumour aggregates: Anonlinear volume-filling chemotactic approach
Luis Almeida ∗ Gissell Estrada-Rodriguez †‡ Lisa Oliver § Diane Peurichard ¶ Alexandre Poulain (cid:107)
Francois Vallette ∗∗ Abstract
Motivated by experimental observations in 3D/organoid cultures derived from glioblastoma,we develop a mathematical model where tumour aggregate formation is obtained as the re-sult of nutrient-limited cell proliferation coupled with chemotaxis-based cell movement. Theintroduction of a chemotherapeutic treatment induces mechanical changes at the cell level, withcells undergoing a transition from rigid bodies to semi-elastic entities. We analyse the influ-ence of these individual mechanical changes on the properties of the aggregates obtained atthe population level by introducing a nonlinear volume-filling chemotactic system of partialdifferential equations. The elastic properties of the cells are taken into account through theso-called squeezing probability , which allows us to change the packing capacity of the aggregates,depending on the concentration of the treatment in the extracellular microenvironment. Weexplore two scenarios for the effect of the treatment: firstly, the treatment acts only on themechanical properties of the cells and, secondly, we assume it also prevents cell proliferation.A linear stability analysis enables us to study the ability of the system to create patterns. Weprovide numerical simulations in 1D and 2D that illustrate the shrinking of the aggregates dueto the presence of the treatment.
Cell migration in the extracellular microenvironment (ECM) and the organisation of cells in responseto chemical and mechanical cues are successfully studied using continuum descriptions based ondifferential equations [1, 42]. In a continuous setting, the chemotactic behaviour of cells ( i.e. theability to move along a chemical gradient) is often modelled using a Keller-Segel system of equations[27]. This model was originally proposed for pattern formation in bacterial populations but turnedout to be pertinent to describe a wide variety of self-organisation behaviours [11, 6, 60, 28, 57]. ∗ Sorbonne Universit´e, Inria, CNRS, Laboratoire Jacques-Louis Lions, F-75005 Paris. email:[email protected] † Sorbonne Universit´e, Laboratoire Jacques-Louis Lions, F-75005 Paris. email: [email protected] ‡ Corresponding author § Centre de Recherche en Canc´erologie et Immunologie Nantes-Angers, UMR 1232, Universit´e de Nantes. email:[email protected] ¶ Sorbonne Universit´e, Inria, CNRS, Laboratoire Jacques-Louis Lions, F-75005 Paris. email: [email protected] (cid:107)
Sorbonne Universit´e, Laboratoire Jacques-Louis Lions, F-75005 Paris. email: [email protected] ∗∗ Centre de Recherche en Canc´erologie et Immunologie Nantes-Angers, UMR 1232, Universit´e de Nantes. email:[email protected] a r X i v : . [ n li n . PS ] J u l ifferent variations of the Keller-Segel model have been adopted in order to better understand theway cells aggregate [8, 4, 40, 47, 16].Cancer cells have been shown to respond to chemical and mechanical signals from components ofthe tumour micro-environment (TME) and vice versa , cells in the TME acquire a more tumourigenicphenotype, which would sustain tumour growth. The interaction of tumour cells with the TME hasbeen the subject of recent biological surveys [22, 18]. Many in vitro (and ex vivo ) experimentshave shown that cells that are cultured on ECM often have a tendency to form aggregate patternsthat are dependent on the particular cell lines and physical properties of the media [20]. The exactconsequences of the dynamic interplay between heterogeneous cellular entities and their response toalterations in the TME have not yet been elucidated. In particular, little is known about the roleof mechanics in the spatial organization of the tumour spheroids. Biological evidence presented in[33, 29, 12] suggests that the formation of aggregates in glioma cells can be explained through achemotaxis process, rather than, e.g. , cell-cell adhesion. In this paper, we follow the chemotacticapproach to explain the formation of glioma aggregates. For the case of breast cancer cells, a recentreport [9] has proposed a chemotaxis-based explanation for spheroid formation based on theoreticalanalysis and numerical simulations of the Turing instabilities of such systems.Inspired by experimental observations in 3D/organoid cultures derived from freshly operatedGlioblastoma (GBM), which reproduce in vivo behaviours as described in [37] (see Section 2 formore details), in this paper we explore a simple setting where GBM aggregate formation is the resultof nutrient-limited cell proliferation coupled with a chemotaxis-based cell movement. A chemother-apeutic treatment introduced after the formation of the aggregates induces mechanical changes atthe cell level. We study the influence of these individual mechanical changes on the characteristicsof the aggregates obtained at the population level.GBM are solid tumours characterised by intra- and inter-tumoural heterogeneity and resistanceto conventional treatments that result in a poor prognosis [31]. They are the most common andaggressive primary brain tumour in adults. Standard treatments include surgical resection (whenpossible), combined with radiotherapy and chemotherapy using the DNA alkylating agent Temo-zolomide (TMZ) [49]. In fact, the overall survival of treated patients is about 15 months versus 3months without treatment, with fewer than 5% of patients surviving longer than 5 years [48].One reason behind this relative therapeutic failure is the poor response of GBM tumours tothis chemotherapeutic treatment due to their plasticity. Several studies have looked at the geneticcompounds of TMZ-resistant cells focusing on the genes responsible for DNA mismatch repair protein[51], while other studies focused on spatial and temporal variations in signalling pathways, whichlead to functional and phenotypic changes in GBM [38]. The communication between the tumourcells and the TME as well as the properties of the ECM have a large impact on tumour evolutionand invasion, as shown in recent studies [15, 7, 59]. From a biological and medical perspective,it is difficult to investigate the connections between clinically observable glioma behaviour andthe underlying molecular and cellular processes. The challenge is to integrate the theoretical andempirical acquired knowledge to better understand the mechanisms and factors that contribute toGBM resistance to treatment. In this context, mathematical models provide useful tools towardsidentifying dependencies between different phenomena and how they are affected by the differenttherapeutic strategies.Much effort has been dedicated to the modelling of GBM formation and invasion of the surround-ing tissue, as well as to improving diagnosis and treatment. The exhaustive review [2] discussesdifferent modelling approaches as well as some of the main mechanisms that are observed in GBMformation and invasion.In this work, we explore two scenarios: the case where the treatment only acts on the mechanicalproperties of the cells, and the case where it also prevents cell proliferation [30]. We adopt a macro-2copic approach where cells are represented by their macroscopic density and are supposed to movein the environment via chemotaxis, i.e. towards zones of high concentrations of a chemoattractantthat is produced by the tumour cells. Moreover, cell proliferation is assumed to depend on the localconcentration of nutrients available. Finally, we suppose that when the treatment - represented byits continuous concentration - is introduced, it diffuses in the environment and is naturally consumedby the cells.Under these hypotheses, which are motivated by the experimental results discussed in Section 2,we obtain a nonlinear volume-filling Keller-Segel model for the cell density, coupled with reactiondiffusion equations for the chemoattractant and treatment concentrations. Moreover, we provide alinear stability analysis that enables us to study the ability of the system to generate patterns, andwe provide numerical simulations in 1D and 2D.The paper is organised as follows. In Section 2 we describe the in vitro experiments and the mainexperimental observations that motivated our model. Section 3 is devoted to the description of themodel for the first part of the experiments (without the treatment, Section 3.1) and the second part,when the treatment is introduced (Section 3.2). In Section 4, we present the stability analysis foreach of these two parts, and Section 5 is devoted to numerical simulations. Section 5.1 presents thenumerical scheme, Section 5.3 presents the results in 1D including a discussion on the comparisonbetween numerical experiments and theoretical predictions of the stability analysis. Section 5.4shows the 2D simulations. Finally, we discuss the results and give some perspectives in Section 6. To address the question of the response of GBM cells to TMZ treatment, we took advantage ofrecently developed 3D biosphere experiments, using GBM patient-derived cultures in a simple 3Dscaffold composed of alginate and gelatin [39]. GBMG5 cells were cultured at a concentration of4 × cells/ml for 14 days until the formation of cell aggregates could be observed, corresponding tothe first part of the experiments ( P1 ). After that, the cultures were treated with 100 µ M TMZ fortwo hours once every week (second part of the experiments, P2 ). Over the 30 days, the proliferationwas determined counting 3 representative samples, and the cell number was determined as follows.The biospheres were dissociated by incubation for 3 min in 100 mM Na-Citrate and the cell numberand cell viability were determined using the Countess optics and image automated cell counter (LifeTechnologies). In addition the aggregates were photographed to analyse their morphology, andthe diameter of the cell aggregates were measured from pictographs using FIJI. To determine thediameters of these cell structures, more than 200 cell aggregates were measured.We show in Figure 1 (I) the mean length (in µ m) of the cell aggregates computed from themicroscopic images without TMZ treatment (round markers), and with TMZ weekly administered(squared markers). Figures 1 A and B show typical microscopic images of the spheroids at day 24,without and with weekly TMZ treatment respectively. In Figure 1 (II), we show the evolution ofthe cell number as a function of time, where we do not observe big changes in cell number with andwithout TMZ, once the carrying capacity of the system is reached.Using clinically relevant concentrations [44] the total number of cells in the biospheres does notseem to be significantly impacted by the TMZ treatment. However, the mean size of the GBMG5cell aggregates decreases when TMZ is introduced weekly as compared to control cultures (Figure 1(I)), suggesting that in the presence of the treatment, GBMG5 cells tend to self-organize into smallerand more compact cell clusters.Another observation supporting a tendency for GBM cells to form more compact structures withhigher intracelluar adhesion under this type of treatment is the increased expression of claudin, a3igure 1: Biological experiments of GBMG5 cells cultured in 3D scaffolds with or without 100 µ MTMZ. (I) Evolution of the mean cell aggregates diameter determined from the microscopic images asfunction of time, without treatment (circle markers) and with a weekly TMZ (square markers). (A)Typical microscopic image on day 24 of control cell aggregates without treatment, (B) microscopicimage at day 24 with 100 µ M of TMZ administrated weekly for 2 hours. (II) Evolution of the totalcell number in the biospheres as function of time, without treatment (circle markers) and with weeklyTMZ (square markers).Figure 2: Claudin expression marking the tight junctions between the cells.4arker of tight junctions formed between cells (Figure 2). In Figure 2 we have GBMG5 cells thatwere cultured for 14 days until the formation of cell aggregates could be observed, and the cultureswere treated with different doses of TMZ until day 23. Furthermore, cell aggregates were fixed thenlabeled with an anti-claudin antibody. The degree of staining was determined in a double-blindexperiment. As one can observe from the images of Figure 2 (left), the cell aggregates seem to besmaller and more compact when TMZ is present compared to the control group (top left figure),and these cultures are associated with higher levels of claudin (see Figure 2 (right)).Inspired by these observations and results presented in [33, 29, 12], we propose a general modelfor the chemotaxis-driven formation of cell clusters and shrinking of the aggregates via the action ofa non-cytotoxic drug whose effect would be characterised by a change in the mechanical properties ofcells that would just modify their packing properties. Moreover, we also consider a second scenariowhere the treatment would also affect cell proliferation and compare the results obtained.
Motivated by the experiments described in Section 2, we assume that glioma cells have a chemotacticbehavior, i.e. they move in response to some signaling chemical (chemoattractant), which is secretedby themselves and diffuses in the environment. The chemotactic movement of cell populations playsa fundamental role for example in gastrulation [17], during embryonic development; it directs themovement of immune cells to sites of infection and it is crucial to understand tumour cell invasion [14]and cancer development [60]. Motivated by these applications, chemotaxis and related phenomenahave received significant attention in the theoretical community, see the reviews [23, 24].We suppose that in normal conditions (first part of the experiments described in Section 2,where there is no treatment), tumour cells proliferate and move via chemotaxis as described before.We suppose that cell proliferation is limited by the nutrients available in the environment, i.e. cellproliferation is only active as long as the local density does not exceed a given threshold correspondingto the carrying capacity of the environment. Moreover, in order to take into account the finite sizeof cells and volume limitations, cell motion is only allowed in locations where the local cell densityis much smaller than another threshold value corresponding to the tight packing state. In normalconditions, cells are supposed to behave as rigid bodies. In stressed conditions however, (second partof the experiments described in Section 2) we suppose that cells respond to the chemotherapeuticstress (induced by the presence of the treatment) by changing their mechanical properties andbehaving as a semi-elastic material.These hypotheses are modelled via a system of partial differential equations (PDE) which corre-sponds to a volume-filling chemotaxis equation [43] for the first part (to describe the self-organizationof cells into aggregates), and a “semi-elastic” volume-filling chemotaxis approach [56] for the secondpart, when the treatment is introduced.For convenience we denote the density of the population of cancer cells in the first part of theexperiments by u ( x , t ) and that of the second part by w ( x , t ). Here x ∈ Ω ⊂ R where Ω is abounded domain. The time t ∈ [0 , T ] where T = T + T represents the total time correspondingto the first and second parts. The main difference between these two populations is the change inthe elastic properties of the cells due to the presence of the treatment. If the concentration of thetreatment is zero, u ( x , t ) = w ( x , t ). Cells follow a biased random walk according to the distributionof the chemoattractant of concentration c ( x , t ) that is secreted by the cells. We start by detailing thedifferent components of the mathematical models corresponding to the two parts of the experimentsdescribed in Section 2. 5 ogistic growth model for cell proliferation As previously described, in order to take intoaccount the nutrient-limited population growth, cell proliferation is modelled by a logistic growthprocess. At the population level, we consider a source term f ( u ) in the PDE for the evolution of thecell density which depends nonlinearly on the local cell density u and reads: f ( u ) = ru (cid:16) − uu max (cid:17) . (1)Here, r > u max corresponds to the maximum density of the population,also referred to as the carrying capacity of the environment. Alternative cell kinetics could be con-sidered. For example, we can assume that the proliferation is also mediated by the chemoattractantconcentration such that f ( u, c ) = ruc (1 − u/u max ) [43]. In this paper we focus on the case given by(1). Chemoattractant dynamics
As described before, we suppose that cell aggregates spontaneouslyemerge as the result of a self-organization phenomenon of chemotaxis type. To this aim, we supposethat the cells themselves produce the signaling chemical (chemoattractant) that drives their motion.The chemical secreted is therefore supposed to be continuously produced by the cells at rate α > d >
0. It is further assumedthat the chemical has a finite lifetime and degrades at rate β >
0. The evolution of the chemicalconcentration c ( x , t ) is therefore given by the following reaction-diffusion equation ∂ t c = d ∆ c + αu − βc , (2)where u is the cell density. Treatment dynamics
We suppose that the treatment is introduced at the beginning of P2 . Thistreatment is supposed to diffuse in the environment with diffusion coefficient d , and to be consumedby the cells at rate δ . This is modelled by a reaction-diffusion equation for the drug concentration M ( x , t ): ∂ t M = d ∆ M − δw , where w ( x , t ) represents the cell density corresponding to the second part of the experiments. Weconsider different initial conditions for the drug: either uniformly distributed in the simulationdomain, or introduced in the center as a very steep Gaussian function (see Section 5). The classical Keller-Segel system of equations [27] describes how cells move along the gradient tolocal maxima of the chemoattractant. At the same time, this chemical, which is produced by thecells, promotes aggregation leading to the so-called “overcrowding scenarios”, and eventually, the celldensity may blow up in finite time (see the comprehensive reviews [23, 50] and references therein).In this paper, in order to take into account volume limitations and the finite cell size, we considera modified version of the Keller-Segel model called the volume-filling approach for cell motion. Thisapproach was introduced in [43], where the authors provided a detailed derivation in one dimensionas well as a comprehensive numerical study of the model.6 u(x,t) q ( u ( x ,t )) M =1 M =2 M =3 Figure 3: Squeezing probability for different values of γ M for ¯ u = 1.Moreover, we introduce the so-called squeezing probability , which describes the probability thata cell finds an empty space at a neighbouring location before moving (see [56]). It takes the form q ( u ( x , t )) = (cid:40) − (cid:16) u ¯ u (cid:17) γ M , for 0 ≤ u ≤ ¯ u , , otherwise , (3)where u ( x , t ) is the cell density, γ M ≥ squeezing parameter , M ≡ M ( x , t ) ≥ u is the crowding capacity which corresponds to the tightpacking state of the cells. The function q ( u ) satisfies the following properties, q (¯ u ) = 0 , < q ( u ) ≤ , and q (cid:48) ( u ) ≤ ≤ u < ¯ u . (4)Moreover, | q (cid:48) ( u ) | is bounded and q (cid:48)(cid:48) ( u ) ≤ γ M is chosen to be γ M = (¯ γ − M + 1 , (5)where ¯ γ is a positive constant, M = 0 when there is no drug in the environment (part P1 of theexperiments), and M ≡ M ( x , t ) when the drug is introduced (part P2 , described bellow). Suchchoice of γ M enables to take into account different forms of squeezing probability, correspondingto different mechanical behaviors of the cells. In Figure 3, we plot the squeezing probabilities asa function of the cell density, corresponding to different values of γ M . We see that when there isno treatment present in the environment ( M = 0, γ M = 1, blue curve of Figure 3), the squeezingprobability decreases linearly with the local cell density, corresponding to cells modelled as solidparticles.However for larger values of γ M (when the treatment is present, M > γ M >
1, red andyellow curves in Figure 3), the squeezing probability becomes a nonlinear function of the cell density,modelling the fact that in the presence of a drug, cells change their mechanical state to behave assemi-elastic entities that can squeeze into empty spaces. We refer to [55] for more details on the linkbetween the cells elastic properties and the squeezing probability.For convenience of notation, we distinguish the case without treatment ( M = 0) for P1 anddefine q ( u ( x , t )) = 1 − u ¯ u . (6)7he Keller-Segel model built with this specific squeezing probability (6) has been widely studied inthe literature, from modelling [43, 55], analytic [58, 21, 32, 16] and numerical [26] perspectives. Complete PDE system for the first part P1
Taking into account all the previous ingredients,the complete PDE system for part P1 (when no treatment is present in the environment) reads: ∂ t u = ∇ · ( d D ( u ) ∇ u − χ u φ ( u ) ∇ c ) + f ( u ) ,∂ t c = d ∆ c + αu − βc , (7)where the first equation describes the evolution of the cell density u and the second is the reaction-diffusion equation for the chemoattractant, previously introduced. The equation for u describesthe volume-filling chemotactic motion associated with the squeezing probability q . This equationhas been obtained as the hydrodynamic limit of the continuous space-time biased random walkmodel that corresponds to the squeezing probability q (see Appendix A for more details on thederivation). In the hydrodynamic limit, the cell density evolves according to a nonlinear transportdiffusion equation with source term, which corresponds to a volume filling Keller-Segel model withlogistic growth. As shown in Appendix A, the density-dependent diffusion coefficient D ( u ) and thechemotactic sensitivity φ ( u ) relate to q via D ( u ) = q − q (cid:48) u , φ ( u ) = q u . For this part P1 , where q is given by (6), these coefficients take the values D = 1 and φ ( u ) = u (1 − u/ ¯ u ). In equation (7), d , χ u , α, β are all positive parameters (see Appendix A for the computationof d and χ u ), and f ( u ) is the logistic growth source term previously mentioned and defined by (1).The PDE system is supplemented with the following zero-flux boundary conditions( d D ( u ) ∇ u − χ u φ ( u ) ∇ c ) · η = 0 , d ∇ c · η = 0 , (8)where η is the outer unit normal at ∂ Ω. The initial conditions are given by u ( x ,
0) = u , c ( x ,
0) = c . (9) We now describe the dynamics of the cell population when the drug is introduced (part P2 of theexperiments described in Section 2). As previously mentioned in Section 2 and motivated by theobservations in [39], where the treatment TMZ does not seem to induce cell death, we suppose thatthe drug only affects the elastic properties of the cells. For a cell density w ( x , t ) the squeezingprobability of part P2 is q ( w ( x , t ) , M ) = 1 − (cid:16) w ¯ u (cid:17) γ M . (10)Note that we have supposed that the crowding capacity ¯ u which corresponds to the tight packingstate remains unchanged from P1 to P2 . This corresponds to the hypothesis that the treatmentdoes not modify the volume of the cells but only change their elastic properties.The complete PDE system corresponding to the second part of the experiments therefore reads: ∂ t w = ∇ · ( d D ( w, M ) ∇ w − χ w φ ( w, M ) ∇ c ) + f ( w ) ,∂ t c = d ∆ c + αw − βc ,∂ t M = d ∆ M − δw , (11)8here f ( w ) is again the logistic growth source term given by (1), and the first equation has beenderived using the squeezing probability q defined by (10). Note that the carrying capacity u max remains unchanged in the two parts of the experiments: we have supposed here that the drug doesnot interact with the nutrients. Analogous to (7), the movement of the cells is described by achemotactic system where, in this case, the diffusion and chemosensitive coefficients depend also onthe concentration of the treatment. These modified coefficients will lead to changes in the size of theaggregates as shown in the numerical experiments in Section 5. The evolution of the concentration c is the same as in (7) where, in this case, the chemoattractant is produced by the new population w . Including proliferation also in this second part allows us to assess the effect of the treatment inthe population at earlier times, while the population of cells is still growing and aggregates are stillforming.Different initial conditions for the cell density and concentration of the treatment will be con-sidered as described in Section 5. Each initial condition for P2 will correspond to a density profilesolution of the system P1 at a given time (solution of (16)), i.e. w ( x ,
0) = u ( x , T ) for some given T .We will consider cases where the treatment is introduced on already formed and stable aggregates(steady state of (16)), but also cases where it is introduced at earlier times (during the formation ofthe aggregates, see Section 5). The initial condition for the concentration of the drug is consideredto be either homogeneously distributed in the simulation domain or introduced in the center. Remark 3.1.
In both parts of the experiments, we assume that the crowding capacity ¯ u is largerthan the carrying capacity u max . The carrying capacity is defined as the maximum population densitybeyond which there are not enough nutrients to support growth, while the crowding capacity describesthe maximum density in an aggregate depending on the space available. The system for part P1 (7) without source term is well known in the literature as the volume-filling Keller-Segel model, for which emergence of patterns has been characterised and is now welldocumented. Pattern formation refers to the phenomenon by which, after varying a bifurcation pa-rameter, the spatially homogeneous steady state loses stability and inhomogeneous solutions appear.In the following, we investigate in which parameter region we can expect instability of homogeneoussolutions, corresponding to the formation of patterns. The linear stability analysis followed here isclassical and follows the lines of [43, 56, 35].We first recall the two systems associated with the dynamics described in parts P1 and P2 ofthe model. Using the fact that in the first part of the experiments, the squeezing probability ischosen to be q ( u ) = 1 − u ¯ u , we have ∂ t u = ∇ · ( d ∇ u − χ u φ ( u ) ∇ c ) + ru (cid:16) − uu max (cid:17) ,∂ t c = d ∆ c + αu − βc , (12)where D = 1 , and φ ( u ) = u (cid:0) − u ¯ u (cid:1) . (13)For part P2 , when the squeezing probability function is given by (10), the system writes ∂ t w = ∇ · ( d D ( w, M ) ∇ w − χ w φ ( w, M ) ∇ c ) + ˜ rw (cid:16) − wu max (cid:17) ,∂ t c = d ∆ c + αw − βc ,∂ t M = d ∆ M − δw , (14)9ith diffusion and chemotactic coefficients given by D ( w, M ) = 1 + ( γ M − (cid:16) w ¯ u (cid:17) γ M and φ ( w, M ) = w (cid:16) − (cid:16) w ¯ u (cid:17) γ M (cid:17) . (15) To get a deeper insight to the behaviour of the system we introduce the characteristic values of thephysical quantities appearing in the models. Denoting by X and T the macroscopic units of spaceand time, respectively, such that ¯ x = x X , ¯ t = tT , then we choose(¯ x , ¯ t ) = (cid:32)(cid:114) βd x , βd d t (cid:33) . Using these new variables, the dimensionless systems write for part P1 ∂ t u = ∇ · ( ∇ u − Aφ ( u ) ∇ c ) + r u (cid:18) − uu max (cid:19) ,ζ∂ t c = ∆ c + u − c . (16)Similarly, we obtain for P2 θ∂ t w = ∇ · ( D ( w, M ) ∇ w − Bφ ( w, M ) ∇ c ) + ˜ r n (cid:16) − wu max (cid:17) ,ζ∂ t c = ∆ c + n − c ,m∂ t M = ∆ M − δ w , (17)where A = χ u d , r = d rd β , ζ = d d , θ = d d , B = χ w d , ˜ r = d ˜ rd β , m = d d , δ = δd . (18)The parameters ζ and m are assumed to be small since the chemoattractant and the chemother-apeutic treatment diffuse faster than the cells. On the other hand, θ (cid:39) u and w are assumed to be diffusing at similar rates. In the following, we state the linearstability for both systems in separate sections. We first consider the system (16), which can be re-written in a more general form as ∂ t u = ∇ · ( ∇ u − Aφ ( u ) ∇ c ) + f ( u ) ,ζ∂ t c = ∆ c + g ( u, c ) , (19)where φ ( u ) is given in (13), f ( u ) = r u (1 − uu max ) and g ( u, c ) = u − c . This system is subject touniformly distributed initial conditions and zero-flux boundary conditions as in (8).The main result in this section is the following theorem, which gives the conditions for patternformation for the system (19). 10 heorem 4.1. Consider ( u ∗ , c ∗ ) a spatially homogeneous steady state. Then pattern formationfor the system (19) with zero flux boundary conditions (8) and initial data (9) is observed if thefollowing conditions are satisfied, f ∗ u + ζ − g ∗ c < , f ∗ u g ∗ c > , ζ − g ∗ c + f ∗ u − ζ − g ∗ u Aφ ( u ∗ ) > ,g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ ) > (cid:112) f ∗ u g ∗ c . (20) The critical chemosensitivity is given by A c = 2 √ r + 1 + r u max (cid:0) − u max ¯ u (cid:1) , (21) and for A > A c patterns can be expected. The wavemodes k are in the interval defined by k = − m − (cid:112) m − f ∗ u g ∗ c < k < k = − m + (cid:112) m − f ∗ u g ∗ c , (22) where m = − ( g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ )) .Proof. See Appendix B for the proof of this theorem.
We now consider the system (17) which corresponds to the second part of the experiments, whenthe treatment is introduced. The parameter range where patterns are observed is summarised inthe following theorem.
Theorem 4.2.
Consider ( w ∗ , c ∗ , M ∗ ) a spatially homogeneous steady state. Also, consider (17)with zero flux boundary conditions given by ( d D ( w, M ) ∇ w − χ w φ ( w, M ) ∇ c ) · η = 0 , d ∇ c · η = 0 , d ∇ M · η = 0 , and initial conditions w ( x ,
0) = u ( x , T ) ,c ( x ,
0) of P2 is equal to c ( x , T ) from P1 ,M ( x ,
0) = (cid:40) constant in Ω , or Ce ( x − x σ , (23) where C is the amplitude, x is the centre of the Gaussian and σ is the width. Then, the criticalchemosensitivity is given by B c = 2 (cid:112) ¯ r D ( u max , M s ) + D ( u max , M s ) + ˜ r u max (cid:16) − (cid:16) u max ¯ u (cid:17) γ Ms (cid:17) , where D ( u max , M s ) = 1 + ( γ M s − (cid:16) u max ¯ u (cid:17) γ Ms . Patterns can be expected if
B > B c and the wavemodes k are in the interval defined by k = − ¯ m − (cid:112) ¯ m − D ( u max , M s )( f ∗ w g ∗ c )2 D ( u max , M s ) < k < k = − ¯ m + (cid:112) ¯ m − D ( u max , M s )( f ∗ w g ∗ c )2 D ( u max , M s ) , for ¯ m = − ( D ( u max , M s ) g ∗ c + g ∗ w Bφ ( u max , M s ) + f ∗ w ) . roof. The proof of this theorem can also be found in Appendix B.
Remark 4.1.
For the case of dimensions, we can rewrite the systems (16) and (17) using polarcoordinates ( ρ, θ ) where we use the transformation x = ρ sin θ , y = ρ cos θ and the Laplace operatoris now given by ∆ p = R ∂∂ρ (cid:16) ρ ∂∂ρ (cid:17) + ρ ∂ ∂θ where R is the radius of the domain. The eigenvalueproblem (45) is now written as − ∆ p ψ k = k ψ k with boundary conditions ∂ψ k /∂ρ = 0 at ρ = R .The eigenfunctions are obtained by separation of variables and are given by ψ k ( x, y ) = R ( ρ )Θ( θ ) .Here Θ( θ ) = e isθ = A cos( sθ ) + B sin( sθ ) for some s ∈ Z . The radial part R ( ρ ) is given in termsof Bessel functions R ( ρ ) = J s ( kρ ) (see [46]) where k = c s,p R and c s,p denotes the the p th zeroderivative of J s , which is a first kind Bessel function of order m . Finally we can write ψ s,pk ( ρ, θ ) = J (cid:0) c s,p R ρ (cid:1) ( A cos( sθ ) + B sin( sθ )) . The stability analysis reveals that several competing effects control the system’s ability to createpatterns (aggregates). The criteria obtained both in P1 or P2 show that the chemotactic sensitivitymust be large enough to compensate the smoothing effect of the diffusion term and of the logisticgrowth. On the other hand, one can observe from the bifurcation formulae that the ratio u max ¯ u (carrying capacity vs density of the tight packing state) plays an important role in the emergenceof patterns: larger values lead to more aggregated states. These results show that the logisticgrowth term has an intrinsic smoothing property, i.e. , it tends to force the density to equate thecarrying capacity, while the chemotactic term acts as a attractive force and creates zones of higherdensity (recall that u max < ¯ u ). The aggregates are an expression of a balance in between these twocompeting effects, which are completely characterised by the stability criterion. A,B max = 1/argmax k { Re(l(k)) } = 1 = 2 = 3 = 5 = 50 (a) A,B max = 1/argmax k { Re(l(k)) } = 1 = 2 = 3 = 5 = 50 (b) Figure 4: Wavenumbers for different values of γ M when (a) r = 0 . u max = u = 0 .
5; and for (b) r = 0 . u max = u = 0 . k c .In order to give more insights about the size of the emerging patterns, we show in Figure 4 thevalues of the inverse of the maximal wavenumber as function of the chemotactic sensitivity (denotedby A for part P1 and B for part P2 ), for different values of the exponent γ M . We recall that γ M = 1in P1 (blue curve) where cells act as rigid bodies and γ M > P2 ). Figure 4a shows the results for growth rate r = 0 . u max = 0 .
5, Figure 4b for r = 0 .
05 (slower growth) and u max = 0 . u = 1. Themaximal wavenumber corresponds to the most unstable mode, i.e. the perturbed wave that will12row the fastest. Therefore, the inverse of this maximal wavenumber is directly related to the sizeof the emerging patterns. As one can observe, an increase in the chemotactic sensitivity parametercorrelates with a decrease in the observed pattern size, suggesting that the aggregates are smaller:larger chemotactic sensitivity leads to more aggregated clusters. Moreover, the aggregate size alsodecreases as cells pass from rigid bodies to semi-elastic entities (when γ M increases). This is dueto the fact that for larger values of γ M , cells are more easily deformed and can aggregate moreefficiently than when they behave as rigid spheres.When we increase the ratio u max ¯ u (compare Figure 4a and 4b), we observe that the critical valueof the chemosensitivity above which patterns are generated is larger than for smaller ratios u max ¯ u .These results highlight once again the smoothing effect of the logistic growth: When the cell tightpacking density is unchanged, decreasing the carrying capacity of the environment enhances celldeath in the aggregates formed by chemotaxis, where cells try to reach the tight packing state. Inthis case, the critical chemosensitivity value must be large enough to compensate for the cell deathinduced by the logistic growth. Moreover, we observe that larger ratios u max ¯ u induce less influence ofthe parameter γ M . The cell aggregation abilities are mainly driven by the chemosensitivity intensityand not so much by the cell mechanical properties for large values of u max ¯ u . In addition to the analytic results obtained in Section 4, we present numerical simulations for thetwo problems (16) and (17). This allows us to investigate the behavior of the solution of the modelsfor different scenarios and range of parameters. It is well-known that a standard discretization ofthe Keller-Segel models can lead to nonphysical solutions due to the convective term. Here, we focuson a numerical method that preserves the non-negativity of the cell density using the upwind finiteelement method described in [10] for the simulation of the Cahn-Hilliard equation.The calculation of the chemotactic coefficient follows the lines of [3]. Indeed, the finite volumescheme proposed in [3] is identical to the numerical method presented in this paper in dimensionone. However, in higher dimensions, and since we also use a finite element method, the numericalscheme presented in [3] differs from the one in this section as detailed in Appendix C.
The domain Ω is discretized by a mesh T h constituted of N el piecewise linear finite elements. Wedenote the total number of nodes by N h and the standard P finite element space by S h . Thestandard L inner product is indicated by ( · , · ) and the lumped scalar product defined in AppendixC by ( · , · ) h . Given N T ∈ N (cid:63) , let ∆ t := T /N T be the time-step where T is the time corresponding tothe end of the simulation. Let t n := n ∆ t , n = 0 , . . . , N T − ∂u h ∂t ≈ u n +1 h − u nh ∆ t . We define u nh ( x ) := N h (cid:88) j =1 u nj χ j ( x ) , and c nh ( x ) := N h (cid:88) j =1 c nj χ j ( x ) , the finite element approximations of the cell density u and the concentration of chemoattractant c where { u nj } j =1 ,...,N h and { c nj } j =1 ,...,N h are unknowns and { χ j } j =1 ,...,N h is the finite element basis.Then, the finite element problem associated with the system (16) reads as follows.13or each n = 0 , . . . , N T −
1, find { u n +1 h , c n +1 h } in S h × S h , such that for all χ ∈ S h (cid:18) u n +1 h − u nh ∆ t , χ (cid:19) h + (cid:0) D ( u nh ) ∇ u n +1 h , ∇ χ (cid:1) =( A ( φ upw ( u nh )) ∇ c nh , ∇ χ ) + r (cid:18) u nh (cid:18) − u nh u max (cid:19) , χ (cid:19) h , (24) ζ (cid:18) c n +1 h − c nh ∆ t , χ (cid:19) h = − (cid:0) ∇ c n +1 h , ∇ χ (cid:1) + (cid:0) u n +1 h − c n +1 h , χ (cid:1) h . (25)The function φ upw is defined by (76) in Appendix C. The finite element scheme associated with thesystem (17) including the effect of the treatment is the following θ (cid:18) w n +1 h − w nh ∆ t , χ (cid:19) h + (cid:0) D ( w nh , M nh ) ∇ w n +1 h , ∇ χ (cid:1) = (cid:16) B (cid:16) φ upw2 ( w nh , M nh ) (cid:17) ∇ c nh , ∇ χ (cid:17) + ˜ r (cid:18) w nh (cid:16) − w nh u max (cid:17) , χ (cid:19) h , (26) ζ (cid:18) c n +1 h − c nh ∆ t , χ (cid:19) h = − (cid:0) ∇ c n +1 h , ∇ χ (cid:1) + (cid:0) w n +1 h − c n +1 h , χ (cid:1) h , (27) m (cid:18) M n +1 h − M nh ∆ t , χ (cid:19) h = − (cid:0) ∇ M n +1 h , ∇ χ (cid:1) − δ (cid:0) M n +1 h , χ (cid:1) h . (28)The method to compute the coefficient φ upw and the diffusion coefficient D ( w nh , M nh ) is the same asfor the problem (24)–(25) (see Appendix C). Here, we comment on the choice of the model parameters that we will use for the numerical simu-lations and how they relate to experimental known data. In [19, 13], the proliferation rate for welloxygenated glioma cells in vitro r was shown to lie between 0 . − . As the proliferationrate relies significantly on the nutrient, also smaller value seems to be biologically admissible in realconditions and following [13] we choose r = 0 . − and r = 0 . − (corresponding to thenon-dimensionalised parameter r = 0 .
05 and r = 0 . d and consumption rate β for the chemoattractant are linked to biological measurementsof the oxygen diffusion in human brain which were estimated in [53, 19] to d = 86 . day − and β = 8640 day − . However, we chose to consider slightly lower values d = 8 .
64 mm day − and β = 864 day − owing to the large variability of this parameter according to the type of tissue (see[52]). For such values and using the scaling of Section 4.1, one unit of time of our model correspondsto approximately 0 . . χ u of glioma cells in response to chemoattractant concentration, the choiceof this parameter is driven by the stability analysis and we find that the interesting regimes areobtained for a dimensionless chemosensitivity in between 7 and 70, corresponding to a chemotac-tic coefficient χ u ∈ [0 . ,
6] mm day − . Moreover, as no measurements for glioma cells diffusioncoefficient are available in the literature, the parameter d is arbitrarily chosen to be 100 timessmaller than the chemoattratant diffusion speed and we choose d = d = 0 .
086 mm day − , i.e thenon-dimensionalised parameter ζ = d d = 0 .
01. 14 .3 Numerical results for a one dimensional case
For all numerical computations we choose the packing capacity ¯ u = 1. We consider different prolifera-tion rates r = 0 . , .
05 and different initial conditions and carrying capacities u max = u = 0 . , . ζ = m = 0 .
01 and θ = 1 since we assume that thechemoattractant c and the treatment diffuse much faster than the cells, while the motility of thecells is not affected by the treatment, so d ≈ d . The initial condition for the cell density u is as-sumed to be randomly distributed in space. Similarly, we can also define the initial chemoattractantconcentration c .In this section we start by solving the systems (16) and (17) on the interval [0 , L ] with homoge-neous non-flux boundary conditions using the method described in Section 5.1. In Appendix D weinvestigate the effect of the size of the domain as well as the effect of the parameters A and B onthe formation and evolution of patterns. Moreover, using (17) we study the effect of the treatmentusing the solution of (16) at the final time T as initial condition. We explore the case when weintroduce the treatment at earlier stages of the formation of the aggregates.We consider two different scenarios for the evolution of the concentration of the treatment. First,we assume that the treatment diffuses very fast in the whole domain so that the concentration ishomogeneous from time T = 0. The other case we consider, which is closer to real experiments,starts with a high concentration of the drug in the centre of the domain and this concentrationdiffuses over time according to the third equation in (17). Comparison with the linear stability analysis
In order to quantify the aggregate sizes andcompare it to the ones predicted by the stability analysis, we use the Fourier transform of thenumerical solution and extract the frequency that corresponds to the maximal Fourier mode. Forthe sake of this analysis, periodic boundary conditions are therefore considered. To this aim, wecompute the discrete Fourier transform F [ u ]( x , t ) = ˆ u ( λ, t ) and define k max = arg max λ ( | ˆ u ( λ ) | ) = arg max λ (cid:16)(cid:112) Re(ˆ u ( λ )) + Im(ˆ u ( λ )) (cid:17) , which corresponds to the frequency of the largest Fourier mode. The inverse ( k max ) − of thismaximal frequency relates to the pattern size. This maximal frequency of the Fourier transformof the solution is expected to correspond to the maximal wavenumber predicted by the stabilityanalysis. We show in Figure 5 the values of ( k max ) − computed from the numerical solution (bluedotted line) compared to the predictions of the stability analysis (red curve), as function of thechemosensitivity, for γ M = 1 (left figure) and γ M = 5 (right figure). As one can observe, we obtain avery good agreement between the numerical values and the predictions of the stability analysis, andwe recover the critical value of the chemosensitivity parameter above which the system generatespatterns. Introduction of the treatment on already-formed aggregates
In this part, we aim to studythe influence of the treatment on already formed aggregates. For this, we let the system run in P1 (without treatment, M = 0 , γ M = 1) until time T = 200 and introduce the treatment uniformly inthe domain ( M = 1, γ M = 5).In Figure 6 left and middle panels, we choose values of the chemosensitivity very close to thecritical values corresponding to k c , where the wavenumbers are very different for the cases γ M = 1and γ M = 5 as we see in Figure 4a. In Figure 6, the blue curves describe the formation of aggregatesfor a time T = 200 without the treatment, while the cells are proliferating with rate r = 0 .
1. We15
10 20 30 40 50 A / k m a x NumericalAnalytic B / k m a x NumericalAnalytic
Figure 5: Comparison of the wavelength obtained analytically, using (55) and (56), and numerically,using the Fourier transform of the solution for (a) γ M = 1 and (b) γ M = 5.consider two different scenarios when introducing the treatment: either cells stop proliferating (redcurves), or they continue with the same rate as before r = 0 . A and B close to the instability threshold ( A = B = 7, left figure), we observe that the aggregates become steeper and the density in each aggregatereaches the packing capacity ¯ u = 1. This clearly leads to more compact aggregates as a result ofthe nonlinearity introduced in (11) by the parameter γ M . The main physical difference betweenchanging the parameter γ M and changing the chemosensitivity coefficients A or B is the following.By changing A or B depending on the concentration of the treatment, we are enhancing aggregationover diffusion, essentially we are changing the motility of cells. By changing γ M the motility, as wellas the elastic properties of the cells in the aggregates are affected. When we introduce the treatmentwhile cells keep proliferating, aggregates tend to merge together since the density is growing, as wesee in Figure 6 middle panel.It is noteworthy that for large values of the chemosensitivity parameter ( A = B = 70, right panelin Figure 6), the treatment does not impact the aggregate dynamics. In this case, cell aggregationis mainly driven by the chemotactic term and the cell mechanical properties have little influence.These observations are in agreement with the stability analysis, which shows that the parameter γ M has more influence when the chemotactic sensitivity is close to the instability threshold. C e ll den s i t y Length of the domain M =1, r =0.1 M =5, r =0 M =5, r =0.1 0 20 4000.20.40.60.81 Figure 6: Aggregation pattern when A = B for A = 7 (left), A = 12 (middle) and A = 70 (right).Blue curve: solution at T = 200 without treatment, red curves: solution at T = T + 200 when thetreatment is introduced uniformly with γ M = 5 and proliferation is stopped, yellow curves: solutionwith the same parameters as the red curves when proliferation goes on with r = 0 . T = 200, we now aim to study the effect of the treatment introduced at different times in theaggregation process. Introduction of the treatment at different times
Here, we consider the case when the treat-ment is introduced at different times in the aggregation process. As before, we let the system runin P1 (without drug, M = 0 , γ M = 1) until time t = T and introduce the treatment uniformlyin the domain ( M = 1, γ M = 5). We consider the cases when the treatment has the ability tostop proliferation, and when the treatment only acts on the cell mechanical properties. In Figure 7,we show the results at time T = T + 200 (red curve), when the treatment is introduced at times T = 50 (left plots), T = 100 (middle plots) and T = 300 (right plots). Figures 7a show the resultswhen the treatment stops proliferation, Figures 7b when the treatment only changes the mechanicalproperties of the cells. For each, the blue curves are the density profiles before introducing the treat-ment. As one can observe, in Figures 7a and 7b, introducing the treatment at different times of theaggregation process have a major impact on the size of the aggregated patterns formed at a lattertime. Introducing the treatment at an earlier time ( T = 50, left figures) enables to obtain smalleraggregates compared to when the treatment is introduced on already formed aggregates ( T = 300,right figures). This effect is more visible when the treatment has the double effect of blocking cellproliferation and changing the elasticity (compare red curves in Figures 7a and 7b). In this case,the earlier the treatment is introduced, the smaller the aggregated patterns. When the treatmentstops proliferation as well as the cell mechanical properties and is introduced at later times (rightpanel of Figure 7a) we recover the observation of the real systems, where the treatment induces ashrinking of the aggregate and favors the formation of more compact cell structures. This effectis not observed when proliferation is active with the treatment, (right panel of Figure 7b) whereaggregates are merging and they are larger than before the treatment introduction. This suggestsindeed that the treatment has the double effect of blocking cell proliferation as well as changing thecell mechanical properties. The model suggests that introducing the treatment at earlier times of thetumour development could enable us to control the size and separation of the tumour aggregates. Introduction of the treatment in the middle of the domain
Finally, we aim to study thecase when the treatment is introduced in the center of the domain and diffuses in the environment.Here we assume that the treatment is not consumed or escapes the domain, therefore δ = 0 in(17). In Figure 8, we show the density profiles of the solution before introducing the treatment(blue curves), and when the treatment is present (red curves), at times T = 0 (left), T = 5(center) and T = 30 (right). The distribution of the treatment follows a Gaussian of the form M ( x ,
0) = Ce ( x − x σ , where C = 40 is the amplitude, x is the center of the Gaussian and σ = 0 .
20 4000.51 C e ll den s i t y T =50 Length of the domain T =100 T =300 M =1 M =5 (a) C e ll den s i t y T =50 Length of the domain T =100 T =300 M =1 M =5 (b) Figure 7: Introduction of the the treatment at different times T = 50 , ,
300 when u = 0 . A = B = 12. The blue curves gives the initial condition for the part P2 of the experiment with thetreatment, represented by the red curve. (a) Without proliferation. (b) With proliferation, r = 0 . T = T + 200. C e ll den s i t y T =0 Length of the domain T =5 T =30 =1= (M)TMZ Figure 8: Evolution of the pattern when the treatment is introduced in the center of the domainwith amplitude 40. The blue curves correspond to the solution without treatment at T = 200 andthe red curves are the solution with the treatment at T = T + 5 and T = T + 30.18 .4 Numerical results for a two dimensional case For the 2D simulations we consider that Ω is a disk of radius R which can be defined as Ω = { ( x, y ) ∈ R : x + y < R } where the boundary is given by ∂ Ω = { ( x, y ) ∈ R : x + y = R } .The proliferation rate is chosen to be r = 0 .
05 and the initial homogeneous density as well as thecarrying capacity are set to u max = u = 0 . u max = u = 0 .
5. The other parameters can befound at the beginning of Section 5.3.In Figure 9 we show the formation of aggregates for different values of A without the treatment,for u = 0 . r = 0 .
05. We observe that for A = 10 (Figure 9a) we do not have patterns, inagreement with the analytic results obtained in Figure 7b since this value of A is less than A c ≈ . u = u max = 0 .
5, the patterns change shape. We observe a transition from spot-like patterns inFigure 9 to maze-like structures in Figure 10. This behaviour has been widely studied experimentally[41], numerically [34, 36] and more recently, also including a volume-filling approach [43]. (a) (b) (c)
Figure 9: Formation of aggregates at T = 200 when u = 0 . r = 0 .
05 and (a) A = 10, (b) A = 20and (c) A = 70. (a) (b) (c) Figure 10: Formation of aggregates at T = 200 without the treatment and when u max = u = 0 . r = 0 .
05 and (a) A = 7, (b) A = 12 and (c) A = 50.Analogous to the one dimensional case, we consider two different initial conditions for the treat-ment: (i) we first include the treatment uniformly in the domain with M = 5, and (ii) we introducethe treatment in the middle of the domain and let it diffuses in space. In these simulations, thetreatment is supposed to block proliferation as well as changing the cell mechanical properties.In order to compare the change in size of the aggregates before and after the treatment, wecompute the difference between the solution of the first part of the experiments u ( x , T ) coming19rom (16), when the aggregates are formed (at T = 200), and the solution w ( x , t ) of (17), oncethe treatment has been inserted (at time T = 200 + T ). In Figure 11, we show the results fordifferent values of the chemosensitivity parameter A = B , when the aggregates have been exposedto a uniform distribution of the treatment for a time T = 30. We notice that the aggregates withoutthe treatment are bigger and therefore the treatment induces a shrinking in the size of the pattern, asdescribed in Section 3. When we compare the results in Figure 11a and 11b for different values of A and B we observe that the effect of the treatment is stronger when the value of the chemosensitivityparameter B is closer to its critical value (see Figure 4b). This is in accordance with the results ofthe stability analysis. The cell mechanical properties (controlled by the parameter γ M ) have lessinfluence on the cell cluster sizes when the chemosensitivity parameter A = B is increased (seeFigure 4b). (a) (b) Figure 11: Difference between the solutions u ( x , − w ( x , A = B = 20 and (b) A = B = 70. Here γ M = 11, r = 0 .
05 and u max = u = 0 . (a) (b) (c) Figure 12: Difference between the solutions when the initial concentration of the treatment is aGaussian function centered in the domain. (a) T = 210, (b) T = 230 (c) T = 300 for T = 200, r = 0 .
05 and u = 0 . δ = 0. InFigure 12 we show the evolution of the difference between the two solutions u ( x , T ) − w ( x , t ), where T = 200 is the time at which the treatment is introduced and T is the duration of the treatment.We explore different times T = 10 , , (a) (b) (c)(d) (e) (f) Figure 13: Effect of the treatment at different times while the aggregates are forming. Top row:Cell aggregates for A = 20 without treatment at (a) T = 30, (b) T = 50 and (c) T = 100. Bottomrow: Cell aggregates at time T = T + 200, when the treatment has been introduced at times (d) T = 30, (e) T = 50 and (f) T = 100.Finally, we also study the effect of the treatment at earlier stages of the formation of the aggre-gates. Figure 13 shows the different patterns obtained during the formation of the aggregates atdifferent times, top row, and the corresponding effect of the treatment, bottom row. For example,introducing the treatment at T = 50 leads to a significant reduction of the size of the pattern witha reasonably low concentration of the treatment. Identifying this specific time in real patients couldmake the treatment much more effective and reduce the spread of the cancer cells. In this paper we propose a mechanism for the effect of certain treatments on tumours formed by achemotaxis type self-organization phenomenon. Inspired by the experiments concerning the action ofTMZ on Glioblastoma cells mentioned in Section 2, we considered the particular case of treatmentsthat do not act as cytotoxic agents but rather induce stress in the environment, which may inducechanges in the mechanical properties of individual tumour cells by making them pass from rigidbodies to semi-elastic entities. We explored two scenarios: a first one where only cell’s plasticity isimpacted by the treatment, and a second one where the treatment has a double effect of preventingcell proliferation as well as changing cell mechanics.Under these hypotheses, we obtained a modified version of the Keller-Segel model, known as thenonlinear volume-filling approach for cell motion, where the cell mechanical properties are taken intoaccount in the form of the so-called squeezing probability. In the nonlinear volume-filling Keller-Segel21odel, this squeezing probability function could be related to the amplitude of the transport termtowards zones of high chemoattractant density (chemosensitivity), as well as with the (nonlinear)diffusion coefficient.By performing a linear stability analysis, we are able to detect and characterise the parameterranges for which the homogeneous distribution is unstable, i.e. the ranges for which patterns appearas the result of the dynamics. We show that the emergence of patterns without treatment ( i.e. when cells act as rigid bodies) is driven by a fine interplay between the chemotactic sensitivity,which tends to aggregate the cells, and the diffusion, together with the logistic growth, which tendto smoothen the solution. We are able to compute the critical chemosensitivity value above whichthe system self-organizes into aggregates, and characterise the size of the aggregates as a functionof the model parameters.Under treatment, i.e. when cells behave as semi-elastic entities, we show that the critical value ofthe chemosensitivity above which patterns emerge is smaller than that in control cultures, showingthat as cells become more elastic, they tend to aggregate more easily than when they behave as rigidentities.We are able to completely characterise the size of the patterns and show that semi-elastic cellscreate smaller aggregates than rigid entities for the same value of the chemosensitivity. These resultssuggest that the mechanical properties of individual cells have a huge impact on the shape and sizeof the aggregated patterns at the population level.Moreover, we show that the ratio between the tight packing cell density and the carrying capacityof the TME plays a major role in the size and shape of the obtained patterns. For large values ofthis ratio, the aggregation abilities of the system are essentially driven by the chemotactic transportterm while the individual cell mechanical behaviour has little impact on the shape and size of thepatterns. However, for smaller values of this ratio, i.e. when the tight packing density is closer tothe carrying capacity of the environment, cell mechanics has a huge influence on the behavior of thepopulation.These results are confirmed by numerical experiments in 1D and 2D for which, given an initialperturbation of the homogeneous cell distribution, we observe the emergence of cell aggregatesand we recover the critical values of the chemosensitivity predicted by the stability analysis. Weobtain a very good correspondence between the simulations and the theoretical predictions, for theappearance of patterns as well as their size.By performing simulations of the whole system, we recover the experimental observations: in-troducing the treatment (TMZ in the experiments mentioned in Section 2) on already-formed ag-gregates leads to the quick formation, of more compact patterns. As the treatment diffuses in thedomain (changing locally the cell mechanical properties as it goes), it sharpens the border of thecell aggregates and leads to denser and well-separated clusters.While the border sharpening of the clusters is independent on whether proliferation is activatedor not during treatment, the shrinking of the clusters is more clear when the treatment has thedouble effect of changing the cell mechanical properties as well as blocking cell proliferation. Indeed,when proliferation is still active in the presence of the treatment, we observe the merging of existingclusters and this results in aggregates being larger than before treatment. These results suggest apossible mechanism for the shrinking of the aggregates observed under the experimental conditionsdescribed in Section 2: TMZ might not only stop cell proliferation, but might also generate a stressin the environment to which cells respond to by changing their mechanical state.While alterations of mechanical properties, around or inside the tumour, are common in solidtumours including GBM, the question of the nature and the regulation of cancer cells throughmechano-sensitive pathways are largely unknown. Recently, in a Drosophila model, glioma pro-gression has been associated to a regulatory loop mediated by the mechano-sensitive ion channel22iezo1 and tissue stiffness [12]. A direct perspective of these works consists in verifying the poten-tial mechano-sensing effect of TMZ proposed in the present model, through direct measures in realsystems by studying the mechanical properties of individual GBM cells which have been exposed toTMZ treatment. The targeting of mechano-sensitive pathways after TMZ treatment may providenew therapeutic angles in GBM and in more general settings.Moreover, it would be interesting to identify other clinical settings where the effects of thetreatment are similar to those of TMZ in GBM, and to check if the effects are due to changes in thetumor cell properties corresponding to the general hypothesis of the model constructed in this work.In the future, better quantitative comparison with experiments will allow for systematic choiceof parameters and validation of the mechanisms we propose here. From a biological point of view, anatural sequel of this work consists in studying the coupled effect of TMZ and irradiation. Indeed,even if TMZ alone seems not to suffice to decrease the tumour mass, the coupling of TMZ treatmentwith irradiation has been shown to have more efficient effects than irradiation alone [48, 5, 54]. Itwould also be interesting to introduce a second treatment with cytotoxic effects in this model, tostudy whether the mechanical changes of individual cells induced by TMZ could explain the betterresponse of the system to irradiation treatments.23
Derivation of the general model
The modified Keller-Segel system (7), including the squeezing probability (3) was derived in onedimension in [43, 56] from a continuous-time and discrete-space random walk. In this appendix wegive a more general (and formal) derivation in R n starting from a kinetic equation, analogous to[40, 45]. For simplicity, we do not include cell proliferation in this derivation.Let us consider a mesoscopic density h ( x , t, v ), where x ∈ R n and v ∈ V = { x ∈ R n : | x | = 1 } ,which evolution is described by the following kinetic equation ∂ t h + v · ∇ h = − ξ ( c ) h ( x , t, v ) + q ( σ ( x , t )) (cid:90) V T ( v , v (cid:48) , c ) h ( x , t, v (cid:48) ) d v (cid:48) . (29)Equation (29) describes a so-called velocity jump model where, individuals moving in an almoststraight line, described by the left hand side, switch velocities after stopping. The right hand sideof (29) describes the density of cells that are stopping with a frequency ξ ( c ) = (cid:82) V T ( v (cid:48) , v , c ) d v (cid:48) ,where the operator T ( v , v (cid:48) , c ) gives the probability of a velocity jump from v (cid:48) to v . Note that thecells change velocity also depending on the chemoattractant concentration c ( x , t ).The second term in the right hand side represents the individuals that start a new trajectorywith a different velocity v (cid:48) . Note that this term is multiplied by the squeezing probability q sincethe change into a new velocity is also determined by the probability of finding a neighbouring spaceavailable.Now we introduce a parabolic scaling ( x , t ) (cid:55)→ (¯ x /ε, ¯ t/ε ) into (29), where the bar denotes thenew variables and ε (cid:28)
1. Letting T ( v , v (cid:48) , c ) be a small perturbation of a random turning event, T ( v , v (cid:48) , c ) = T + εT T ( v (cid:48) , ∇ c ), we have, after dropping the bars, ε ∂ t h ε + ε v · ∇ h ε = q | V | T u ε ( x , t ) + εqT (cid:90) V T ( v (cid:48) , ∇ c ) h ε ( x , t, v (cid:48) ) d v (cid:48) − h ε ( x , t, v ) T | V | − εT h ε ( x , t, v ) | V | T ( v , ∇ c ) , (30)where u ε ( x , t ) = 1 | V | (cid:90) V h ε ( x , t, v ) d v . (31)Dividing by ε we get ε∂ t h ε + v · ∇ h ε = T ε (cid:16) qu ε ( x , t ) | V | − h ε ( x , t, v ) | V | (cid:17) + T (cid:16) q (cid:90) V T ( v (cid:48) , ∇ c ) h ε ( x , t, v (cid:48) ) d v (cid:48) − h ε ( x , t, v ) | V | T ( v , ∇ c ) (cid:17) . (32)The leading order terms in (32), when ε →
0, give h ( x , t, v ) = qu ( x , t ). Integrating with respect to v ∈ V in (30) we obtain a macroscopic conservation equation ∂ t u + ∇ · j = 0 , (33)where j ε = ε | V | (cid:82) V v h ε ( x , t, v ) d v is the mean direction of the cells.Finally, we have to obtain j ε in terms of the macroscopic density u ε . For that, we multiply (30)by v and integrate in V to get ε | V | ∂ t j ε + ∇ · (cid:90) V v ⊗ v h ε d v = − T | V | j ε − | V | T (cid:90) V v T ( v , ∇ c ) h ε d v . (34)24etting h ε ( x , t, v ) = qu + εh ⊥ ( x , t, v ) + O ( ε ) we obtain, in the limit as ε → j = − (cid:82) V v ⊗ v d v T | V | ∇ ( qu ) − (cid:82) V v T ( v , ∇ c ) d v | V | uq . (35)Substituting (35) into the conservation equation (33) we obtain ∂ t u − ∇ · (cid:16) d ∇ ( qu ) − χ ( c ) uq (cid:17) = 0 , (36)where d = (cid:82) V v ⊗ v d v T | V | and χ ( c ) = − (cid:82) V v T ( v , ∇ c ) d v | V | . (37)Considering the squeezing probability function q ( u ) in (3) and using the chain rule for differen-tiation, ∇ ( qu ) = ( q − q (cid:48) u ) ∇ u where q (cid:48) = d q d u we finally have, assuming χ ( c ) = χ u ∇ c , ∂ t u = ∇ · ( d ( q − q (cid:48) u ) ∇ u − χ u uq ∇ c )= ∇ · ( d D ( u ) ∇ u − χ u φ ( u ) ∇ c ) . (38) Remark A.1.
With the choice (3) for the squeezing probability, the diffusion of the cells is enhanced (cid:16) d D ( u )d u > (cid:17) . This means that the elastic collisions of the cells may increase, locally, the randommotion component. B Stability analysis
First part of the experiments
We first observe that the homogeneous distributions u ( x , t ) = u ∗ and c ( x , t ) = c ∗ are steady-states solutions of system (16) for u ∗ and c ∗ such that f ( u ∗ ) = 0 and g ( u ∗ , c ∗ ) = 0. In order to study their stability, we consider the system without spatial variations ∂ t u = f ( u ) , ζ∂ t c = g ( u, c ) , (39)and linearize the solution at ( u ∗ , c ∗ ). We obtain ∂ t σ = Gσ , where σ = (cid:18) u − u ∗ c − c ∗ (cid:19) and G = (cid:18) f ∗ u ζ − g ∗ u ζ − g ∗ c (cid:19) , (40)where the quantities f ∗ u , g ∗ u and g ∗ c are the linearization slopes of f and g : f ∗ u = f (cid:48) ( u ∗ ) , g ∗ u = ∂ u g ( u ∗ , c ∗ ) , g ∗ c = ∂ c g ( u ∗ , c ∗ ). The steady state is linearly stable if tr( G ) < G ) >
0, whichimposes the following constraints on the kinetic functions f ( u ) and g ( u, c ), f ∗ u + ζ − g ∗ c < f ∗ u g ∗ c > . (41)Note that in our case, f ∗ u = − r , g ∗ u = 1 , g ∗ c = − i.e. the ability of the system to create patterns, we introduce a smallparameter ε (cid:28) u = u ∗ + ε ˜ u ( x , t ) , c = c ∗ + ε ˜ c ( x , t ) . (42)25e substitute (42) into (19) and, computing the first order terms with respect to ε and neglectinghigher order terms, the linearized system reads ∂ t u = ∆ u − Aφ ( u ∗ )∆ c + uf ∗ u + cf ∗ c ,ζ∂ t c = ∆ c + ug ∗ u + cg ∗ c , (43)where φ ( u ∗ ) = u ∗ q ( u ∗ ). We now look for perturbations of the form u ( x , t ) = (cid:88) k a k ( t ) ψ k ( x ) and c ( x , t ) = (cid:88) k b k ( t ) ψ k ( x ) , (44)where ( ψ k ) k ≥ is an orthonormal basis of L (Ω) and satisfies the following spatial eigenvalue problem − ∆ ψ k = k ψ k , ∂ψ k ∂η = 0 . (45)Then, the linearized system (43) can be written as ∂ t ( a k ) = − a k k + Aφ ( u ∗ ) b k k + a k f ∗ u + b k f ∗ c ,ζ∂ t ( b k ) = − b k k + a k g ∗ u + b k g ∗ c , (46)where k is the spatial eigenfunction, also called the wavenumber and 1 /k is proportional to thewavelength ω . In matrix form we can write (46) as ∂ t X k ( t ) = P k ( t ) X k ( t ) where X k = (cid:18) a k b k (cid:19) , P k = (cid:18) − k + f ∗ u Aφ ( u ∗ ) k + f ∗ c ζ − g ∗ u ζ − ( − k + g ∗ c ) (cid:19) . (47) Remark B.1.
Since the solutions of the eigenvalue problem (45) are simply sines and cosines, the“size” of various spatial patterns is measured by the wavelength of the trigonometric functions. Forexample, in one dimension when < x < L , ψ ∝ cos( nπx/L ) and the wavelength is ω = 1 /k = L/nπ ,where n ∈ Z . If the matrix P k has eigenvalues with positive real part, then the homogeneous steady state( u ∗ , c ∗ ) is unstable, resulting in pattern formation. The characteristic polynomial related to (47) isgiven by (cid:96) + a ( k ) (cid:96) + b ( k ) = 0 where a ( k ) = (1 + ζ − ) k − ( f ∗ u + ζ − g ∗ c ) , (48) b ( k ) = ζ − k − ζ − ( g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ )) k + ζ − f ∗ u g ∗ c . (49)The eigenvalues (cid:96) determine the temporal growth of the eigenmodes, and we require R e ( (cid:96) ( k )) > k (cid:54) = 0since we already guaranteed that the steady state is stable in the absence of spatial perturbations, i.e. R e ( (cid:96) ( k = 0)) < a ( k ) >
0, hence the instability can only occur if b ( k ) < k so that the characteristic polynomial has one positive and one negative root.This implies k − ( g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ )) k + f ∗ u g ∗ c < . (50)We also know from (41) that f ∗ u g ∗ c >
0, then a necessary (but not sufficient) condition for b ( k ) < g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ ) > . emark B.2. The bifurcation between spatially stable and unstable modes occurs when the criticalexpression b min ( k ) = 0 is satisfied. Moreover, to satisfy (50) the minimum b min must be negative [35]. Differentiation with respectto k in (49) leads to b min ( k ) = − ( g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ )) f ∗ u g ∗ c . (51)Hence, the condition b min < g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ ) > (cid:112) ( f ∗ u g ∗ c ) . (52)To summarise, we have obtained the following conditions for the generation of spatial patternsfor the chemotaxis system (19), f ∗ u + ζ − g ∗ c < , f ∗ u g ∗ c > , ζ − g ∗ c + f ∗ u − ζ − g ∗ u Aφ ( u ∗ ) > ,g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ ) > (cid:112) f ∗ u g ∗ c . (53)From the analysis in this section, and using the particular forms of φ ( u ), f and g as in (19),it is easy to see that the spatially homogeneous steady states are (0 ,
0) and ( u max , u max ). Wecan check that (0 ,
0) is an unstable steady state, therefore we only work with ( u max , u max ) which,on the contrary, is stable. The first and second properties in (53) are immediately satisfied, i.e. , − ( r + ζ − ) < r /ζ >
0, respectively. Finally, we have to check that the third and fourthconditions are satisfied as well. We have that − − r + Aφ ( u ∗ ) > √ r . (54)Therefore, (54) is a necessary condition for pattern formation for the original system (12). Consid-ering A as a bifurcation parameter, we can obtain a critical value A c , so that we observe patternformation if A > A c . From (54) we get A c = 2 √ r + 1 + r u max (cid:0) − u max ¯ u (cid:1) . (55)The corresponding critical wavenumber k c is obtained from (51) using (54) as follows, k c = g ∗ c + f ∗ u + g ∗ u A c φ ( u ∗ )2 = (cid:112) f ∗ u g ∗ c = √ r . (56)This means that, within the unstable range, R e ( (cid:96) ( k )) > k c . The range of linearly unstable modes k < k < k is obtained from b ( k ) = 0, k = − m − (cid:112) m − f ∗ u g ∗ c < k < k = − m + (cid:112) m − f ∗ u g ∗ c , (57)where m = − ( g ∗ c + f ∗ u + g ∗ u Aφ ( u ∗ )). Second part of the experiments
Following the same steps as before we linearise the system(17) to get ∂ t w = D ( w ∗ , M ∗ )∆ w − Bφ ( w ∗ , M ∗ )∆ c + wf ∗ w ,ζ∂ t c = ∆ c + wg ∗ w + cg ∗ c ,m∂ t M = ∆ M − δ w , (58)27here D ( w ∗ , M ∗ ) = 1 + ( γ ∗ M − (cid:16) w ∗ ¯ u (cid:17) γ ∗ M , φ ( w ∗ , M ∗ ) = w (cid:16) − (cid:16) w ∗ ¯ u (cid:17) γ ∗ M (cid:17) , and γ ∗ M is given by (5) evaluated at M ∗ . As in (44), we let w ( x , t ) = (cid:88) k a k ( t ) ψ k ( x ) , c ( x , t ) = (cid:88) k b k ( t ) ψ k ( x ) , M ( x , t ) = (cid:88) k c k ( t ) ψ k ( x ) , (59)where ψ k ( x ) satisfies (45) and we obtain a system ∂ t X k ( t ) = P k ( t ) X k ( t ) where X k = a k b k c k , P k = − D ( w ∗ , M ∗ ) k + f ∗ w Bφ ( w ∗ , M ∗ ) k ζ g ∗ w ζ ( − k + g ∗ c ) 0 − δ m − m k . (60)Similar to the previous section, the characteristic polynomial is given by a ( k ) (cid:96) + b ( k ) (cid:96) + c ( k ) (cid:96) + d ( k ) = 0 where a ( k ) = − b ( k ) = − (cid:16) D + 1 ζ + 1 m (cid:17) k + f ∗ w + g ∗ c ζ , (61) c ( k ) = − (cid:16) D ζ + D m + 1 mζ (cid:17) k + (cid:16) g ∗ c D ζ + g ∗ c ζm + f ∗ w ζ + f ∗ w m + g ∗ w Bφ ζ (cid:17) k − f ∗ w g ∗ c ζ , (62) d ( k ) = − D mζ k + 1 mζ ( D g ∗ c + g ∗ w Bφ + f ∗ w ) k + 1 mζ ( − f ∗ w g ∗ c ) k . (63)In general, the stability analysis for this cubic polynomial will require the Ruth–Hurwitz stabilitycriterion [25] which states that the steady state is unstable if the coefficients of a ( k ) (cid:96) + b ( k ) (cid:96) + c ( k ) (cid:96) + d ( k ) = 0 satisfy the condition1( a ( k )) (cid:0) b ( k ) c ( k ) − a ( k ) d ( k ) (cid:1) < . However, from (60) we observe that one of the eigenvalues of the matrix P k is given by (cid:96) = − k m < (cid:18) − D ( w ∗ , M ∗ ) k + f ∗ w Bφ ( w ∗ , M ∗ ) k ζ g ∗ w ζ ( − k + g ∗ c ) (cid:19) , (64)following the same analysis as for the case without TMZ.The characteristic polynomial (cid:96) + ¯ a ( k ) (cid:96) + ¯ b ( k ) = 0 related to (64) has coefficients¯ a ( k ) = (cid:16) D ( w ∗ , M ∗ ) + 1 ζ (cid:17) k − f ∗ w − g ∗ c ζ , (65)¯ b ( k ) = D ( w ∗ , M ∗ ) ζ k − (cid:16) D ( w ∗ , M ∗ ) g ∗ c ζ + Bφ ( w ∗ , M ∗ ) g ∗ n ζ + f ∗ w ζ (cid:17) k + f ∗ w g ∗ c ζ . (66)28or the steady state to be unstable we require, as before, that R e ( (cid:96) ( k )) >
0. Since ¯ a ( k ) > b ( k ) <
0. Computing d¯ b ( k )d k = 0 from (66) we obtain k = D ( w ∗ , M ∗ ) g ∗ c + g ∗ w Bφ ( w ∗ , M ∗ ) + f ∗ w D ( w ∗ , M ∗ ) . (67)Hence from the condition ¯ b min ( k ) < D ( w ∗ , M ∗ ) g ∗ c + g ∗ w Bφ ( w ∗ , M ∗ ) + f ∗ w > (cid:112) D ( w ∗ , M ∗ ) f ∗ w g ∗ c . (68)The spatially homogeneous steady state is ( w ∗ , c ∗ , M ∗ ) = ( u max , u max , M s ), where M s = | Ω | − (cid:82) Ω M ( x ,
0) d x . Therefore, from (68) we obtain a critical constant B c so that for any B > B c we observe pattern formation. This critical constant is given by B c = 2 (cid:112) ¯ r D ( u max , M s ) + D ( u max , M s ) + ˜ r u max (cid:16) − (cid:16) u max ¯ u (cid:17) γ Ms (cid:17) , (69)where D ( u max , M s ) = 1 + ( γ M s − (cid:16) u max ¯ u (cid:17) γ Ms . (70)The corresponding critical wavemode is given by k c = D ( u max , M s ) g ∗ c + g ∗ w B c φ ( u max , M s ) + f ∗ w D ( u max , M s ) = (cid:112) D ( u max , M s )( f ∗ w g ∗ c ) D ( u max , M s ) . (71)Finally, the unstable modes are k < k c , where from ¯ b ( k ) = 0 we get k = − ¯ m − (cid:112) ¯ m − D ( u max , M s )( f ∗ w g ∗ c )2 D ( u max , M s ) < k < k = − ¯ m + (cid:112) ¯ m − D ( u max , M s )( f ∗ w g ∗ c )2 D ( u max , M s ) , (72)for ¯ m = − ( D ( u max , M s ) g ∗ c + g ∗ w Bφ ( u max , M s ) + f ∗ w ). C Description of the numerics
We denote H (Ω) = W , (Ω) the usual Sobolev space and the standard L inner product is denotedby ( · , · ). Let T h , h >
0, be a quasi-uniform mesh of the domain Ω consisting of N disjoint piecewiselinear mesh elements K such that the discretized domain Ω h = (cid:83) K ∈T h K . Let h K := diam( K ) and h = max K h K and for d = 2, we choose linear triangular elements. In addition, we assume that themesh is acute i.e. for d = 2, each angle of the triangles can not exceed π . We must stress that for d = 2 since the domain Ω is circular, a small error of approximation is committed using Ω h . Weconsider the standard finite element space associated with T h S h := { χ ∈ C (Ω) : χ (cid:12)(cid:12) K ∈ P ( K ) , ∀ K ∈ T h } ⊂ H (Ω) , (73)where P ( K ) denotes the space of first-order polynomials on K . Let N h be the total number of nodesof T h , J h the set of nodes and { x j } j =1 ,...,N h their coordinates. We call { χ j } j =1 ,...,N h the standard29agrangian basis functions associated with the spatial mesh. We denote by π h : C (Ω) → S h thestandard Lagrangian interpolation operator We also need the lumped scalar product to define theproblem ( η , η ) h = (cid:90) Ω π h ( η ( x ) η ( x )) d x ≡ (cid:88) x j ∈ J h (1 , χ j ) η ( x j ) η ( x j ) , η , η ∈ C (Ω) . We define the standard mass and stiffness finite element matrices as G and K , where G ij = (cid:90) Ω χ i χ j d x, for i, j = 1 , . . . , N h ,K ij = (cid:90) Ω ∇ χ i ∇ χ j d x, for i, j = 1 , . . . , N h . In the following finite element approximation of the Keller-Segel problem, the mass matrix is lumped, i.e. the matrix becomes diagonal with each term being the row-sum of the corresponding row of thestandard mass matrix, G l,ii := N h (cid:88) j =1 G ij , for i = 1 , . . . , N h . In order to describe how the chemotactic coefficient φ uwp is computed, let us rewrite the discreteequation (24) into its matrix form( G L + ∆ tK D ) u n +1 = G L u n + ∆ tK φ c n + ∆ tG L g n , where u n and c n are the vectors of coefficients which are the unknowns of the problem and g n is avector defined by (cid:2) g n (cid:3) i = (cid:18) u nh (cid:18) − u nh u max (cid:19)(cid:19) ( x i ) , for i = 1 , . . . , N h . We define the finite element matrices associated with the diffusion K D and the advection K φ K D,ij = (cid:90) Ω D ( u nh ) ∇ χ i ∇ χ j d x for i, j = 1 , . . . , N h , (74) K φ,ij = (cid:90) Ω φ upw ( u nh ( x i ) , u nh ( x j )) ∇ χ i ∇ χ j d x for i, j = 1 , . . . , N h . (75)In (74), the integral is computed using Gauss quadrature to deal with a potential choice of nonlinearfunctional for D ( u nh ). The exactness of the quadrature is obtained using the adequate number ofGauss points since D ( u nh ) is a polynomial of order γ M + 1.The chemotactic coefficient is computed using an upwing approach. For each element and de-pending on the direction of the gradient of the chemoattractant we have φ upw ( u nh ( x i ) , u nh ( x j )) = u nh ( x j ) (cid:16) − (cid:16) u nh ( x i ) u (cid:17)(cid:17) , if c nh ( x j ) − c nh ( x i ) < ,u nh ( x i ) (cid:16) − (cid:16) u nh ( x j ) u (cid:17)(cid:17) , otherwise . (76)Therefore, the chemotactic coefficient is chosen as function of the sign of the difference of chemoat-tractant between nodes connected by an edge. The property of non-negativity of the cell densitysatisfied by our numerical scheme can be proved using similar arguments as in [10].30 One dimensional numerical results
Influence of the domain size
The unstable wavenumbers are discrete values, k = nπ/L , thatsatisfy the relation (57) from Section 4.2. The wavemode n determines the number of aggregatesdepending on the length of the domain. For A = 7 and the parameters specified at the beginningof this section we have 0 . < nπ/L < .
4. As shown in Figure 14, as we increase the length of thedomain, the number of aggregates also increases. When the domain is large, as in Figure 14c, weobserve that some aggregates are merging together while others are emerging, i.e. , they are formedfrom a zone of low cell density. This process is called coarsening [43] and is not observed in a smalldomain such as in Figure 14a.
Length of the domain (L) C e ll den s i t y (a) Length of the domain (L) C e ll den s i t y (b) Length of the domain (L) C e ll den s i t y (c) Figure 14: Relationship between the wavenumber k and the length of the domain L for A = 7.For the remaining simulations in this section we fix the length of the domain to L = 40, and allsimulations are performed with carrying capacity u max = 0 . Comparison of the results with the stability analysis predictions
Here, we study theinfluence of the chemosensitivity parameter A on the pattern dynamics and size of the aggregates,in the presence or absence of TMZ. In order to compare the solutions to the predictions of the stabilityanalysis, the initial condition is a small perturbation around the homogeneous distribution u = 0 . r = 0 .
1. For such parameters, using theresults of Section 4, the critical value of the chemosensitivity parameter without the treatment (in P1 , when M = 0 and therefore γ M = 1 ) is A c ≈ .
92 and with the treatment uniformly distributed(for γ M = 5) is B c ≈ . i.e. without any treatment, we show in Figure 15 theformation of patterns at different times, for A = 7 (close to the instability threshold, Figure 15a), A = 50 (Figure 15b) and A = 150 (Figure 15c). We observe here the process of merging andemerging patterns through time.As predicted by the stability analysis, larger values of the chemotactic sensitivity A favor theemergence of smaller aggregates (compare Figures 15a and 15c for A = 7 and A = 150, respectively).This is due to the fact that for larger chemotactic sensitivity, cells are more attracted to zones ofhigh concentration of chemoattractant. The creation of patterns instead of the expansion of ahomogeneous cell distribution is due to an instability which results from a positive feedback loopbetween the production of the chemical by the cells on one hand, and their attraction to high densityzones of this chemical on the other. The chemotactic sensitivity A must be large enough to triggerthis instability, in order to compensate the competing effects of diffusion and of the logistic growthterm, which on the contrary, tends to regulate the local cell density to the carrying capacity of theenvironment u max , and therefore induces cell death inside the aggregated patterns for which the31
20 4000.20.40.60.81 C e ll den s i t y T=100
Length of the domain
T=300
T=500 (a) C e ll den s i t y T=100
Length of the domain
T=300
T=500 (b) C e ll den s i t y T=100
Length of the domain
T=300
T=500 (c)
Figure 15: Formation of the pattern without the treatment for (a) A = 7, (b) A = 50 and (c) A = 150 when r = 0 . u max = u = 0 .
5. 32ensity is above u max .When the drug is introduced uniformly in the domain starting from a homogeneous distributionof cells ( i.e. for M = 1, γ M = 5) at time t = 0, we show in Figure 16 the formation of patterns atdifferent times, for chemosensitivity B = 5 (close to the instability threshold, Figure 16a), B = 30(Figure 16b) and B = 150 (Figure 16c). Note that in this case we also let cells to proliferate withrate r = 0 . C e ll den s i t y T=100
Length of the domain
T=300
T=500 (a) C e ll den s i t y T=100
Length of the domain
T=300
T=500 (b) C e ll den s i t y T=50
Length of the domain
T=500
T=1000 (c)
Figure 16: Formation of the pattern with treatment ( M = 1 , γ M = 5), included at t = 0 when u = 0 . r = 0 . B = 5, (b) B = 30 and (c) B = 150.As one can see in Figure 16, we first observe again that increasing the chemosensitivity parameter B results in the formation of smaller cell aggregates (compare Figures 16a and 16b). Very close tothe instability threshold (Figure 16a), the system converges quickly to one aggregate, while for largervalues of B (Figure 16c) a large number of well-separated small aggregates arises. These clustersmerge in time to form bigger clusters as for the case without the treatment. Moreover, comparingFigures. 15c and 16c, we clearly observe that varying the mechanical state of cells ( i.e. passing from γ M = 1 to γ M = 5), leads to a change in the cell’s aggregate size. When cells are more elastic, theytend to create smaller aggregates than when they behave as rigid spheres.33 cknowledgments The authors would like to thank K. J. Painter and H. Gimperlein for fruitful discussions and helpfulsuggestions. This research was funded by the Plan Cancer/ITMO HTE (tumor heterogeneity andecosystem) program (MoGlimaging).
References [1] A. Agosti, S. Marchesi, G. Scita, and P. Ciarletta. Modelling cancer cell budding in-vitro asa self-organised, non-equilibrium growth process.
Journal of Theoretical Biology , 492:110203,2020.[2] J. Alfonso, K. Talkenberger, M. Seifert, B. Klink, A. Hawkins-Daarud, K. Swanson,H. Hatzikirou, and A. Deutsch. The biology and mathematical modelling of glioma invasion: areview.
Journal of the Royal Society Interface , 14(136):20170490, 2017.[3] L. Almeida, F. Bubba, B. Perthame, and C. Pouchol. Energy and implicit discretization of theFokker-Planck and Keller-Segel type equations.
Network and Heterogeneous Media , 14(1):23–41,2019.[4] W. Alt. Biased random walk models for chemotaxis and related diffusion approximations.
Journal of mathematical biology , 9(2):147–177, 1980.[5] L. Barazzuol, R. Jena, N. G. Burnet, J. C. Jeynes, M. J. Merchant, K. J. Kirkby, and N. F.Kirkby. In vitro evaluation of combined temozolomide and radiotherapy using x rays and high-linear energy transfer radiation for glioblastoma.
Radiation research , 177(5):651–662, 2012.[6] E. Ben-Jacob, I. Cohen, and H. Levine. Cooperative self-organization of microorganisms.
Ad-vances in Physics , 49(4):395–554, 2000.[7] M. Brandao, T. Simon, G. Critchley, and G. Giamas. Astrocytes, the rising stars of the glioblas-toma microenvironment.
Glia , 67(5):779–790, 2019.[8] F. Bubba, T. Lorenzi, and F. R. Macfarlane. From a discrete model of chemotaxis withvolume-filling to a generalized Patlak–Keller–Segel model.
Proceedings of the Royal SocietyA , 476(2237):20190871, 2020.[9] F. Bubba, C. Pouchol, N. Ferrand, G. Vidal, L. Almeida, B. Perthame, and M. Sabbah. Achemotaxis-based explanation of spheroid formation in 3d cultures of breast cancer cells.
Journalof Theoretical Biology , 479:73–80, 2019.[10] F. Bubba and A. Poulain. A positivity-preserving finite element scheme for the relaxed cahn-hilliard equation with single-well potential and degenerate mobility.
ArXiv , page 20, 2019.[11] M. A. Chaplain and G. Lolas. Mathematical modelling of cancer invasion of tissue: dynamicheterogeneity.
Networks & Heterogeneous Media , 1(3):399, 2006.[12] X. Chen, S. Wanggou, A. Bodalia, M. Zhu, W. Dong, J. J. Fan, W. C. Yin, H.-K. Min, M. Hu,D. Draghici, et al. A feedforward mechanism mediated by mechanosensitive ion channel piezo1and tissue mechanics promotes glioma aggression.
Neuron , 100(4):799–815, 2018.3413] M. C. Colombo, C. Giverso, E. Faggiano, C. Boffano, F. Acerbi, and P. Ciarletta. Towards thepersonalized treatment of glioblastoma: integrating patient-specific clinical data in a continuousmechanical model.
PLoS One , 10(7):e0132887, 2015.[14] J. Condeelis, R. H. Singer, and J. E. Segall. The great escape: when cancer cells hijack the genesfor chemotaxis and motility.
Annual Review of Cell and Developmental Biology , 21:695–718,2005.[15] S. L. Di Jia, D. Li, H. Xue, D. Yang, and Y. Liu. Mining tcga database for genes of prognosticvalue in glioblastoma microenvironment.
Aging (Albany NY) , 10(4):592, 2018.[16] Y. Dolak and C. Schmeiser. The Keller–Segel model with logistic sensitivity function and smalldiffusivity.
SIAM Journal on Applied Mathematics , 66(1):286–308, 2005.[17] D. Dormann and C. J. Weijer. Chemotactic cell movement during dictyostelium developmentand gastrulation.
Current Opinion in Genetics & Development , 16(4):367–373, 2006.[18] Y. A. Fouad and C. Aanei. Revisiting the hallmarks of cancer.
American Journal of CancerResearch , 7(5):1016, 2017.[19] H. B. Frieboes, J. S. Lowengrub, S. Wise, X. Zheng, P. Macklin, E. L. Bearer, and V. Cristini.Computer simulation of glioma growth and morphology.
Neuroimage , 37:S59–S70, 2007.[20] P. Friedl. Prespecification and plasticity: shifting mechanisms of cell migration.
Current opinionin cell biology , 16(1):14–23, 2004.[21] Y. Han, Z. Li, J. Tao, and M. Ma. Pattern formation for a volume-filling chemotaxis modelwith logistic growth.
Journal of Mathematical Analysis and Applications , 448(2):885–907, 2017.[22] D. Hanahan and R. A. Weinberg. Hallmarks of cancer: the next generation.
Cell , 144(5):646–674, 2011.[23] T. Hillen and K. J. Painter. A user’s guide to pde models for chemotaxis.
Journal of Mathe-matical Biology , 58(1-2):183, 2009.[24] D. Horstmann. The Keller-Segel model in chemotaxis and its consequences.
I, Jahresber.Deutsch. Math.-Verein , pages 103–165, 2004.[25] A. Hurwitz et al. On the conditions under which an equation has only roots with negative realparts.
Selected Papers on Mathematical Trends in Control Theory , 65:273–284, 1964.[26] M. Ibrahim and M. Saad. On the efficacy of a control volume finite element method for thecapture of patterns for a volume-filling chemotaxis model.
Computers & Mathematics withApplications , 68(9):1032–1051, 2014.[27] E. F. Keller and L. A. Segel. Model for chemotaxis.
Journal of Theoretical Biology , 30(2):225–234, 1971.[28] J. S. Kennedy and D. Marsh. Pheromone-regulated anemotaxis in flying moths.
Science ,184(4140):999–1001, 1974.[29] Y. Kim and S. Kumar. Cd44-mediated adhesion to hyaluronic acid contributes to mechanosens-ing and invasive motility.
Molecular Cancer Research , 12(10):1416–1429, 2014.3530] S. Y. Lee. Temozolomide resistance in glioblastoma multiforme.
Genes & Diseases , 3(3):198–210, 2016.[31] D. N. Louis, A. Perry, G. Reifenberger, A. Von Deimling, D. Figarella-Branger, W. K. Cavenee,H. Ohgaki, O. D. Wiestler, P. Kleihues, and D. W. Ellison. The 2016 world health organizationclassification of tumors of the central nervous system: a summary.
Acta Neuropathologica ,131(6):803–820, 2016.[32] M. Ma, C. Ou, and Z.-A. Wang. Stationary solutions of a volume-filling chemotaxis model withlogistic growth and their stability.
SIAM Journal on Applied Mathematics , 72(3):740–766, 2012.[33] G. D. Maurer, D. P. Brucker, and J. P. Steinbach. Loss of cell-matrix contact increases hypoxia-inducible factor-dependent transcriptional activity in glioma cells.
Biochemical and biophysicalresearch communications , 515(1):77–84, 2019.[34] H. Meinhardt. Models for positional signalling with application to the dorsoventral patterningof insects and segregation into different cell types.
Development , 107(Supplement):169–180,1989.[35] J. Murray.
Mathematical biology II: spatial models and biomedical applications . Springer NewYork, 2001.[36] B. Nagorcka and J. Mooney. From stripes to spots: prepatterns which can be produced in theskin by a reaction-diffusion system.
Mathematical Medicine and Biology: A Journal of the IMA ,9(4):249–267, 1992.[37] K. Oizel, C. Chauvin, L. Oliver, C. Gratas, F. Geraldo, U. Jarry, E. Scotet, M. Rabe, M.-C. Alves-Guerra, R. Teusan, et al. Efficient mitochondrial glutamine targeting prevails overglioblastoma metabolic plasticity.
Clinical Cancer Research , 23(20):6292–6304, 2017.[38] L. Oliver, L. Lalier, C. Salaud, D. Heymann, P. F. Carton, and F. Vallette. Drug resistance inglioblastoma: Are persisters the key to therapy?
Under review , 2020.[39] L. Oliver, A. ´Alvarez Arenas, C. Salaud, et al. A simple 3d cell culture method for studyingthe interactions between primary glioblastoma cells and tumor-activated stromal cells.
Underreview , 2020.[40] H. G. Othmer and T. Hillen. The diffusion limit of transport equations ii: Chemotaxis equations.
SIAM Journal on Applied Mathematics , 62(4):1222–1250, 2002.[41] Q. Ouyang and H. L. Swinney. Transition from a uniform state to hexagonal and striped turingpatterns.
Nature , 352(6336):610–612, 1991.[42] K. Painter. Modelling cell migration strategies in the extracellular matrix.
Journal of Mathe-matical Biology , 58(4-5):511, 2009.[43] K. J. Painter and T. Hillen. Volume-filling and quorum-sensing in models for chemosensitivemovement.
Canadian Applied Mathematics Quarterly , 10(4):501–543, 2002.[44] W. Roos, L. Batista, S. Naumann, W. Wick, M. Weller, C. Menck, and B. Kaina. Apoptosisin malignant glioma cells triggered by the temozolomide-induced dna lesion o 6-methylguanine.
Oncogene , 26(2):186–197, 2007. 3645] J. Saragosti, V. Calvez, N. Bournaveas, A. Buguin, P. Silberzan, and B. Perthame. Mathe-matical description of bacterial traveling pulses.
PLoS Computational Biology , 6(8):e1000890,2010.[46] W. Sarfaraz and A. Madzvamuse. Domain-dependent stability analysis of a reaction–diffusionmodel on compact circular geometries.
International Journal of Bifurcation and Chaos ,28(08):1830024, 2018.[47] A. Stevens. The derivation of chemotaxis equations as limit dynamics of moderately interactingstochastic many-particle systems.
SIAM Journal on Applied Mathematics , 61(1):183–212, 2000.[48] R. Stupp, M. E. Hegi, W. P. Mason, M. J. Van Den Bent, M. J. Taphoorn, R. C. Janzer, S. K.Ludwin, A. Allgeier, B. Fisher, K. Belanger, et al. Effects of radiotherapy with concomitant andadjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomisedphase iii study: 5-year analysis of the eortc-ncic trial.
The Lancet Oncology , 10(5):459–466,2009.[49] R. Stupp, W. P. Mason, M. J. Van Den Bent, M. Weller, B. Fisher, M. J. Taphoorn, K. Belanger,A. A. Brandes, C. Marosi, U. Bogdahn, et al. Radiotherapy plus concomitant and adjuvanttemozolomide for glioblastoma.
New England Journal of Medicine , 352(10):987–996, 2005.[50] T. Suzuki.
Chemotaxis, reaction, network . World Scientific Publishing Co. Pte. Ltd., Hacken-sack, NJ, 2018. Mathematics for self-organization.[51] I. Talhaoui, B. T. Matkarimov, T. Tchenio, D. O. Zharkov, and M. K. Saparbaev. Aberrantbase excision repair pathway of oxidatively damaged dna: implications for degenerative diseases.
Free Radical Biology and Medicine , 107:266–277, 2017.[52] L. Tang, A. L. Van De Ven, D. Guo, V. Andasari, V. Cristini, K. C. Li, and X. Zhou. Com-putational modeling of 3d tumor growth and angiogenesis for chemotherapy evaluation.
PloSone , 9(1):e83962, 2014.[53] A. Toma, L. R. C. Castillo, T. A Schuetz, S. Becker, A. Mang, A. R´egnier-Vigouroux, andT. M Buzug. A validated mathematical model of tumour-immune interactions for glioblastoma.
Current Medical Imaging , 9(2):145–153, 2013.[54] J. van Rijn, J. J. Heimans, J. van den Berg, P. van der Valk, and B. J. Slotman. Survival ofhuman glioma cells treated with various combination of temozolomide and x-rays.
InternationalJournal of Radiation Oncology* Biology* Physics , 47(3):779–784, 2000.[55] Z. Wang. On chemotaxis models with cell population interactions.
Mathematical Modelling ofNatural Phenomena , 5(3):173–190, 2010.[56] Z. Wang and T. Hillen. Classical solutions and pattern formation for a volume filling chemotaxismodel.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 17(3):037108, 2007.[57] S. Ward. Chemotaxis by the nematode caenorhabditis elegans: identification of attractants andanalysis of the response by use of mutants.
Proceedings of the National Academy of Sciences ,70(3):817–821, 1973.[58] D. Wrzosek. Volume filling effect in modelling chemotaxis.
Mathematical Modelling of NaturalPhenomena , 5(1):123–147, 2010. 3759] R. W¨urth, A. Bajetto, J. K. Harrison, F. Barbieri, and T. Florio. Cxcl12 modulation of cxcr4and cxcr7 activity in human glioblastoma stem-like cells and regulation of the tumor microen-vironment.
Frontiers in Cellular Neuroscience , 8:144, 2014.[60] H. Yamaguchi, J. Wyckoff, and J. Condeelis. Cell migration in tumors.