Kinetic models with non-local sensing determining cell polarization and speed according to independent cues
KKinetic models with non-local sensing determining cellpolarization and speed according to independent cues
Nadia Loy ∗ Luigi Preziosi † Abstract
Cells move by run and tumble, a kind of dynamics in which the cell alternates runs overstraight lines and re-orientations. This erratic motion may be influenced by external factors,like chemicals, nutrients, the extra-cellular matrix, in the sense that the cell measures theexternal field and elaborates the signal eventually adapting its dynamics.We propose a kinetic transport equation implementing a velocity-jump process in which thetransition probability takes into account a double bias, which acts, respectively, on the choiceof the direction of motion and of the speed. The double bias depends on two different non-local sensing cues coming from the external environment. We analyze how the size of thecell and the way of sensing the environment with respect to the variation of the externalfields affect the cell population dynamics by recovering an appropriate macroscopic limit anddirectly integrating the kinetic transport equation. A comparison between the solutions ofthe transport equation and of the proper macroscopic limit is also performed.
It is well known that a classical migration mode of bacteria consists in alternating runs or swimsover straigth lines and tumbles (Berg, 1983). Also eukariotic cells alternate persistent crawlingsalong a polarization axis with reorientation phases in which they re-polarize as a result of severalexternal stimuli. These can be represented by nutrients, chemical factors, or noxious substances(Block et al., 1983), and also by environmental cues, such as density, stiffness, and structure of theextra-cellular matrix (ECM). In addition, the presence of other cells along the way can act in a twofold-way, either attractive due to the interaction of adhesion molecules (e.g, cadherin complexes)expressed on the cellular membrane, or repulsive when the region is becoming too overcrowded.Cells measure all these signals by transmembrane receptors located on their protrusions thatcan extend up to several cell diameters. The captured chemical or physical signals activate in turntransduction pathways downstream that lead to the cell response. This includes in particular thepolarization of the cell with the formation of a “head” and a “tail” and the activation of adhesionmolecules and traction forces leading to motion (Adler, 1966). In the framework of kinetic models,in the present paper we will focus on how the environmental sensing over a finite radius can betranslated in the choice of the direction of motion and of the speed of the cell. In fact, in order todo that we develop a kinetic model in which the distribution function depends on cell speed andorientation, respectively a scalar and a unit vector, in addition to the usual dependence on spaceand time. In the tumbling phase, then, the new orientation and subsequent speed are chosen asa result of a sensing over a finite neighbourhood of the cell, giving the kinetic model a non-localcharacter.Focusing on mesoscopic models of cell migration and referring to Hillen and Painter (2008) foran extensive and more general review on PDE models of chemotaxis, we can observe that the run ∗ 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] ) a r X i v : . [ q - b i o . CB ] J u l nd tumble dynamics can be modeled by a stochastic process called velocity-jump process thatwas introduced by Stroock (1974). In this work, the author derived a linear transport equationfrom the velocity-jump process. Such an equation describes the movement of a single-particledistribution function like in the Boltzmann equation (Cercignani, 1987). The main elements ofsuch a process are the tumbling frequency, the mean speed, and the transition probability thatdescribes the probability of choosing a certain velocity after re-orientation. The mean speed, meanruntime (which is the inverse of the frequency), and the tumbling probability may be measuredfrom individual patterns of members of the population. Another advantage of such models is thatthey admit finite propagation speed. Alt (1980) and Othmer et al. (1988) introduced the biasinduced by an external stimuli in the linear transport equation. Dickinson (2000) generalized themodel to a transport equation in an anisotropic environment and applied to an environment witha gradient of a stimulus.Hillen (2006) proposed a Boltzmann-like model taking into account of the interaction betweencells and ECM, also including the degradation of ECM fibers as a function of the angle between theECM fiber and the cell velocity. In particular, the author recovered macroscopic limits (diffusiveor hyperbolic) according to the structure of the matrix. Painter (2008) used a similar model todescribe ECM remodelling resulting from the migration of cells in a heterogeneous environmentby integrating numerically the kinetic equation.An extension of this Boltzmann-like model taking also into account of the interaction betweencells and between cells and ECM is proposed by Chauviere et al. (2007a) and Chauviere et al.(2007b). A further extension is suggested by Chauviere and Preziosi (2010) where the operatorsdescribing interactions between cells and between cells and ECM include the bias due to an externalcue. Surulescu and coworkers (Engwer et al., 2017, 2015a,b, 2016; Stinner et al., 2015) insteaddiscussed in more details cell-ECM adhesion phenomena, and, focusing on tumour growth, mainlygliomas, also included proliferation and therapies.Many of the papers above determine for the proposed mesoscopic models the related macro-scopic, usually diffusive, limits and implement them numerically. The diffusive limit of transportequations is deeply described by Hillen and Othmer (2000). More in detail, Othmer and Hillen(2002) investigated various forms and orders of the bias which may be included in the transitionprobability. However, in some regimes a hyperbolic limit is more suitable (for example in presenceof formation of networks, as studied by Chalub et al. (2004), Hwang et al. (2005), Filbet et al.(2005)).Coming to the sensing strategies and then to the determination of cell re-orientation andspeed, they have been considered in several models. In position jump processes, the transitionprobabilities model the strategy of sensing and may typically be (i) local, if the information at thepresent position is considered, (ii) neighbour based, if the information at the target jump site istaken into account, (iii) gradient based, if the local difference in information between the targetand the local site is considered, (iv) local average, if the average of the information between thepresent and target site as done by Othmer and Stevens (2001) and by Painter and Sherratt (2003).In general, a way of including cell sensing is to consider a non-local average of the external fields.For what concerns bacteria and cells dynamics, Othmer and Hillen (2002) and Hillen et al. (2007)introduced a finite sampling radius and they defined a non-local gradient as the average of theexternal field on a surface which represents the membrane of the cell. In particular, the Authorsderived a macroscopic diffusive model with a non-local gradient both from a position jump processand from a velocity jump process by postulating that the non-local sensing is a bias of higher order.The notion of non-local sensing was also used for cell adhesion and haptotaxis by Armstronget al. (2006) yielding a macroscopic integro-differential equation. This model was recently derivedfrom a space jump process by Buttensch¨on et al. (2018). Similar equations were proposed in 2Dset-ups by Colombi et al. (2017), Colombi et al. (2015), and applied also to model crowd dynamicsand traffic flow by Tosin and Frasca (2011).Eftimie (2012) and Carrillo et al. (2015) proposed non-local kinetic models including repulsion,alignement and attraction. They distinguish cells going along the two directions in the one-dimensional set-up with the possibility of switching between the two directions with a constantspeed. The one-dimensional kinetic model takes then a discrete structure that is then integrated2irectly.As already stated, in this paper we want to include the sensing strategy in a velocity-jumpprocess and to model the tactic and kinetic response of a cell as a consequence of the non-localsensing of the environment. In particular, we focus on the fact that different fields can influencerespectively cell direction and speed. Because of these characteristics we name the model as adouble-bias non-local kinetic model. To be specific, a chemical factor may influence the orientationof a cell (taxis), while the cell population density or the matrix density may influence the speed ofthe cell once the direction is chosen. For instance, a volume filling effect can hamper the motion ina certain direction due to cell overcrowding. Similarly, fiber networks with too narrow pores canhamper or even completely inhibit migration, a phenomenon called physical limit of migration byte Boekhorst et al. (2016) and Wolf et al. (2013).Therefore, cell response is related to the choice of the transition probability evaluated in anon-local way over the neighborhood of the cell. Several examples are given, some in order toshow the link with existing models (Buttensch¨on et al., 2018; Hillen et al., 2007; Painter et al.,2010, 2015; Painter and Hillen, 2002), others to point out the novelties of the proposed model.In particular, from the numerical point of view, we simulate directly the kinetic equation takinginto account of different phenomena, e . g . , chemotaxis, ECM steryc hindrance, adhesion inducedaggregation, durotaxis.We also discuss how the choices made on the transition probability and on the size of thesampling volume determine the type of macroscopic limit that can be performed. Then, wecompare the outcome of the numerical integration of the transport equation with the propermacroscopic limits and the related analytical solution that can be obtained for the stationarycase. This shows a satisfactory agreement when the hypothesis necessary to perform the limitsare satisfied.The plan of the article is as follows. After a section introducing the general and preliminaryconcepts, Section 3 introduces the general structure of the double-bias non-local kinetic model.Then, in Section 4 the macroscopic limits are discussed, distinguishing when it is possible toperform a parabolic scaling or a hyperbolic scaling according to the modelling assumptions. InSection 5 we start with a simpler single-bias model in which there is no cue affecting the speed,that is determined through a given distribution function. Several sensing strategies determiningcell orientation are discussed, starting from local responses to non-local averages and strategiescomparing the level of chemicals in different points of the neighbourhood, also including thepossible dependence of the encounter rate from taxis factors. In Section 6 the double-bias modelis discussed. Again, for clarity we first focus on the simpler case in which there is no cue influencingcell orientation, so that it is randomly chosen. We finally give an example for the complete model,already presented in its general form in Section 3. The statistical description of the cell population is performed through the distribution density p = p ( t, x , v p ) which is parametrized by the time t >
0, the position x ∈ Ω ⊆ R d and thevelocity pair v p = ( v, ˆ v ) ∈ V p = [0 , U ] × S d − , being S d − the unit sphere boundary in R d and U the maximal cell speed, so that ˆ v is the unit vector and v is the speed, i . e . , the modulus ofthe velocity vector v = v ˆ v . The choice of representing the distribution function p depending onthe direction and on the speed of the velocity, instead of the velocity vector v , lies in the needof separating the mechanisms governing cell polarization and motility, for instance in response ofchemotaxis and in presence of other influencing factors (Alt, 1980). The mesoscopic model consistsin the transport equation for the cell distribution: ∂p∂t ( t, x , v p ) + v · ∇ p ( t, x , v p ) = J [ p ]( t, x , v p ) (1)where the operator ∇ denotes the spatial gradient. The term J [ p ]( t, x , v p ), named turning opera-tor , is an integral operator that describes the change in velocity which is not due to the free-particle3ransport. It may describe the classical run and tumble behaviors, contact guidance phenomena,or cell-cell interactions. For the moment we will consider the classical run and tumble, e . g . , ran-dom re-orientations, which, however, may be biased by external cues. Therefore, our turningoperator will be the implementation of a velocity-jump process in a kinetic transport equation asintroduced by Alt (1980).A macroscopic description for the cell population can be classically recovered through thedefinition of moments of the distribution function p :- the cell number density ρ ( t, x ) = (cid:90) V p p ( t, x , v p ) d v p ; (2)- the cell mean velocity U ( t, x ) = 1 ρ ( t, x ) (cid:90) V p p ( t, x , v p ) v d v p ; (3)- the cell variance-covariance matrix P ( t, x ) = (cid:90) V p ( v − U ) ⊗ ( v − U ) p ( t, x , v p ) d v p ; (4)- the cell speed variance E ( t, x ) = (cid:90) V p | v − U | p ( t, x , v p ) d v p . (5)We remark that, because of the definition of V p , the integrals over the velocity space are (cid:90) V p f ( v, ˆ v ) d v p = (cid:90) U (cid:90) S d − f ( v, ˆ v ) d ˆ v dv. (6)In particular, if the dependence of f on v and ˆ v can be factorized, i . e . f ( v, ˆ v ) = f ( v ) f (ˆ v ), wehave that (cid:90) V p f ( v, ˆ v ) d v p = (cid:90) U f ( v ) dv (cid:90) S d − f (ˆ v ) d ˆ v . (7)The integral over S d − , which we remind to be the boundary of the unit sphere, is not a surfaceintegral. It has to be interpreted as (cid:90) S d − f (ˆ v ) d ˆ v = (cid:90) π f (cos( θ ) , sin( θ )) dθ in 2D and similarly in 3D.Computing the moments of the transport equation (1) allows to obtain evolution equationsfor the moments above. Macroscopic limits can then be achieved by classical procedures. Inparticular, the diffusive limit of transport equations with velocity jump processes is deeply treatedby Hillen and Othmer (2000), Othmer and Hillen (2002), and Hillen (2006), where the hyperbolicscaling is treated as well. The interest here is focused on when, according to the hypothesis doneon the transition probability and on the sensing radius, it is possible to perform either a parabolicor a hyperbolic scaling, leading to the related diffusive and hyperbolic limits. In taxis processes cells are capable of detecting and measuring external signals through membranereceptors located along cell protrusions that can extend over a finite radius. Information of bothmechanical and chemical origin are then transduced and act as control factors for the dynamicsof cells. Therefore, in the turning operator of the kinetic model we will include the evaluation ofmean fields in a neighborhood of the re-orientation position. In particular, we will consider two4ifferent control factors S and S (cid:48) which are positive quantities defined on Ω. They respectivelybias cell polarization (and therefore orientation) and speed, once cell orientation is set. Hence, wewill have a transition probability and a frequency of tumbling which depend on both S and S (cid:48) .The general form of the turning operator which implements a velocity jump processes is J [ p ]( x , v p ) = G [ p ]( x , v p ) − L [ p ]( x , v p ) , (8)where the gain term is G [ p ]( x , v p ) = (cid:90) V p µ ( x , v (cid:48) p ) T [ S , S (cid:48) ]( x , v p | v (cid:48) p ) p ( t, x , v (cid:48) p ) d v (cid:48) p , (9)the loss term is L [ p ]( x , v p ) = (cid:90) V p µ ( x , v p ) T [ S , S (cid:48) ]( x , v p (cid:48)(cid:48) | v p ) p ( t, x , v p ) d v p (cid:48)(cid:48) , (10)and v (cid:48) p is the pre-turning velocity of the gain term and v p (cid:48)(cid:48) is the post-turning velocity of the lossterm. The so-called turning kernel T [ S , S (cid:48) ]( x , v p | v (cid:48) p ) is the probability for a cell in x of choosingthe velocity v p after a re-orientation biased by the external fields S and S (cid:48) given the pre-turningvelocity v (cid:48) p . Being a transition probability, it satisfies (cid:90) V p T [ S , S (cid:48) ]( x , v p (cid:48)(cid:48) | v p ) d v p (cid:48)(cid:48) = 1 , ∀ x ∈ Ω , ∀ v p ∈ V p . (11)Also the turning frequency µ may depend on the external signal and/or on its gradient, and on themicroscopic orientation, e . g . through its polarization, typically through v · ∇S (see for examplethe work by Othmer and Hillen (2002)). In fact, it is proved that the speed and the turning ratesof individuals depend not only on the magnitude of an external signal S but also on its temporaland spatial variations as highlighted by Berg (1983), Fisher et al. (1989), Koshland (1981), Solland Wessels (1998).As done by Chauviere et al. (2007a) and Chauviere et al. (2007b), in the following, wewill assume that cells retain no memory of their velocity prior to the re-orientation, i . e . , T = T [ S , S (cid:48) ]( x , v p ). The independence from the pre-tumbling velocity lies in the fact that the choiceof the new velocity is linked to the slow interaction process also related to cell ruffling and sensingwhich is responsible for the biased re-orientation. However, the assumption might be restrictivein some cases, as it does not include, for instance, the case in which the sensing region depends onthe incoming velocity through a polarization-dependent expression of transmembrane proteins. Italso excludes persistence effects in which the re-orientation direction depends on the pre-tumblingpolarization of the cell. Under such an assumption, the operator (8) simplifies to J [ p ]( x , v p ) = (cid:90) V p µ ( x , v (cid:48) p ) p ( t, x , v (cid:48) p ) d v (cid:48) p T [ S , S (cid:48) ]( x , v p ) − µ ( x , v p ) p ( t, x , v p ) , (12)where (cid:82) V p µ ( x , v (cid:48) p ) p ( t, x , v (cid:48) p ) d v (cid:48) p represents the fraction of re-orientating cells per time unit in x .If also the frequency µ ( x ) does not depend on the microscopic velocity, then we have that J [ p ]( t, x , v p ) = µ ( x ) (cid:16) ρ ( t, x ) T [ S , S (cid:48) ]( x , v p ) − p ( t, x , v p ) (cid:17) . (13)However, in the following, we will also consider a ’weak’ dependence of the frequency on theincoming velocity and in particular on cell polarization.We observe that the distribution function nullifying the operator in (13) is simply p ( t, x , v p ) = ρ ( t, x ) T [ S , S (cid:48) ]( x , v p ) , (14)which, in general, is time dependent, but it may also be a stationary equilibrium state.5s a consequence of (11), for any function T [ S , S (cid:48) ] the operator J [ p ] satisfies the property (cid:90) V p J [ p ]( t, x , v p ) d v p = 0 , (15)that expresses mass conservation. On the other hand, the evaluation of the first moment of J [ p ]with respect to the velocity variable yields (cid:90) V p J [ p ]( t, x , v p ) v d v p = µ ( x ) ρ ( t, x ) (cid:0) U S , S(cid:48) ( x ) − U ( t, x ) (cid:1) (16)that indicates no preservation of momentum in general. Similarly, also energy is not preserved.In (16) the vector U is the mean velocity of the cells defined by (3), while U S , S(cid:48) is a macroscopicvelocity due to the variations of S , S (cid:48) , defined by U S , S(cid:48) ( x ) = (cid:90) V p T [ S , S (cid:48) ]( x , v p ) v d v p . (17)It is the mean outgoing velocity after a re-orientation biased by S and S (cid:48) . Another importantquantity is the variance-covariance matrix of T , which is defined as D S , S(cid:48) ( t, x ) = (cid:90) V p ( v − U S , S(cid:48) ) ⊗ ( v − U S , S(cid:48) ) T [ S , S (cid:48) ]( x , v p ) d v p . (18)Bisi et al. (2008) and Pettersson (2004) deal with an operator which has a structure similar to(13). They prove that, provided that the probability distribution has initially a finite mass andenergy and non-absorbing boundary conditions hold, the function (14), which makes (13) vanish, isa stable asymptotic equilibrium state. It is stationary, which means that ρ is stationary. We shallconsider no-flux boundary conditions ( i . e . (59)) which are non-absorbing boundary conditions.Moreover, U S , S(cid:48) and D S , S(cid:48) are the mean velocity and variance-covariance tensor of the cells at theequilibrium. T From the physiological point of view, the decision process determining the new cell velocity canbe split in two parts. First the cell decides where to go and polarizes, forming the so-called cell“head” and “tail”. Then, it starts to move in that direction with a certain speed. These twoprocesses are mainly independent and might be influenced by different chemical and mechanicalcues. For instance, environment microstructure, availability of intracellular space and of adhesionsites, expression of integrins, activity of motor proteins, calcium influx, availability of ATP allinfluence cell speed, while gradients of free and bound chemical factors, of ECM components andECM stiffness all influence cell polarization. A cell can even be polarized but unable to move. This,for instance, occurs if the cell is treated with latrunculin as shown by Devreotes and Janetopoulos(2003), or because of overcrowding due to the presence of too many cells in the direction where itwould like to move, giving rise to the so-called volume filling effect (see the work by Painter andHillen (2002) and references therein), or because of a too dense ECM, giving rise to the so-calledphysical limit of migration (Wolf et al., 2013; Arduino and Preziosi, 2015; Giverso et al., 2017;Scianna and Preziosi, 2013, 2014; Scianna et al., 2013).As already stated, the choice of the new velocity is a result of a sensing activity of the cellin its neighborhood performed to monitor the quantity of two generally different fields S and S (cid:48) ,respectively determining cell polarization and speed. This information is then weighted by twofunctions γ S : R + (cid:55)−→ R + and γ (cid:48) S(cid:48) : R + (cid:55)−→ R + , both with compact support that is related to the finite size of the influencing neighborhood relatedrespectively to the fields S and S (cid:48) around the position x and weighting the information accordingto the distance from the cell center. They may be6 Dirac deltas, e . g . , γ S ( λ ) = δ ( λ − R S ), if the cells only measure the information perceived ona surface of given radius R S , • Heaviside functions, e . g . , γ S ( λ ) = H ( R S − λ ) if the cells explore the whole volume of thesphere centered in x with radius R S and weight the information uniformly, or • decreasing functions of the distance λ from the cell center x taking for instance into accountthat the probability of making longer protrusion decreases with the distance, so the sensingof closer regions is more accurate.The distance R S is, therefore, a measure of the capability of sensing and detecting of the cell. Itmay be the measure of cell protrusions, or it may be bounded by Uµ which is the mean linear tracttravelled between two consecutive re-orientations. Berg and Purcell (1977) have shown that thesampling volume in which the signal is significantly detected depends on how rapidly the receptorson the cell’s membrane process the signal. The same considerations are valid for γ S(cid:48) with a radius R S(cid:48) which in general, might be different from R S .The above considerations justify the following splitting T [ S , S (cid:48) ]( x , v p ) = c ( x ) (cid:90) R + γ S ( λ ) T ˆ v λ [ S ]( x ) dλ (cid:90) R + γ S(cid:48) ( λ (cid:48) ) ψ ( x , v |S (cid:48) ( x + λ (cid:48) ˆ v )) dλ (cid:48) . (19)The quantity T ˆ v λ [ S ]( x ) is a functional which acts on S and describes the way the cell measuresthe quantity S around x along the direction ˆ v and, therefore, the bias intensity in the directionˆ v . In some cases, we shall write for instance T ˆ v λ [ S ]( x ) = b (cid:0) S ( x + λ ˆ v ) (cid:1) , (20)where b is a quantity depending on the value of S in a point at a distance λ from x along thedirection ˆ v . So, if b is an increasing function of S and the signal is stronger in the direction ˆ v ,then there will be a higher probability for the cell to move along ˆ v than along − ˆ v . We shall alsoconsider other forms for T ˆ v λ [ S ]( x ) as well.The quantity ψ ( x , v |S (cid:48) ( x + λ (cid:48) ˆ v )) is a conditional probability distribution function of the cell speedgiven the value of S (cid:48) in a neighborhood of x in the direction ˆ v that is then weighted by γ S(cid:48) . Again,we shall discuss this further in Section 6. Finally, the factor c ( x ) is a normalizing factor which isgiven by c ( x ) = 1 (cid:90) S d − Γ (cid:48) (cid:90) R + γ S ( λ ) T ˆ v λ [ S ]( x ) dλ d ˆ v (21)where Γ (cid:48) = (cid:90) R + γ S(cid:48) ( λ (cid:48) ) dλ (cid:48) , so that (11) is satisfied. In order to highlight the driving macroscopic phenomenon which may be diffusion or convection,we will discuss the appropriate macroscopic limits of the kinetic equations. A way to do that isto perform, respectively, a parabolic scaling and then a diffusive limit, or a hyperbolic scaling andthen a hydrodynamic limit. In these approaches it is assumed that there is a small parameter (cid:15) (cid:28) τ = (cid:15) t, ξ = (cid:15) x ,τ = (cid:15)t, ξ = (cid:15) x . Following the work by Othmer and Hillen (2002), we also assume that T can be expanded as T [ S , S (cid:48) ]( ξ , v p ) = T [ S , S (cid:48) ] ( ξ , v p ) + (cid:15)T [ S , S (cid:48) ] ( ξ , v p ) + O ( (cid:15) ) (22)7his means that there are different orders of bias as will be better illustrated in the examples tofollow. If we assume that µ = O (1), we may denote J and J the related operators defined by T and T respectively, and we assume that (cid:90) V p T [ S , S (cid:48) ] ( ξ , v p ) d v p = 1 , (23a)and (cid:90) V p T [ S , S (cid:48) ] i ( ξ , v p ) d v p = 0 ∀ i ≥ . (23b)Coherently, the velocity and tensor of the second moments of the transition probability may bewritten as U S , S(cid:48) = U S , S(cid:48) + (cid:15) U S , S(cid:48) + O ( (cid:15) ) (24)and D S , S(cid:48) = D S , S(cid:48) + (cid:15) D S , S(cid:48) + O ( (cid:15) ) (25)where U i S , S(cid:48) = (cid:90) V p T [ S , S (cid:48) ] i v d v p and D i S , S(cid:48) = (cid:90) V p T [ S , S (cid:48) ] i ( v − U i S , S(cid:48) ) ⊗ ( v − U i S , S(cid:48) ) d v p for i ≥
1. In both limits we consider an expansion of the distribution function p in the form p = p + (cid:15)p + O ( (cid:15) ) , (26)where the function p is the equilibrium function given by the solution to the leading order partof the turning operator J [ p ]( v p ) = 0 . (27)As J [ p ]( v p ) = µ ( ρ T [ S , S (cid:48) ] ( v p ) − p ( v p )), Eq.(27) is satisfied if and only if p ( v p ) = ρ T [ S , S (cid:48) ] ( v p ) . (28) The diffusive limit of velocity jump processes is deeply treated by Hillen and Othmer (2000)and Othmer and Hillen (2002). In particular, in (Othmer and Hillen, 2002) macroscopic modelsfor chemosensitive movement are systematically derived for different orders of (cid:15) of the turningprobability. In the same spirit as in the works above, because of the conservation of mass, we havethat • all the mass is in p , i . e . , ρ = ρ, ρ i = 0 ∀ i ≥ , (29)where ρ i = (cid:90) V p p i d v ; • (cid:90) V p p i v d v p = 0 ∀ i ≥ i . e ., T = T ( v p | v (cid:48) p ) and it is required that (cid:90) V p T ( v p | v (cid:48) p ) d v (cid:48) p = 1 . (31)Thanks to this property and to other regularity assumptions on T Hillen and Othmer (2000)show that the function identically equal to 1 is an eigenfunction of the turning operator witheigenvalue 0 that is a simple eigenvalue. In the present work there is no dependence on thepre-tumbling velocity and (31) does not hold true. Furthermore, because of the structure of theturning operator, the density p making the leading part of the turning operator vanish (28) is aneigenfunction of the velocity jump operator and it depends on the velocity. Moreover, as all thefunctions nullifying the turning operators are proportional to T up to a multiplicative constant,we have that the eigenvalue 0 is simple and the subspace of L ( V p ) defined by (cid:104) T [ S , S (cid:48) ] ( ξ , v p ) (cid:105) is the kernel of the turning operator and it has dimension equal to one. Therefore, the appropriatescalar product allowing to determine uniquely the solvability condition is the one proposed byHillen (2006) (cid:104) f, g (cid:105) = (cid:90) V p f gT − d v p . (32)Whence the equilibrium function p belongs to the subspace generated by the eigenfunction of J with eigenvalue 0, i . e . , (cid:104) T [ S , S (cid:48) ] ( ξ , v p ) (cid:105) . This means that the equilibrium function has meanand the variance-covariance tensor equal to U S , S(cid:48) and D S , S(cid:48) respectively. For i ≥
1, the functions p i ∈ (cid:104) T [ S , S (cid:48) ] ( ξ , v p ) (cid:105) ⊥ ⊆ L ( V p ) for a.e. ξ , where (cid:104) T [ S , S (cid:48) ] ( ξ , v p ) (cid:105) ⊥ is the orthogonal subspaceto the subspace generated by T [ S , S (cid:48) ] ( ξ , v p ). Thanks to the scalar product (32) we have that g ∈ (cid:104) T [ S , S (cid:48) ] ( ξ , v p ) (cid:105) ⊥ if and only if (cid:82) V p g ( v p ) d v p = 0. This justifies the assumptions (29).The dependence of p on v p has also consequences on the conditions for the isotropy of T andof the diffusion tensor are not the same as the ones considered by Hillen and Othmer (2000) andOthmer and Hillen (2002). The case in which T does not depend on the pre-tumbling velocity istreated by Hillen (2006) where it is assumed that the turning probability T is of order zero in (cid:15) .Therefore, we illustrate the diffusive limit procedure in the case in which the turning probabilitydoes not depend on the pre-tumbling velocity and it can be expressed as in (22), and we alsoillustrate the necessary conditions which are needed in order to perform it.The discriminating factor is whether T is such that U S , S(cid:48) = (cid:90) V p T [ S , S (cid:48) ] v d v p = for a.e. ξ , (33) i . e . , the leading order macroscopic velocity vanishes, or not. In this case we are led to perform aparabolic scaling, i . e . , τ = (cid:15) t , ξ = (cid:15) x , (cid:15) (cid:28) (cid:15) ∂p∂τ ( τ, ξ , v p ) + (cid:15) v · ∇ p ( τ, ξ , v p ) = J [ p ] + (cid:15) J [ p ] + O ( (cid:15) ) . (34)We suppose that µ ∼ O (1). From now on, to simplify the notation, we will drop the dependencieson τ and ξ . Comparing terms of equal order in (cid:15) , we then obtain the following system of equations:In (cid:15) : J [ p ]( v p ) ≡ µ (cid:16) ρ T [ S , S (cid:48) ] ( v p ) − p ( v p ) (cid:17) = 0 (35)In (cid:15) : ∇ · (cid:0) p ( v p ) v (cid:1) = J [ p ]( v p ) + J [ p ]( v p ) == µ (cid:0) ρ T [ S , S (cid:48) ] ( v p ) − p ( v p ) (cid:1) + µρ T [ S , S (cid:48) ] ( v p ) (36)9n (cid:15) : ∂∂τ p ( v p ) + ∇ · (cid:0) p ( v p ) v (cid:1) = J [ p ]( v p ) + J [ p ]( v p ) + J [ p ]( v p ) . (37)As already discussed, Eq. (35) implies (28), that is the equilibrium state of order zero. From Eq.(36), by inverting J , we get p ( v p ) = − µ ∇ · (cid:0) v p (cid:1) + ρ T [ S , S (cid:48) ] ( v p ) . (38)The solvability condition for inverting J is that p ∈ (cid:104) T [ S , S (cid:48) ] ( v p ) (cid:105) ⊥ . Because of (32) thismeans that (cid:90) V p (cid:20) −∇ · µ (cid:0) v p (cid:1) + ρ T [ S , S (cid:48) ] ( v p ) (cid:21) d v p = 0 for a.e. ξ . (39)Equation (39) is satisfied because (cid:90) V p v p ( v p ) d v p = ρ (cid:90) V p T [ S , S (cid:48) ] ( v p ) v d v p = 0 for a.e. ξ , (40)that is (33), and (cid:90) V p T [ S , S (cid:48) ] ( v p ) d v p = 0 for a.e. ξ , (41)because of (23b).The evolution equation for ρ (= ρ ) is obtained from the solvability condition at O ( (cid:15) ), i . e . , theintegral on V p of Eq. (37) which is given by ∂ρ∂τ + ∇ · (cid:16) ρ U S , S(cid:48) (cid:17) = ∇ · (cid:18) µ ∇ · (cid:0) D S , S(cid:48) ρ (cid:1)(cid:19) (42)where we used (30). The term U S , S(cid:48) is the chemotactic velocity, which is present if there is a biasof order (cid:15) as in the work by Othmer and Hillen (2002). In the following sections, we will show thatunder appropriate conditions of sensing of the cells, microscopic rules may allow to recover thechemotaxis equation in which the chemotactic velocity is proportional to the gradient of S andthe chemotactic sensitivity depends on the sensitivity function, which captures the cell capabilityof sensing and detecting the signal. We now generalize the above procedure in the case it is possible to write µ ( x , v p ) = µ ( x ) + (cid:15)µ ( x , v p ) + O ( (cid:15) ). In this case the rescaled transport equation reads (cid:15) ∂p∂τ ( v p ) + (cid:15) v · ∇ p ( v p ) = µ (cid:16) ρ (cid:0) T [ S , S (cid:48) ] ( v p ) + (cid:15)T [ S , S (cid:48) ] ( v p ) (cid:1) − p ( v p ) (cid:17) + (cid:90) V p (cid:15)µ ( v (cid:48) p ) p ( v (cid:48) p ) d v (cid:48) p (cid:16) T [ S , S (cid:48) ] ( v p ) + (cid:15)T [ S , S (cid:48) ] ( v p ) (cid:17) − (cid:15)µ ( v p ) p ( v p )+ µ (cid:15) ρT [ S , S (cid:48) ] ( v p ) + (cid:15) T [ S , S (cid:48) ] ( v p ) (cid:90) V p µ ( v (cid:48) p ) p ( v (cid:48) p ) d v (cid:48) p − (cid:15) µ ( v p ) p ( v p ) + O ( (cid:15) ) , (43)where also p has to be intended as expanded in terms of (cid:15) . Comparing terms of equal order in (cid:15) ,we obtain the same equation as before in (cid:15) , while at the first order ∇ · (cid:0) p ( v p ) v (cid:1) = µ (cid:0) ρ T [ S , S (cid:48) ] ( v p ) − p ( v p ) (cid:1) + µ ρ T [ S , S (cid:48) ] ( v p )+ (cid:90) V p µ ( v (cid:48) p ) p ( v (cid:48) p ) d v (cid:48) p T [ S , S (cid:48) ] ( v p ) − µ ( v p ) p ( v p ) , (44)10nd at second order ∂∂τ p ( v p ) + ∇ · (cid:0) p ( v p ) v (cid:1) = µ (cid:0) ρ T [ S , S (cid:48) ] ( v p ) − p ( v p ) (cid:1) + µ ρ T [ S , S (cid:48) ] ( v p ) + (cid:90) V p µ ( v (cid:48) p ) p ( v (cid:48) p ) d v (cid:48) p T [ S , S (cid:48) ] ( v p )+ (cid:90) V p µ ( v (cid:48) p ) p ( v (cid:48) p ) d v (cid:48) p T [ S , S (cid:48) ] ( v p ) − µ ( v p ) p ( v p ) + µ ρ T [ S , S (cid:48) ] ( v p )+ (cid:90) V p µ ( v (cid:48) p ) p ( v (cid:48) p ) d v (cid:48) p T [ S , S (cid:48) ] ( v p ) − µ ( v p ) p ( v p ) . (45)From equation (44), we get p ( v p ) = − µ ∇ · (cid:0) v p ( v p ) (cid:1) + ρ T [ S , S (cid:48) ] ( v p ) + ρ µ T [ S , S (cid:48) ] ( v p ) (cid:0) ¯ µ − µ ( v p ) (cid:1) , (46)where ¯ µ = (cid:90) V p µ ( v (cid:48) p ) T [ S , S (cid:48) ] ( v (cid:48) p ) d v (cid:48) p . (47)In this case the solvability condition for inverting J becomes (cid:90) V p (cid:20) − µ ∇ · (cid:0) v p (cid:1) + ρ T [ S , S (cid:48) ] ( v p ) + ρ µ T [ S , S (cid:48) ] ( v p ) (cid:0) ¯ µ − µ ( v p ) (cid:1)(cid:21) d v p = 0 , (48)that is again satisfied because (cid:90) V p v p ( v p ) d v p = 0 for a.e. ξ , (49) (cid:90) V p T [ S , S (cid:48) ] ( v p ) d v p = 0 for a.e. ξ (50)and (cid:90) V p T [ S , S (cid:48) ] ( v p ) (cid:0) ¯ µ − µ ( v p ) (cid:1) d v p = 0 for a.e. ξ . (51)The first two equations vanish as proved in the preceding section, while recalling the definition(47) and (23a), we have that (51) is trivially satisfied.As in the previous section, the solvability condition at O ( (cid:15) ), i . e . , the integral on V p of (45) givesthe evolution equation for ρ∂ρ∂τ + ∇ · (cid:16) ρ (cid:0) U S , S(cid:48) + U ,µ S , S(cid:48) (cid:1)(cid:17) = ∇ · (cid:18) µ ∇ · (cid:0) D S , S(cid:48) ρ (cid:1)(cid:19) (52)where U S , S(cid:48) = (cid:90) V p T [ S , S (cid:48) ] v d v p and U ,µ S , S(cid:48) = 1 µ (cid:90) V p { T [ S , S (cid:48) ] ( v p )[¯ µ − µ ( v p )] } v d v p = − µ (cid:90) V p T [ S , S (cid:48) ] ( v p ) µ ( v p ) v d v p , (53)where the first term vanished because of (33). So, the chemotactic velocity presents a correctionto the term U S , S(cid:48) given by U ,µ S , S(cid:48) . In the particular case in which µ is independent of v p weobviously recover (42) as U ,µ S , S(cid:48) = . 11 .3 The hyperbolic limit The discriminating factor on whether a diffusive or a hyperbolic limit can be performed lies in theproperties of the chosen turning probability T and specifically on the satisfaction of (33). So, atvariance with the previous section, we here assume that T does not satisfy (33). In this case, ahyperbolic scaling, i . e . τ = (cid:15), ξ = (cid:15) x , (cid:15) (cid:28) (cid:15) ∂p∂τ ( v p ) + (cid:15) v · ∇ p ( v p ) = µ (cid:2) ρ (cid:0) T [ S , S (cid:48) ] ( v p ) + (cid:15)T [ S , S (cid:48) ] ( v p ) (cid:1) − p ( v p ) (cid:3) . (54)As the equilibrium state is the same as before, we consider a Chapman-Enskog expansion of p inthe form p ( v p ) = ρ T [ S , S (cid:48) ] ( v p ) + (cid:15)g + O ( (cid:15) ) . (55)where g ∈ (cid:104) T [ S , S (cid:48) ] ( v p ) (cid:105) ⊥ ⊆ L ( V p ) with the same scalar product as before. Substituting (55) in(54) and integrating on V p the equation at the order (cid:15) , we obtain ∂ρ∂τ + ∇ · (cid:0) ρ U S , S(cid:48) (cid:1) = 0 , (56)thanks to (29), where (with the dependencies) U S , S(cid:48) ( ξ ) = (cid:90) V p T [ S , S (cid:48) ] ( ξ , v p ) v d v p . (57)In this case we have then a drift-driven phenomenon.If we consider a frequency weakly depending on the velocity, the rescaled transport equationmodifies into (cid:15) ∂p∂τ ( v p ) + (cid:15) v · ∇ p ( v p ) = µ (cid:2) ρ (cid:0) T [ S , S (cid:48) ] ( v p ) + (cid:15)T [ S , S (cid:48) ] ( v p ) (cid:1) − p ( v p ) (cid:3) + (cid:15)T [ S , S (cid:48) ] ( v p ) (cid:90) V p µ ( v (cid:48) p ) p ( v (cid:48) p ) d v (cid:48) p − µ ( v p ) p ( v p ) + O ( (cid:15) ) , (58)and, thanks to (29) and (47), we recover, again, (56) with U S , S(cid:48) given by (57) .
Isotropy and anisotropy
As we said at the beginning of the section, the necessary and sufficientconditions for the isotropy of the diffusion tensor are not the same as the ones considered by Hillenand Othmer (2000), as the transition probability does not depend on the incoming velocity. Thisimplies that the mean outgoing velocity for a given incoming velocity is always zero and, therefore,the adjoint persistence is zero. In our case, we have that D S , S(cid:48) is isotropic if and only if T [ S , S (cid:48) ]is isotropic as a function of ˆ v , i . e . , if T [ S , S (cid:48) ] is invariant with respect to rotations of π as afunction of ˆ v . This also implies that U S , S(cid:48) is zero and that the directional persistence is zero.Viceversa, if U S , S(cid:48) is not zero, D S , S(cid:48) is anisotropic.
We shall consider the following biological no-flux condition (Plaza, 2019) (cid:90) V p p ( τ, ξ , v p ) v · n ( ξ ) d v p = 0 , ∀ ξ ∈ ∂ Ω , τ > n ( ξ ) the outer normal to the boundary ∂ Ω in the point ξ , and we remind that v p = ( v, ˆ v ).This choice is motivated by the fact that in experiments in vitro there is no exchange of biologicalmaterial between the system and its environment. Equation (59) implies that there is no normal12iological aspect Section J Macrolimit
Application Figure
Orientation
Local sensing 5.2.1 (80) (83)
Encounter rate dependingon the incoming velocity
Speed
Random polarization 6.1.1 (130) (131) ECM steryc hindrance 9(Non-local speed sensing)
Orientation and Speed
Non-local orientation 6.1.3 (133) (134) Chemotaxis and 11and speed sensing ECM steryc hindranceTable 1: Summary of modelsmass flux across the boundary (Lemou and Mieussens, 2008). At the macroscopic level (59) gives(Plaza, 2019) (cid:16) D S , S(cid:48) ∇ ρ − ρ U S , S(cid:48) (cid:17) · n ( ξ ) = 0 , on ∂ Ω . for the diffusive limit, whilst for the hyperbolic limit the corresponding boundary condition is U ( ξ ) · n = 0 , on ∂ Ω . Two important classes of kinetic boundary conditions which satisfy (59) are the regular reflectionboundary operators and the non-local (in velocity) boundary operators of diffusive types definedby Palcewski (1992) (see also the work by Lods (2005)). Let us denote the boundary operator R [ p ]( τ, ξ , v, ˆ v ) = p ( τ, ξ , v (cid:48) , ˆ v (cid:48) ) . Two main examples of the regular reflection boundary conditions are • bounce back reflection condition p ( τ, ξ , v (cid:48) , ˆ v (cid:48) ) = p ( τ, ξ , v, − ˆ v ) if ˆ v · n ≥ p ( τ, ξ , v (cid:48) , ˆ v (cid:48) ) = p ( τ, ξ , v, ˆ v ) elsewhere. • specular reflection boundary condition p ( τ, ξ , v (cid:48) , ˆ v (cid:48) ) = p ( τ, ξ , v, ˆ v − (ˆ v · n ) n ) if ˆ v · n ≥ p ( τ, ξ , v (cid:48) , ˆ v (cid:48) ) = p ( τ, ξ , v, ˆ v ) elsewhere.Diffusive boundary conditions are prescribed in the form p ( τ, ξ , v (cid:48) , ˆ v (cid:48) ) = (cid:90) ˆ v ∗ · n ≥ h ( ξ , v p , v ∗ p ) p ( τ, ξ , v ∗ , ˆ v ∗ ) | v ∗ · n | d v ∗ p where h = h ( ξ , v p , v ∗ p ) is a non-negative measurable function called the Gaussian equilibriumfunction satisfying (cid:90) ˆ v ∗ · n ≤ h ( ξ , v p , v ∗ p ) | v ∗ · n | d v ∗ p = 1 .
13 linear combination of a regular reflection with a diffusive boundary operator is called a Maxwell-type boundary operator and reads p ( τ, ξ , v (cid:48) , ˆ v (cid:48) ) = α ( ξ ) p ( τ, ξ , v, V (ˆ v )) + (1 − α ( ξ )) M ( ξ , v, ˆ v ) (cid:90) ˆ v ∗ · n ≥ p ( τ, ξ , v ∗ , ˆ v ∗ ) | v ∗ · n | d v ∗ p , (62)where M ( ξ , v, ˆ v ) is the Maxwellian of the wall and V (ˆ v ) = − ˆ v for the bounce back reflectioncondition and V (ˆ v ) = ˆ v − (ˆ v · n ) n for the specular reflection. We shall consider regular reflectionboundary conditions for the one-dimensional simulations and Maxwell type boundary conditionsfor the two-dimensional simulations (see section 5.3.1). In this section we will specialize the above kinetic models and macrosopic limits to several situa-tions of biological interest, that are listed in Table 1 in order to guide the reader through them.For sake of clarity, for the moment we will forget about the S (cid:48) field, assuming that once the cell ispolarized its speed is determined through a distribution function over the speed ψ ( x , v ). The casein which the distribution function depends on a signaling cue S (cid:48) will be denoted as double biasand will be treated in Section 6. Therefore, ψ ( x , v ) describes the random unbiased probability fora cell at position x of having modulus v after a random re-orientation. In this case we expect afactorized turning probability, as the distributions for the (biased) direction and for the speed areindependent.We will apply the model to several non-local sensing dynamics, comparing the results with otherexisting models. Let us introduce the quantities B [ S ]( x , ˆ v ) = (cid:90) R + γ S ( λ ) T ˆ v λ [ S ]( x ) dλ , (63)and ¯ B [ S ]( x , ˆ v ) = c ( x ) B [ S ]( x , ˆ v ) , (64)where c ( x ) = 1 (cid:90) S d − (cid:90) R + γ S ( λ ) T ˆ v λ [ S ]( x ) dλ d ˆ v is the normalization constant.In this case, the turning probability can then be factorized as T [ S ]( x , v p ) = ψ ( x , v ) ¯ B [ S ]( x , ˆ v ) . (65)We also introduce ¯ U ( x ) = (cid:90) U ψ ( x , v ) v dv , (66)that is the mean speed and D ( x ) = (cid:90) U ψ ( x , v )( v − ¯ U ) dv , (67)that is the variance of the distribution function of the speed. The final expression of the operator J [ p ] then becomes J [ p ]( t, x , v p ) = µ ( x ) (cid:16) ρ ( t, x ) ψ ( x , v ) ¯ B [ S ]( x , ˆ v ) − p ( t, x , v p ) (cid:17) . (68)The distribution which makes (68) vanish corresponds to a macroscopic velocity U S ( x ) = ¯ U ( x ) (cid:90) S d − ¯ B [ S ]( x , ˆ v )ˆ v d ˆ v (69)14nd to a diffusion tensor D S ( x ) = (cid:90) V p ψ ( x , v ) ¯ B [ S ]( x , ˆ v )( v − U S ) ⊗ ( v − U S ) d v p . (70)The expansion (22) of T reflects now on the expansion of ¯ B ¯ B [ S ]( ξ , ˆ v ) = ¯ B [ S ] ( ξ , ˆ v ) + (cid:15) ¯ B [ S ] ( ξ , ˆ v ) , with U i S ( ξ ) = ¯ U ( ξ ) (cid:90) S d − ¯ B [ S ] i ( ξ , ˆ v )ˆ v d ˆ v , (71)and D i S ( ξ ) = (cid:90) V p ψ ( x , v ) ¯ B [ S ] i ( ξ , ˆ v )( v − U i S ) ⊗ ( v − U i S ) d v p . (72)We remark that if (33) holds true ( i . e . , if U S vanishes) D S ( ξ ) = D ( ξ ) (cid:90) S d − ¯ B [ S ]( ξ , ˆ v )ˆ v ⊗ ˆ v d ˆ v . (73)The diffusive limit in the case in which the turning frequency does not depend on the microscopicvelocity, i . e . Eq. (42), now reads ∂ρ∂τ + ∇ · (cid:16) ρ U S (cid:17) = ∇ · (cid:16) µ ∇ · (cid:0) D S ρ (cid:1)(cid:17) , (74)while the diffusive limit in the case in which µ depends on the microscopic velocity, i . e . Eq. (52),now reads ∂ρ∂τ + ∇ · (cid:16) ρ (cid:0) U S + U ,µ S (cid:1)(cid:17) = ∇ · (cid:16) µ ∇ · (cid:0) D S ρ (cid:1)(cid:17) , (75)where U ,µ S = 1 µ (cid:90) V p µ ( v p ) ψ ( v ) ¯ B [ S ] (ˆ v ) v d v p , (76)while the hyperbolic limit is ∂ρ∂τ + ∇ · (cid:0) ρ U S (cid:1) = 0 . (77) The function T λ ˆ v , introduced to account for a bias in the choice of the direction, has to be chosenaccording to the specific taxis process to be considered. In particular, as the cell scouts and detectsthe signal around itself, the functional T λ ˆ v will depend on S through the quantity S ( x + λ ˆ v ) asintroduced by Othmer and Hillen (2002) and later treated by Hillen et al. (2007), with λ = R S . Forus, λ ∈ [0 , R S ], where R S accounts for the size of the sensing neighborhood and the informationis weighted through γ S ( λ ). This quantity will be considered in order to define for every x theprobability of going in direction ˆ v , which is ¯ B [ S ( x )](ˆ v ). According to the characteristics of thisprobability, we will recover a macroscopic model, depending on the value of (cid:90) S d − ¯ B [ S ] ( ξ , ˆ v )ˆ v d ˆ v (78)which discriminates between a diffusion driven or a drift driven phenomena, according to whether(33) is satisfied or not. 15 .2.1 Local sensing In the first example we consider a transition rate that depends only on the information given bythe control factor at that site, i . e ., γ S ( λ ) = δ ( λ − B [ S ] in the transitionprobability is simply ¯ B [ S ]( x , ˆ v ) = 1 | S d − | , (79)where | · | denotes the measure of the set. Then the turning operator reads J [ p ]( t, x , v p ) = µ ( x ) (cid:16) | S d − | ρ ( t, x ) ψ ( x , v ) − p ( t, x , v p ) (cid:17) . (80)This is what we expect, as the measure of S in x does not affect the choice of the direction. Inthis model, S ( x ) may only affect the frequency of turning. As (33) is satisfied and B [ S ] = 0, weperform a diffusive limit. In addition, since (cid:90) S d − ˆ v ⊗ ˆ v d ˆ v = | S d − | d I , (81) D S ( ξ ) is isotropic D S ( ξ ) = D ( ξ ) I , (82)and, if µ ∼ O (1), we have that ∂ρ∂τ ( τ, ξ ) = ∇ · (cid:16) µ ( ξ ) ∇ (cid:0) D ( ξ ) ρ ( τ, ξ ) (cid:1)(cid:17) . (83)Therefore, there is no directional bias which comes from the transition probability, coherently withthe fact that no non-local information is taken into account. In order to consider situations in which organisms are too small to perform a non-local measure-ment, like E.Coli, following the ideas by Block et al. (1983) and Othmer and Hillen (2002), weallow here the encounter rate to weakly depend on a signal S , still working with local models.For sake of simplicity, we assume a time-independent signal, or however changes over a time scalemuch larger than the free-fly time. In fact, E.Coli can measure a pathwise gradient in time, i . e . , D S Dt = ∂∂t S + v · ∇S . Block et al. (1983) show that the movement of E.Coli may be modelled witha Poisson process with turning frequency depending on the temporal gradient of S µ ( x , v p ) = µ r ( x ) exp (cid:20) − f ( S ) (cid:18) ∂∂t S + v · ∇S (cid:19)(cid:21) (84)being µ r the turning rate in absence of the external signal. The turning operator reads J [ p ]( t, x , v p ) = − µ r exp (cid:20) − f ( S ) (cid:18) ∂∂t S + v · ∇S (cid:19)(cid:21) p ( t, x , v p )+ 1 | S d − | ψ ( ξ , v ) (cid:90) V p µ r exp (cid:20) − f ( S ) (cid:18) ∂∂t S + v (cid:48) · ∇S (cid:19)(cid:21) p ( t, x , v (cid:48) p ) d v (cid:48) p . (85)For instance, one can take the function proposed by Lapidus and Schiller (1976) i . e . f ( S ) = C K D ( K D + S ) where K D is the dissociation constant for the attractant. Hence, in the spirit of thework by Othmer and Hillen (2002) by a parabolic scaling, we can expand (84) as µ ( ξ ) = µ r ( ξ ) (cid:16) − (cid:15)f ( S ) v · ∇S + O ( (cid:15) ) (cid:17) , µ ( ξ ) = µ r ( ξ ) , (86)and µ ( ξ , v p ) = − µ r ( ξ ) f ( S ) v · ∇S . (87)We then obtain U S = because B [ S ] = 0 and, recalling (76), U ,µ S = − − µ r µ r f ( S ) 1 | S d − | (cid:90) U ψ ( ξ , v ) vdv (cid:90) S d − ˆ v ⊗ ˆ v d ˆ v ∇S = ¯ U ( ξ ) f ( S ) d ∇S ( ξ )and the macroscopic equation reads ∂ρ∂τ ( τ, ξ ) + ∇ · (cid:16) ρ ( τ, ξ ) ¯ U ( ξ ) f ( S ) d ∇S ( ξ ) (cid:17) = ∇ · (cid:16) µ r ( ξ ) ∇ (cid:0) D ( ξ ) ρ ( τ, ξ ) (cid:1)(cid:17) . (88)that is a chemotaxis model also in the case of local sensing, when an organism is too small formeasuring a gradient, but it is able to measure the temporal gradient of S along its path. We now consider the case in which the choice of the new velocity depends on a suitable averageof the value of the signal S perceived through transmembrane receptors by a cell that extendingits protrusions can scout its neighborhood.We shall consider in this case T ˆ v λ [ S ]( x ) = b (cid:0) S ( x + λ ˆ v ) (cid:1) that models the fact that in order to decide its new direction, the cell measures S in x + λ ˆ v andevaluates an average of a quantity which is linked to this measure ( b ). Therefore, recalling (68)the operator for biased random motion including taxis is J [ p ]( t, x , v p ) = µ ( x ) (cid:32) ρ ( t, x ) ψ ( x , v ) c ( x ) (cid:90) R + b ( S ( x + λ ˆ v )) γ S ( λ ) dλ − p ( t, x , v p ) (cid:33) , (89)if the turning frequency does not depend on the microscopic velocity. In order to understand whatto expect from the integration of (89) that will be performed in Section 5.3, we can look at theasymptotic limit.For a general b and R S , the macroscopic velocity of order zero is U S ( ξ ) = ¯ U ( ξ ) (cid:90) S d − (cid:18)(cid:90) R + b ( S ( ξ + λ ˆ v )) γ S ( λ ) dλ (cid:19) ˆ v d ˆ v (cid:90) S d − (cid:18)(cid:90) R + b ( S ( ξ + λ ˆ v )) γ S ( λ ) dλ (cid:19) d ˆ v . (90)As, in this case, (33) is generally not satisfied, we have a non-zero drift term and a hyperboliclimit needs be performed, leading to the integro-differential equation ∂ρ∂τ + ∇ · (cid:0) ρ U S (cid:1) = 0 . (91)In particular, if b ( S ) = S , this drift term measures the dominant direction in the extracellularfactor in the neighbourhood of the cell as U S ( ξ ) = ¯ U ( ξ ) (cid:90) S d − (cid:18)(cid:90) R + S ( ξ + λ ˆ v )) γ S ( λ ) dλ (cid:19) ˆ v d ˆ v (cid:90) S d − (cid:18)(cid:90) R + S ( ξ + λ ˆ v ) γ S ( λ ) dλ (cid:19) d ˆ v . (92)17owever, it is possible to simplify the turning operator (89), if R S is much smaller than thecharacteristic length l S of variation of S , e . g . , l S = 1 / max |∇S|S . In fact, we can expand b in aTaylor series b ( S ( x + λ ˆ v )) = b ( S ( x )) + λb (cid:48) ( S ( x )) ∇S ( x ) · ˆ v + O ( λ ) . (93)Therefore, if b ( S ( x )) does not vanish, the operator may be approximated by¯ B [ S ]( x , ˆ v ) = 1 | S d − | (cid:20) b (cid:48) ( S ( x )) b ( S ( x )) ∇S ( x ) · ˆ v (cid:21) (94)where Λ = (cid:90) R + γ S ( λ ) λdλ (cid:90) R + γ S ( λ ) dλ . (95)For instance, if only the signals at a membrane having a distance R S from x is taken into account, i . e . γ S ( λ ) = δ ( λ − R S ), then Λ = R S . If instead all the signals between the cell center and itsmembrane are uniformly mediated, i . e . γ S ( λ ) = H ( R S − λ ), then Λ = R S J [ p ]( t, x , v p ) = µ ( x ) (cid:20) ρ ( t, x ) ψ ( x , v ) | S d − | (cid:16) b (cid:48) ( S ( x )) b ( S ( x )) ∇S ( x ) · ˆ v (cid:17) − p ( t, x , v p ) (cid:21) . (96)We remark that the smallness of R S and the approximation (93) localizes the non-local integralmodel (89) into (96). Under these limit assumptions it is possible to perform a parabolic scalingthat gives ¯ B [ S ]( ξ , ˆ v ) = 1 | S d − | (cid:18) (cid:15) Λ b (cid:48) ( S ( ξ )) b ( S ( ξ )) ∇S ( ξ ) · ˆ v (cid:19) . (97)Hence, ¯ B [ S ] ( ξ , ˆ v ) = 1 | S d − | and ¯ B [ S ] ( ξ , ˆ v ) = 1 | S d − | Λ b (cid:48) ( S ( ξ )) b ( S ( ξ )) ∇S ( ξ ) · ˆ v . In this case recalling(74) we obtain ∂ρ∂τ ( τ, ξ ) + ∇ · (cid:16) ρ ( τ, ξ ) χ ( S ( ξ )) ∇S ( ξ ) (cid:17) = ∇ · (cid:16) µ ( ξ ) ∇ (cid:0) D ( ξ ) ρ ( τ, ξ ) (cid:1)(cid:17) , (98)where χ ( S ( ξ )) = Λ ¯ U ( ξ ) d b (cid:48) ( S ( ξ )) b ( S ( ξ )) . (99)So, the chemotactic sensitivity χ takes into account of the kinetic response through ¯ U and of thesensing capability through γ S ( λ ) contained in Λ. More importantly, different signal-dependentsensitivity models can be obtained according to the choice of b . Viceversa, any relation on thechemotactic sensitivity can be obtained by a proper b given by b ( S ) = exp (cid:20) d Λ ¯ U (cid:90) χ ( S ) d S (cid:21) . Trivially, if b is independent of S , one has no chemotaxis. If b is proportional to S (and S (cid:54) = 0,always), then χ ( S ( ξ )) = Λ ¯ U ( ξ ) d S ( ξ ) . (100)Let us now introduce the parameter η = R S l S . (101)18e have seen that if the sensing radius is smaller then the characteristic length of variation of thechemotactic field, i . e . η (cid:28)
1, then (89) can be simplified to (96), whereas this is not possible ifthe sensing radius is larger then l S , i . e . η (cid:29)
1. This leads to different macroscopic limits (98) and(91) respectively. This different macroscopic behavior and the proper choice of the scaling may bejustified thanks to a nondimensionalization argument. We shall rescale the variables as ξ = x l S , ˜ v p = v p U ref , τ = tσ t , where we shall choose U ref = R S ¯ µ being ¯ µ a reference frequency. The time scale σ t can be chosen as a diffusion time T Diff scale ora drift time scale T Drift . In general we may write (Hillen and Othmer, 2000) T Diff = l S D , D = U ref ¯ µ , T Drift = l S U ref The regime is diffusive - and we will choose σ t = T Diff - if the frequency ¯ µ is very large withrespect to the reference time scale σ t , i . e . if we can find a small parameter (cid:15) such that¯ µ = (cid:15) − σ t . The latter is equivalent to l S = O (cid:18) U ref (cid:15) (cid:19) which implies the hierarchy T Drift (cid:28) T Diff . (102)In the present case this is equivalent to η = R S l S < . On the other hand, the macroscopic regime is hyperbolic, and we choose a faster time scale if¯ µ = (cid:15) − σ . In this case the appropriate choice will be σ t = T Drift = l S U ref = ηµ as the hierarchy (5.2.3) does not hold anymore. Hence, the following relation holds η ∼ (cid:15) (cid:29) . With respect to (100), different chemotactic coefficients can be obtained by different choices ofthe bias function b ( S ): for instance, if b ( S ) = k + S n with k >
0, then we have a saturatingchemotactic sensitivity χ ( S ( ξ )) = Λ ¯ U ( ξ ) d S n − ( ξ ) k + S n ( ξ ) . (103)Another example is the receptor binding process. In this case we may consider a saturatingdependence on the signal concentration b (cid:0) S ( x + λ ˆ v ) (cid:1) = S n ( x + λ ˆ v ) k + S n ( x + λ ˆ v ) , (104)19ielding χ ( S ( ξ )) = Λ ¯ U ( ξ ) d kn S ( ξ )[ k + S n ( ξ )] . (105)We notice that in this case, if we choose n = 1 and γ S ( λ ) = δ ( λ − R S ), i . e . , a membrane sensing,we recover the model proposed by Hillen et al. (2007).Generally speaking, b increasing with S will give rise to chemoattraction, while b decreasing with S will lead to a repulsive effect. For instance, if b ( S ( x + λ ˆ v )) = 1 k + S n ( x + λ ˆ v ) , we have (103) but with the opposite sign, i . e . , alignment in the opposite direction with respect tothe gradient of S .At this stage we did not specify whether S represents a diffusing or a matrix bound chemical.Even more, one can deal with durotaxis in a completely analogous way. In this case the signal S is the perceived stiffness of the ECM and the sensing of the mechanical properties of the ECMaround the cell will induce motion toward stiffer region of the ECM if b is an increasing functionof S , as we shall see as last application in Section 5.3. In some cases the signals perceived by the cell in differently localized receptors on its membrane isamplified by a polarization dynamics involving PTEN and the phosforillation of PIP2 into PIP3(Ambrosi et al., 2005). This causes the formation of a “head” and a “tail” in the cell that choosesthe direction accordingly. In order to mimick this phenomenon, we assume here that the turningrate depends on what is sensed in x + λ ˆ v and in x − λ ˆ v , i . e . , T λ ˆ v [ S ( x )] = α + βb (cid:0) S ( x + λ ˆ v ) , S ( x − λ ˆ v ) (cid:1) , α > β > b ∈ ( − ,
1) is a quantity comparing the values S ( x − λ ˆ v ) and S ( x + λ ˆ v ) in a way that keeps T λ ˆ v always positive. Hence, J [ p ]( t, x , v p ) = µ ( x ) (cid:34) ρ ( t, x ) ψ ( x , v ) c ( x ) (cid:32) α Λ + β (cid:90) R + b (cid:0) S ( x + λ ˆ v ) , S ( x − λ ˆ v ) (cid:1) γ S ( λ ) dλ (cid:33) − p ( t, x , v p ) (cid:35) , (107)being Λ = (cid:90) R + γ S ( λ ) dλ . In the attractive case, the cell is most likely to go where there is a largerconcentration of chemical, and therefore we might take b (cid:0) S ( x + λ ˆ v ) , S ( x − λ ˆ v ) (cid:1) = S ( x + λ ˆ v ) − S ( x − λ ˆ v )2 k + S ( x + λ ˆ v ) + S ( x − λ ˆ v ) . (108)On the other hand, in the repulsive case, the cell tends to go where there is a smaller concentrationof chemical, and therefore we might take b (cid:0) S ( x + λ ˆ v ) , S ( x − λ ˆ v ) (cid:1) = S ( x − λ ˆ v ) − S ( x + λ ˆ v )2 k + S ( x + λ ˆ v ) + S ( x − λ ˆ v ) . (109)As done in the previous section, if we can assume that the size of the neighborhood providinginformation through signaling is small, one can perform a Taylor expansion of the function T ˆ v λ in λ = 0 and write T ˆ v λ = α + β λ ∇S ( x ) · ˆ v k + S ( x ) (110)20n the attractive case (the repulsive case is similar with a minus sign, or equivalently a negative β ). The biased transition probability becomes¯ B [ S ]( ξ , ˆ v ) = 1 | S d − | (cid:18) βα Λ ∇S ( x ) · ˆ v k + S ( x ) (cid:19) , (111)that for a logarithmic dependence of b from S has the same structure as (97). In this limit theoperator for biased random motion including taxis becomes J [ p ] = µ (cid:18) ρ ( t, x ) ψ ( x , v ) 1 | S d − | (cid:18) βα Λ ∇S ( x ) · ˆ v k + S ( x ) (cid:19) − p ( t, x , v ) (cid:19) . (112)Equation (111) satisfies (33) and, therefore, a diffusive limit can be performed. As, ¯ B [ S ] ( ξ , ˆ v ) = ψ ( ξ ; v ) 1 | S d − | and ¯ B [ S ] ( ξ , ˆ v ) = βα Λ ∇S ( x ) · ˆ v k + S ( x ) , we get a drift-diffusion equation (98) with χ ( S ) = β ¯ U Λ dα k + S . (113)Even if in the limit of small sensing radii the comparative sensing is almost the same as the non-local sensing introduced in the previous section (in the sense that the diffusive limit are similar),it allows to add some details. In fact, the ratio of the coefficients β and α measures the differentweight of the diffusive and the advective (chemotactic) contributions. Further, we may include in β a dependence on other substances like S itself or auto-inducers to model quorum sensing. Forexample in the case of the cellular slime mold Dictyostelium discoideum, the response to the auto-inducer Netrin-1 may be attractive in high concentrations of cAMP, whilst it may be decreasingin case of low levels of cAMP (Deery and Gomer, 1999). Painter and Hillen (2002) consider aresponse in the form β ( w ) = 1 − ww ∗ , being w ∗ a saturation level. In addition, when the sensingradius is not small the kinetic models are different, reflecting the fact that the mechanisms arefundamentally different, with the comparative sensing giving rise to a stronger response, especiallywhen considering that the reception of the signal can be amplified by feedback mechanisms arisingfrom the activation of proper protein cascades inside the cells, giving rise, for instance to thesegregation of PIP2 and PIP3 in different parts of the cell. We simulate the kinetic model with turning operator given by (89) with b ( S ) = S both in oneand two dimensions. In order to do that, we use the numerical scheme proposed by Filbet andVauchelet (2010) in which a kinetic model for chemotaxis is simulated in two-dimensions using avan Leer scheme for the space transport. We consider a computational domain that will be in the form [ x min , x max ] × [0 , U ] in the onedimensional case and [ x min , x max ] × [ y min , y max ] × V p in the two dimensional case. The computa-tional domain is discretized with a Cartesian mesh X h × V p h , where X h and V p h are defined by(in two dimensions) X h = { x i,j = ( x i , y j ) = ( x min + i ∆ x, y min + j ∆ y ) , ≤ i ≤ n x , ≤ j ≤ n y } V p h = { v l,k = v k (cos θ l , sin θ l ) , θ l = ( l + 1 / θ, ≤ l ≤ n ang − , v k = v + k ∆ v, ≤ k ≤ n v } where ∆ x = x max − x min n y , ∆ y = y max − y min n y , ∆ v = Un v , ∆ θ = 2 πn ang . Denoting by p ni,j,l,k anapproximation of the distribution function p ( t n , x i,j , v p l,k ), where v p l,k = ( v k , ˆ v l ). We introduce21he first order splitting p n +1 / i,j,l,k − p ni,j,l,k ∆ t + v · ∇ x ,h p ni,j,l,k = 0 (transport step) p n +1 i,j,l,k − p n +1 / i,j,l,k ∆ t = µ (cid:16) ρ n +1 i,j T [ S , S (cid:48) ] i,j,l,k − p n +1 i,j,l,k (cid:17) (relaxation step)where h = (∆ t, ∆ x, ∆ y ), v · ∇ x ,h p ni,j,l,k is an approximation of the transport operator v · ∇ p computed with a Van Leer scheme. It is a high resolution monotone, conservative scheme whichis second order if the solution is smooth and first order near the shocks. T [ S , S (cid:48) ] i,j,l,k is thediscretization of the transition probability T . We observe that as the turning operator preservesmass and the turning probability is known and does not depend on p , the relaxation step may beimplicit and we may consider the density at time n + 1. In particular the density is computed byusing a trapezoidal rule ρ ni,j = ∆ v ∆ θ n v (cid:88) k =0 n ang − (cid:88) l =0 p ni,j,l,k . Concerning boundary conditions, in the one dimensional case we consider regular reflecting condi-tions. In one dimension, in a domain [0 , L ], the bounce-back and the specular reflection boundaryconditions coincide, that is p ( t, x = 0 , v ) = p ( t, x = 0 , − v ) and p ( t, x = L, − v ) = p ( t, x = L, v ).We do not consider Maxwell type conditions as only the outgoing speed would be affected. In thetwo dimensional case, the regular reflection is biologically unrealistic, as cells do not bounce backnor they collide with the wall as hard spheres. Therefore, Maxwell type boundary conditions aremore realistic, and we shall consider for the Maxwellian to the wall M ( x , v p ) = T [ S , S (cid:48) ]( x , v p ) , being T the asymptotic equilibrium of the system with this class (no-flux) of boundary conditions. Chemotactic motion
In the first simulation we consider a fixed gaussian distribution of chemoattractant S ( x ) = 0 . (cid:20) ( x − . . (cid:21) , and a constant initial condition for the cell population, as shown as in Figure 1 (a). From Figures1(b), 2, we can observe that cells tend to assume a profile similar to that of the chemoattractant(see also Figure 2). This is not surprising since for a turning rate µ and distribution function ψ ( v )independent of x (and therefore constant ¯ U and D ) in one dimension the stationary solution of(98) is given by ρ ( x ) = C [ b ( S ( x ))] m with m = Λ µ ¯ UD , (114)and the constant C given by the initial mass, i . e . , ρ ( x ) = (cid:90) Ω ρ ( x ) dx [ b ( S ( x ))] m (cid:90) Ω [ b ( S ( x ))] m dx . (115)The maximum of the cell stationary state in response to the chemoattractant depends, as shownin Figure 2, both on the size of the sensing neighborhood and on the kind of sensing. A largersensing radius leads to a stronger motility of cells and, therefore, to a higher peak in the steadystate. This is because a bigger sensing radius allows to scout a bigger neighborhood. Furthermore,the sensing of a volume ( γ S = H ) leads to an average of S over a bigger region, with respect to22 S = δ . This trend is also confirmed by the diffusive limit. In fact, the maximum density dependson the power m defined in (114). Having set µ = 50 and chosen a uniform distribution in [0 , U = 1 and D = 13 , then m = 150Λ with Λ = 0 .
02 for γ S = δ ( λ − R S ) and Λ = 0 .
01 for γ S = H ( λ − R S ). (a) (b)(c) Figure 1: Evolution of a constant macroscopic cell density in presence of a gaussian distributionof chemoattractant ((a) and (b)). The chemotactic sensitivity is b ( S ) = S , ψ ( v ) = 1 U , U = 2, and µ = 50. The sensing radius is R S = 0 . γ S = δ ( λ − R S ). The stationary distribution functionis given in (c), while the stationary density in Fig. 2 (a).Referring to Table 2, when R S = 0 . (cid:28) l S , the collision operator (89) is well approximatedby (96) and we can compare the solution between the kinetic model and Equation (98) with χ ( S ) = Λ ¯ U S . As shown in Figures 2(a),(b), the macroscopic density of the stationary solution ofthe kinetic model with turning operator given by (89) is very close to the analytic solution (114)of the diffusive limit. In addition, one can observe that a larger sensing radius gives rise to astronger sensitivity as reflected by a larger Λ and therefore a larger χ ( S ) in the advection-diffusiveequation (98) and a larger m in (114). The comparison between the two solutions remains quitegood for R S = 0 .
2, which corresponds to η = 0 .
01. Instead, as shown in Figure 2 and Table 2,an R S = 2 is not much smaller than l S , yielding to larger discrepancies. The trends obtained for23 a) (b)(c) Figure 2: Comparison of the stationary density of the kinetic model with the stationary solution(115) of the diffusion limit (98) for (a) γ S = δ ( λ − R S ) and (b) γ S = H ( λ − R S ). In (c) thesolutions for R S = 2 are given showing discrepancies between the density of the solution of thekinetic model and of the diffusive limit (98) obtained when η (cid:28) l S R S γ S Λ η relative L ∞ error17.964 0.02 δ δ δ S = δ and γ S = H are similar. However, the transient time is larger if the sensing function is theHeavyside function.In the second simulation reported in Figure 3, we use a comparative sensing kernel (112). Here R S = 0 .
02 and we consider different values of α and β . We may observe that if α is much largerthan β , then the diffusive dynamics is dominant.Figure 3: Stationary density of the kinetic model with a comparative sensing kernel for R S = 0 . γ S = δ ( λ − R S ).Finally, in the set of simulations shown in Figure 4, we start from a macroscopic gaussiandistribution for cells, that moves due to the perception of the chemoattractant which is distributedlinearly. Cells gradually move to the right until they reach the right boundary where the specularreflection boundary condition, corresponding to a no-flux condition for the macroscopic density,keeps them in the domain. Cell-cell adhesion
Using the same framework, we may also consider cellular adhesion by considering S = ρ as abiasing signal. In this case, taking for instance b ( ρ ) = ρ and γ ρ a Dirac delta, the turning operatorreads J [ p ]( t, x , v p ) = µ ( x ) [ ρ ( t, x ) ψ ( x , v ) c ( x ) ρ ( x + R ρ ˆ v ) − p ( t, x , v p )] , As initial condition we take p (0 , x, v p ) = 0 . U (1 + 0 .
09 sin x ) moving both to the right and to theleft, that corresponds to a small perturbation of the constant initial condition with unitary massthat would stay constant. However, cell-cell adhesion triggers an instability so that the smallperturbation amplifies, reaching a peaked distribution in the center of the domain as shown inFigure 5.In two dimensions we show cell-cell adhesion for a four peaks distribution. Boundary conditionsare of Maxwell-type, with α = 0 . Durotaxis
We now consider the case in which cells are guided by the rigidity of the extracellular matrixmoving towards stiffer regions. We consider a stiffness profile as in Fig. 7 (red line) and b ( S ) = S .The sensing function is again a Dirac delta and the sensing radius in the simulation shown inthe top row of Figure 7 (red line) is R S = 0 .
02. So, the cells within that radius start moving tothe stiffer region. Eventually, a stationary state is reached that depends on the sensing radius.In particular, Figure 7 (bottom row) clearly shows that for a small radius R S = 0 .
02 (whichcorresponds to η = 0 .
02) the analytic solution of the advection-diffusion equation (115) is similarto the macroscopic density of the solution to the kinetic model. If, instead, R S = 2 . t = 0 , . , ,
10 inpresence of a linear chemoattractant distribution S ( x ) = 0 . x shown in red. The sensing radius R S = 0 .
02 and other parameters are the same as in Figure 1.Figure 5: Temporal evolution of the macroscopic density starting from a perturbed constantdistribution under the action of cell-cell adhesion. Here γ ρ = δ ( λ − R ρ ), with R ρ = 1.26 a) t = 0 (b) t = 2 . t = 7 . t = 12 . t = 15 (f) t = 17 . t = 20 (h) t = 50 Figure 6: Temporal evolution of the macroscopic density starting from the four peaks distributionin (a) under the action of cell-cell adhesion. Here γ ρ = δ ( λ − R ρ ), with R ρ = 0 . η = 3), we have that (115) fails in approximating the solution to the kineticequation, while the solution to the hyperbolic macroscopic equation well approximates the profileof the cell density.In Figure 8 we present the corresponding simulations in two dimensions. We may observethat the qualitative behavior is the same. Here we show the results with Maxwell type boundaryconditions ( α = 0 in (62)). We checked that the macroscopic stationary state is the same withregular reflection boundary conditions ( α = 1 in (62)). In the previous section we have assumed that the external signal only affects the choice of theturning direction keeping the distribution function ψ determining the speed always the same and inparticular unaffected by environmental cues or external chemical signals. Instead, in cell motilityone should distinguish polarization mechanisms from the ability to move as they can depend ondifferent factors. In fact, as stated in the introduction, a cell can even be polarized but unableto move. For this reason, in this section we consider a double bias, due to a second factor S (cid:48) evaluated in a non-local way that influences the speed, once the direction ˆ v is chosen accordingto the bias ruled by S . These can be represented not only by different external free and boundchemical factors, or subcellular cues involved in cell motility, but also by cell and ECM density.In order to do that, we allow the distribution function ψ ( x , v ) to depend on the level of thesignal S (cid:48) . In order to stress this dependence, we use the following notations ψ = ψ ( x , v |S (cid:48) ( x + λ (cid:48) ˆ v ))for the distribution function, ¯ v ( x |S (cid:48) ( x + λ (cid:48) ˆ v )) for its mean, and D ( x |S (cid:48) ( x + λ (cid:48) ˆ v )) for its variance,where |S (cid:48) reads as ”given S (cid:48) in a point x + λ (cid:48) ˆ v in the neighborhood of the cell”. Then, the cellpolarized along ˆ v determines its speed averaging the signal over the sensing radius R S(cid:48) through aweight function γ S(cid:48) .In the same spirit as Eq. (63), introducing alsoΨ[ S (cid:48) ]( x , v | ˆ v ) = (cid:90) R + γ S(cid:48) ( λ (cid:48) ) ψ ( x , v |S (cid:48) ( x + λ (cid:48) ˆ v )) dλ (cid:48) , (116)the turning probability in (12) or (13) reads T [ S , S (cid:48) ]( x , v p ) = c ( x ) B [ S ]( x , ˆ v )Ψ[ S (cid:48) ]( x , v | ˆ v ) . (117)27 a) (b) (c)(d) (e) Figure 7: Cell migration under the action of durotaxis. The red line refers to the stiffness of theextracellular matrix. The initial condition is uniformly ditributed in the velocity space and has amacroscopic density which is a gaussian centered in 1 . t = 0 , ,
10 and for R S = 0 .
02 and in the bottom row the stationary statesfor R S = 0 .
02 (d) and R S = 2 . a) (b) T = 0(c) t = 0 . t = 6 (e) t = 22(f) t = 0 . t = 3 (h) t = 22 Figure 8: Evolution of an initial distribution given in (b) in presence of a heterogeneous environ-ment with stiffness given in (a). In (c)-(f) R S = 0 .
02, while in (f)-(h) R S = 2 . c can be factorized as c ( x ) = c ( x ) c ( x ) , with c ( x ) = 1 (cid:90) R + γ S(cid:48) ( λ (cid:48) ) dλ (cid:48) = 1Γ (cid:48) , and c ( x ) = 1 (cid:90) S d − (cid:90) R + γ S ( λ ) T ˆ v λ [ S ]( x ) dλ d ˆ v . In this way, ¯ B [ S ]( x , ˆ v ) = c ( x ) (cid:90) R + γ S ( λ ) T ˆ v λ [ S ]( x ) dλ (118)and ¯Ψ[ S (cid:48) ]( x , v | ˆ v ) = 1Γ (cid:48) (cid:90) R + γ S(cid:48) ( λ (cid:48) ) ψ ( x , v |S (cid:48) ( x + λ (cid:48) ˆ v )) dλ (cid:48) (119)are both distribution functions and the turning probability reads T = ¯ B ¯Ψ.Going back to the general case, the turning operator reduces to J [ p ]( x , v p ) = µ ( x ) (cid:16) ρ ( t, x ) c ( x ) B [ S ]( x , ˆ v )Ψ[ S (cid:48) ]( x , v | ˆ v ) − p ( t, x , v p ) (cid:17) (120)and the distribution function which nullifies the turning operator is p = ρ ( t, x ) c ( x ) B [ S ]( x , ˆ v )Ψ[ S (cid:48) ]( x , v | ˆ v ) . Now, U S , S(cid:48) in Eq. (17) reads U S , S(cid:48) ( x ) = c ( x ) c ( x ) (cid:90) S d − B [ S ](ˆ v ) ¯ U S(cid:48) ( x | ˆ v )ˆ v d ˆ v (121)where ¯ U S(cid:48) ( x | ˆ v ) = c ( x ) (cid:90) R + ¯ v ( x |S (cid:48) ( x + λ (cid:48) ˆ v )) γ S(cid:48) ( λ (cid:48) ) dλ (cid:48) (122)and the diffusion tensor is D S , S(cid:48) ( x ) = c ( x ) (cid:90) V p B [ S ](ˆ v ) Ψ[ S (cid:48) ]( v | ˆ v )( v − U S , S(cid:48) ) ⊗ ( v − U S , S(cid:48) ) d v p . (123)Let us now assume that, like B , Ψ[ S (cid:48) ] may be written, up to re-scaling, asΨ[ S (cid:48) ] = Ψ[ S (cid:48) ] + (cid:15) Ψ[ S (cid:48) ] + O ( (cid:15) ) . The means and variances of Ψ[ S (cid:48) ] and Ψ[ S (cid:48) ] will be, respectively, denoted by ¯ U S(cid:48) , ¯ U S(cid:48) , ¯ D S(cid:48) , ¯ D S(cid:48) .Up to re-scaling, the re-scaled turning probability becomes T [ S , S (cid:48) ]( v p ) = T [ S , S (cid:48) ] ( v p ) + (cid:15)T [ S , S (cid:48) ] ( v p ), where, now T [ S , S (cid:48) ] ( v p ) = c ( x ) B [ S ] (ˆ v )Ψ[ S (cid:48) ] ( v | ˆ v ) (124)and T [ S , S (cid:48) ] ( v p ) = c ( x ) B [ S ] (ˆ v )Ψ[ S (cid:48) ] ( v | ˆ v ) + c ( x ) B [ S ] (ˆ v )Ψ[ S (cid:48) ] ( v | ˆ v ) . (125)30n this case, Eqs.(23a) and (23b) are satisfied if (cid:90) U Ψ[ S (cid:48) ] dv = (cid:90) U Ψ[ S (cid:48) ] dv (cid:90) U Ψ[ S (cid:48) ] i dv = 0 ∀ i ≥ . Therefore, the macroscopic velocity of order 0 is U S , S(cid:48) ( ξ ) = c ( ξ ) c ( ξ ) (cid:90) S d − B [ S ] (ˆ v ) ¯ U S(cid:48) ( ξ | ˆ v )ˆ v d ˆ v (126)while the macroscopic velocity of order 1 is U S , S(cid:48) ( ξ ) = c ( ξ ) (cid:90) S d − (cid:16) B [ S ] (ˆ v )Ψ ( v | ˆ v ) + B [ S ] (ˆ v )Ψ ( v | ˆ v ) (cid:17) ˆ v d ˆ v (127)and the equilibrium diffusion tensor D S , S(cid:48) ( ξ ) = c ( ξ ) (cid:90) V p B [ S ] (ˆ v ) ¯ D S(cid:48) ( ξ | ˆ v )( v − U S , S(cid:48) ) ⊗ ( v − U S , S(cid:48) ) d v p . (128)The diffusion tensor is in general anisotropic, as ¯ D may be not isotropic in ˆ v .If (33) is satisfied, the equilibrium tensor is D S , S(cid:48) ( ξ ) = c ( ξ ) (cid:90) S d − B [ S ] (ˆ v ) ¯ D S(cid:48) ( ξ | ˆ v )ˆ v ⊗ ˆ v d ˆ v . (129)where ¯ D S(cid:48) ( ξ | ˆ v ) = (cid:90) R + D ( ξ |S (cid:48) ( x + λ (cid:48) ˆ v )) γ S(cid:48) ( λ (cid:48) ) dλ (cid:48) . In this case, we get as a diffusive limit Eq. (42), or (52) if the frequency depends weakly on thepolarization. Otherwise, if (33) is not satified, we get (56) as a hyperbolic limit.
First of all, we observe that if γ S(cid:48) ( λ ) = γ S ( λ ) = δ ( λ −
0) are Dirac deltas, then the local modelwith T [ S , S (cid:48) ] = | S d − | ψ ( x , v ) is recovered. Hence U S , S(cid:48) is zero and ¯ U S(cid:48) ( x | ˆ v ) and ¯ D S(cid:48) ( x | ˆ v ) are themean speed and variance of ψ ( x , v ) as defined in Section 5.1.In this section we shall at first consider γ S ( λ ) = δ ( λ − i . e . , when there is no bias in thedecision of the direction of motion. We shall do this in order to analyze the influence of the signal S (cid:48) . In particular, we will see that if, for instance, S (cid:48) is the ECM density then in the limit werecover a model describing the steryc hindrance of motion. Finally, we will analyse models whereboth S and S (cid:48) have a non-local influence in determining cell motion. In order to understand the meaning of the transport model we here focus on the particular casein which there is no bias in the decision of the direction of motion. This does not mean that thesituation is isotropic, because the speed can have different distribution functions in the randomlychosen direction of motion because the distribution of S (cid:48) sensed ahead may be different.In this case the turning operator (117) simplifies to T [ S , S (cid:48) ]( x , v p ) = ¯Ψ[ S (cid:48) ]( x , v | ˆ v ) , (130)where ¯Ψ is defined in (119), the macroscopic velocity of the transition probability is U S(cid:48) = 1 | S d − | Γ (cid:48) (cid:90) S d − ¯ U S(cid:48) ˆ v d ˆ v , D S(cid:48) = (cid:90) V p ¯Ψ[ S (cid:48) ]( x , v | ˆ v )( v − U S(cid:48) ) ⊗ ( v − U S(cid:48) ) d ˆ v . Up to re-scaling, we obviously have that¯Ψ[ S (cid:48) ] = ¯Ψ[ S (cid:48) ] + (cid:15) ¯Ψ[ S (cid:48) ] + O ( (cid:15) ) , and we can consider the following cases: • If ¯Ψ[ S (cid:48) ] ( ξ , v | ˆ v ) is even as a function of ˆ v , i . e . if ¯ U is even as a function of ˆ v , we have that U S (cid:48) ( ξ ) = . We can therefore perform a diffusive limit to get ∂ρ∂τ = ∇ · (cid:20) µ ∇ · (cid:0) D S(cid:48) ρ (cid:1)(cid:21) , if ¯Ψ[ S (cid:48) ] is zero and ∂ρ∂τ + ∇ · (cid:16) ρ U S(cid:48) (cid:17) = ∇ · (cid:20) µ ∇ · (cid:0) D S(cid:48) ρ (cid:1)(cid:21) , (131)otherwise. The tensor D S(cid:48) may be anisotropic, as ¯ U and ¯ D may be anisotropic as functionsof ˆ v . • If ¯Ψ[ S (cid:48) ] ( ξ , v | ˆ v ) is not even as a function of ˆ v , i . e . if ¯ U is not even as a function of ˆ v , we canperform a hyperbolic limit to get ∂ρ∂τ + ∇ · (cid:0) ρ U S(cid:48) (cid:1) = 0 . Therefore, we see that, even if there is no directional sensing, as the speed sensing is conditionedby the direction, there may be a drift driven phenomenon or an anisotropic diffusive one.
We shall now focus on the effect of the density M ( x ) of extracellular matrix. Following theexperimental results by Palecek et al. (1997), Paul et al. (1993), Wolf et al. (2013), we take intoaccount that there is an optimal ECM density or pore cross section for cell migration. In fact, forlower densities cell motion results hampered for the lack of adhesion sites that the cell can use toexert traction. On the other hand, for higher densities cell motion is hampered as well by excessiveadhesion and lack of space among the fibers. Actually, it is known that there is a maximum value ofECM density, or better a minimum value of pore area cross section, above which the fiber networkis too tight to be penetrated (the so-called physical limit of migration (Wolf et al., 2013; Sciannaand Preziosi, 2013; Giverso et al., 2017). It is more controversial whether there is a minimumvalue of ECM density denoted in the following by M below which cells are unable to migratebecause of the absence of ECM fibers on which to crawl. The existence of these physical limits ofmigration will be examined in more detail in a future work.One can then take such experimentally determined functions as an indication of the dependenceof the mean ¯ v of the distribution function ψ from the density M of extracellular matrix. In thefollowing simulation, to mimick such a behaviour we assume that¯ v ( M ) = v Max MM Max exp (cid:20) M max − MM max (cid:21) , (132)and choose ψ to be ψ ( v ) = 4 v ¯ v e − v/ ¯ v .
32e consider an initial condition which is uniformly distributed in the velocity space and withmacroscopic density a gaussian centered in x = 1 . M , M ] = [0 . ,
1] shown in Fig. 9 (a) and (b)(green line). We take here γ M ( λ (cid:48) ) = H M ( R M − λ (cid:48) ) as sensing function, meaning that the speed is determined by uniformlyaveraging the density of ECM in the direction ˆ v looking ahead up to a distance R M that in thesimulations is 0 . M max < M = 0 . v max M M max exp [ − ( M − M max ) /M max ].Hence, they are strongly slowed down by the ECM. On the other hand, if M max > M = 1 (as inthe simulations reported in Fig. 9 (b)) cell mean speed is on the increasing branch of (132), i.e.,always larger than the value given above. So, they progressively invade the whole spatial domainand crawl across of the ECM. (a) (b)(c) (d) Figure 9: Initial condition (cyan line) and density at t = 200 (full line) for M max = 1 /
10 (a)and M max = 1 / . R M = 0 .
04 and γ S is a Heavyside function. The corresponding final distribution functions at t = 200 are given in (c) and (d).In two dimensions, if a matrix density is distributed as in Fig. 10 (a), we observe that cellstend to go where the value of the matrix density is 0 . ∗ M max , i . e . where the speed is the largest.Upon reaching this zone with maximum mean speed (see Fig.10(e)) diffusion looks anisotropicbecause cells tend to follow the region with optimal ECM density staying away from the regionswhere matrix density gets closer to M max where motility is hampered. We shall now consider a real double bias, which means that there is a non-local sensing of thefield S for the directional bias, i . e . , γ S ( λ ) = δ ( λ − R S ) or γ S ( λ ) = H ( R S − λ ). In addition we take33 a) (b) t = 0 (c) t = 0 . t = 1 (e) t = 10 (f) t = 60 Figure 10: Matrix density, with M max = 1 (a). The sensing radius is R M = 0 .
25 and γ S is a Diracdelta. The time evolution of the two dimensional macroscopic density is given in (b)-(f). T ˆ v λ [ S ]( x ) = b (cid:0) S ( x + λ ˆ v ) (cid:1) . Then the turning operator reads J [ p ]( t, x , v p ) = µ ( x ) (cid:16) c ( x )Γ (cid:48) Ψ( x , v |S (cid:48) , ˆ v ) (cid:90) R + b ( S ( x + λ ˆ v )) γ S ( λ ) dλ − p ( t, x , v p ) (cid:17) (133)As (33) is hardly satisfied, the macroscopic limit is advective ∂ρ∂τ + ∇ · (cid:34) ρ (cid:90) S d − (cid:32) ¯ U S(cid:48) ( ξ | ˆ v ) (cid:90) R + b ( S ( ξ + λ ˆ v )) γ S ( λ ) dλ (cid:33) ˆ v d ˆ v (cid:35) = 0 . (134)On the other hand, if R S can be considered small, the turning operator reads J [ p ]( t, x , v p ) = µ ( x ) (cid:20) | S d − | Γ (cid:48) Ψ[ S (cid:48) ]( x , v | ˆ v ) (cid:18) b (cid:48) b ∇S · ˆ v (cid:19) − p ( t, x , v p ) (cid:21) . In this case T [ S , S (cid:48) ] = 1 | S d − | Γ (cid:48) Ψ[ S (cid:48) ] and T [ S , S (cid:48) ] = 1 | S d − | Γ (cid:48) (cid:18) Ψ[ S (cid:48) ] Λ b (cid:48) b ∇S · ˆ v + Ψ[ S (cid:48) ] (cid:19) Condition (33) is satisfied if Ψ[ S (cid:48) ] is even as a function of ˆ v for all x and the diffusive limit is(42) with advective velocity U S , S(cid:48) ( x ) = 1 | S d − | Γ (cid:48) (cid:18)(cid:90) S d − ¯ U S(cid:48) ˆ v d ˆ v + Λ b (cid:48) b T ∇S (cid:19) , T = (cid:90) S d − ¯ U S(cid:48) ˆ v ⊗ ˆ v d ˆ v . (135)It may be anisotropic as, in general, the mean speed is anisotropic.If Ψ[ S (cid:48) ] is not even as a function of ˆ v , we must perform a hyperbolic limit that gives ∂ρ∂τ + ∇ · (cid:0) ρ U S , S(cid:48) (cid:1) = 0 with U S , S(cid:48) = (cid:90) S d − ¯ U S(cid:48) ( x | ˆ v )ˆ v d ˆ v . (136)This shows that even if the directional sensing gives rise to an even distribution function of orderzero in ˆ v and, then, to a diffusive behavior, the speed sensing may change this tendency and wemay have a drift driven dynamics.As an application we shall now consider the case of cell migration in the extra-cellular matrixunder the action of a chemoattractant. We then integrate the transport equation (1) with theoperator (133) with b ( S ) = S , the mean ¯ v given by (132). We take here the same sensing functionas before, i . e ., γ M ( λ (cid:48) ) = H M ( R M − λ (cid:48) ). Also γ S ( λ ) = H ( R S − λ ).In Figure 11, the initial gaussian distribution of cells, then, starts travelling to the right underthe action of the linearly distibuted chemoattractant. However, at variance with what happendin Figure 4, when cells encounter the denser area of ECM they slow down. If M max is large( M max = 1 / . R S = 0 . M max is small ( M max = 1 / R S = 1), cells are strongly hampered by the ECM as shownin the second row of Figure 8. The kinetic model developed in this paper is based on several observations, all related to a par-ticular feature of the model. Specifically,1. in order to move cells react to external stimuli by polarizing and deciding their directionof motion, forming a “head” and “tail”. Then they try to move in that direction formingadhesion sites and polimerizing actin at the head and removing adhesion sites at the tail. Itis therefore convenient to write the kinetic model with the orientation and speed, separately,as independent variables;2. Both orientation and speed depend on cell sensing of their environment over a finite radiusrelated, for instance, to the length of their protrusion. This gives the kinetic model a non-local character;3. Orientation and speed can be related to different chemical and mechanical cues, with theformer influencing mainly orientation and the latter mainly speed (though this is not alwaysthe case). This induces the inclusion of different cues determining the changes of orientationand speed in the turning operator;4. While chemical cues are measured by the cells extending their protrusions, many mechanicalcues are related to the fact that the nucleus, which represents the stiffest organelle of thecell, has difficulty in passing through narrow pores in the ECM or through densely packedcells. So, the sensing radii involved in the identification of the orientation and of the speedmight be different with the former being larger than the latter.A particular attention was devoted to the identification of the proper macroscopic limit accord-ing to the properties of the turning operator. However, all simulations integrate the kinetic model.In some cases, the solution is compared with the solution of the related macroscopic equations.35 a) (b)(c) t = 5 (d) t = 15 (e) t = 150(f) t = 5 (g) t = 15 (h) t = 200(i) (j) (k) Figure 11: Evolution of the macroscopic cell density in response to a chemoattractant and an het-erogeneous distribution of extracellular matrix. (a): initial macroscopic density; (b): distributionfunction. (c,d,e): M max = 1 / , R S = 1. (f,g,h): M max = 1 / . , R S = 0 .
04. (i,j,k): evolution ofthe distribution function corresponding respectively to (f,g,h).36he model was in fact applied to several sensing strategies, including when the signal is perceivedlocally, averaged over a region, or perceived at the cell “head” and “tail” and then compared.Several applications were also considered as chemotaxis (or equivalently haptotaxis), durotaxis,cell-cell adhesion. The effect of the presence of an heterogeneous ECM was also considered inabsence of any cue affecting cell orientation, then yielding a random polarization, and in presenceof a chemical cue.We have left to a coming paper the case in which the density of cell or of ECM is so large thatthey represent physical limits for migration. This introduces a dependence of the sensing radiusdetermining the speed that depends on the density of cells or of ECM. In fact, for instance, if thecharacteristic pore size of the ECM is very small, then cell can protrude its cytoskeleton across thedense ECM and scout for chemical signals and polarize, but the cell cannot penetrate the ECMbecause the nucleus is stuck. This holds true even in the case in which the layer of dense ECMis very thin, like in basal membrane or in cell lining, e.g., endothelial or mesothelial lining. Thecell would have no difficulties in moving after passing through the membrane but cannot cross it.One should then address the question of deducing a kinetic model capable of dealing with suchsituations that are of great interest because of their relationship with the spread of metastasis.37 knowledgements
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-2022and Compagnia San Paolo that finances NL PhD scholarship.
References
Adler, J. (1966). Chemotaxis in bacteria.
Science , 153(3737):708–116.Alt, W. (1980). Biased random walk models for chemotaxis and related diffusion approximations.
Journal of Mathematical Biology , 9(2):147–177.Ambrosi, D., Gamba, A., and Serini, G. (2005). Cell directional persistence and chemotaxis invascular morphogenesis.
Bulletin of Mathematical Biology , 67(1):195–195.Arduino, A. and Preziosi, L. (2015). A multiphase model of tumour segregation in situ by aheterogeneous extracellular matrix.
International Journal of Non-Linear Mechanics , 75:22–30.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.Berg, H. and Purcell, E. (1977). Physics of chemoreception.
Biophysical Journal , 20(2):193 – 219.Berg, H. C. (1983).
Random Walks in Biology . Princeton University Press, rev - revised edition.Bisi, M., Carrillo, J. A., and Lods, B. (2008). Equilibrium solution to the inelastic boltzmannequation driven by a particle bath.
Journal of Statistical Physics , 133(5):841–870.Block, S. M., Segall, J. E., and Berg, H. C. (1983). Adaptation kinetics in bacterial chemotaxis.
Journal of Bacteriology , 154(1):312–323.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.Cercignani, C. (1987).
The Boltzmann Equation and its Applications . Springer, New York.Chalub, F. A. C. C., Markowich, P. A., Perthame, B., and Schmeiser, C. (2004). Kinetic modelsfor chemotaxis and their drift-diffusion limits.
Monatshefte f¨ur Mathematik , 142(1):123–141.Chauviere, A., Hillen, T., and Preziosi, L. (2007a). Modeling cell movement in anisotropic andheterogeneous network tissues.
Networks & Heterogeneous Media , 2(2):333–357.Chauviere, A., Hillen, T., and Preziosi, L. (2007b). Modeling the motion of a cell population inthe extracellular matrix.
Discrete and Continuous Dynamical Systems B , 2007(Supplementalvolume):250–259.Chauviere, A. and Preziosi, L. (2010).
A Mathematical Framework to Model Migration of a CellPopulation in the Extracellular Matrix .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.Colombi, A., Scianna, M., and Tosin, A. (2015). Differentiated cell behavior: a multiscale approachusing measure theory.
Journal of Mathematical Biology , 71:1049–1079.38eery, W. J. and Gomer, R. H. (1999). A putative receptor mediating cell-density sensing indic-tyostelium.
Journal of Biological Chemistry , 274(48):34476–34482.Devreotes, P. and Janetopoulos, C. (2003). Eukaryotic chemotaxis: Distinctions between direc-tional sensing and polarization.
Journal of Biological Chemistry , 278(23):20445–20448.Dickinson, R. B. (2000). A generalized transport model for biased cell migration in an anisotropicenvironment.
Journal of Mathematical Biology , 40(2):97–135.Eftimie, R. (2012). Hyperbolic and kinetic models for self-organized biological aggregations andmovement: a brief review.
Journal of Mathematical Biology , 65(1):35–75.Engwer, C., Hillen, T., Knappitsch, M., and Surulescu, C. (2015a). Glioma follow white mattertracts: a multiscale dti-based model.
Journal of Mathematical Biology , 71(3):551–582.Engwer, C., Hunt, A., and Surulescu, C. (2015b). Effective equations for anisotropic glioma spreadwith proliferation: A multiscale approach and comparisons with previous settings.
MathematicalMedicine and Biology : a Journal of the IMA , 33:435–459.Engwer, C., Knappitsch, M., and Christina, S. (2016). A multiscale model for glioma spreadincluding cell-tissue interactions and proliferation.
Mathematical Biosciences & Engineering ,13:443–460.Engwer, C., Stinner, C., and Surulescu, C. (2017). On a structured multiscale model for acid-mediated tumor invasion: The effects of adhesion and proliferation.
Mathematical Models andMethods in Applied Sciences , 27:1355–1390.Filbet, F., Laurencot, P., and Perthame, B. (2005). Derivation of hyperbolic models for chemosen-sitive movement.
Journal of Mathematical Biology , 50:189–207.Filbet, F. and Vauchelet, N. (2010). Numerical simulation of a kinetic model for chemotaxis.
Kinetic and Related Models , 3:B348–B366.Fisher, R., Merkl, R., and G¨unther, G. (1989). Quantitative analysis of cell motility and chemo-taxis in dictyostelium discoideum by using an image processing system and a novel chemotaxischamber providing stationary chemical gradients.
The Journal of Cell Biology , 108:973–84.Giverso, C., Arduino, A., , and Preziosi, L. (2017). How nucleus mechanics and ecm microstructureinfluence the invasion of single cells and multicellular aggregates.
Bulletin of MathematicalBiology , 80:1–29.Hillen, T. (2006). M5 mesoscopic and macroscopic models for mesenchymal motion.
Journal ofMathematical Biology , 53:585–616.Hillen, T. and Othmer, H. G. (2000). The diffusion limit of transport equations derived fromvelocity-jump processes.
SIAM Journal of Applied Mathematics , 61:751–775.Hillen, T. and Painter, K. J. (2008). A user’s guide to pde models for chemotaxis.
Journal ofMathematical Biology , 58(1):183–217.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.Hwang, H., Kang, K., and Stevens, A. (2005). Drift-diffusion limits of kinetic models for chemo-taxis: A generalization.
Discrete and Continuous Dynamical Systems - B , 5(2):319–334.Koshland, D. E. (1981). Bacterial chemotaxis as a model behavioral system.
The Quarterly Reviewof Biology , 56(4):473–474. 39apidus, R. and Schiller, R. (1976). Model for the chemotactic response of a bacterial population.
Biophysical Journal , 16(7):779 – 789.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.Lods, B. (2005). Semigroup generation propertiesof streaming operators with noncontractiveboundary conditions.
Mathematical and Computer Modelling , 42:1441–1462.Othmer, H. and Hillen, T. (2002). The diffusion limit of transport equations ii: Chemotaxisequations.
SIAM Journal of Applied Mathematics , 62:1222–1250.Othmer, H. and Stevens, A. (2001). Aggregation, blowup, and collapse: The abc’s of taxis inreinforced random walks.
SIAM Journal on Applied Mathematics , 57:1044–1081.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.
Can. Appl. Math. Q. , 10:501–543.Painter, K. J. (2008). Modelling cell migration strategies in the extracellular matrix.
Journal ofMathematical Biology , 58(4):511.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.Painter, K. J. and Sherratt, J. A. (2003). Modelling the movement of interacting cell populations.
Journal of Theoretical Biology , 225(3):327 – 339.Palcewski, A. (1992).
Velocity averaging for boundary value problems . Series On Advances InMathematics For Applied Sciences. World Scientific Publishing Company.Palecek, S., Loftus, J., H. Ginsberg, M., A. Lauffenburger, D., and Horwitz, A. (1997). Integrin-ligand binding properties govern cell migration speed through cell-substratum adhesiveness.
Nature , 385:537–40.Paul, A. D., John, A. S., John, A. Q., Steven, M. A., and Douglas, A. L. (1993). Maximal migrationof human smooth muscle cells on fibronectin and type iv collagen occurs at an intermediateattachment strength.
The Journal of Cell Biology , 122:729 – 737.Pettersson, R. (2004). On solutions to the linear boltzmann equation for granular gases.
TransportTheory and Statistical Physics , 33(5-7):527–543.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 .Scianna, M. and Preziosi, L. (2013). Modeling the influence of nucleus elasticity on cell invasionin fiber networks and microchannels.
Journal of Theoretical Biology , 317:394 – 406.Scianna, M. and Preziosi, L. (2014). A cellular potts model for the mmp-dependent and-independent cancer cell migration in matrix microtracks of different dimensions.
ComputationalMechanics , 53:485–497. 40cianna, M., Preziosi, L., and Wolf, K. (2013). A cellular potts model simulating cell migrationon and in matrix environments.
Mathematical Biosciences and Engineering : MBE , 10:235–61.Soll, D. and Wessels, D. (1998).
Motion Analysis of Living Cells . Techniques in Modern BiomedicalMicroscopy. Wiley.Stinner, C., Surulescu, C., and Meral, G. (2015). A multiscale model for ph-tactic invasion withtime-varying carrying capacities.
IMA Journal of Applied Mathematics (Institute of Mathemat-ics and Its Applications) , 80:1300–1321.Stroock, 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.te Boekhorst, V., Preziosi, L., and Friedl, P. (2016). Plasticity of cell migration in vivo and insilico.
Annual Review of Cell and Developmental Biology , 32:491–526.Tosin, A. and Frasca, P. (2011). Existence and approximation of probability measure solutions tomodels of collective behaviors.
Networks & Heterogeneous Media , 6(1):561–596.Wolf, K., Te Lindert, M., Vortmeyer (Krause, M., Alexander, S., te Riet, J., L Willis, A., M Hoff-man, R., Figdor, C., J Weiss, S., and Friedl, P. (2013). Physical limits of cell migration: Controlby ecm space and nuclear deformation and tuning by proteolysis and traction force.