Stability of a non-local kinetic model for cell migration with density dependent orientation bias
SStability of a non-local kinetic model for cell migration withdensity dependent orientation bias
Nadia Loy ∗ Luigi Preziosi †‡ January 23, 2020
Abstract
The aim of the article is to study the stability of a non-local kinetic model proposed byLoy and Preziosi (2019a). We split the population in two subgroups and perform a linearstability analysis. We show that pattern formation results from modulation of one non-dimensional parameter that depends on the tumbling frequency, the sensing radius, the meanspeed in a given direction, the uniform configuration density and the tactic response to the celldensity. Numerical simulations show that our linear stability analysis predicts quite preciselythe ranges of parameters determining instability and pattern formation. We also extend thestability analysis in the case of different mean speeds in different directions. In this case, forparameter values leading to instability travelling wave patterns develop.
Keywords : Kinetic model, Non-local interactions, Stability, Cell migration.
Cells decide where to go by sensing the surrounding environment, extending their protrusionsover a distance that can reach several cell diameters. In order to describe such a non-local actionseveral mathematical models have been recently proposed in the literature.Othmer and Hillen (2002) and Hillen et al. (2007) introduced a finite sampling radius anddefined a non-local gradient as the average of the external field on a surface which represents themembrane of the cell. In particular, the Authors derived a macroscopic diffusive model with a non-local gradient both from a position jump process and from a velocity jump process by postulatingthat the non-local sensing is a bias of higher order.With the aim of modelling cell-cell adhesion and haptotaxis Armstrong et al. (2006) proposeda macroscopic integro-differential equation where the integral over a finite radius is in chargeof describing non-local sensing. Recently, Buttensch¨on et al. (2018) derived this model from aspace jump process, while Buttensch¨on (2018), Buttensch¨on and Hillen (2020) studied the steadystates and bifurcations of the macroscopic adhesion model and their stability. Further studiesconcerning these models on bounded domains are proposed in (Buttensch¨on and Hillen, 2019).Other macroscopic models describing cell migration with non-local measures of the environmentwere proposed by Painter et al. (2010, 2015) and Painter and Hillen (2002). Schmeiser and Nouri(2017) considered a kinetic model with velocity jumps biased towards the chemical concentration ∗ Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca degli Abruzzi24, 10129 Torino, Italy, and Department of Mathematics “G. Peano”, Via Carlo Alberto 10 ,10123 Torino, Italy( [email protected] ) † Department of Mathematical Sciences “G. L. Lagrange”, Dipartimento di Eccellenza 2018-2022, Politecnico diTorino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy ( [email protected] ) ‡ Corresponding author: [email protected] a r X i v : . [ q - b i o . CB ] J a n radient. Similar equations were also proposed in 2D set-ups by Colombi et al. (2015, 2017), andapplied to model crowd dynamics and traffic flow for instance by Tosin and Frasca (2011).Eftimie et al. (2007b) proposed a non-local kinetic models including repulsion, alignement andattraction. They distinguish cells going along the two directions of a one-dimensional set-up withthe possibility of switching between the two directions keeping the same speed (they thus havetwo velocities). The one-dimensional kinetic model takes then a discrete velocity structure. Itsintegration shows a wide variety of patterns. The linear stability analysis of the model is thenperformed in Eftimie et al. (2007a). Carrillo et al. (2015) derived the macroscopic limits of thismodel. In Eftimie et al. (2017) and Bitsouni and Eftimie (2018) the model was applied to model,respectively, tumour dynamics and cell polarisation in heterogeneous cancer cell populations. Thelinear stability analysis of a kinetic chemotaxis equation coupled with a macroscopic diffusionequation for the chemical dynamics is presented by Perthame and Yasuda (2018).Loy and Preziosi (2019a) proposed a non-local kinetic model with double-bias on the basis ofthe observation that different fields can influence respectively cell polarization and speed. Then theturning operator is characterized by the presence of two integrals, one determining the probabilityof polarizing in a certain direction after averaging the sensing of a chemical or mechanical cueover a finite neighborhood, and the other setting the probability of moving with a certain speedin the chosen direction after averaging the sensing of another chemical or mechanical cue overa possibly different finite radius. The model was then modified in Loy and Preziosi (2019b)to take into account of physical limits of migrations, that might hamper the real possibility ofcells of measuring beyond physical barriers or to move in certain regions because of physicallyconstraining situations, such as too dense extracellular matrix or cell overcowding. In such papersit was shown that the presence of density dependent cues, such as cell-cell adhesion and volumefilling, may generate instabilities and the formation of patterns. In fact, on the one hand cells maybe attracted due to the mutual interaction of transmembrane adhesion molecules ( e . g . , cadherincomplexes). On the other hand, they may want to stay away from overcrowded areas. It seemsthat, depending on the sensing kernel and on other modelling parameters, both the wish to staytogether and to stay away might lead to instability.In this paper we study the stability of the homogeneous configuration when density dependentcues influence cell polarization. The case in which the density distribution also affects cell speedwill be treated in a following paper.With this aim in mind, in Section 2 we briefly recall the non-local kinetic model, then focalizingto the case of orientational biases. Restricting to the one-dimensional case, in Section 3 the linearstability analysis is performed, first for a general speed distribution function and then in theparticular case of a Dirac delta. At this stage the sensing kernel is still general and it is provedthat if it is a non-increasing function of the distance, then staying away strategies are alwaysstable. On the contrary, Section 3.1 shows that a localized sensing kernel, i . e . a Dirac delta,might lead to instability of a finite wavelength if for instance, cell speed is sufficiently small, or theturning rate or the sensing radius is sufficiently high, or in the case of volume filling effects, if celldensity is sufficiently high. Section 3.2 and Section 3.3 then respectively focus on a homogeneousand a decreasing sensing kernel over a finite range, specifically a Heaviside function and a rampfunction. It these cases, as well as for the localized sensing kernel, cell-cell adhesion strategies leadto long wave instabilities. Section 4 reports some simulations of the cases above, while Section5 discusses the case of an asymmetric speed distribution function in the two directions. In theunstable case this leads to the formation of travelling waves instabilities similar to those reportedby Eftimie et al. (2007a), Eftimie et al. (2007b), Carrillo et al. (2015), Eftimie (2012). A finalsection draws some conclusions pointing out possible developments. In the model introduced by Loy and Preziosi (2019a) the cell population is described at a meso-scopic level by the distribution density p = p ( t, x , v, ˆ v ) parametrized by the time t >
0, the position x ∈ Ω ⊆ R d , the speed v ∈ R + and the polarization direction ˆ v ∈ S d − where S d − is the unit2phere boundary in R d . We remark that the distribution function p depends separately on velocitymodulus and direction, instead of the velocity vector v = v ˆ v . This is due to the need of separatingthe subcellular mechanisms governing cell polarization and motility. In fact, cells respond both totactic factors affecting the choice of the direction, and to kinetic factors, typically of mechanicalorigins, influencing cell speed.The mesoscopic model consists in the transport equation for the cell distribution ∂p∂t ( t, x , v, ˆ v ) + v · ∇ p ( t, x , v, ˆ v ) = J [ p ]( t, x , v, ˆ v ) (1)where the operator ∇ denotes the spatial gradient then coupled with proper initial and boundaryconditions. In particular, we shall consider no-flux boundary conditions that are defined as (Plaza,2019) (cid:90) R + (cid:90) S d − p ( t, x , v, ˆ v ) v · n ( x ) d ˆ v dv = 0 , ∀ x ∈ ∂ Ω , t > n ( x ) the outer normal to the boundary ∂ Ω in the point x . Equation (2) implies that thereis no mass flux across the boundary (Lemou and Mieussens, 2008).A macroscopic description for the cell population can be classically recovered through thedefinition of moments of the distribution function p . For instance, the cell number density will begiven by ρ ( t, x ) = (cid:90) S d − (cid:90) R + p ( t, x , v, ˆ v ) dv d ˆ v , (3)and the cell mean velocity by U ( t, x ) = 1 ρ ( t, x ) (cid:90) S d − (cid:90) R + p ( t, x , v, ˆ v ) v dv d ˆ v . (4)The term J [ p ]( t, x , v, ˆ v ), named turning operator , is an integral operator that describes thechange in velocity which is not due to free-particle transport. It may describe the classical runand tumble behaviors, random re-orientations, which, however, may be biased by external cues.In the present case, the turning operator will be the implementation of a velocity-jump process ina kinetic transport equation as introduced by Stroock (1974) and then by Othmer et al. (1988).The turning operator that we are going to consider is J [ p ]( t, x , v, ˆ v ) = µ ( x ) (cid:16) ρ ( t, x ) T ( x , v, ˆ v ) − p ( t, x , v, ˆ v ) (cid:17) , (5)which is obtained assuming that cells retain no memory of their velocity prior to the re-orientation,and also turning rates do not depend on the orientation of the individual cell.Following Loy and Preziosi (2019a), we consider here a transition probability that depends onthe non-local sensing of the macroscopic density of the cell population in a neighborhood of thecell, that writes as T [ ρ ]( x , v, ˆ v ) = c ( t, x ) (cid:90) R + γ R ( λ ) b (cid:0) ρ ( t, x + λ ˆ v ) (cid:1) dλ ψ ( v | ˆ v ) , (6)where c ( t, x ) is a normalization constant, that is c − ( t, x ) = (cid:90) S d − (cid:34)(cid:90) R + b ( ρ ( t, x + λ ˆ v )) γ R ( λ ) dλ (cid:35) ψ ( v | ˆ v ) d ˆ v , so that the integral of T over the velocity space is one. In this way, T is a mass preservingtransition probability.The term b describes the response of the cells to the tactic cue, in this case the cell density ρ itself, around x along the direction ˆ v and, therefore, the bias intensity in the direction ˆ v while3he sensing kernel γ R ( λ ) weights the collected density signal with respect to the distance λ from x . In particular, γ R , that we shall assume to be a positive valued L ( R + ) function, has a compactsupport in [0 , R ] where R is the maximum extension of cell protrusions, determining the furthestpoints cells can reach to measure the external signals. So, integrals over R + are actually integralsover the finite interval [0 , R ]. Specifically, if • γ R ( λ ) = δ ( λ − R ) is a Dirac delta, then cells only measure the information perceived on aspherical surface of given radius R , • if γ R ( λ ) = H ( R − λ ) is a Heaviside function, then cells explore the whole volume of thesphere centered in x with radius R and weight the information uniformly, or • if γ R ( λ ) is a decreasing function of λ , then closer information play a bigger role with re-spect to farther ones, taking for instance into account that the probability of making longerprotrusions decreases with the distance, so the sensing of closer regions is more accurate.In conclusion, the integral in (6) is such that if the signal is stronger in the direction ˆ v , thenthere will be a higher probability for the cell to move along ˆ v than along − ˆ v . Other functionalforms for the b term could be considered as discussed by Loy and Preziosi (2019a).The function ψ = ψ ( v | ˆ v ) is the density distribution of the speeds that we assume to dependon the direction ˆ v ∈ S d − . In fact, Loy and Preziosi (2019a) introduce a density distribution ψ describing the probability of having a certain speed in a given direction given the non local sensingof a tactic external cue in that direction. Therefore, ψ also depends on the spatial variable throughthe kinetic cue. For simplicity we will drop this extra dependency.In the following the mean of the probability density function ψ will be denoted by V and itsvariance by s . As ψ depends on the direction ˆ v , V and s will depend on the direction ˆ v as well. In order to perform a stability analysis of the uniform configuration, we will assume that µ isconstant and for the moment we will start assuming that ψ = ψ ( v | ˆ v ) does not depend on ˆ v . Wewill treat the case in which ψ = ψ ( v | ˆ v ) in Section 5. Hence, we consider the equation ∂p∂t ( t, x , v, ˆ v )+ v ˆ v ·∇ p ( t, x , v, ˆ v ) = µ (cid:34) c ( t, x ) ρ ( t, x ) (cid:90) R + γ R ( λ ) b ( ρ ( t, x + λ ˆ v )) dλ ψ ( v ) − p ( t, x , v, ˆ v ) (cid:35) . (7)We will also carry out the analysis in the one-dimensional case and call e and − e the two possibledirections characterizing the one dimensional problem. The population is then splitted into twosubgroups p + and p − corresponding to the groups of cells respectively going to the right and tothe left, so that p ( t, x, v, ˆ v ) = p + ( t, x, v ) δ (ˆ v − e ) + p − ( t, x, v ) δ (ˆ v + e ) . (8)Coherently, we define ρ ( t, x ) = ρ + ( t, x ) + ρ − ( t, x ) with ρ ± ( t, x ) = (cid:90) R + p ± ( t, x, v ) dv . The system of equations satisfied by p + and p − is ∂p + ∂t ( t, x, v ) + v ∂p + ∂x ( t, x, v ) = µ [ ρ ( t, x ) T + [ ρ ]( v ) − p + ( t, x, v )] ∂p − ∂t ( t, x, v ) − v ∂p − ∂x ( t, x, v ) = µ [ ρ ( t, x ) T − [ ρ ]( v ) − p − ( t, x, v )] , (9)4here T + [ ρ ]( v ) = (cid:90) R + γ R ( λ ) b ( ρ ( t, x + λ )) dλ (cid:90) R + γ R ( λ ) [ b ( ρ ( t, x + λ )) + b ( ρ ( t, x − λ ))] dλ ψ ( v ) (10)and T − [ ρ ]( v ) = (cid:90) R + γ R ( λ ) b ( ρ ( t, x − λ )) dλ (cid:90) R + γ R ( λ ) [ b ( ρ ( t, x + λ )) + b ( ρ ( t, x − λ ))] dλ .ψ ( v ) (11)It can be proved that the local asymptotic equilibrium states of (9) are (Loy and Preziosi, 2019a) p + ∞ = ρ ∞ T + [ ρ ∞ ]( v ) , p −∞ = ρ ∞ T − [ ρ ∞ ]( v ) . (12)They are homogeneous if and only if p + ∞ ( v ) = p −∞ ( v ) = ρ ∞ ψ ( v )2 . (13)Hence, once ψ is chosen, all the possible spatially homogeneous and stationary solutions aredetermined by the choice of ρ ∞ . We explicitly notice that the macroscopic densities will be ρ + ∞ = ρ −∞ = ρ ∞ p + ( t, x, v ) = p + ∞ ( v ) + ˆ p + ( t, x, v ) , p − ( t, x, v ) = p −∞ ( v ) + ˆ p − ( t, x, v ) . Hence, ρ ( t, x ) = ρ ∞ + (cid:90) R + (cid:90) S d − (cid:2) ˆ p + ( t, x, v ) δ (ˆ v − e ) + ˆ p − ( t, x, v ) δ (ˆ v + e ) (cid:3) d ˆ v dv = ρ ∞ + ˆ ρ + ( t, x ) + ˆ ρ − ( t, x )where ˆ ρ ± ( t, x ) = (cid:90) R + ˆ p ± ( t, x, v ) dv and we define ˆ ρ = ˆ ρ + + ˆ ρ − . Thus, the system (9) becomes ∂ ˆ p + ∂t ( t, x, v ) + v ∂ ˆ p + ∂x ( t, x, v ) = µ (cid:2)(cid:0) ρ ∞ + ˆ ρ (cid:1) T + [ ρ ∞ + ˆ ρ ]( v ) − p + ∞ ( v ) − ˆ p + ( t, x, v ) (cid:3) ,∂ ˆ p − ∂t ( t, x , v ) − v ∂ ˆ p + ∂x ( t, x, v ) = µ (cid:2)(cid:0) ρ ∞ + ˆ ρ (cid:1) T − [ ρ ∞ + ˆ ρ ]( v ) − p + ∞ ( v ) − ˆ p − ( t, x, v ) (cid:3) . We now need to linearize the transition probabilities T ± [ ρ ∞ + ˆ ρ ]( v ) where, for instance, T + [ ρ ∞ + ˆ ρ ]( v ) = (cid:90) R + b ( ρ ∞ + ˆ ρ ( t, x + λ )) γ R ( λ ) dλ (cid:90) R + [ b ( ρ ∞ + ˆ ρ ( t, x + λ )) + b ( ρ ∞ + ˆ ρ ( t, x − λ ))] γ R ( λ ) dλ ψ ( v ) . Assuming ˆ ρ small (and then neglecting perturbation of higher order), we may perform a Taylorexpansion and write T + [ ρ ∞ + ˆ ρ ]( v ) ≈ (cid:90) R + [ b ( ρ ∞ ) + b (cid:48) ( ρ ∞ )ˆ ρ ( x + λ )] γ R ( λ ) dλ (cid:90) R + [ b ( ρ ∞ ) + b (cid:48) ( ρ ∞ )ˆ ρ ( x + λ ) + b ( ρ ∞ ) + b (cid:48) ( ρ ∞ )ˆ ρ ( x − λ )] γ R ( λ ) dλ ψ ( v ) ≈ b (cid:48) ( ρ ∞ ) b ( ρ ∞ )Γ R (cid:90) R + ˆ ρ ( t, x + λ ) γ R ( λ ) dλ b (cid:48) ( ρ ∞ ) b ( ρ ∞ )Γ R (cid:90) R + ˆ ρ ( t, x + λ ) + ˆ ρ ( t, x − λ )2 γ R ( λ ) dλ ψ ( v )2 (14)5eing Γ R = (cid:90) R + γ R ( λ ) dλ . Expanding the denominator, we eventually have T + [ ρ ∞ + ˆ ρ ]( v ) ≈ (cid:34) b (cid:48) ( ρ ∞ ) b ( ρ ∞ )Γ R (cid:90) R + [ˆ ρ ( t, x + λ ) − ˆ ρ ( t, x − λ )]2 γ R ( λ ) dλ (cid:35) ψ ( v )2 . (15)Analogously T − [ ρ ∞ + ˆ ρ ]( v ) ≈ (cid:34) − b (cid:48) ( ρ ∞ ) b ( ρ ∞ )Γ R (cid:90) R + [ˆ ρ ( t, x + λ ) − ˆ ρ ( t, x − λ )]2 γ R ( λ ) dλ (cid:35) ψ ( v )2 . (16)Therefore, the right hand sides of the system (3) become (cid:0) ρ ∞ + ˆ ρ ( t, x ) (cid:1) T ± [ ρ ∞ + ˆ ρ ]( v ) − p ±∞ ( v ) − ˆ p ± ( t, x, v ) ≈ =0 (cid:122) (cid:125)(cid:124) (cid:123) ρ ∞ ψ ( v )2 − p ±∞ ( v ) + (cid:34) ˆ ρ ( t, x ) ± ρ ∞ b (cid:48) ( ρ ∞ ) b ( ρ ∞ )Γ R (cid:90) R + [ˆ ρ ( t, x + λ ) − ˆ ρ ( t, x − λ )]2 γ R ( λ ) dλ (cid:35) ψ ( v )2 − ˆ p ± ( v ) . Let us now consider perturbations in the formˆ p ± ( t, x, v ) = g ± ( v ) e ikx + σt where g ± have densities defined as ρ g ± = (cid:82) R + g ± ( v ) dv and, then, ˆ ρ ± = ρ g ± e ikx + σt .Substitution in (3) leads to σg ± ± ikvg ± + µg ± = µ (cid:34) ρ g ± ρ ∞ b (cid:48) ( ρ ∞ ) b ( ρ ∞ )Γ R ρ g (cid:90) R + e ikλ − e − ikλ γ R ( λ ) dλ (cid:35) ψ ( v )2 , (17)where we set ρ g = ρ + g + ρ − g . Now, as e ikλ − e − ikλ = 2 i sin( kλ ), definingˆ γ R ( k ) = 1Γ R (cid:90) R + sin( kλ ) γ R ( λ ) dλ , (18)the normalized unilateral sine transform of γ R , the system (17) rewrites as ( σ + ikv + µ ) g + = µ [1 + i B ˆ γ R ( k )] ψ ( v )2 (cid:90) R + (cid:2) g + ( v ) + g − ( v ) (cid:3) dv, ( σ − ikv + µ ) g − = µ [1 − i B ˆ γ R ( k )] ψ ( v )2 (cid:90) R + (cid:2) g + ( v ) + g − ( v ) (cid:3) dv, (19)where B = ρ ∞ b (cid:48) ( ρ ∞ ) b ( ρ ∞ ) . (20)We remark that the dimensionless number B can be either negative or positive, according to thefact that b is a decreasing or an increasing function of ρ . From the phenomenological point of viewthe former case ( B < b (cid:48) ( ρ ∞ ) <
0) corresponds to a predisposition of cells at a density ρ ∞ to re-polarize toward regions that have lower cell densities and move away from crowded areas,the latter case ( B > b (cid:48) ( ρ ∞ ) >
0) corresponds to a predisposition of cells to re-orient towardregions with a density higher than ρ ∞ , e . g . a sort of adhesion-like behaviour due to the fact thatcells wants to stay together.We observe that by summing the two equations in (19) we readily have (cid:0) σ + µ (cid:1)(cid:0) g + + g − (cid:1) + ikv (cid:0) g + − g − (cid:1) = µψ ( v ) (cid:90) R + (cid:2) g + ( v ) + g − ( v ) (cid:3) dv.
6f we then integrate over R + , by defining ρ g U g = (cid:90) R + (cid:0) g + − g − (cid:1) v dv the mean momentum of theperturbation, we get ( σ + µ ) ρ g + ikρ g U g = µρ g , and, then, the mean speed of the perturbation is U g = i σk . Looking for unstable situations (so that it is sure that we are not dividing by zero) the system(19) can be written as g + = µ i B ˆ γ R ( k ) σ + ikv + µ ψ ( v )2 (cid:90) R + (cid:2) g + ( v ) + g − ( v ) (cid:3) dv,g − = µ − i B ˆ γ R ( k ) σ − ikv + µ ψ ( v )2 (cid:90) R + (cid:2) g + ( v ) + g − ( v ) (cid:3) dv. (21)If we integrate the two equations over R + and sum them, we obtain the following solvabilitycondition (for non trivial solutions) µ (cid:90) R + (cid:20) i B ˆ γ R ( k ) σ + ikv + µ + 1 − i B ˆ γ R ( k ) σ − ikv + µ (cid:21) ψ ( v ) dv = 1 , (22)or µ (cid:90) R + σ + µ + B ˆ γ R ( k ) kv ( σ + µ ) + k v ψ ( v ) dv = 1 . (23)In order to give an analytical discussion of the result, let us take, as an example, ψ ( v ) = δ ( v − V ).In this case the integral condition (23) takes the algebraic form µ σ + µ + B ˆ γ R ( k ) kV ( σ + µ ) + k V = 1 , (24)so that the dispersion relation reads σ + µσ + k V − µ B kV ˆ γ R ( k ) = 0 . The most dangerous eigevalue is then σ = − µ + √ ∆2 , with ∆ = µ − k V + 4 µ B kV ˆ γ R ( k ) . (25)Therefore, the instability condition (cid:60) e ( σ ) >
0, given by ∆ >
0, is B ˆ γ R ( k ) k > Vµ . (26)that, introducing the dimensionless number V = VµR , (27)may be rewritten as B ˆ γ R ( k ) Rk > V . (28)Therefore, one can conclude that the following proposition holds.7 roposition 3.1. If B ≤ and γ R ( λ ) > is non increasing, then the homogeneous configurationis always stable.If B > and γ R ( λ ) > is such that λγ R ( λ ) ∈ L ( R + ) , then long waves k ≈ are unstable forsufficiently small values of V = Vµ B .Proof. In order to prove the statement it is enough to observe that if the sensing kernel γ R ( λ ) isnon increasing then ˆ γ R ( k ) > ∀ k >
0, vanishing only in the trivial cases k = 0 or in the limitof constant sensing kernel in R + . Then condition (28) is never satisfied, being its r.h.s. strictlypositive.On the other hand, under the stated integrability conditionslim k → ˆ γ R ( k ) k = (cid:90) ∞ λγ R ( λ ) dλ (cid:90) ∞ γ R ( λ ) dλ > . Hence, if B >
0, then condition (28) is satisfied for sufficiently small ratios Vµ B leading to instability.We observe that the statement also holds for the relevant case of a sensing kernel with compactsupport, for instance for the Heaviside and ramp kernels that will be respectively consideredin Sections 3.2 and 3.3. On the other hand, it does not hold for the localized sensing kernel γ R ( λ ) = δ ( λ − R ) that will be considered in the following section. In fact, for such kernelsinstability is possible also for B < Vµ B that is unstable at leastto long waves.We conclude this general part of the analysis by observing that in unstable situations, themaximum growth rate σ max = σ ( k max ) can be identified by evaluating the stationary points of(25), that are obtained when1Γ R k max R (cid:90) R + [sin( k max λ ) + k max λ cos( k max λ )] γ R ( λ ) dλ = 2 V b , (29)where V b = V / B . (30) We assume now that also the sensing function is a Dirac delta, i . e . γ R = δ ( λ − R ), so that thetwo equations satisfied by p + and p − (9) specialize as ∂p + ∂t ( t, x, v ) + v ∂p + ∂x ( t, x, v ) = µ (cid:20) ρ ( t, x ) b ( ρ ( t, x + R )) b ( ρ ( t, x + R )) + b ( ρ ( t, x − R )) ψ ( v ) − p + ( t, x, v ) (cid:21) ,∂p − ∂t ( t, x, v ) − v ∂p − ∂x ( t, x, v ) = µ (cid:20) ρ ( t, x ) b ( ρ ( t, x − R )) b ( ρ ( t, x + R ) + b ( ρ ( t, x − R )) ψ ( v ) − p − ( t, x, v ) (cid:21) . (31)In this case, Γ R = 1 and ˆ γ R ( k ) = sin( kR ), so that the criterium (28) now reads B sin( kR ) kR > V . (32)From (29) the maximum growth rate is obtained for k max such thatsin( k max R ) k max R + cos( k max R ) = 2 V b . (33)To discuss the stability properties it is now useful to distinguish two cases according to thesign of b (cid:48) , which means the sign of B or of V b . 8 a) b (cid:48) ( ρ ∞ ) < b (cid:48) ( ρ ∞ ) > Figure 1: Stability diagram for a localized sensing kernel γ R = δ ( λ − R ). The red dashed linedelimits the unstable region, in (a) when b (cid:48) ( ρ ∞ ) < b (cid:48) ( ρ ∞ ) >
0. The blue dottedlines evidentiate the dimensionless wave numbers k max R with local maxima of the growth rate(given by (33)) in the unstable regime. The lowest curves correspond to the most dangerous wavenumbers in both cases. b (cid:48) ( ρ ∞ ) < b corresponds to a pre-disposition of cells to re-orient toward regions that are less crowded than the uniform stationarysolution ρ ∞ .In this case, as B <
0, the instability condition (32) becomessin( kR ) kR < V b < , (34)and the critical condition, i . e . the first value for whichsin( kR ) kR = V b , (35)is given by V b,cr = min x> sin xx = − m ≈ − . b ( ρ ) = (cid:16) − ρρ th (cid:17) + , then for ρ ∞ < ρ th , V b = −V (cid:16) ρ th ρ ∞ − (cid:17) and there is instabilityif ρ ∞ > ρ th V / ( V + m ).The critical wave number is then obtained when (cid:20) sin( kR ) kR (cid:21) (cid:48) k = k cr = 0 ⇐⇒ tan( k cr R ) = k cr R, (36)that is when k cr R = ¯ απ with ¯ α ≈ .
43. Therefore, the critical wave length Λ cr = 2 πk cr is such thatΛ cr = 2¯ α R ≈ . R. If V b < − m the system (3) is always stable, while if V b ∈ ( − m,
0) there are unstable wave numbers.Using (35), Figure 1a identifies for any |V b | < |V b,cr | ≈ .
22 the range of unstable wave numbersand the instability region is the one to the left of the red dashed curve. In this unstable region local9aximum growth rates are represented by the blue curves, with the longest waves (correspondingto the lowest curve) being the most unstable ones.We then have instability of a finite wavelength if for instance, cell speed is sufficiently small,or the turning rate or the sensing radius is sufficiently high, or in the case of volume filling effects,if cell density if sufficiently high b (cid:48) ( ρ ∞ ) > b corresponds to a predispo-sition of cells to re-orient toward regions that are more crowded than in the uniform configuration ρ ∞ . This might be for instance related to an adhesion-like behaviour for cells that want to staytogether. An example of this case is b ( ρ ) = ρ , leading to B = 1 and V b = VRµ .In this case the instability condition becomessin( kR ) kR > V b > . (37)We may trivially observe that, as sin xx <
1, if V b > V b ∈ [0 ,
1) wave numbers to the left of the red dashed curves represented inFigure 1(b) are unstable. Again the blue lines represent local maxima for the growth rates withthe lowest curve corresponding to the most unstable dimensionless wave number.
The same stability analysis can be performed for a sensing function γ R that is a Heaviside function, i . e . γ R ( λ ) = H ( R − λ ). In this case, ˆ γ R ( k ) = 1 − cos( kR ) kR . (38)Contrary to the localized kernel, for the uniform kernel the first statement of the Proposition holdsand the uniform configuration with ρ = ρ ∞ is stable when B ≤
0, that is cells at that densityprefer to avoid overcrowding.On the other hand, if B > − cos( kR ) k R > V b . (39)Therefore, there are unstable waves if V b < and, recalling (29), the maximum growth rate isachieved for k = k max such that sin( k max R ) k max R = 2 V b . (40)Referring to Figure 2(a), one then has instability in the region of the V b − kR plane to the leftof the red dashed curve with the most unstable wave number identified by the lowest blue curve. A sensing function γ R that decreases with the distance from the present position of the cell meansthat the cell gives more importance to the local information than to distant ones. An example isgiven by the ramp function γ R ( λ ) = (cid:18) − λR (cid:19) + . where ( f ) + is the positive part of f . In this case,ˆ γ R ( k ) = 2 kR − sin( kR ) k R , (41)10 a) (b) Figure 2: Stability diagram for a uniform (a) and a ramp (b) sensing kernel. The blue dotted linerepresents k max R given respectively by (a) (40) and (b) (43). The unstable region is the one tothe left of the red dashed line, i . e . the values of k max R also satisfy (39) in (a) and (42) in (b).Again, independently of the specific form of the decreasing kernel, the uniform configuration with ρ = ρ ∞ is stable when B ≤
0, that is, if cells at that density prefer to avoid overcrowding.On the other hand, if B > kR − sin( kR ) k R > V b . (42)Therefore, there are unstable waves if V b < and, recalling (29), the maximum growth rate isachieved for k = k max such thatsin( k max R ) − k max R cos( k max R ) k max R = V b . (43)In Figure 2(b), again one then has instability in the region of the V b − kR plane to the left ofthe red dashed curve with the most unstable wave identified by the lowest blue curve. In this section we show some numerical tests in order to show the consistency of our linear stabilityanalysis. We simulate equation (7) with specular reflective or periodic boundary conditions, bothsatisfying Eq. (2). Specular reflective boudary conditions in one dimension are such that p − ( t, L, v ) = p + ( t, L, − v ) , p + ( t, , v ) = p − ( t, , − v ) . (44)We perform a first order splitting for the relaxation and transport step that we perform usinga Van Leer scheme. For further details, we address the reader to (Loy and Preziosi, 2019a) and(Filbet and Vauchelet, 2010). We remark that in our case the ψ ( v | ± e ) is a Gaussian with mean V ± and variance s . In particular we consider s = 10 − , so that ψ is close in sense of measureto a Dirac delta and the perturbation of the dispersion relation is of the same order of s . In order to mimick the dynamics of cells that are more likely to re-orient where there are less cellsand tend to avoid overcrowded areas, corresponding to the case b (cid:48) ( ρ ) < b ( ρ ) = (cid:18) − ρρ th (cid:19) + . a) (b)(c) Figure 3: Temporal evolution of ρ ( t, x ) from ρ ( x ) = 0 . . πx/ V b ≈ − . V b ≈ − . Ω ρ ( x ) < ρ th . Otherwise, physicalconstraint effects should be taken into account as in Loy and Preziosi (2019b). Specifically, weconsider a Dirac delta sensing function γ and we set ρ th = 0 . ρ ( x ) =0 . . πx/ ρ ∞ ≈ . B ≈ − . V b is about − .
22. Having set V = 0 .
25 and R = 0 .
2, in Figure 3(a) we use µ = 15, so that V = 1 /
12 andeverywhere |V b | is above the critical value leading to instability. In fact, its value is V b ≈ − . V = 0 . R = 0 . µ = 5, so that V = 1 / V b ≈ − . ρ = ρ ∞ ≈ . |V b | closer to zero and a practical coincidence forvalues above 0.06.We recall that in this case the first statement of the Proposition does not hold while in thecase of smoother non increasing kernels (e.g., the Heaviside or ramp kernels treated below), wealways have stability. We consider now the case in which cells prefer to stay together and reorient toward regions withmore crowded areas. Specfically, we take b ( ρ ) = ρ so that b (cid:48) ( ρ ) = 1 is positive. Therefore, inthis case B = 1 regardless of the density distribution and then V b = V . Having set V = 0 .
25 and R = 0 .
04, changes in V b correspond to changes in µ .In Figure 4 we show linear stability and instability in the case of cell-cell adhesion and perfectreflective boundary conditions. The sensing function is a Dirac delta in the top row and a Heavisidefunction in the bottom row.Recalling Figure 1, in Figure 4(a), µ = 5 .
5, so that V b ≈ . > µ = 8, so that V b ≈ . < µ = 180, so that V b ≈ . > .
5, corresponds to astable condition, while in Figure 4(d), we µ = 260, so that V b ≈ . < . V = 0 . ψ We now consider Eq. (1) with (5) and (6) allowing the speed probability distribution to dependon ˆ v , i.e. ψ = ψ ( v | ˆ v ). In the one-dimensional case it means that cells that go to the right and tothe left have different speed distributions.Following the same procedure as in Section 3, we obtain again Eq.(9) and the local asymptoticequilibria (12), but with with T + and T − defined as in (10) and (11) with ψ ( v ) respectivelysubstituted by ψ ( v | e ) ≡ ψ + ( v ) and ψ ( v | − e ) ≡ ψ − ( v )The local asymptotic equilibria are stationary and homogeneous if and only if p + ∞ ( v ) = ρ ∞ ψ + ( v )2 , p −∞ ( v ) = ρ ∞ ψ − ( v )2 , which, on the contrary of the symmetric case, are no longer equal.Carrying out the same computations as in Section 3, it is convenient to keep Eq.(23) in thefollowing modified form µ (cid:90) R + (cid:20) i B ˆ γ R ( k )( σ + µ ) + ikv ψ + ( v ) + 1 − i B ˆ γ R ( k )( σ + µ ) − ikv ψ − ( v ) (cid:21) dv = 1 , (45)13 a) localized kernel (b) localized kernel(c) uniform kernel (d) uniform kernel Figure 4: Evolution of the density distribution in the adhesion case starting from the initialcondition ρ ( x ) = 0 . . πx/ ρ ∞ ≈ . V b correspond to stable cases,while in the right column V correspond to the unstable case. Specifically, in (a) V b ≈ . V b ≈ . V b ≈ . V b ≈ . a) localized lernel (b) uniform kernel(c) decreasing kernel Figure 5: Comparison of unstable evolutions for the same value of V = 0 . V =0 . , R = 0 . , µ = 100) starting form the initial condition ρ ( x ) = 0 . . πx/ ρ ∞ ≈ . γ R ( λ ) = δ ( λ − R ), in (b) γ R ( λ ) = H ( R − λ ) and in (c) γ R = (cid:0) − λR (cid:1) + .15eing ˆ γ R ( k ) defined as in (18).In order to understand the stability behaviour we will consider ψ ± ( v ) = δ ( v − V ± ). This allowsto integrate (45) to get µ (cid:20) σ + µ + k B ˆ γ R ( k ) V + + V − ik V + − V − (cid:21) = ( σ + µ ) + k V + V − + ik V + − V − σ + µ ) , or σ + σ (cid:2) µ + ik ( V + − V − ) (cid:3) + k V + V − − B µ ˆ γ R ( k ) k V + + V − ikµ V + − V − . The dispersion relation then reads σ = − µ − ik ( V + − V − ) + (cid:112) µ − k ( V + + V − ) + 2 B µ ˆ γ R ( k ) k ( V + + V − )2 , (46)where we notice that the speed difference ( V + − V − ) affects the imaginary part while the realpart is affected by the mean speed ( V + + V − ) /
2. Furthermore, the fact that V + (cid:54) = V − impliesthat the wave frequencies σ are always complex, hence we expect moving patterns.According to the type of kernels one then has γ R ( λ ) = δ ( λ − R ) = ⇒ B sin( kR ) kR > ¯ V ,γ R ( λ ) = H ( R − λ ) = ⇒ B − cos( kR ) k R > ¯ V ,γ R ( λ ) = (cid:0) − λR (cid:1) + = ⇒ B kR − sin( kR ) k R > ¯ V , (47)where ¯ V = V + + V − µR . (48)that reflect the conditions (32), (39), and (42) found in the symmetric case.In Figure 6 we present some tests in the asymmetric case both with perfectly reflective boundaryconditions and periodic boundary conditions in the case of a Dirac delta sensing function. In (a)-(c) we consider adhesion, i.e. b ( ρ ) = ρ in an unstable situation. In particular, Figure 6(c) showsthat the imaginary part of the eigenvalue is always stricly positive and in fact both in (a) and (b)there is a moving pattern. In (d)-(f) we have unstable configurations in the case of volume fillingFigure 6(f) shows again that as V + (cid:54) = V − there are complex eigenvalues and therefore we havea moving pattern in the direction − e that is what we expect as V + < V − . In Figure 6(g) wepresent the test with the same value of V b as in Figure 6(d)-(e) but with V + = V − . As V + = V − in Figure 6(g), the pattern is symmetric and there are real eigenvalues (see Figure (i)). In Figure6(g) the boundary conditions are periodic, but, as V + = V − , we would obtain the same resultwith specular reflection boundary conditions. In Figure 6(h) we instead consider adhesion and wehave a stable configuration as V b >
1, even if we have an anisotropic setting as V + (cid:54) = V − : thetransient shows an asymmetric behavior, but the solution goes to a stationary homogeneous case. We analyzed the stability properties of a non-local kinetic equation implementing a velocity jumpprocess in which the transition probability models a directional response to a non-local evaluationof the macroscopic cell density. We identified a stability condition (28) that depends on two non-dimensional parametres: B that takes into account of the directional response through a non-localmeasure of the cell density weighted by a sensing kernel b ( ρ ), and V , that takes into account of themotility properties of cells, specifically their mean speed, sensing radius and tumbling frequency.We proved that if B is negative, corresponding to a response that tends to avoid crowding, then16 a) Perfect relfection BCs (b) Periodic BCs (c)(d) Perfect relfection BCs (e) Periodic BCs (f)(g) Periodic BCs (h) Periodic BCs (i) Figure 6: Evolution in the asymmetric case for a localised sensing kernel. The initial conditionis always ρ ( x ) = 0 . .
01 sin( πx/ ρ ∞ ≈ . µ = 3 , R =0 . , V + = 1 , V − = 0 . , V ≈ b . µ = 200 , R = 0 . , V + = 0 . , V − =0 . , ρ th = 1 , V b ≈ − . µ = 200 , R = 0 . , V + = V − = 0 . , ρ th = 1 , V b ≈− . µ = 1 , R = 0 . , V + = 0 . , V − = 1 , V b ≈ . V / |B| smallenough. On the other hand, if B is positive, corresponding to an adhesion-like behavior, thenthe homogeneous solution is unstable to long wave for values of V / B that are sufficiently small.Considering that V = VµR , this means that for fixed B instability occurs in stiff regimes, e.g. alarge tumbling frequency leads to instability.On the other hand, if V is fixed, the stiffness of the response b ( ρ ∞ ) determines the transitionfrom stability to instability. In the case of increasing b at ρ ∞ , cells tend to go towards zones thatare more crowded. If b (cid:48) ( ρ ∞ ) is sufficiently high, there is linear instability, for all the analyzedsensing functions. If b (cid:48) ( ρ ∞ ) is negative, cells tend to go towards regions that are less crowded. Inthis case instability can be triggered for example in the case in which the sensing function is aDirac delta, i . e . when cells do not evaluate the density in between their position and the sensingposition x ± R .Numerical simulations show that the linear stability analysis predicts pattern formation (orstability) quite sharply and it is able to catch the characteristic wavelengths of the pattern.We restricted to the case in which cell density only affects the direction of motion, e.g., cellsturn away if they sense an overcrowded area. Of course, cell density can also affect their speed.This case is currently under study. Actually, the most interesting case is the one in which celldensity has an ambivalent effect. For instance, cells are attracted by the other cells because ofcell-cell adhesion but at the same time they want to stay away from overcrowding and they do notwant to overcome a certain threshold density (volume filling). Acknowledgements
This work was partially supported by Istituto Nazionale di Alta Matematica, Ministry of Educa-tion, Universities and Research, through the MIUR grant Dipartimenti di Eccellenza 2018-2022,Project no. E11G18000350001, and the Scientific Reseach Programmes of Relevant National In-terest project n. 2017KL4EF3. NL also acknowledges Compagnia di San Paolo that funds herPh.D. scholarship.
References
Armstrong, N. J., Painter, K. J., and Sherratt, J. A. (2006). A continuum approach to modellingcell-cell adhesion.
Journal of Theoretical Biology , 243(1):98–113.Bitsouni, V. and Eftimie, R. (2018). Non-local parabolic and hyperbolic models for cell polarisationin heterogeneous cancer cell populations.
Bulletin of Mathematical Biology , 80(10):2600–2632.Buttensch¨on, A. (2018).
Integro-partial differential equation models for cell-cell adhesion and itsapplication . PhD thesis, University of Alberta.Buttensch¨on, A. and Hillen, T. (2019). Non-local adhesion models for microorganisms on boundeddomains. arXiv:1903.06635.Buttensch¨on, A. and Hillen, T. (2020). Non-local cell adhesion models: Steady states and bifur-cations. arXiv:2001.00286.Buttensch¨on, A., Hillen, T., Gerisch, A., and Painter, K. J. (2018). A space-jump derivationfor non-local models of cell-cell adhesion and non-local chemotaxis.
Journal of MathematicalBiology , 76(1):429–456.Carrillo, J., Hoffmann, F., and Eftimie, R. (2015). Non-local kinetic and macroscopic models forself-organised animal aggregations.
Kinetic & Related Models , 8:413–441.Colombi, A., Scianna, M., and Preziosi, L. (2017). Coherent modelling switch between pointwiseand distributed representations of cell aggregates.
Journal of Mathematical Biology , 74(4):783–808. 18olombi, A., Scianna, M., and Tosin, A. (2015). Differentiated cell behavior: a multiscale approachusing measure theory.
Journal of Mathematical Biology , 71:1049–1079.Eftimie, R. (2012). Hyperbolic and kinetic models for self-organized biological aggregations andmovement: a brief review.
Journal of Mathematical Biology , 65(1):35–75.Eftimie, R., de Vries, G., A Lewis, M., and Lutscher, F. (2007a). Modeling group formation andactivity patterns in self-organizing collectives of individuals.
Bulletin of Mathematical Biology ,69:1537–65.Eftimie, R., Perez, M., and Buono, P.-L. (2017). Pattern formation in a nonlocal mathematicalmodel for the multiple roles of the tgf- β pathway in tumour dynamics. Mathematical Biosciences ,289:96 – 115.Eftimie, R., Vries, G., and Lewis, M. (2007b). Complex spatial group patterns result from differentanimal communication mechanisms.
Proceedings of the National Academy of Sciences of theUnited States of America , 104:6974–9.Filbet, F. and Vauchelet, N. (2010). Numerical simulation of a kinetic model for chemotaxis.
Kinetic and Related Models , 3:B348–B366.Hillen, T., Painter, K. J., and Schmeiser, C. (2007). Global existence for chemotaxis with finitesampling radius.
Discrete & Continuous Dynamical Systems - B , 7(1):125–144.Lemou, M. and Mieussens, L. (2008). A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit.
SIAM Journal on ScientificComputing , 31(1):334–368.Loy, N. and Preziosi, L. (2019a). Kinetic models with non-local sensing determining cell polar-ization and speed according to independent cues.
Journal of Mathematical Biology-accepted(arXiv:1906.11039v3) .Loy, N. and Preziosi, L. (2019b). Modelling physical limits of migration by a kinetic model withnon-local sensing. Preprint ( arXiv:1908.08325v1 ).Othmer, H. and Hillen, T. (2002). The diffusion limit of transport equations ii: Chemotaxisequations.
SIAM Journal of Applied Mathematics , 62:1222–1250.Othmer, H. G., Dunbar, S. R., and Alt, W. (1988). Models of dispersal in biological systems.
Journal of Mathematical Biology , 26(3):263–298.Painter, J. K. and Hillen, T. (2002). Volume-filling and quorum-sensing in models for chemosen-sitive movement.
Canadian Applied Mathematics Quarterly , 10:501–543.Painter, K. J., Armstrong, N. J., and Sherratt, J. A. (2010). The impact of adhesion on cellularinvasion processes in cancer and development.
Journal of Theoretical Biology , 264(3):1057–1067.Painter, K. J., Bloomfield, M. J., Sherratt, J. A., and Gerisch, A. (2015). A nonlocal model forcontact attraction and repulsion in heterogeneous cell populations.
Bulletin of MathematicalBiology , 77:1132–1165.Perthame, B. and Yasuda, S. (2018). Stiff-response-induced instability for chemotactic bacteriaand flux-limited keller–segel equation.
Nonlinearity , 31(9):4065–4089.Plaza, R. G. (2019). Derivation of a bacterial nutrient-taxis system with doubly degenerate cross-diffusion as the parabolic limit of a velocity-jump process.
Journal of Mathematical Biology ,78:1681–1711.Schmeiser, C. and Nouri, A. (2017). Aggregated steady states of a kinetic model for chemotaxis.
Kinetic and Related Models , 10(1):313 – 327. 19troock, D. W. (1974). Some stochastic processes which arise from a model of the motion of abacterium.
Zeitschrift f¨ur Wahrscheinlichkeitstheorie und Verwandte Gebiete , 28(4):305–315.Tosin, A. and Frasca, P. (2011). Existence and approximation of probability measure solutions tomodels of collective behaviors.