Multi-cue kinetic model with non-local sensing for cell migration on a fibers network with chemotaxis
MMulti-cue kinetic model with non-local sensing for cellmigration on a fibers network with chemotaxis
Martina Conte ∗ Nadia Loy †‡ July 14, 2020
Abstract
Cells perform directed motion in response to external stimuli that they detect by sensingthe environment with their membrane protrusions. In particular, several biochemical and bio-physical cues give rise to tactic migration in the direction of their specific targets. This definesa multi-cue environment in which cells have to sort and combine different, and potentiallycompetitive, stimuli. We propose a non-local kinetic model for cell migration in presence oftwo external factors both influencing cell polarization: contact guidance and chemotaxis. Wepropose two different sensing strategies and we analyze the two resulting models by recoveringthe appropriate macroscopic limit in different regimes, in order to see how the size of the cell,with respect to the variation of both external fields, influences the overall behavior. Moreover,we integrate numerically the kinetic transport equation in a two-dimensional setting in orderto investigate qualitatively various scenarios.
Keyword.
Kinetic equations, multiscale modeling, multi-cue, non-local, hydrodynamic limit,cell migration, contact guidance, chemotaxis
AMS subject classifications.
Cell migration is a fundamental mechanism in a huge variety of processes, such as embryogenesis,wound healing, angiogenesis, immune response and tumor stroma formation and metastasis.During such processes, cells sense the environment and respond to external factors that inducea certain direction of motion towards specific targets ( taxis ): this results in a persistent migrationin a certain preferential direction. The guidance cues leading to directed migration may be bio-chemical or biophysical. Biochemical cues can be, for example, soluble factors or growth factorsthat give rise to chemotaxis , which involves a mono-directional stimulus. Other cues generat-ing mono-directional stimuli include, for instance, bound ligands to the substratum that induce haptotaxis , durotaxis , that involves migration towards regions with an increasing stiffness of theECM, electrotaxis , also known as galvanotaxis , that prescribes a directed motion guided by anelectric field or current, or phototaxis , referring to the movement oriented by a stimulus of light[37]. Important biophysical cues are some of the properties of the extracellular matrix (ECM),first among all the alignment of collagen fibers and its stiffness. In particular, the fiber align-ment is shown to stimulate contact guidance [25, 24]. Contact guidance is a key mechanism in anumber of in vivo situations in which cells tend to migrate crawling on the fibers, thus following ∗ BCAM - Basque Center for Applied Mathematics, Alameda de Mazarredo, 14, 48009 Bilbao, Spain( [email protected] ) † 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] ) ‡ Corresponding author: [email protected] a r X i v : . [ q - b i o . CB ] J u l he directions imposed by the network structure of the ECM. This is a bi-directional cue, as, ifthe fibers network is not polarized, there is no preferential sense of migration along them. Forexample, during wound healing fibroblasts migrate efficiently along collagen or fibronectin fibersin connective tissues; in cancer spread and metastasis formation, cancer cells migrate through thestromal tissue and are thus facilitated to reach blood and lymphatic vessels [58, 52, 53].In many processes there are several directional cues that may induce different simultaneousstimuli. While the cell response to each of them has been largely studied, from both an intracellularand a migrative point of view, cell responses to a multi-cue environment are much less understood.The fundamental issue is the way cells rank, integrate or hierarchize multiple cues, in particularwhen these give conflicting stimuli, because, for example, they are not co-aligned [54]. Somestudies have shown that there may be competition or cooperation between different stimuli inthe directional response of a cell in a multi-cue environment. Considering the angle between therelative orientation of the directional cues, in the mono-directional case they compete when thisangle is π , whereas they collaborate when this angle is 0. Bi-directional cues, such as contactguidance, compete when the angle is π/
2. Then, many intermediate scenarios may happen andguidance stimuli submit or prevail according to other factors, among all their average concentrationand intensity, that relates to the steepness of the gradient for taxis processes and to the degreeof alignment for contact guidance. In particular, regarding the external environment, the averagevalue of the directional cue (fiber density, molecule concentration, etc.) and the steepness ofthe gradient, or the degree of fiber alignment, are fundamental parameters that can be quantified.While, for cell migration, the angle between the polarization direction and the preferential directionimposed by the guidance cue can be measured, as well as the displacement, the mean squareddisplacement and the persistence time [17]. However, in general, when cues are aligned, a simpleadditive mechanism is not what governs multi-cue migration [37], even if it is weighted by theaverage cue concentrations or intensities.In the framework of kinetic models, in the present paper we will focus on how the environmentalsensing of two different stimuli over a finite radius can influence the choice of the direction ofmotion of a cell. In particular, we combine chemotaxis, a mono-directional biochemical cue, withcontact guidance, defining the new orientation of the cells as a result of the sensing of the twocues over a finite neighborhood, that gives a non-local character to the model. In particular,the combination of chemotaxis and contact-guidance happens in vivo in a variety of situations,for example in wound healing and in breast cancer. In wound healing, fibers guide cells towardsthe provisional clot, whilst in breast cancer cells follow the aligned fibers at the tumor-stromainterface for migrating out of the primary tumor. Chemotaxis accelerates and enhances theseprocesses [37, 8, 52, 53]. Therefore, a deep understanding of multi-cue migrational responses is akey step for the comprehension of both physiologic and pathologic processes, but also for buildingengineered tissues, as their structure is realized for guiding cell migration in a focused way [37].There are not many experimental studies concerning chemotaxis and contact guidance, as wellas other combinations of directional guidances cues [37]. One of the main reasons is the difficulty indesigning environments for controlling multiple directional cues, in particular soluble factors andaligned fibers and fibrous materials. For example, in one of the first works studying in vitro contactguidance of neutrophil leukocytes on fibrils of collagen [62], it is shown that migration is moreefficient in the direction of alignment, instead of in the perpendicular direction; in the presenceof chemotaxis, obtained by adding a chemoattractant, they observe that these cues cooperate orcompete in dependence on their relative orientation. In particular, the chemotactic response islower for cells trying to cross fibers in the perpendicular direction. In [8], it is shown that alignmentalong the fibers is greater in presence of a co-aligned chemoattractant. In [41], the authors studyhow multiple uniformly distributed cues quantitatively regulate random cell migration. One of thelatest works concerning the competition between chemotaxis and contact guidance shows that lesscontractile cells are dominated by chemotaxis, while contact guidance might dominate in morecontractile cells [55]. This suggests that, as amoeboid cells are less contractile, while mesenchymalcells are more contractile, and there may be a switching between amoeboid and mesenchymalmigration, perhaps there can also be a switching between the dominance of chemotaxis (amoeboidmigration) and contact guidance (mesenchymal migration) [63]. One of the most interesting 2D2latforms, allowing to study contact guidance and chemotaxis, was proposed in [60], in which theauthors demonstrated an additive effect of chemical gradients and fiber alignment by measuringthe persistence time; they also observed that cells were directed by fiber alignment and there wasno effect of the chemical gradient when fibers were aligned perpendicular to it. A similar settingwas also used for studying the dependence of contact guidance on the cell cycle [51]. However, Inthe case of different multi-directional cues, totally different scenarios may happen, e . g . in [54] it isshown that for contact guidance and electrotaxis in the cornea, electrotaxis wins when competingwith the direction of alignment of the fibers.There is a huge variety of mathematical models concerning cell migration. They range frommicroscopic models (also called individuals based models), that describe migration at the celllevel, up to macroscopic ones, that describe collective cell-migration at a tissue level. There aremany examples of individual based models regarding chemotaxis ([16, 26] and references therein)and migration on the ECM [13, 57, 56]. Concerning macroscopic models, first among all thefamous Keller and Segel model is a drift-diffusion model postulated at the macroscopic level[32]. Many efforts were made in order to encompass the defects of the Keller and Segel model,as well as for deriving it from lower scale models (see [33, 30, 43, 44, 3] and references therein).Between microscopic and macroscopic models there are mesoscopic models that are an intermediaterepresentative scale, as they include microscopic dynamics and describe the statistical distributionof the individuals. They also allow, for instance in the case of kinetic theory, to recover theappropriate macroscopic regime which inherit some details of the microscopic dynamics, thusgiving more significance to some of the parameters [43]. Some examples are [14, 9, 19]. The twomajor models for contact guidance at the mesoscopic level were proposed in [27] and [18], bothlocal models in the physical space. Concerning multiple cues, not many models exist. In [34], theauthors propose a macroscopic drift-diffusion model derived from a space jump process in whichthey include the response to multiple chemicals. A recent review for macroscopic PDEs includingmultiple-taxis has been proposed in [35]. In [61], the authors propose one of the first models forboth contact guidance and chemotaxis, derived from a microscopic dynamics description. In arecent work [1], the authors propose a microscopic stochastic model for studying contact guidanceand add chemotaxis in order to study migration at the tumor-stroma interface for classifying TACS(tumor associated collagen signature). In [10], a kinetic model for cell-cell interactions on a fibersnetwork in presence of a tactic cue is considered. In [39, 40], the authors propose a non-localkinetic model with a double biasing cue: the first one affecting the choice of the direction andthe second one affecting the speed, including, through the non-locality, the sensing of macroscopicquantities performed by the cell, that depends on the cell size, i . e . , on its maximum protrusionlength.As already stated, in this paper we want to include chemotaxis and contact guidance as di-rectional cues guiding cell polarization. In particular, we analyze two possible sensing strategiesthat a cell could apply for exploring the neighborhood around, and that determine the choicefor the transition probability for the transport model. The cell can measure the guidance cuesindependently, and, then, choose the new orientation using the collected information, eventuallyweighted in different ways. Otherwise, it can measure the two directional stimuli, weighting themequally, and assuming a conditioning of one cue on the other. Therefore, cell response is related tothe choice of the sensing strategy, and the macroscopic overall effect of the two cues would also beaffected. Moreover, we shall consider for the first time a non-local sensing of the fibers distributiondefined at a mesoscopic level; this allows for many intermediate scenarios in the analysis aboutthe collaborative or competitive effect of the cues. For a better understanding, we discuss howthe choices made on the transition probability, together with the size of the sampling volume andthe characteristics of the two cues determine the macroscopic behavior. Specifically, in section2, we shall present the mathematical framework, while in section 3 we shall introduce the twoclasses of models, that describe the different strategies for the sensing of a double cue, along withthe corresponding macroscopic limits in various regimes, depending on the cell size and on thevariability of the external cues. In section 4, some numerical simulations of the kinetic models willbe presented for investigating qualitatively various scenarios in a two-dimensional setting.3 Mathematical framework
The cell population will be described at a mesoscopic level through the distribution density p = p ( t, x , v, ˆ v ) that, for every time t > x ∈ Ω ⊆ R d , gives the statistical distributionof the speeds v ∈ [0 , U ], where U is the maximal speed a cell can achieve, and of the polarizationdirections ˆ v ∈ S d − , being S d − the unit sphere boundary in R d . The velocity vector, thus, willbe given by v = v ˆ v .Then, a macroscopic description for the cell population can be classically recovered throughthe definition of moments of the distribution function p . In particular, we recover the cell numberdensity ρ ( t, x ) ρ ( t, x ) = (cid:90) S d − (cid:90) U p ( t, x , v, ˆ v ) dv d ˆ v (1)the momentum ρ ( t, x ) U ( t, x ) = (cid:90) S d − (cid:90) U v p ( t, x , v, ˆ v ) dv d ˆ v (2)the cell mean velocity U ( t, x ) = 1 ρ ( t, x ) (cid:90) S d − (cid:90) U v p ( t, x , v, ˆ v ) dv d ˆ v (3)and the energy tensor D ( t, x ) = (cid:90) S d − (cid:90) U ( v − U ) ⊗ ( v − U ) p ( t, x , v, ˆ v ) dv d ˆ v . (4)The mesoscopic model consists in the transport equation for the cell distribution ∂p∂t ( t, x , v, ˆ v ) + v · ∇ p ( t, x , v, ˆ v ) = J [ p ]( t, x , v, ˆ v ) (5)where the operator ∇ denotes the spatial gradient, so that the term v ·∇ p takes into account the freeparticle transport. The term J [ p ]( t, x , v, ˆ v ) is the turning operator that describes the scattering ofthe microscopic velocity in direction and speed. This is related to the typical microscopic dynamicsof the cell, that is the run and tumble [7, 4]. The run and tumble prescribes an alternation ofruns over straight lines and re-orientations: the choice of the new direction may be random or itmay be biased by the presence of external factors, that may attract or repel the cell as well asincrease the time spent in a run. The run and tumble is classically modeled by a scattering of themicroscopic velocity called velocity jump process [59], characterized by a turning frequency µ anda transition probability T . The general form of the turning operator which implements a velocityjump process at a kinetic level is given by J [ p ]( x , v, ˆ v ) = µ ( x ) (cid:90) S d − (cid:90) U (cid:104) T ( x , v, ˆ v | v (cid:48) , ˆ v (cid:48) ) p ( t, x , v (cid:48) , ˆ v (cid:48) ) − T ( x , v (cid:48) , ˆ v (cid:48) | v, ˆ v ) p ( t, x , v, ˆ v ) (cid:105) dv (cid:48) d ˆ v (cid:48) (6)where we assumed that the turning frequency does not depend on the microscopic velocity. Thetransition probability T ( x , v, ˆ v | v (cid:48) , ˆ v (cid:48) ) is also called turning kernel and it is a conditional probabilitysatisfying, ∀ x ∈ Ω, (cid:90) S d − (cid:90) U T ( x , v, ˆ v | v (cid:48) , ˆ v (cid:48) ) dvd ˆ v = 1 , ∀ v (cid:48) ∈ [0 , U ] , ˆ v (cid:48) ∈ S d − . (7)Thanks to this property, the operator (6) reads J [ p ]( t, x , v, ˆ v ) = µ ( x ) (cid:32)(cid:90) S d − (cid:90) U T ( x , v, ˆ v | v (cid:48) , ˆ v (cid:48) ) p ( t, x , v (cid:48) , ˆ v (cid:48) ) dv (cid:48) d ˆ v (cid:48) − p ( t, x , v, ˆ v ) (cid:33) . T ( x , v, ˆ v | v (cid:48) , ˆ v (cid:48) ) = T ( x , v, ˆ v ) (8)as classically done in the pioneering work concerning kinetic equations for velocity jump processes[59, 45, 27]. This assumption, along with the assumption on the turning frequency, is due to thefact that we shall consider directional cues which are sensed non-locally, and, therefore, the mostrelevant aspect will be the measured preferential direction instead than the incoming velocity. Thelatter (8) allows to write the turning operator as J [ p ]( t, x , v, ˆ v ) = µ ( x ) (cid:16) ρ ( t, x ) T ( x , v, ˆ v ) − p ( t, x , v, ˆ v ) (cid:17) . (9)The mean macroscopic velocity after a tumble is given by the average of T U T ( x ) = (cid:90) S d − (cid:90) U v T ( x , v, ˆ v ) dv d ˆ v (10)and the diffusion tensor by the variance-covariance matrix D T ( x ) = (cid:90) S d − (cid:90) U T ( x , v, ˆ v )( v − U T ) ⊗ ( v − U T ) dv d ˆ v . (11)Arguing as in [49, 6], we can prove a linear version of the classical H-Theorem for the linearBoltzmann equation (5)-(9) with p = p (0 , x , v, ˆ v ) ∈ L (Ω × [0 , U ] × S d − ). In particular theMaxwellian M ( x , v, ˆ v ) = ρ ∞ ( x ) T ( x , v, ˆ v ) , making the turning operator vanish, is the local asymptotic stable equilibrium of the system. Asalready remarked by [39], this implies that T is the local asymptotic equilibrium steady state ofthe system. Therefore U T and D T are the mean velocity and diffusion tensor of the cell populationat equilibrium. Since we are going to consider two-dimensional bounded domains without loss of cells and nocells coming in, we shall assume conservation of mass. Therefore, we will require that the chosenboundary condition is no-flux [50] (cid:90) S d − (cid:90) U p ( t, x , v, ˆ v )ˆ v · n ( x ) dv d ˆ v = 0 , ∀ x ∈ ∂ Ω , t > , (12)being n ( x ) the outward normal to the boundary ∂ Ω in the point x . This class of boundaryconditions is part of the wider class of non-absorbing boundary conditions. Denoting the boundaryoperator as R [ p ]( t, x , v, ˆ v ) = p ( t, x , v (cid:48) , ˆ v (cid:48) ) , there are two important classes of kinetic boundary conditions which satisfy (12): the regularreflection boundary operators and the non-local (in velocity) boundary operators of diffusive type.We address the reader to the works [48] and [38] for the definition of these boundary operators.In the present work, we shall consider specular reflection boundary conditions p ( t, x , v (cid:48) , ˆ v (cid:48) ) = p (cid:18) t, x , v, ˆ v − v · n ) n | ˆ v − v · n ) n | (cid:19) , n · ˆ v ≤ , (13)that means that cells are reflected with an angle of π/ .3 Macroscopic limits In order to investigate the overall trend of the system, the macroscopic behavior is typicallyanalyzed. By integrating Eq. (5) with (9) on S d − × [0 , U ], thanks to Eq. (7), we have that ∂ t ρ ( t, x ) + ∇ · ( ρ ( t, x ) U ( t, x )) = 0 , i . e . , the mass is conserved pointwise and in the entire domain, because of no-flux boundaryconditions (after integration on Ω). If we multiply Eq. (5) with (9) by v ˆ v , and we then integratethe result on S d − × [0 , U ], we see that the momentum is not conserved ∂ t ρ ( t, x ) U ( t, x ) + ∇ · ( ρ ( t, x ) D T ( t, x )) = µ ( x ) ( ρ ( t, x ) U T ( x ) − ρ ( t, x ) U ( t, x )) . We can observe that, if we multiply the transport equations by increasing orders n of power of v and, then, we integrate on the velocity space, we obtain a non-closed system of macroscopicequations, since the equations describing the evolution of n th moment of p contain the ( n + 1) th moment. Therefore, we need some procedures to obtain a closed evolution equation (or system ofequations) for the macroscopic quantities. In particular, we are interested in the evolution of ρ ( t, x )in the emerging regime of the system. Therefore, we shall consider a diffusive or a hydrodynamicscaling of the transport equation (5) with (9), resulting from a proper non-dimensionalization of thesystem. Diffusive and hydrodynamic limits for transport equations with velocity jump processeshave been widely treated in [29, 43, 27, 39, 2, 23]. Formally, we introduce a small parameter (cid:15) (cid:28) ξ = (cid:15) x , (14)being ξ the macroscopic spatial variable. According to the other characteristic quantities of thesystem of study, the macroscopic time scale τ will be τ = (cid:15) t, (15)that is the parabolic scaling representing a diffusion dominated phenomenon, or τ = (cid:15)t, (16)that is the hyperbolic scaling that represents a drift driven phenomenon. Up to the spatial scaling(14), we have that the transition probability may be expanded as T ( ξ , v, ˆ v ) = T ( ξ , v, ˆ v ) + (cid:15)T ( ξ , v, ˆ v ) + O ( (cid:15) ) . Therefore, the corresponding means and diffusion tensors will be given by U iT ( ξ ) = (cid:90) S d − (cid:90) U T i ( ξ , v, ˆ v ) v dvd ˆ v (17)and D iT ( ξ ) = (cid:90) S d − (cid:90) U T i ( ξ , v, ˆ v )( v − U iT ) ⊗ ( v − U iT ) dv d ˆ v . (18)Considering a Hilbert expansion of the distribution function pp = p + (cid:15)p + O ( (cid:15) ) , (19)if there is conservation of mass, we have that all the mass is in p [29], i . e . , ρ = ρ, ρ i = 0 ∀ i ≥ , (20)where ρ i = (cid:90) S d − (cid:90) U p i dv d ˆ v . Furthermore, for performing the diffusive limit we shall assumethat (cid:90) S d − (cid:90) U p i v dv d ˆ v = 0 ∀ i ≥ i . e . , forchoosing τ = (cid:15) t ) is U T = 0 , (21)meaning that the leading order of the drift vanishes, which is coherent with the fact that the timescale τ = (cid:15) t is chosen because the phenomenon macroscopically is diffusion-driven. The diffusivelimit procedure prescribes to re-scale (5)-(9) with (14)-(15) and to insert (19) in the re-scaledequation. By comparing equal order of (cid:15) , we obtain the macroscopic diffusive limit, given by(dropping the dependencies) ∂∂τ ρ + ∇ · (cid:0) U T ρ (cid:1) = ∇ · (cid:20) µ ∇ · (cid:0) D T ρ (cid:1)(cid:21) , (22)being D T ( ξ ) = (cid:90) S d − (cid:90) U T ( ξ , v, ˆ v ) v ⊗ v dvd ˆ v the diffusion motility tensor. Equation (22) is a diffusion-advection equation, where U T is thedrift velocity of first order. If (21) does not hold, a hyperbolic scaling is required, that gives ∂∂τ ρ + ∇ · (cid:0) ρ U T (cid:1) = 0 . (23)This is an advection equation modeling a drift driven phenomenon. We address the reader to [39]for further details.Concerning the boundary conditions, at the macroscopic level (12) gives [50] (cid:16) D T ∇ ρ − ρ U T (cid:17) · n = 0 , on ∂ Ω , for the diffusive limit, whilst for the hyperbolic limit the corresponding boundary condition is U T · n = 0 , on ∂ Ω . In this section, we shall introduce the transition probability modeling a decision process of a cell inpresence of a double directional guidance cue: a fibrous ECM and a chemoattractant. In particular,we shall consider amoeboid cells [63] moving by contact guidance without proteolysis: cells hitthe fiber and then move along the direction of the fiber itself. It has been shown experimentally,for example in the case of glioma cancer cells [31], that randomly disposed fibers imply isotropicdiffusion of cells, while aligned fibers cause anisotropic diffusion of cells along the preferentialdirection of the fibers themselves. The first transport model for contact guidance was proposedby [27], further studied and developed by [46, 10, 11] and applied to the study of glioma by[47, 22, 21, 15, 20]. The model proposed by [27] prescribes a distribution of fibers on the space ofdirections, given by the unit sphere in R n , q = q ( x , ˆ v ) , x ∈ Ω , ˆ v ∈ S d − (24)that satisfiesQ1: q ( x , ˆ v ) > , ∀ x ∈ Ω , ˆ v ∈ S d − Q2: (cid:90) S d − q ( x , ˆ v ) d ˆ v = 1 , ∀ x ∈ ΩQ3: q ( x , ˆ v ) = q ( x , − ˆ v ) , ∀ x ∈ Ω , ˆ v ∈ S d − , 7here the last condition means that we are considering a non-polarized network of fibers, so thatcells are able to go in both senses in every direction. Being, then, q ( x , ˆ v ) a probability density, wecan define the mean direction of the fibers E q ( x ) = (cid:90) S d − q ( x , ˆ v ) ˆ v d ˆ v , (25)and the diffusion tensor of the fibers, given by the variance-covariance matrix of q D q ( x ) = (cid:90) S d − q ( x , ˆ v ) (ˆ v − E q ) ⊗ (ˆ v − E q ) d ˆ v . (26)As we consider a non polarized fibers network, we have that E q ( x ) = 0 , (27)meaning that there is no mean direction in the dynamics. The tensor (26) is symmetric andpositive definite, when q is a regular probability distribution, and, thus, it is diagonalizable. Eacheigenvalue represents the diffusivity in the direction of the corresponding eigenvector, meaningthat, if the eigenvalues are equal, there is isotropic diffusion, while, if they are different, there is apreferential direction of motion, i . e . anisotropy. Therefore, the model introduced in [27], as shownin [46], allows to reproduce isotropic/anisotropic diffusion on a non-polarized fibers network.Concerning chemotaxis, we shall consider a chemoattractant in the region Ω defined by astrictly positive definite function S = S ( x ) : Ω (cid:55)−→ R + . (28)We consider that the sensing performed by the cells is non-local, as they may extend theirprotrusions, through which they sense the environment, up to several cell diameters [5]. Themaximum length R of a protrusion is called sensing radius and it has been first introduced in[43] for modeling a non-local gradient of a chemical and, then, used in a number of works (see[12] for a review and references therein) for describing the sensing of macroscopic quantities. Inparticular, in [39] and, later, in [40] the authors propose a double bias model, in which two cuesare sensed non-locally and they affect cell polarization and speed. In the present work we shalldrop the sensing of a cue that affects the speed, that will be unbiased, and we will extend themodel proposed in [39] to a double sensing of cues affecting the polarization of the cell.Therefore, in the model both S and q will be sensed non-locally by a cell that, starting fromits position x , extends its protrusions in every direction ˆ v ∈ S d − up to the distance R , given bythe sensing radius. In particular, assuming a non-local sensing of the fibers network will allow toreproduce a wider range of migration strategies, that a cell can perform in order to cleverly reachthe chemoattractant, with respect to a local sensing. Therefore, we shall consider the quantities S ( x + λ ˆ v ) , q ( x + λ ˆ v , ˆ v ) , ∀ x ∈ Ω , ∀ ˆ v ∈ S d − , λ ≤ R. Of course, next to the border of the domain Ω, we shall always consider λ such that x + λ ˆ v ∈ Ω.In order to analyze qualitatively the impact of the non-locality at the macroscopic level, westudy, as previously done in [39, 40], the impact of the directional cues S and q with respect tothe size of the cell, that is related to its sensing radius R . Thus, we introduce the characteristiclength of variation of S as l S := 1max x ∈ Ω |∇S|S . (29)It allows to approximate S ( x + λ ˆ v ) with a positive quantity S ( x + λ ˆ v ) ∼ S ( x ) + λ ∇S · ˆ v ≥ ∀ λ ≤ R if R < l S (30)where we neglected higher order terms in λ . Beside the above defined characteristic length ofvariation of the chemoattractant l S , we define an analogue quantity for the fibers distribution. Wechoose l q := 1max x ∈ Ω max ˆ v ∈ S d − |∇ q · ˆ v | q . (31)8n this case, we can approximate q ( x + λ ˆ v , ˆ v ) with a positive quantity q ( x + λ ˆ v , ˆ v ) ∼ q ( x , ˆ v ) + λ ∇ q · ˆ v ≥ ∀ λ < R if R < l q . (32)In particular, this definition of l q takes into account the variation of directionality of the fibersin space, that is what actually influences the cell orientation, more than spatial variation of thedensity of the extracellular matrix. We analyze the possible scenarios depending on the relationbetween R , l S and l q .In analogy to [39], let us now introduce the parameters η q := Rl q (33)and η S := Rl S , (34)that quantify the capability of measuring of the cell with respect to the characteristic lengths ofvariation of the sensed guidance cues q and S . In particular, η i < , i = q, S , means that thesensing radius is smaller than the characteristic length of variation of q ( S , respectively) and theidea is that a single instantaneous sensing of the cell is not capable of catching the total spatialvariability of q ( S , respectively), while if η i > , i = q, S , the sensing radius is large enough inorder to capture the spatial variability of q ( S , respectively). If we consider the two cues separately,in the first case we expect that the sensing of q ( S , respectively) induces a diffusive behavior, whilein the second scenario the overall behavior induced by q ( S , respectively) is drift-driven.As we are considering the two guidance cues simultaneously affecting cell polarization, we nowtake into account for limit cases: i ) η q , η S (cid:29) ii ) η q , η S (cid:28) iii ) η S (cid:28) , η q (cid:29) iv ) η S (cid:29) , η q (cid:28) i ), a Taylor expansion cannot be used, since there is no guarantee that the first orderapproximations are positive, as well as in case iii ) and iv ) for q and S , respectively.In order to quantify the relative contribution of chemotaxis to contact guidance, we mayintroduce the parameter η = η q η S (35)that is larger than 1 if contact guidance prevails, whilst it is smaller then 1 if chemotaxis isstronger. Due to (33) and (34), we have that, despite its definition, η does not depend on the sizeand sensing capability of the cell, as η = η q η S = l S l q . In particular, if l S is larger than l q , i . e . η > q is steeper than the one of S , thus enhancing a stronger effect ofcontact guidance on the dynamics. We may also observe that in case iii ) we have always that η > iv ) we always have η < i . e . contact guidance is weaker then chemotaxis.We shall propose two different transition probabilities describing two different sensing strate-gies: in the first model the sensings of q and S are independent, while in the second model aunique sensing is performed. In the first model, we shall introduce a transition probability that isthe product of two different independent sensings T [ q, S ]( x , v, ˆ v ) = c ( x ) (cid:90) R + γ S ( λ ) S ( x + λ ˆ v ) dλ (cid:90) R + γ q ( λ ) q ( x + λ ˆ v , ˆ v ) dλ ψ ( v ) . (36)In this case the cell located in position x measures along the direction ˆ v the field S ( x + λ ˆ v ) weightedby γ S , and, independently, the quantity q ( x + λ ˆ v , ˆ v ), weighted by γ q . The sensing functions γ S γ q have compact support in [0 , R ] and they may be Dirac deltas centered in R , if the cell onlymeasures the guidance cues on its membrane (only on x + R ˆ v for every ˆ v ), or Heaviside functionsif the cell measures and gives the same weight to q and S from x to x + R ˆ v in every direction.Formally the transition probability might be seen as the product of the independent probabilitiesof q and S , i . e . T [ q, S ] = ˆ T [ q ] ˆ T [ S ].The second model prescribes a simultaneous averaging of the guidance cues S and q , i . e . , T [ q, S ]( x , v, ˆ v ) = c ( x ) (cid:90) R + γ ( λ ) S ( x + λ ˆ v ) q ( x + λ ˆ v , ˆ v ) dλ ψ ( v ) . (37)This transition probability describes a cells in position x that measures in the direction ˆ v the twoquantities S ( x + λ ˆ v ) and q ( x + λ ˆ v ), weighting both with γ , that is a sensing function. Formally,as the two sensing are not independent and, therefore, factorized, we have a conditioning of S given q and viceversa, i . e . , T [ q, S ] = ˜ T [ S| q ] ˜ T [ q ] = ˜ T [ q |S ] ˜ T [ S ].In (36) and (37), c ( x ) is a normalization coefficient. Moreover the probability density ψ is thedistribution of the speeds on the interval [0 , U ] and satisfies (cid:90) U ψ ( v ) dv = 1 . We introduce its mean speed ¯ U = (cid:90) U v ψ ( v ) dv (38)and the second moment D = (cid:90) U v ψ ( v ) dv , (39)such that the variance of ψ is given by σ ψ = 12 ( D − ¯ U ).We shall refer to the transport model (5)-(9) with (36) as non-local independent sensing model ,in which the cell averages the two cues independently according to two different sensing functions γ q , γ S . On the other hand, the transport model (5)-(9) with (37) is defined as non-local dependentsensing model , describing cells that sense the two cues at the same time and average them with aunique sensing kernel γ . In the next sections we shall analyze the macroscopic limits for the twomodels in the scenarios i ) − iv ) and we shall compare the two models. We first consider the non-local independent sensing case (5)-(9) with (36). We recall the expressionof the transition probability T [ q, S ]( x , v, ˆ v ) = c ( x ) (cid:90) R + γ S ( λ ) S ( x + λ ˆ v ) dλ (cid:90) R + γ q ( λ ) q ( x + λ ˆ v , ˆ v ) dλ ψ ( v ) . The average of T , that will be the equilibrium velocity of the cell population, is given by U T ( x ) = c ( x ) ¯ U (cid:90) S d − ˆ v (cid:32)(cid:90) R + γ S ( λ ) S ( x + λ ˆ v ) dλ (cid:90) R + γ q ( λ ) q ( x + λ ˆ v , ˆ v ) dλ (cid:33) d ˆ v . (40) Case i ) In this case, we shall choose (cid:15) = min (cid:26) η q , η S (cid:27) .
10s a consequence of the fact that T cannot be expanded in powers of (cid:15) after re-scaling with (14),we have that U T = U T given by (40). Therefore, we have to perform a hyperbolic scaling thatleads to the following macroscopic equation for the cells macroscopic density: ∂∂τ ρ ( τ, ξ ) + ∇ · ( ρ ( τ, ξ ) U T ( ξ )) = 0 , (41)with U T ( ξ ) given by the re-scaling of (40) with (14). Case ii ) In this case, we can expand both S ( x + λ ˆ v ) and q ( x + λ ˆ v , ˆ v ) and consider the approx-imations (30) and (32) for λ < min { l q , l S } . Therefore, we approximate the transition probabilityby substituting (30) and (32) in (36), and, thus, we obtain the following approximation for theturning kernel T [ q, S ], that reads T [ q, S ]( x , v, ˆ v ) = c ( x ) (cid:104) Γ S Γ q S ( x ) q ( x , ˆ v ) + Γ S Γ q S ( x ) ∇ q · ˆ v + Γ S Γ q q ( x , ˆ v ) ∇S · ˆ v (cid:105) ψ ( v ) (42)where we neglected higher orders terms in λ . In the latter c ( x ) = 1 S ( x ) Γ S Γ q and Γ S i := (cid:90) R + λ i γ S ( λ ) dλ i = 0 , qi := (cid:90) R + λ i γ q ( λ ) dλ i = 0 , . The quantities Γ q , Γ S are the weighted (by γ q , γ S ) measures of the sensed linear tracts in everydirection, whilst Γ q , Γ S are the averages of γ q , γ S on [0 , R ].We can, then, introduce the small parameter (cid:15) = min { η q , η S } (43)and re-scale the space variable as ξ = (cid:15) x , getting T [ q, S ]( ξ , v, ˆ v ) = q ( ξ , ˆ v ) ψ ( v ) , (44)meaning that the equilibrium is determined by the fibers distribution, and T [ q, S ]( ξ , v, ˆ v ) = (cid:20) Γ q ∇ q · ˆ v + Γ S q ( ξ , ˆ v ) ∇SS ( ξ ) · ˆ v (cid:21) ψ ( v )where Γ S := Γ S Γ S , Γ q := Γ q Γ q . Because of (27) and (44), we have that U T ( ξ ) = 0, meaning that we are in a diffusive regime,and the diffusive limits leads to the advection-diffusion equation (22). The explicit form for thezero-order macroscopic diffusion tensor is D T ( ξ ) = D (cid:90) S d − q ( ξ , ˆ v )ˆ v ⊗ ˆ v d ˆ v = D D q ( ξ ) , (45)and for the macroscopic first-order velocity is U T ( ξ ) = ¯ U (cid:90) S d − (cid:18) Γ q ∇ q · ˆ v + Γ S ∇SS ( ξ ) · ˆ v q ( ξ , ˆ v ) (cid:19) ˆ v d ˆ v = ¯ U Γ q (cid:90) S d − ( ∇ q · ˆ v ) ˆ v d ˆ v + ¯ U Γ S ∇SS (cid:90) S d − ˆ v ⊗ ˆ v q ( ξ , ˆ v ) d ˆ v = ¯ U (cid:20) Γ q ∇ · D q + Γ S D q ∇SS (cid:21) . (46)11herefore, the diffusion-advection equation (22) reads (dropping the dependencies) ∂∂τ ρ + ∇ · (cid:2)(cid:0) χ S D q ∇S + χ q ∇ · D q (cid:1) ρ (cid:3) = ∇ · (cid:20) µ ∇ · (cid:0) D D q ρ (cid:1)(cid:21) , (47)where χ S ( ξ ) := ¯ U Γ S S ( ξ ) , χ q := ¯ U Γ q (48)are the sensitivities. The diffusion represented by the motility tensor of the cells (45) only dependson the fibers distribution, while the advective term has two contributions differently weighted bythe sensitivities (48). We remark that, in this regime, we obtain the same macroscopic behaviorpostulated by Keller and Segel [32], with the logarithmic chemotactic sensitivity χ S given in(48). The term D q ∇S depends on both the fibers distribution and the chemotactic field; it nevervanishes if ∇S is not the null vector, since it may be proved that D q is invertible. In the caseof randomly disposed fibers, corresponding to the isotropic case, i . e . , when D q is proportional tothe identity matrix, then D q ∇S is parallel to ∇S , that, thus, represents the anisotropy direction.On the other hand, when D q is anisotropic, if ∇S is not parallel to the eigenvector correspondingto the highest eigenvalue of D q , then the migration does not follow the dominant direction of thefibers, but rather its projection on ∇S . Moreover, the second contribution in the drift term, i . e . , ∇ · D q , is a measure of the velocity field induced by the spatial variation of the distribution ofthe fiber directions, that determines the microscopic velocities of the cells. This term vanishes ifthe fibers distribution is homogeneous in space. Therefore, if q is homogeneous in space, even incase of competing cues, i . e . , E q ⊥ ∇S , in general the advective term U T does not vanish, whilein case of cooperating cues, i . e ., ∇S is an eigenvector of D q with eigenvalue D ∇S , migration is indirection ∇S with a kinetic factor χ S D ∇S . In intermediate scenarios, migration happens in theprojection D q ∇S , but, if q is not homogeneous, the dynamics is more complex and, even in caseof cooperation, we cannot conclude anything about additivity effects. Case iii ) In this case, we can only expand with Taylor series the chemoattractant, as in (30),and the turning kernel (36) may be approximated as T [ q, S ]( x , v, ˆ v ) = c ( x ) (cid:104) S ( x ) Γ S (cid:90) R + γ q ( λ ) q ( x + λ ˆ v , ˆ v ) dλ + Γ S ( ∇S · ˆ v ) (cid:90) R + γ q ( λ ) q ( x + λ ˆ v , ˆ v ) dλ (cid:105) ψ ( v )(49)where we neglected higher order terms in λ . Here, the normalization coefficient reduces to c ( x ) = 1Γ S Γ q S ( x ) . In this case we may choose (cid:15) = min (cid:26) η q , η S (cid:27) , and, re-scaling the space variable as (14), we get T [ q, S ]( ξ , v, ˆ v ) = 1Γ q (cid:90) R + γ q ( λ ) q ( ξ + λ ˆ v , ˆ v ) dλ ψ ( v ) (50)and T [ q, S ]( ξ , v, ˆ v ) = Γ S Γ q (cid:18) ∇SS · ˆ v (cid:19) (cid:90) R + γ q ( λ ) q ( ξ + λ ˆ v , ˆ v ) dλ ψ ( v ) . Equation (50) indicates that the equilibrium distribution is a non-local average of the fibers dis-tribution according to the sensing kernel γ q and normalized by the measure of the sensed lineartract Γ q over the direction ˆ v . Its average is U T ( ξ ) = ¯ U Γ q (cid:90) R + γ q ( λ ) E q ( ξ + λ ˆ v ) dλ ξ + λ ˆ v ∈ Ω and (27) holds true. Therefore, we perform the diffusive limit thatleads to (22) with D T ( ξ ) = D (cid:90) S d − q (cid:90) R + γ q ( λ ) q ( ξ + λ ˆ v , ˆ v ) dλ ˆ v ⊗ ˆ v d ˆ v . Let us now define D λq ( ξ ) = (cid:90) S d − q ( ξ + λ ˆ v , ˆ v ) ˆ v ⊗ ˆ v d ˆ v , (51)that, for each point ξ , is the diffusion tensor of the fibers on a circle of radius λ , and¯ D q = 1Γ q (cid:90) R + γ q ( λ ) D λq dλ , (52)that is a weighted diffusion tensor of the fibers in the whole neighborhood sensed by the cells, sothat D T ( ξ ) = D ¯ D q ( ξ ) (53)and U T ( ξ ) = ¯ U c ( ξ ) (cid:90) S d − (cid:32) Γ S ( ∇S · ˆ v ) (cid:90) R + γ q ( λ ) q ( ξ + λ ˆ v , ˆ v ) dλ (cid:33) ˆ v d ˆ v = ¯ U c ( ξ ) Γ S ∇S (cid:90) R + γ q ( λ ) (cid:90) S d − ˆ v ⊗ ˆ v q ( ξ + λ ˆ v , ˆ v ) d ˆ v dλ == ¯ U Γ S ¯ D q ( ξ ) ∇SS ( ξ ) = χ S ( ξ ) ¯ D q ( ξ ) ∇S . (54)We have defined the chemotactic sensitivity as χ S ( ξ ) := ¯ U Γ S S ( ξ ) , that is a function of the chemical alone, as it is the cue inducing a diffusive behavior. Here, theadvection velocity is related to a non-local average of the diffusion tensor of the fibers ¯ D q pro-jected on ∇S , and it cannot be decomposed into two contributions because of the large size of thecell with respect to the spatial variability of the fibers distribution. Therefore, in this case theadditivity effect of the two cues is not evident and the possible scenarios are many more. Remark
If we consider γ q = δ ( λ −
0) we obtain a local sensing of fibers. Without chemotaxis wewould have the classical model for contact guidance [27], that gives rise, at the macroscopic level,to a fully anisotropic diffusive equation. The presence of a non-local chemoattractant, even when
R < l S , gives rise to a drift correction term proportional to D q ∇S . Case iv ) The last case allows only for the Taylor expansion of the distribution function q , as in(32). Therefore, the turning kernel may be approximated as T [ q, S ]( x , v, ˆ v ) = (cid:104) c ( x ) Γ q q ( x , ˆ v ) (cid:90) R + γ S ( λ ) S ( x + λ ˆ v ) dλ + c ( x )Γ q ( ∇ q · ˆ v ) (cid:90) R + γ S ( λ ) S ( x + λ ˆ v ) dλ (cid:105) ψ ( v )(55)where c ( x ) − := 2 (cid:90) S d − Γ q q ( x , ˆ v ) (cid:90) R + γ S ( λ ) S ( x + λ ˆ v ) dλ d ˆ v c ( x ) − := 2 (cid:90) S d − Γ q ( ∇ q · ˆ v ) (cid:90) R + γ S ( λ ) S ( x + λ ˆ v ) dλ d ˆ v , both different from zero. In this case we may choose (cid:15) = min (cid:26) η S , η q (cid:27) and, by re-scaling (55) with (14), we get T [ q, S ] = T [ q, S ]. Hence U T ( ξ ) does not vanish in Ω, asit is given by U T ( ξ ) = ¯ U Γ q c ( ξ ) (cid:90) S d − ˆ v q ( ξ , ˆ v ) (cid:90) R + γ S ( λ ) S ( ξ + λ ˆ v ) dλ d ˆ v + ¯ U Γ q c ( ξ ) (cid:90) S d − ˆ v ⊗ ˆ v ∇ q (cid:90) R + γ S ( λ ) S ( ξ + λ ˆ v ) dλ d ˆ v , (56)and the macroscopic equation is given by (23). The mean velocity (56) is a linear combination ofa non-local measure of the chemoattractant S over the fibers network and a non-local measure of S weighted by the directional average of the spatial variability of the fiber direction. Remark
If we consider a local sensing for the chemoattractant, i . e . γ S = δ ( λ − ∇ · D q , and the measure of S does not affect thechoice of the direction. In this case, if ∇ q vanishes, the model reduces to a fully anisotropicdiffusive equation [27]. Concerning the non-local dependent sensing case (5)-(9) with (37), we recall the expression of thetransition probability T [ q, S ]( x , v, ˆ v ) = c ( x ) (cid:90) R + γ ( λ ) S ( x + λ ˆ v ) q ( x + λ ˆ v , ˆ v ) dλ ψ ( v ) , with c ( x ) := (cid:90) S d − (cid:90) R + γ ( λ ) S ( x + λ ˆ v ) q ( x + λ ˆ v , ˆ v ) dλ . The macroscopic velocity is here given by U T ( x ) = c ( x ) ¯ U (cid:90) S d − ˆ v (cid:90) R + γ ( λ ) S ( x + λ ˆ v ) q ( x + λ ˆ v , ˆ v ) dλ d ˆ v . (57)The macroscopic limits can be performed as in the previous section and the choice of the parameter (cid:15) will be the same for the cases i ) − iv ), since it does not depend on the kind of model (independentor dependent sensing), but only on η S and η q . Case i ) In this case we cannot consider the expansions (32) and (30), and, thus, we cannotexpand the turning kernel, whose non vanishing average is given by (57). Therefore, we performa hyperbolic limit leading to (23) with macroscopic velocity (57).14 ase ii ) When, instead, the maximum sensing radius R is smaller than both the characteristiclengths, we may consider the positive expansions (32) and (30) and substitute them in (37).Neglecting the higher order terms in λ , we get the approximation T [ q, S ]( x , v, ˆ v ) = c ( x ) (cid:104) S ( x ) Γ q ( x , ˆ v ) + S ( x ) Γ ∇ q · ˆ v + Γ q ( x , ˆ v ) ∇S · ˆ v (cid:105) ψ ( v ) (58)with c ( x ) = 1 S ( x ) Γ and Γ i := (cid:90) R λ i γ ( λ ) dλ , i = 0 , . Re-scaling the space variable as in (14), we find T [ q, S ]( ξ , v, ˆ v ) = q ( ξ , ˆ v ) ψ ( v )and T [ q, S ]( ξ , v, ˆ v ) = Γ (cid:104) ∇ q · ˆ v + q ( ξ , ˆ v ) ∇SS · ˆ v (cid:105) ψ ( v )with Γ := Γ Γ . Therefore, U T ( ξ ) = 0, because of (27), and we can perform a diffusive scaling that leads to thezero-order macroscopic diffusion tensor D T ( ξ ) = D D q ( ξ ) , (59)and to the macroscopic first-order velocity U T ( ξ ) = ¯ U Γ ∇ · D q ( ξ ) + ¯ U Γ D q ( ξ ) ∇SS . (60)The macroscopic advection-diffusion equation (22) now reads (dropping the dependencies) ∂∂τ ρ + ∇ · (cid:20) χ (cid:18) ∇ · D q + D q ∇SS (cid:19) ρ (cid:21) = ∇ · (cid:20) µ ∇ · (cid:0) D D q ρ (cid:1)(cid:21) (61)where χ := ¯ U Γ . Similar considerations to the case ii ) of the non-local independent sensing model may be done,except that there is a unique sensitivity χ that weights equally the two contributions to theadvection term (60). Case iii ) In this case, we expand only the chemoattractant S ( x + λ ˆ v ), as in (30), and the turningkernel (37) can be approximated as T [ q, S ]( x , v, ˆ v ) = c ( x ) (cid:104) S ( x ) (cid:90) R + γ ( λ ) q ( x + λ ˆ v , ˆ v ) dλ + ( ∇S · ˆ v ) (cid:90) R + λ γ ( λ ) q ( x + λ ˆ v , ˆ v ) dλ (cid:105) ψ ( v )(62)with c ( x ) := 1Γ S ( x ) . Re-scaling the space variable as in (14), we find T [ q, S ]( ξ , v, ˆ v ) = 1Γ (cid:90) R + γ ( λ ) q ( ξ + λ ˆ v , ˆ v ) dλ ψ ( v ) , T [ q, S ]( ξ , v, ˆ v ) = 1Γ (cid:18) ∇SS · ˆ v (cid:19) (cid:90) R + λ γ ( λ ) q ( ξ + λ ˆ v , ˆ v ) dλ ψ ( v ) . The macroscopic velocity of zero order is then U T ( ξ ) = ¯ U Γ (cid:90) S d − (cid:90) R + γ ( λ ) q ( ξ + λ ˆ v , ˆ v ) dλ ˆ v d ˆ v , (63)and, again, it vanishes because of ξ + λ ˆ v ∈ Ω and (27). Therefore, the macroscopic diffusion-advection equation is given by (22) with D T ( ξ ) = D Γ (cid:90) R + D λq ( ξ ) γ ( λ ) dλ = D ¯ D q (64)and U T ( ξ ) = ¯ U Γ (cid:90) R + λ D λq ( ξ ) γ ( λ ) dλ ∇SS ( ξ ) = ¯ U ¯ D q ( ξ ) ∇SS ( ξ ) , (65)where we defined ¯ D q ( ξ ) = 1Γ (cid:90) R + λ D λq ( ξ ) γ ( λ ) dλ (66)as an average of the weighted diffusion tensor of the fibers in the whole neighborhood sensed bythe cells, differently form the case iii ) of the non-local independent model. Case iv ) In this case, again, we can only consider the positive approximation (32), and thetransition probability rewrites as T [ q, S ]( x , v, ˆ v ) = (cid:104) c ( x ) q ( x , ˆ v ) (cid:90) R + γ ( λ ) S ( x + λ ˆ v ) dλ + c ( x ) ∇ q · ˆ v (cid:90) R + λ γ ( λ ) S ( x + λ ˆ v ) dλ (cid:105) ψ ( v )(67)where c ( x ) − := 2 (cid:90) S d − q ( x , ˆ v ) (cid:90) R + γ ( λ ) S ( x + λ ˆ v ) dλ d ˆ v and c ( x ) − := 2 (cid:90) S d − ( ∇ q · ˆ v ) (cid:90) R + λ γ ( λ ) S ( x + λ ˆ v ) dλ d ˆ v , both different from zero. As before, by re-scaling (67) with (14), we get T [ q, S ] = T [ q, S ] and wehave that the average velocity U T = U T (cid:54) = 0. In particular, it is given by U T ( ξ ) : = ¯ Uc ( ξ ) (cid:90) S d − ˆ v q ( ξ , ˆ v ) (cid:90) R + γ ( λ ) S ( ξ + λ ˆ v ) dλ d ˆ v + ¯ Uc ( ξ ) (cid:90) S d − ˆ v ⊗ ˆ v ∇ q ( ξ , ˆ v ) (cid:90) R + λ γ ( λ ) S ( ξ + λ ˆ v ) dλ d ˆ v (68)and, thus, we perform a hyperbolic limit leading to (23). The mean velocity (68) is a linearcombination of a non-local measure of the chemoattractant S over the fibers network and a non-local average of S weighted by the directional average of the spatial variability of the fiber direction.16 ase non-local independent sensing (5)-(9)-(36) non-local dependent sensing (5)-(9)-(37) i ) drift dominated drift dominated U T = c ¯ U (cid:90) S d − ˆ v (cid:90) R γ S ( λ ) S ( ξ + λ ˆ v ) dλ (cid:90) R γ q ( λ ) q ( ξ + λ ˆ v , ˆ v ) dλd ˆ v U T = c ¯ U (cid:90) S d − ˆ v (cid:90) R γ ( λ ) S ( ξ + λ ˆ v ) q ( ξ + λ ˆ v , ˆ v ) dλ d ˆ v ii ) drift-diffusion drift-diffusion D T = D D q D T = D D q U T = ¯ U (cid:20) Γ q ∇ · D q + Γ S D q ∇SS (cid:21) U T = ¯ U Γ (cid:20) ∇ · D q + D q ∇SS (cid:21) iii ) drift-diffusion drift-diffusion D T = D ¯ D q D T = D ¯ D q U T = ¯ U Γ S ¯ D q ∇SS U T = ¯ U ¯ D q ∇SS iv ) drift dominated drift dominated U T = ¯ U Γ q c (cid:90) S d − ˆ v q (cid:90) R γ S ( λ ) S ( ξ + λ ˆ v ) dλ d ˆ v U T := ¯ Uc (cid:90) S d − ˆ v q (cid:90) R γ ( λ ) S ( ξ + λ ˆ v ) dλ d ˆ v + ¯ U Γ q c (cid:90) S d − ˆ v ⊗ ˆ v ∇ q (cid:90) R γ S ( λ ) S ( ξ + λ ˆ v ) dλ d ˆ v + ¯ Uc (cid:90) S d − ˆ v ⊗ ˆ v ∇ q (cid:90) R λ γ ( λ ) S ( ξ + λ ˆ v ) dλ d ˆ v Table 1: Summary of the models (dropping the local dependencies in ξ ). We can observe that, if γ q = γ S = γ = δ ( λ − R ), the two non-local transport models for independentand dependent sensing are the same, while, if the sensing kernels are not dirac deltas (even if γ q = γ S = γ ), the transport models are always different. Instead, at the macroscopic level, withany choice of the sensing functions the models coincide only in case ii ). In this case, in fact, themacroscopic limits are different only if γ q (cid:54) = γ S , while in the cases iii ) and iv ) they are different ifthe sensing kernel are not dirac deltas (even if γ S = γ q = γ ). The relevant difference concerns themacroscopic transport velocities (see (54) and (65) for the case iii , and (56) and (68) for the case iv ). In fact, in the cases iii ) and iv ), for the non-local dependent sensing model, as only one cueis considered non-locally and both cues are averaged with the same sensing function γ , we have aweighted average on λ of the non-local quantities, that results in the weighted averages (65) andthe second term of (68). These remarks are summarized in Table 2. γ q = γ S = γ = δ γ q = γ S = γ (cid:54) = δ γ q (cid:54) = γ S Meso models (5)-(9)-(36) and (5)-(9)-(37) = (cid:54) = (cid:54) =Macro models case i ) = (cid:54) = (cid:54) =Macro models case ii ) = = (cid:54) =Macro models case iii ) = (cid:54) = (cid:54) =Macro models case iv ) = (cid:54) = (cid:54) = Table 2: Summary of the comparison of the models for different choices of the sensing functions.= indicates the cases in which the models coincide, while (cid:54) = the ones in which the models aredifferent. 17
Numerical simulations
We shall now propose two-dimensional numerical simulations in order to illustrate the behavior ofthe kinetic transport models for non-local independent sensing and non-local dependent sensing.In particular, we shall integrate numerically the transport equation as in [39] and, then, we shallcompute the macroscopic density (1). Concerning the fibers network, a classical used distributionis the Von Mises distribution [42]˜ q ( x , ˆ v ) = 12 πI ( k ( x )) e k ( x ) u ( x ) · ˆ v where I ν ( k ) is the modified Bessel function of first kind of order ν and u ( x ) = (cos( θ q ( x )) , sin( θ q ( x ))) . It can be proved that E ˜ q ( x ) = u ( x ) [28], and, therefore, θ q ( x ) is the mean direction in thespace [0 , π ) of the fibers located at point x . As we are dealing with cell migrating on a non-polarized network of fibers, we shall consider the symmetric version, namely the Bimodal VonMises distribution q ( x , ˆ v ) = 14 πI ( k ( x )) (cid:16) e k ( x ) u ( x ) · ˆ v + e − k ( x ) u ( x ) · ˆ v (cid:17) , that also satisfies Q3; its variance is [28] D q ( x ) = 12 (cid:18) − I ( k ) I ( k ) (cid:19) I + I ( k ) I ( k ) u ⊗ u , where I is the identity tensor in R × , while k and u are functions of x . Moreover, the variancein the space [0 , π ) is the scalar D q ( x ) = 12 (cid:90) π q ( θ − θ q ) dθ = (cid:18) − I ( k ) I ( k ) (cid:19) that represents the degree of alignment of the fibers at point x . As a first example, we shall present the particular case in which the sensing of q is local. Thisillustrates the effect of a second directional cue when dealing with a cell population migrating bycontact guidance and evaluating the local alignment of the fibers over a non-polarized network.Formally, we are dealing with (36) in which γ q = δ ( λ − q = { x = ( x, y ) ∈ Ω s.t. x ≤ x ≤ x } (69)with x = 1 . x = 3 . θ q = π/
2. In particular, for ( x, y ) ∈ Ω q , k ( x, y ) = 700, such that D q = 5 · − . In the rest ofthe domain Ω − Ω q fibers are uniformly distributed. The chemoattractant has a Gaussian profile S ( x, y ) = m S (cid:112) πσ S e − (( x, y ) − ( x S , y S )) σ S . (70)In particular, in Test 1 (see Fig. 1) we choose ( x S , y S ) = (4 , , m S = 10 , σ S = 0 .
1. The initialcondition for the cell population is a Gaussian ρ ( x, y ) = r e − (( x, y ) − ( x , y )) σ (71)18ith r = 0 . σ = 0 .
1. In this first test, the initial condition for the cell population is centeredin ( x , y ) = (2 . , . i . e . , the center of the region Ω q (see Fig. 1a). Without chemoattractant,because of the presence of highly aligned fibers, we would expect that cells diffuse anisotropicallyin the preferential direction of the fibers ± π/
2, forming the well known ellipsis [46], that representscells moving with the same probability along direction π/ − π/
2. In the present case, due tothe presence of a chemoattractant, the symmetry is broken, and, even if q describes a non-polarizedfibers network, there is a preferential sense of motion (see Fig. 1d-1f). In particular, cells migratealong the fibers in the direction identified by θ q = π/
2, corresponding to the preferential senseimposed by the presence of the chemoattractant in the upper-right corner of the domain Ω. Giventhis directional setting, the cell population dynamics is also greatly affected by the strength of thechemoattractant, that depends on m S and σ S , the degree of the alignment D q , that depends on k ( x, y ), and by the sensing radius R . Another important aspect is the sensing function γ S , thatinfluences the transient dynamics and, especially, the relaxation time. This appears to be doublein the case of a Heaviside function, since the kernel γ S doubles when computed with a Heavisidefunction instead of a Dirac delta (see also [39]). (a) Initial cell distribution (b)
Initial average polarization (c)
Center of mass: trajectory. (d) t=1.25 (e) t=3.75 (f) t=12.5 (g) t=1.25 (h) t=3.75 (i) t=12.5
Figure 1:
Test 1
Evolution of the initial distribution given in (a) for the case of local q andnon-local chemoattractant S with sensing function γ S = δ ( λ − R ). In (b), S is a Gaussian centredin (4 ,
4) and with m S = 10 and σ S = 0 .
1. The sensing radius of the cells is set to R = 0 . t = 1. Figs. (d)-(f): evolution of the macroscopic density. Figs. (g)-(i): polarizations of thecells. 19e also analyzed the average polarization of the cells at every position x , that is given bythe momentum (2). The microscopic directions of cells are initially randomly distributed andthey start from a vanishing initial speed (see Fig. 1b). Then, they start to align along thefibers and to migrate upward in the direction individuated by the angle π/
2, since cells sense thechemoattractant (see Figs. 1g-1h). Eventually when cells reach the level y = 4, the microscopicdirections polarize towards the chemoattractant (see Fig. 1i). The center of mass plotted in Fig.1c stays in the region Ω q during the migration of cells along the fibers bundle in Ω q , and it movesout of Ω q only when it reaches y = 4. The black dots are plotted every ∆ t = 1 and it is clearthat the highest acceleration happens when cells are on the bundle of fibers, while they are sloweddown when they start to move out of the fibers stripe Ω q . As a second test, we present both the non-local independent sensing model and the non-localdependent sensing model. We shall now consider a non-local sensing of the distribution of fibers.In particular, we assume fibers distributed similarly to the previous test, i . e . , fibers shall be highlyaligned in Ω q given, this time, by x = 2 . x = 2 . x, y ) ∈ Ω q , k ( x, y ) = 100, that corresponds to D q = 0 . θ q ( x, y ) = π/
2. In the region Ω − Ω q fibers areuniformly distributed. The initial condition of the cell population is (71) with in ( x , y ) = (1 , . m S = 10 and σ S = 0 .
05. Weshall compare the dynamics of the cells in four settings:1. local fiber distribution and non-local chemoattractant, as in Test 1, i . e . , (36) with γ q = δ ( λ −
0) and γ S = δ ( λ − R );2. non-local sensing with a Dirac Delta for both q and S ; this corresponds to both (36) and(37) with γ q = γ S = γ = δ ( λ − R );3. non-local independent sensing with Heaviside sensing functions for both S and q , i . e . , (36)with γ q = γ S = H ( R − λ );4. non-local dependent sensing for q and S , dealing with (37) and γ = H ( R − λ ).Results of these simulations are shown in Fig. 2. We can observe that, in the 1-4 settings, cells startfrom (1 , . S , they crossthe aligned fibers region Ω q and climb up this region in the direction π/
2. Eventually, in all thecases, cells reach the chemoattractant, but the dynamics, as well as the transient time, is influencedby the different sensing kernels, even though the differences are not extremely appreciable, andby the local or non-local sensing strategy. Although settings 3 and 4 in Fig. 2, that are related tothe case of independent and dependent cues, respectively, do not show very strong differences, incase 3 (see Figs. 2k-2n) the tendency of going in both the direction π/
2, determined by the fibers,and π/
4, determined by the chemoattractant, appears more marked because of the independentsensing. In contrast, this behavior results the least evident in the case in which cells deal with alocal sensing of the fibers (setting 1), resulting also in a general slow down of the dynamics. i ) − iv ) We now present a comparison of the macroscopic behaviors of the cells, depending on the relationbetween R , l S and l q , i . e . , we compare the cases i ) , ii ) , iii ) and iv ). In particular, we shall dothis for the non-local independent sensing model with γ q = γ S = H ( R − λ ), as this is the casein which the transport model is different from the dependent sensing model. Additionally, theindependence of the two sensings allows to visualize more efficiently the two distinct directionaleffects (contact guidance and chemotaxis). 20 a) Initial condition for cells (b)
Initial fiber distribution (c) t=1.25 (d) t=3.75 (e) t=5 (f) t=6.25 (g) t=1.25 (h) t=3.75 (i) t=5 (j) t=6.25 (k) t=1.25 (l) t=3.75 (m) t=5 (n) t=6.25 (o) t=1.25 (p) t=3.75 (q) t=5 (r) t=6.25
Figure 2:
Test 2
Time evolution of the initial distribution given in Fig. 2a in the four settings1-4. The sensing radius of the cells is R = 0 . m S =10 , σ S = 0 .
05 and ( x S , y S ) = (4 , q and non-localchemoattractant, γ S = δ ( λ − R ). Setting 2 is represented in Figs. (g)-(j): non-local q and S withsensing functions γ q = γ S = δ ( λ − R ). Setting 3 is represented in Figs. (k)-(n): non-local q and S ,independent sensing with γ q = γ S = H ( R − λ ). Setting 4 is represented in Figs. (o)-(r): non-local q and S , dependent sensing with γ = H ( R − λ ).21e shall consider the turning kernel describing contact guidance lead by a q with mean direction θ q ( x, y ) = 3 π/ ∀ ( x, y ) ∈ Ω and coefficient k ( x, y ), modulating the strength of the alignment, givenby a gaussian distribution k ( x, y ) = m k e − (( x, y ) − ( x k , y k )) σ k (72)where ( x k , y k ) = (2 . , .
5) and σ k = 0 .
15 (Fig. 3d). This mimics the situation of fibers morealigned in the central circular region and uniformly disposed in the rest of the domain. We shallconsider different values of m k in order to obtain different values of l q : m k = 10 corresponds to l q ≈ .
031 and m k = 100 corresponds to l q ≈ . l q for aBimodal Von Mises distribution of fibers q are given in Appendix A. The chemoattractant is (70)with ( x S , y S ) = (4 . , .
5) and m S = 10. In the simulations, we shall consider three different valuesfor the variance of the chemoattractant σ S in order to obtain different values of l S : σ S = 0 .
05 thatcorresponds to l S = 0 .
002 in Fig. 3a, σ S = 0 .
25 that corresponds to l S = 0 .
055 in Fig. 3b and σ S = 1 . l S = 0 .
25 in Fig. 3c. The initial distribution of cells for all the testspresented in Figs. 4, 5, 6, 7 and 8 is given by (71) with ( x , y ) = (1 . , . r = 0 . σ = 0 . l S l q R Case η Figure0 .
002 0 . . i ) < .
25 0 . . i ) (cid:29) .
055 0 .
031 0 . ii ) > .
25 0 . . iii ) (cid:29) .
002 0 .
031 0 . iv ) < Table 3: Summary of the simulations presented in Test 3.In Fig. 4, we consider the case in which η S , η q (cid:29) i . e . , we are dealing with case i ). Themacroscopic behavior is strongly hyperbolic with macroscopic velocity given by (40). In fact, inFig. 4 we can observe that the behavior is not diffusive and the cluster of cells is quite com-pact. Moreover, when cells reach the region in which fibers are strongly aligned in the direction3 π/ π/ S with σ S = 1 . l S = 0 .
25 (see Fig. 3c).Concerning the fibers, we have m k = 100, so that l q ≈ . R = 0 . i ), but the behavior is different with respect to the previoussimulation in Fig. 4. The chemoattractant in Fig. 3c, in fact, is spread over the whole domainand, actually, the quantity l S is almost 10 times the l S considered in Fig. 3a and used for thesimulation in Fig. 4. Even though we are still in a strongly hyperbolic case and cells are guidedby the strong drift (40), as R is slightly larger then l S and l S is large, the cell cluster diffuses a bitmore in the domain. When it reaches the region of strongly aligned fibers, it starts to surroundthat region (see Figs. 5a-5c), but, as η S = 2 . O (1), some cells, that do not surround the region,are slowed down and partially tend to align along the fibers. In Fig. 5d, for instance, we have ahigh density of cells both in the strongly aligned fiber region and in the region of high density ofchemoattractant. Eventually, cells manage to overcome the area of highly aligned fibers and theytend to converge to the chemoattractant profile (see Figs. 5e-5f). Now, the the overall dynamicsis greatly affected by the fibers and, in fact, η (cid:29) ii ), since the sensing radius R = 0 .
02 is smaller than both l S = 0 .
055 and l q ≈ . a) Chemoattractant S with σ S = 0 . (b) Chemoattractant S with σ S = 0 . (c) Chemoattractant S with σ S = 1 . (d) Fibers distribution
Figure 3:
Test 3
Three different chemoattractants used for comparing models i ) − iv ). Thechemoattractant profile is given by (70) with m S = 10 and (a) σ S = 0 .
05, corresponding to l S = 0 . σ S = 0 .
25, corresponding to l S = 0 . σ S = 1 .
8, corresponding to l S = 0 .
25. The fibers distribution in sketched in (d).of the system is described by the diffusion-advection equation (47) with macroscopic velocity (46).Actually, in Fig. 6 we can observe a highly diffusive behavior, as the macroscopic density of cellshas invaded almost the half of the domain before even starting to be influenced by the fibers. Ifwe compare the same time step in Figs. 6b and 5b, we see that the cells are in both cases reachingthe fibers and feeling the region in which fibers are aligned the most. However, in Fig. 5b the cellcluster is much more compact than in Fig. 6b, where, instead, cells already occupied half of thedomain, because of diffusion, and we have high density of cells both closely to the strongly alignedfiber region and around the initial position. Therefore, cells start surrounding the central regionof strongly aligned fibers, because they already sense the chemoattractant, and, once overcomethis area, they tend to the chemoattractant profile (see Figs. 6c-6f). In particular, in the transienttime, cells accumulate the most at the sides of the region with highly aligned fibers. In this specificsetting, η > iii ), since the sensing radius R = 0 . l S = 0 .
25 but it is larger then l q ≈ . η S <
1, we have that the chemoattractant induces a strong diffusivity, but being η q >
1, the alignment of fibers strongly affects the dynamics (see Figs. 7c-7d). Comparing, inaddition, Figs. 6b and 7b, we have now that the highest cell concentration is in the mean fiberdirection θ q = 3 π/ η (cid:29) R = 0 .
02 smaller than l q ≈ . l S = 0 . a) t=1.25 (b) t=1.875 (c) t=2.5 (d) t=3.75 (e) t=5 (f) t=6.25 Figure 4:
Test 3
Case i ) with non-local q and S , sensed with an independent sensing through thekernels γ q = γ S = H ( R − λ ). S is given in Fig. 3a with m S = 10 and σ S = 0 .
05, so that l S = 0 . q has a space dependent parameter k given by (72) with m k = 100, so that l q ≈ . R = 0 . (a) t=2.5 (b) t=5 (c) t=10 (d) t=15 (e) t=22.5 (f) t=27.5 Figure 5:
Test 3
Case i ) with non-local q and S , independent and sensing with γ q = γ S = H ( R − λ ). S is given in Fig. 3c, that corresponds to l S = 0 .
25, while for the fiber distribution m k = 100, sothat l q ≈ . R = 0 . η q is smaller than 1, and they start movingin a region with randomly disposed fibers (see Fig. 8a). Then, they mainly follow the preferentialdirection π/ a) t=2.5 (b) t=5 (c) t=7.5 (d) t=10 (e) t=15 (f) t=20 Figure 6:
Test 3
Case ii ) with non-local q and S , independent and sensing with γ q = γ S = H ( R − λ ). S is given in Fig. 3b, that corresponds to l S = 0 . m k = 10, so that l q ≈ . R = 0 . (a) t=2.5 (b) t=5 (c) t=10 (d) t=20 (e) t=30 (f) t=60 Figure 7:
Test 3
Case iii ) with non-local q and S , independent and with sensing function γ q = γ S = H ( R − λ ). S is given in Figure 3c, so that l S = 0 .
25, while for the fiber distribution m k = 100,corresponding to l q ≈ . R = 0 . η S (cid:29)
1. Here chemotaxis is slightly dominating thedynamics and, in fact, η <
1. 25 a) t=1.25 (b) t=2.5 (c) t=5 (d) t=7.5 (e) t=10 (f) t=15 Figure 8:
Test 3
Case iv ) with non-local q and S , independent sensing with γ q = γ S = H ( R − λ ). S is given in Fig. 3a, that corresponds to l S = 0 . m k = 10, so that l q ≈ . R = 0 . We now consider a domain Ω divided in several regions, each of them characterized by a differentaverage direction of the fibers. In particular, we shall do this in the case of independent sensingmodel with γ q = γ S = H ( R − λ ), as for Test 3; the independence of the two sensings, in fact, allowsto visualize more efficiently the two distinct directional effects. As first scenario, we shall considerthe domain schematized in Fig. 9a; in each subdomain we have k ( x, y ) = 50, that correspondsto D q = 0 . r = 0 .
1, while the chemoattractant has a gaussian profile (70) centered in ( x S , y S ) = (4 , m S = 10 and σ S = 0 .
5, as shown in Fig. 9b. We observe that cells do not migrate collectivelytowards the chemoattractant, but they divide into two main separated clusters (see Figs. 9f -9h): in fact, although the sensing radius R = 0 . k ( x, y ) = 50. The initial condition of the cell population is (71) with ( x , y ) =(4 , .
5) and r = 0 .
1, while the chemoattractant has a gaussian profile (70) centered in ( x S , y S ) =(2 , .
5) with m S = 10 and σ S = 0 .
05, as shown in Fig. 10c and 10b, respectively. We observe thatcells do not migrate directly towards the chemoattractant, as they sense the heterogeneous fibrousenvironment and, consequently, adapt their migration to it. In particular, cells that are able toreach and sense the isotropic subdomain where the fibers are uniformly distributed (defined by1 ≤ x ≤ ≤ y ≤ ≤ x ≤ ≤ y ≤
2, they follow the direction offiber alignment, that is π/
4, perpendicular to the favorable direction imposed by S . However,the sensing radius R = 0 . π/ ≤ y ≤ π/ ≤ y ≤
4, to reach the chemoattractant.26 a) Fibers distribution. (b)
Chemoattractant S . (c) Initial condition for the cell. (d) t=0.04 (e) t=1.6 (f) t=2.904 (g) t=18.4 (h) t=44 (i) t=67.2
Figure 9:
Test 4
Migration of cells in an heterogenous domain as illustrated in (a). The sensingradius of the cells is R = 0 .
8. The chemoattractant (b) is (70) with m S = 10 and σ S = 0 .
5. Theinitial cell profile (c) evolves in time as illustrated in (d)-(i).
We have proposed a kinetic model for describing cell migration in a multi-cue environment. Inparticular, in the same spirit as [39], we have considered that cells, as they can extend protrusionsup to several cell diameters, perform a non-local sensing of the environment up to a distance R (named the sensing radius) from its nucleus. In the present model, there are two guidance cuesaffecting the polarization, and, therefore, the direction of motion of the cells: contact guidance,that is a bi-directional cue, and a chemical gradient, that is a mono-directional cue. We remarkthat for the first time in this work a non-local sensing in the physical space of the mesoscopicdistribution of fibers is considered. In particular, we introduced two classes of models: in the firstone, the cells perform an independent sensing of the fibers and of the chemical in its neighborhood,while in the second class of models the cells average the chemical and the fibers with the samesensing kernel.In the two cases, a particular attention was devoted to the identification of the proper macro-scopic limit according to the properties of the turning operator. We detected two parameters, η q and η S , that measure the relation between the cell sensing radius and the characteristic lengthsof variation − l S and l q − of the two cues, and discriminate between a diffusion-driven regimewith an advective correction and a drift-driven regime. In particular, when the sensing radiusdoes not exceed the characteristic length of the chemoattractant, the bi-directional nature of the27 a) Fibers distribution. (b)
Chemoattractant S . (c) Cells initial condition. (d) t=0.5 (e) t=1 (f) t=1.5(g) t=2.5 (h) t=3.5 (i) t=4.5
Figure 10:
Test 4
Migration of cell in an heterogenous domain as illustrated in (a). The sensingradius of the cells is R = 0 .
7. The chemoattractant (b) is (70) with m S = 10 and σ S = 0 .
05. Theinitial cell profile (c) evolves in time as illustrated in (d)-(i).fibers allows for a diffusive regime; otherwise the hyperbolic scaling leads to macroscopic drift.A common feature in the different cases is the dependency of the macroscopic velocity on boththe fibers network and the chemoattractant. This aspect enhances the non-trivial influence ofcontact guidance on the cell drift, although we considered a non polarized fibers network. Thisinterdependence is in accordance with the model proposed in [61]. Moreover, in absence of achemoattractant, this impact on the drift term could persist for spatial heterogenous fiber distri-butions. This is in accordance to what is observed in [27] and it represents a step forward withrespect to [61], in which the drift is a function of contact guidance only through to the presenceof a chemical gradient, i . e . , without chemoattractant there will be no drift.The numerical simulations of the transport model pointed out the main features characterizingthe two classes of models and the possible scenarios that they are able to capture. We observedthat the presence of two cues influencing cell polarization, even when the fibers are sensed locally,ensures a preferential sense of motion for cells laying on regions of highly aligned non-orientedfibers. Test 3 allowed to show the importance of deriving the macroscopic equations from an un-derlying microscopic dynamics and in the appropriate regime: a directly postulated drift-diffusionequation would not capture the exact dynamics in all the possible regimes. The competitive orcollaborative effects of the cues depend, in a first instance, on the angle between their relativeorientations, i . e . , the direction of fiber alignment θ q and the gradient of the chemoattractant.28oreover, especially for the cases of competitive cues, determining which one is the dominantcue depends on their relative strengths, in terms of both concentration and intensity (degree ofalignment of the fiber k ( x ) or steepness of the chemoattractive gradient). We introduced theparameter η = l S /l q that, independently on the cell size or its sensing capability, quantifies therelative contribution of guidance to chemotaxis and provide a first separation between the casesof fiber-dominating and chemotaxis-dominating dynamics ( η (cid:29) η (cid:28)
1, respectively). Thepresented framework also allows for the direct calculation of parameters that can be used to quan-tify directed cell migration and to set its efficiency, like, for instance, mean square displacement,persistence time, directional persistence and mean speed [45].Additionally, the non locality brings an further level of detail to the model, allowing to obtaindifferent macroscopic behaviour depending on the characteristics of the two sensing. In fact, wedid not observe strong differences between the independent and the dependent sensing models,when we assume in the former the same sensing kernel for fibers and chemoattractant, i . e . , when γ q = γ S . However, if there are biological observations sustaining the possibility that a cell mightimplement different strategies for sensing the underlying fibers network and the chemoattractant,it would be possible to use the proposed model, in its independent sensing version, to investigatethis scenario and to compare the possible outcomes of this sensing approach with the case of aunique and common sensing strategy.Potentially, the case of competitive cues, combined with the non-local aspect of the model,could lead to interesting further analysis. As observed in the last numerical tests, the combinationof heterogenous landscapes of fiber with chemoattractive agents show how the cell density candivide and cross the domain using different migration strategies. This leads to natural questionsabout the deeper mechanisms leading the competition between the two cues, considering, forinstance, the possible role of cell adhesion in recovering collective migration.We remark that, even if simulations were performed in a two dimensional setting, the trans-port model (and its macroscopic limits, as a consequence) is formulated in a general d-dimensionalsetting. Hence, a possible future development is to perform simulations in the three dimensionalcase, that would be much more realistic for mimicking in-vivo migration of cells in the extracellularmatrix. Moreover, the model that we proposed may be adapted to describe other directional cuesthat might describe, among others, haptotactic, durotactic or electrotactic mechanisms. Further-more, in the same spirit as in [40] we could enrich this model with a non-constant sensing-radius,as it may vary according to the spatial and directional variability of the external guidance cues.Lastly, this study was restricted to the case in which the cues affect only cell polarization, consid-ering a uniform distribution of the speeds. However, similarly to what is done in [39, 40], it maybe modified to model a multi-cue environment in which one of the signals affects also the speedof the cells. A Estimation of l q Let us consider the fiber density distribution q ( x , ˆ v ) defined by a bimodal Von Mises Fisher q ( x , ˆ v ) = 14 πI ( k ( x )) (cid:16) e k ( x ) u · ˆ v + e − k ( x ) u · ˆ v (cid:17) , where k ( x ) ∈ C (Ω) and I ν ( k ( x )) denotes the modified Bessel function of first kind of order ν .We now want to give an estimation for the range of variability of the characteristic length l q ,defined as: l q := 1max x ∈ Ω max ˆ v ∈ S d − |∇ q · ˆ v | q . ∂I ∂k = I ( k ) I ( k ) , we have that ∇ q = (cid:0) e k ( x ) u · ˆ v − e − k ( x ) u · ˆ v (cid:1) πI ( k ( x )) ∇ k ( u · ˆ v ) − (cid:0) e k ( x ) u · ˆ v + e − k ( x ) u · ˆ v (cid:1) πI ( k ( x )) ∂I ∂k ∇ k == (cid:0) e k ( x ) u · ˆ v − e − k ( x ) u · ˆ v (cid:1) πI ( k ( x )) ∇ k ( u · ˆ v ) − (cid:0) e k ( x ) u · ˆ v + e − k ( x ) u · ˆ v (cid:1) πI ( k ( x )) I ( k ( x )) I ( k ( x )) ∇ k Since q ( x , ˆ v ) >
0, we have: ∇ q · ˆ v q = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) e k ( x ) u · ˆ v − e − k ( x ) u · ˆ v (cid:1)(cid:0) e k ( x ) u · ˆ v + e − k ( x ) u · ˆ v (cid:1) ( u · ˆ v ) − I ( k ( x )) I ( k ( x )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ||∇ k || cos( ∇ k · ˆ v )where || · || denotes the L -norm and we use the fact that || ˆ v || = 1. Therefore, (cid:12)(cid:12)(cid:12)(cid:12) ∇ q · ˆ v q (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) e k ( x ) u · ˆ v − e − k ( x ) u · ˆ v (cid:1)(cid:0) e k ( x ) u · ˆ v + e − k ( x ) u · ˆ v (cid:1) ( u · ˆ v ) − I ( k ( x )) I ( k ( x )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ||∇ k || | cos( ∇ k · ˆ v ) | Recalling that | a − b | ≤ | a | + | b | , − ≤ (cid:0) e k ( x ) u · ˆ v − e − k ( x ) u · ˆ v (cid:1)(cid:0) e k ( x ) u · ˆ v + e − k ( x ) u · ˆ v (cid:1) ≤ | cos ( · ) | ≤
1, we get (cid:12)(cid:12)(cid:12)(cid:12) ∇ q · ˆ v q (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:18) (cid:12)(cid:12)(cid:12)(cid:12) I ( k ( x )) I ( k ( x )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ||∇ k || . Considering Eq. (1.12) in [36] for ν = 1, we obtain that (cid:12)(cid:12)(cid:12)(cid:12) I I (cid:12)(cid:12)(cid:12)(cid:12) <
1, and, therefore, (cid:12)(cid:12)(cid:12)(cid:12) ∇ q · ˆ v q (cid:12)(cid:12)(cid:12)(cid:12) < ||∇ k || that implies max x ∈ Ω max ˆ v ∈ S d − (cid:12)(cid:12)(cid:12)(cid:12) ∇ q · ˆ v q (cid:12)(cid:12)(cid:12)(cid:12) < x ∈ Ω ||∇ k || . This translates into l q ≥
12 max x ∈ Ω ||∇ k || . (73)In particular, if there exists x such that ∇ k ( x ) · ˆ v = 1 and, at the same time, also satisfies ∇ k ( x ) (cid:107) u , then (73) is true with the equal sign. In particular, for the symmetry of (72) and (70)we shall consider l q ≈
12 max x ∈ Ω ||∇ k || . Acknowledgments
The authors would like to thank Prof. Luigi Preziosi for fruitful discus-sions and valuable comments. This work was partially supported by Istituto Nazionale di AltaMatematica, Ministry of Education, Universities and Research, through the MIUR grant Dipar-timento di Eccellenza 2018-2022, Project no. E11G18000350001, and the Scientific Reseach Pro-grammes of Relevant National Interest project n. 2017KL4EF3. NL also acknowledges Compagniadi San Paolo. This research was also partially supported by the Basque Government through theBERC 2018- 2021 program and by the Spanish State Research Agency through BCAM SeveroOchoa excellence accreditation SEV-2017-0718. MC has received funding from the EuropeanUnions Horizon 2020 research and innovation programme under the Marie Skodowska- Curie grantagreement No. 713673. The project that gave rise to these results received the support of a fellow-ship from la Caixa Foundation (ID 100010434). The fellowship code is LCF/BQ/IN17/11620056.30 eferences [1] Y. Azimzade, A. A. Saberi, and M. Sahimi. Regulation of migration of chemotactic tumorcells by the spatial distribution of collagen fiber orientation.
Phys. Rev. E , 99:062414, 2019.[2] N. Bellomo, A. Bellouquid, J. Nieto, and J. Soler. Multicellular biological growing sys-tems: Hyperbolic limits towards macroscopic description.
Math. Mod. Meth. Appl. S. ,17(supp01):1675–1692, 2007.[3] N. Bellomo, A. Bellouquid, Y. Tao, and M. Winkler. Toward a mathematical theory ofkeller–segel models of pattern formation in biological tissues.
Math. Mod. Meth. Appl. S. ,25(09):1663–1763, 2015.[4] H. C. Berg.
Random Walks in Biology . Princeton University Press, revised edition, 1983.[5] H. C. Berg and E. M. Purcell. Physics of chemoreception.
Biophys. J. , 20(2):193–219, 1977.[6] M. Bisi, J. A. Carrillo, and B. Lods. Equilibrium solution to the inelastic boltzmann equationdriven by a particle bath.
J. Stat. Phys. , 133(5):841–870, 2008.[7] S. M. Block, J. E. Segall, and H. C. Berg. Adaptation kinetics in bacterial chemotaxis.
J.Bacteriol. Res. , 154(1):312–323, 1983.[8] B. A. Bromberek, P. A. J. Enever, D. I. Shreiber, M. D. Caldwell, and R. T. Tranquillo.Macrophages influence a competition of contact guidance and chemotaxis for fibroblast align-ment in a fibrin gel coculture assay.
Exp. Cell Res. , 275(2):230–242, 2002.[9] F. A. C. C. Chalub, P. A. Markowich, B. Perthame, and C. Schmeiser. Kinetic models forchemotaxis and their drift-diffusion limits.
Monatsh. Math. , 142(1):123–141, 2004.[10] A. Chauviere, T. Hillen, and L. Preziosi. Modeling cell movement in anisotropic and hetero-geneous network tissues.
Netw. Heterog. Media , 2(2):333–351, 2007.[11] A. Chauviere, T. Hillen, and L. Preziosi. Modeling the motion of a cell population in theextracellular matrix.
Discrete Cont. Dyn.-B , 2007(Supplemental volume):250–259, 2007.[12] L. Chen, K. J. Painter, C. Surulescu, and A. Zhigun. Mathematical models for cell migration:a nonlocal perspective. arXiv preprint arXiv:1911.05200 , 2019.[13] A. Colombi, M. Scianna, and L. Preziosi. Coherent modelling switch between pointwise anddistributed representations of cell aggregates.
J. Math. Biol. , 74(4):783–808, 2017.[14] A. Colombi, M. Scianna, and A. Tosin. Differentiated cell behavior: a multiscale approachusing measure theory.
J. Math. Biol. , 71:1049–1079, 2015.[15] M. Conte, L. Gerardo-Giorda, and M. Groppi. Glioma invasion and its interplay with nervoustissue and therapy: A multiscale model.
J. Theo. Biol. , 486:110088, 2020.[16] E. Di Costanzo, M. Menci, E. Messina, R. Natalini, and A. Vecchio. A hybrid model ofcollective motion of discrete particles under alignment and continuum chemotaxis.
DiscreteCont. Dyn.-B , 25:443–472, 2020.[17] R. Dickinson and R. T. Tranquillo. Stochastic model of biased cell migration based on bindingfluctuations of adhesion receptors.
J. Math. Biol. , 19:563–600, 1991.[18] R. B. Dickinson. A generalized transport model for biased cell migration in an anisotropicenvironment.
J. Math. Biol. , 40(2):97–135, 2000.[19] R. Eftimie. Hyperbolic and kinetic models for self-organized biological aggregations andmovement: a brief review.
J. Math. Biol. , 65(1):35–75, 2012.3120] C. Engwer, T. Hillen, M. Knappitsch, and C. Surulescu. Glioma follow white matter tracts:a multiscale dti-based model.
J. Math. Biol. , 71(3):551–582, 2015.[21] C. Engwer, M. Knappitsch, and C. Surulescu. A multiscale model for glioma spread includingcell-tissue interactions and proliferation.
Math. Biosci. Eng. , 13:443–460, 2016.[22] C. Engwer, C. Stinner, and C. Surulescu. On a structured multiscale model for acid-mediatedtumor invasion: The effects of adhesion and proliferation.
Math. Mod. Meth. Appl. S. ,27:1355–1390, 2017.[23] F. Filbet, P. Lauren¸cot, and B. Perthame. Derivation of hyperbolic models for chemosensitivemovement.
J. Math. Biol. , 50(2):189–207, 2005.[24] P. Friedl. Prespecification and plasticity: shifting mechanisms of cell migration.
Curr. Opin.Cell Biol. , 16:1423, 2004.[25] P. Friedl and E.-B. Brocker. The biology of cell locomotion within three dimensional extra-cellular matrix.
Cell Mol Life Sci. , 57:41–64, 2000.[26] R. Gininait, R. E. Baker, P. M. Kulesa, and P. K. Maini. Modelling collective cell migration:neural crest as a model paradigm.
J. Math. Biol. , 80:481–504, 2019.[27] T. Hillen. M5 mesoscopic and macroscopic models for mesenchymal motion.
J. Math. Biol. ,53(4):585–616, 2006.[28] T. Hillen, A. Murtha, K. J. Painter, and A. Swan. Moments of the von mises and fischerdistributions and applications.
Math. Biosci. Eng. , 14(3):673–694, 2017.[29] T. Hillen and H. G. Othmer. The diffusion limit of transport equations derived from velocity-jump processes.
SIAM J. Appl. Math. , 61:751–775, 2000.[30] T. Hillen and K. J. Painter. A user’s guide to pde models for chemotaxis.
J. Math. Biol. ,58(1):183–217, 2008.[31] J. Johnson, M. O. Nowicki, C. H. Lee, E. A. Chiocca, M. S. Viapiano, S. E. Lawler, and J. JLannutti. Quantitative analysis of complex glioma cell migration on electrospun polycapro-lactone using time-lapse microscopy.
Tissue Eng. Part C-Me , 15(4):531–540, 2009.[32] E. F. Keller and L. A. Segel. Initiation of slime mold aggregation viewed as an instability.
J.Theo. Biol. , 26(3):399–415, 1970.[33] P. J. Kevin. Mathematical models for chemotaxis and their applications in self-organisationphenomena.
J. Theor. Biol. , 481:162–182, 2019.[34] P. J. Kevin, P. K. Maini, and H. G. Othmer. Development and applications of a model forcellular response to multiple chemotactic cues.
J. Math. Biol. , 41(4):285–314, 2000.[35] N. Kolbe, N. Sfakianakis, C. Stinner, C. Surulescu, and J. Lenz. Modeling multiple taxis:tumor invasion with phenotypic heterogeneity, haptotaxis, and unilateral interspecies repel-lence. arXiv preprint arXiv:2005.01444 , 2020.[36] A. Laforgia and P. Natalini. Some inequalities for modified bessel functions.
J. Inequal. Appl. ,2010(1):253035, 2010.[37] L. Lara and I. Schneider. Directed cell migration in multi-cue environments.
Integr. Biol. ,5(11):1306–1323, 2013.[38] B. Lods. Semigroup generation propertiesof streaming operators with noncontractive bound-ary conditions.
Math. Comput. Model. , 42:1441–1462, 2005.3239] N. Loy and L. Preziosi. Kinetic models with non-local sensing determining cell polarizationand speed according to independent cues.
J. Math. Biol. , 80:373–421, 2019.[40] N. Loy and L. Preziosi. Modelling physical limits of migration by a kinetic model withnon-local sensing.
J. Math. Biol. , 2019. In Press.[41] G. Maheshwari, A. Wells, L. G. Griffith, and D. A. Lauffenburger. Biophysical integrationof effects of epidermal growth factor and fibronectin on fibroblast migration.
Biophys. J. ,76(5):2814–2823, 1999.[42] K. V. Mardia and P. E. Jupp.
Directional statistics , volume 494. John Wiley & Sons, 2009.[43] H. Othmer and T. Hillen. The diffusion limit of transport equations ii: Chemotaxis equations.
SIAM J. Appl. Math. , 62:1222–1250, 2002.[44] H. Othmer and A. Stevens. Aggregation, blowup, and collapse: The ABC’s of taxis inreinforced random walks.
SIAM J. Appl. Math. , 57:1044–1081, 2001.[45] H. G. Othmer, S. R. Dunbar, and W. Alt. Models of dispersal in biological systems.
J. Math.Biol. , 26(3):263–298, 1988.[46] K. J. Painter. Modelling cell migration strategies in the extracellular matrix.
J. Math. Biol. ,58(4):511–543, 2008.[47] K. J. Painter and T. Hillen.
Transport and anisotropic diffusion models for movement inoriented habitats , volume 2071, pages 177–222. Lect. Notes Math., Springer - verlag -, 2013.[48] A. Palcewski.
Velocity averaging for boundary value problems , pages 1–284. Ser. Adv. Math.Appl. Sci. World Scientific Publishing Company, 1992.[49] R. Pettersson. On solutions to the Linear Boltzmann equation for granular gases.
TransportTheor. Stat. , 33(5-7):527–543, 2004.[50] R. G. Plaza. Derivation of a bacterial nutrient-taxis system with doubly degenerate cross-diffusion as the parabolic limit of a velocity-jump process.
J. Math. Biol. , 78(6):1681–1711,2019.[51] K. E. Pourfarhangi, E. Hoz, A. Cohen, and B. Gligorijevic. Contact guidance is cell cycle-dependent.
APL Bioeng. , 2:031904, 2018.[52] P. P. Provenzano, K. W. Eliceiri, J. M. Campbell, and et al. Collagen reorganization at thetumor-stromal interface facilitates local invasion.
BMC Med. , 4(1):38, 2006.[53] P. P. Provenzano, K. W. Eliceiri, and P. J. Keely. Shining new light on 3d cell motility andthe metastatic process.
Trends Cell Biol. , 19(11):638–648, 2009.[54] A. M. Rajnicek, L. E. Foubister, and C. D. McCaig. Prioritising guidance cues: Direc-tional migration induced by substratum contours and electrical gradients is controlled by arho/cdc42 switch.
Dev. Biol. , 312(1):448–460, 2007.[55] S. W. Rhee, A. M. Taylor, C. H. Tu, D. H. Cribbs, C. Cotman, and N. Li Jeon. Patternedcell culture inside microfluidic devices.
Lab Chip , 51:102–107, 2005.[56] D. Schlter, I. Ramis-Conde, and M. Chaplain. Computational modeling of single-cell migra-tion: The leading role of extracellular matrix fibers.
Biophys. J. , 103:1141–51, 2012.[57] M. Scianna, L. Preziosi, and K. Wolf. A cellular potts model simulating cell migration onand in matrix environments.
Math. Biosci. Eng. , 10:235–261, 2013.[58] P. Steeg. Targeting metastasis.
Nat. Rev. Cancer. , 16:201–218, 2016.3359] D. W. Stroock. Some stochastic processes which arise from a model of the motion of abacterium.
Z. Wahrscheinlichkeit , 28(4):305–315, 1974.[60] H. Sundararaghavan, R. Saunders, D. Hammer, and J. Burdick. Fiber alignment directs cellmotility over chemotactic gradients.
Biotechnol. Bioeng. , 110(4):1249–1254, 2013.[61] M. A. Wagle and R. T. Tranquillo. A self-consistent cell flux expression for simultaneouschemotaxis and contact guidance in tissues.
J. Math. Biol. , 41(4):315–330, 2000.[62] P. C. Wilkinson and J. M. Lackie. The influence of contact guidance on chemotaxis of humanneutrophil leukocytes.
Exp. Cell Res. , 145(2):255–264, 1983.[63] K. Wolf, I. Mazo, H. Leung, K. Engelke, U. H. von Andrian, E. I. Deryugina, A. Y. Stron-gin, E.-B. Br¨ocker, and P. Friedl. Compensation mechanism in tumor cell migration: mes-enchymalamoeboid transition after blocking of pericellular proteolysis.