Modelling physical limits of migration by a kinetic model with non-local sensing
MModelling physical limits of migration by a kinetic modelwith non-local sensing
Nadia Loy ∗ Luigi Preziosi †‡ August 23, 2019
Abstract
Migrating cells choose their preferential direction of motion in response to different signalsand stimuli sensed by spanning their external environment. However, the presence of densefibrous regions, lack of proper substrate, and cell overcrowding may hamper cells from mov-ing in certain directions or even from sensing beyond regions that practically act like physicalbarriers. We extend the non-local kinetic model proposed by Loy and Preziosi (2019) to in-clude situations in which the sensing radius is not constant, but depends on position, sensingdirection and time as the behaviour of the cell might be determined on the basis of informa-tion collected before reaching physically limiting configurations. We analyse how the actualpossible sensing of the environment influences the dynamics by recovering the appropriatemacroscopic limits and by integrating numerically the kinetic transport equation.
Keywords : Biased cell migration, Extracellular matrix, Taxis, Physical limit of migra-tion.
During their motion, cells sense the external environment thanks to their protrusions which mayextend up to several cell diameters. The captured chemical or mechanical signals activate trans-duction pathways inside the cell leading to the cell response which consists in ( i ) the formation ofa “head” and a “tail”, ( ii ) in the triggering of actin polimerization at the front edge and depolar-ization at the rear of the cell, and ( iii ) in the activation of adhesion molecules and traction forcesleading eventually to motion (Abercrombie and Heaysman, 1953; Adler, 1966; Block et al., 1983).The above steps can be somewhat distinguished in polarization and mobility mechanisms. Forinstance, Devreotes and Janetopoulos (2003) showed that D. Discoideum cells polarize in responseto cAMP even when treated with inhibitors of the cytoskeleton, such as latrunculin A, that inhibitcell motion.After external stimuli determine the preferential direction of a cell, in addition to internalcauses, other environmental cues may promote or hamper the movement in that direction, such asthe cell density. For instance, the presence of other cells influences cell migration in a two-fold way.On one hand cells may be attracted due to the mutual interaction of transmembrane adhesionmolecules ( e . g . , cadherin complexes). On the other hand, cells may stay away from too crowded ∗ Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca degli Abruzzi24, 10129 Torino, Italy, and Department of Mathematics “G. Peano”, Via Carlo Alberto 10 ,10123 Torino, Italy( [email protected] ) † Department of Mathematical Sciences “G. L. Lagrange”, Dipartimento di Eccellenza 2018-2022, Politecnico diTorino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy ( [email protected] ) ‡ Corresponding author: [email protected] a r X i v : . [ q - b i o . CB ] A ug egions or just lean on their boundaries. Another important migration determinant is the extra-cellular matrix (ECM), i . e . the network of macro-molecules (such as proteoglycans, collagene,fibronectin, and elastin) representing one of the main non-cellular component of all tissues andorgans. In fact, it provides cells with a physical scaffold. Its density, stiffness, and microstructurehighly influence the behaviour of cells and in particular their migration mode (te Boekhorst et al.,2016). Again, the presence of ECM is necessary to form focal adhesion sites through the activationof integrins that are used by the cell to exert traction forces, but if it is too dense it might representa steric obstacle to cell motion.In particular, experiments (see, for instance, Schoumacher et al. (2010); Shankar et al. (2010))show that while the cytoplasm is very flexible and able to accommodate nearly any pore size(including 1 µm gaps in collagen gels and 0 . µm pores of polycarbonate membranes), the cellnucleus is five- to ten-fold stiffer than the surrounding cytoplasm and, with a typical diameter of3 − µm , might be larger than the ECM fiber spacing (Davidson et al., 2014). During MMP-independent migration ( i . e . , when the proteolytic machinery is inhibited) and in spite of cellcytoplasm protrusions into the ECM trying to pull the nucleus inside, the stiff nucleus may beunable to squeeze through narrow pores, setting a critical pore size below which MMP-inhibitedcells remain trapped (Wolf et al., 2013). This phenomenon is named physical limit of migrationand it has been recently studied from the modelling point of view by Scianna et al. (2013); Sciannaand Preziosi (2013, 2014) using cellular Potts models and by Arduino and Preziosi (2015); Giversoet al. (2014, 2017) using continuum mechanics methods.The aim of this article is to include such effects extending the kinetic model developed byLoy and Preziosi (2019), where possibly independent cues are sensed by cells in a non-local wayand used to determine respectively cell polarization and speed. The transport kinetic equationimplements then a velocity jump process that describes the movement performed by the cell as analternation of runs over straight lines and re-orientations (also called tumbles or turnings) (Stroock,1974), also considering the bias induced by external stimuli, as done by Alt (1980) and Othmeret al. (1988). Such an equation describes the evolution of a single-particle density distribution likein the Boltzmann equation (Cercignani, 1987). The main elements of the velocity-jump processare the tumbling frequency, the mean speed, and the transition (also called turning or tumbling)probability that describes the probability of choosing a certain velocity after re-orientation. Themean speed, mean runtime (which is the inverse of the frequency), and the tumbling probabilitymay be measured from individual patterns of members of the population.In the literature many models considered different sensing strategies and their relation to thedetermination of cell re-orientation and speed. Focusing on non-local aspects, tipically in positionjump processes the transition probabilities depend on the acquisition of information sensed at acertain location ( e . g . at the target or at the present location) or by averaging the signal over acertain neighborhood 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 externalfields. Othmer and Hillen (2002) and Hillen et al. (2007) introduced a finite sensing radius anddefined a non-local gradient as the average of the external field on a surface which representsthe membrane of the cell. This idea was also used for cell adhesion and haptotaxis by Armstronget al. (2006) yielding a macroscopic integro-differential equation. Buttensch¨on et al. (2018) recentlyderived this type of models from a space jump process. Other macroscopic models describing cellmigration with non-local measures of the environment were proposed by Painter et al. (2010,2015); Painter and Hillen (2002). Schmeiser and Nouri (2017) considered a kinetic model withvelocity jumps biased towards the chemical concentration gradient. Similar equations were alsoproposed in 2D set-ups by Colombi et al. (2015, 2017), and applied to model crowd dynamics andtraffic flow for instance by Tosin and Frasca (2011). Eftimie (2012) and Carrillo et al. (2015) alsoproposed a non-local kinetic models including repulsion, alignement and attraction and used it formodelling tumor dynamics (Eftimie et al., 2017) and cell polarisation in heterogeneous cancer cellpopulations (Bitsouni and Eftimie, 2018).In Loy and Preziosi (2019) the transition probability of the transport equation models theprobability of choosing a certain velocity direction and speed according to an environmental sens-ing over a finite neighbourhood of the cell, giving the kinetic model a non-local character. In2articular, a double bias is considered, as cells perform a double sensing both of a tactic cue in-fluencing the probability of polarizing in a certain direction and of an environmental cue, usuallyof mechanical origins, affecting its speed in the polarization direction. For instance, in presenceof a chemoattractant giving a preferential direction of motion, a volume filling effect can hamperthe motion in a certain direction due to cell overcrowding.In the present work, then, the focus is on how to include in kinetic models the effect of physicallimits of migration. These include the dependence of the sensing radius on physical characteristicsof the environment or on the omotypic or heterotypic cell distributions, that might even hamperthe real possibility of cells of measuring in a certain region. In particular, this means that theinterval with valuable information to determine cell motion, depends on space, time and directionof sensing. We also discuss how limiting the sampling volume determines the type of macroscopiclimit that can be performed.With this aim in mind, the plan of the article is then the following. In Section 2 the kineticmodelling framework is introduced and then specialised to the case of physical limits of migrationin Section 3. In particular, the concept of limitation of sensing radius and its effect on motility isqualitatively described. Section 4 focuses on deducing the appropriated macroscopic limit of thekinetic model showing that in presence of physical limits of migration the macroscopic speed doesnot vanish in general and therefore the appropriate limit is hyperbolic. Section 5 starts analysingthe model from a simpler case in which there is no cue determining a preferred polarization inspace, e . g . , no chemotaxis, but only mechanical cues influencing cell speed. Simulations of thekinetic model are performed focusing on volume filling effects in Section 5.1 and on cell-ECMinteractions in Section 5.2. A particular attention is paid to pointing out what is the effect ofdifferent modelling assumptions and how some choices might lead to non-physical results. Then,in Section 6 chemotaxis and cell-cell adhesion are added as mechanisms influencing cell polarizationin order to consider the simultaneous presence of physical limits of migration and cell polarization.A final section draws some conclusions pointing out to other open issues. Let us describe the cell population through the distribution density p = p ( t, x , v, ˆ v ) parametrizedby the time t >
0, the position x ∈ Ω ⊆ R d , the speed v ∈ [0 , U ], where U is the maximal speed acell can achieve, and the polarization direction ˆ v ∈ S d − where S d − is the unit sphere boundaryin R d . The choice of representing the distribution function p depending on velocity modulus anddirection, instead of the velocity vector v = v ˆ v , lies in the need of separating the mechanismsgoverning cell polarization and motility, for instance in response of chemotaxis and in presence ofother factors, typically of mechanical origins, influencing cell speed.The mesoscopic model consists in the transport equation for the cell distribution ∂p∂t ( t, x , v, ˆ v ) + v · ∇ p ( t, x , v, ˆ v ) = J [ p ]( t, x , v, ˆ v ) (1)where the operator ∇ denotes the spatial gradient. The term J [ p ]( t, x , v, ˆ v ), named turningoperator , is an integral operator that describes the change in velocity which is not due to free-particle transport. It may describe the classical run and tumble behaviors, contact guidancephenomena, or cell-cell interactions. For the moment we will consider the classical run and tumble, e . g . , random re-orientations, which, however, may be biased by external cues. Therefore, ourturning operator will be the implementation of a velocity-jump process in a kinetic transportequation as introduced by Stroock (1974) and then by Othmer et al. (1988).Defining V p = [0 , U ] × S d − , a macroscopic description for the cell population can be classicallyrecovered through the definition of moments of the distribution function p as follows- the cell number density ρ ( t, x ) = (cid:90) V p p ( t, x , v, ˆ v ) dv d ˆ v ; (2)3 the cell mean velocity U ( t, x ) = 1 ρ ( t, x ) (cid:90) V p p ( t, x , v, ˆ v ) v dv d ˆ v ; (3)- the cell variance-covariance matrix P ( t, x ) = (cid:90) V p p ( t, x , v, ˆ v )[ v − U ( t, x )] ⊗ [ v − U ( t, x )] dv d ˆ v ; (4)- the cell speed variance E ( t, x ) = (cid:90) V p p ( t, x , v, ˆ v ) | v − U ( t, x ) | dv d ˆ v . (5)We remark that, because of the definition of V p , the integrals over the velocity space area simple double integral. In particular, if the dependence of f on v and ˆ v can be factorized, i . e . f ( t, x , v, ˆ v ) = f ( t, x , v ) f ( t, x , ˆ v ), we have that (cid:90) V p f ( t, x , v, ˆ v ) dv d ˆ v = (cid:90) U f ( t, x , v ) dv (cid:90) S d − f ( t, x , ˆ v ) d ˆ v , (6)where the integral over the boundary of the unit sphere S d − is not a surface integral, but has tobe interpreted as (cid:90) S d − f ( t, x , ˆ v ) d ˆ v = (cid:90) π f ( t, x , cos θ, sin θ ) dθ in 2D and similarly in 3D.The general form of the turning operator which implements a velocity jump processes is J [ p ]( t, x , v, ˆ v ) = (cid:90) V p µ ( x , (cid:48) v, (cid:48) ˆ v ) T ( x , v, ˆ v | (cid:48) v, (cid:48) ˆ v ) p ( t, x , (cid:48) v, (cid:48) ˆ v ) d (cid:48) v d (cid:48) ˆ v − (cid:90) V p µ ( x , v, ˆ v ) T ( x , v (cid:48) , ˆ v (cid:48) | v, ˆ v ) p ( t, x , v, ˆ v ) dv (cid:48) d ˆ v (cid:48) , (7)where (cid:48) v = (cid:48) v (cid:48) ˆ v is the pre-turning velocity of the gain term and v (cid:48) = v (cid:48) ˆ v (cid:48) is the post-turningvelocity of the loss term. The so-called turning kernel T ( x , v, ˆ v | (cid:48) v, (cid:48) ˆ v ) is the probability for a cellin x of re-orienting along ˆ v and moving with speed v given the pre-turning polarization direction (cid:48) ˆ v and speed (cid:48) v . Being a transition probability, it satisfies (cid:90) V p T ( x , v (cid:48) , ˆ v (cid:48) | v, ˆ v ) dv (cid:48) d ˆ v (cid:48) = 1 , ∀ x ∈ Ω , ∀ ( v, ˆ v ) ∈ V p , (8)which allows to simplify (7) to J [ p ]( t, x , v, ˆ v ) = (cid:90) V p µ ( x , (cid:48) v, (cid:48) ˆ v ) T ( x , v, ˆ v | (cid:48) v, (cid:48) ˆ v ) p ( t, x , (cid:48) v, (cid:48) ˆ v ) d (cid:48) v d (cid:48) ˆ v − µ ( x , v, ˆ v ) p ( t, x , v, ˆ v ) . (9)As done by Stroock (1974); Hillen (2006); Chauviere et al. (2007a,b), in the following we willassume that cells retain no memory of their velocity prior to the re-orientation, i . e . , T = T ( x , v, ˆ v ).The independence from the pre-tumbling velocity lies in the fact that the choice of the new velocityis linked to the slow interaction process also related to cell ruffling and sensing which is responsiblefor the biased re-orientation. However, the assumption might be restrictive in some cases, as itexcludes, for instance, persistence effects in which the re-orientation direction depends on the pre-tumbling polarization of the cell and the case in which the sensing region depends on the incomingvelocity through a polarization-dependent expression of transmembrane receptors.4ssuming also that the frequency µ does not depend on the microscopic velocity allows tosimplify considerably (9) in J [ p ]( t, x , v, ˆ v ) = µ ( x ) (cid:16) ρ ( t, x ) T ( x , v, ˆ v ) − p ( t, x , v, ˆ v ) (cid:17) . (10)Following Loy and Preziosi (2019), we consider here the fact that in taxis processes cells arecapable of detecting and measuring external signals through membrane receptors located alongcell protrusions that can extend over a finite radius. Information are then transduced and act ascontrol factors for the dynamics of cells. Therefore, in the turning operator of the kinetic modelwe will include the evaluation of mean fields in a neighbourhood of the re-orientation position.In particular, we will consider a control factor S determining cell polarization (and thereforeorientation) and a control factor S (cid:48) determining cell speed in that specific orientation direction.Hence, we will have a transition probability that depends on both S and S (cid:48) that can be writtenas T [ S , S (cid:48) ]( x , v, ˆ v ) = 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) , (11)where c ( x ) is a normalization constant, so that according to (8) the integral of T over the velocityspace is one. In this way, T is a mass preserving transition probability that takes into accounttwo different external fields S and S (cid:48) in order to choose the post-tumbling velocity. Specifically,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 weighted by the sensing kernel γ S ( λ ) where λ measures the distance from x . In particular, γ S hasa compact support in [0 , R max S ] where R max S is the maximum extension of cell protrusions, i . e . thefurthest point cells can reach to measure the external signals. In all the applications consideredhere we take T ˆ v λ [ S ]( x ) = b (cid:0) S ( x + λ ˆ v ) (cid:1) , (12)but other functionals could be considered as discussed by Loy and Preziosi (2019).The density distribution of the modulus ψ = ψ ( x , v |S (cid:48) ( x + λ (cid:48) ˆ v )) describes the choice of thespeed v ∈ [0 , U ] for cells located in x if oriented along ˆ v according to the value of another externalfield S (cid:48) measured along this direction and weighted by another function γ S(cid:48) with compact supportin [0 , R max
S(cid:48) ]. In fact, the notation |S (cid:48) reads as ”given S (cid:48) in a point x + λ (cid:48) ˆ v in the neighbourhood ofthe cell”. In the following the mean of ψ , that is normalized to 1 i . e ., (cid:82) U ψdv = 1 will be denotedby ¯ v ( x |S (cid:48) ( x + λ (cid:48) ˆ v )) and its variance by D ( x |S (cid:48) ( x + λ (cid:48) ˆ v )).Figure 1: Cell extending protrusions toward the top left corner but unable to move becausethe nucleus (yellow region on the bottom right corner) can not pass through the pores of theextracellular matrix (by courtesy of P. Friedl and K. Wolf, under permission).5 Modelling Physical Limits of Migration
The main novelty of this paper with respect to Loy and Preziosi (2019) is to consider the existenceof physical limits of migration that hamper the cell from sensing or moving beyond a physicalbarrier. To be specific, for instance, we want to deal with volume filling effects, that occur whencells ahead are too packed for the coming cell to overcome the crowd, or the presence of regionswith high density of extracellular matrix (ECM), or better, pores so small that the cell nucleuscan not pass through them (see, for instance, the work by Wolf et al. (2013)). This means that R max S and R max S(cid:48) can not be reached, or even if they are, the decision is taken on the basis ofinformation collected in the space preceding the overcrowded region. Hence, the sensing supportof the weight functions γ S and γ S(cid:48) can be reduced. We mean, for instance, that a cell that isencountering on its way a physical barrier, e . g . , a basal membrane, might be unable to squeezethrough the fine net, but it can still extend cytoplasmic protrusions beyond the dense area (seeFigure 2). Even if the cell is getting information, generally denoted by M , that there would bespace to move beyond the barrier, then this does not influence the possibility to pass through it,because the barrier constitutes an impediment to go further. However, we observe that in orderto test the physical limit of migration we allow the cell to poke a bit in it up to a small depth ∆.For instance, the nucleus in contact with a very small pore in the ECM, still tries to squeeze a bitin, or at least leaning on the entrance of the pore, part of the nucleus will still be in the pore. Thisalso corresponds to the limited, but not absent, freedom cells have also in constrained situations,like in the case presented in Video 1 of Wolf et al. (2013).In principle, the limitation on the sensing radius can apply to both the polarization signal S and the speed signal S (cid:48) , though it is easier to find examples of the latter case. In mathematicalterms we then assume that there might be a cue M ( t, x ) related to the sensing of S (cid:48) characterizedby a threshold value M th representing the physical limit of migration. We will sometimes callit a mechanical signal, though it might be more general than that, like, for instance, the lackof adhesion sites in a certain region to allow cell traction and then migration in a certain area.Other examples regard volume filling when S (cid:48) and M are equal to the density of cells ρ and ECMhindrance effects where they are related to the matrix density M , or, better, to its characteristicpore cross section.We then define R M ( t, x , ˆ v ) = R max S(cid:48) if M ( t, x + λ ˆ v ) ≤ M th , ∀ λ ∈ [0 , R max S(cid:48) ) , inf { λ ∈ [0 , R max S(cid:48) ) : M ( t, x + λ ˆ v ) > M th } , otherwise , (13)and R MS(cid:48) ( t, x , ˆ v ) = min { R M ( t, x , ˆ v ) + ∆ , R max S(cid:48) } . (14)Figure 2 highlights the values of R M (gray circles) and R MS(cid:48) (black squares) related to severalpossible landscapes of M . For instance, the landscape corresponding to the top full line, the cellis already in a region with M > M th and so it can hardly move, R M = 0 and R MS(cid:48) = ∆. The sameoccurs for the second full line from the top. In this case the cell is very close to the border of thebarrier, but still inside it. Conversely, the third full line from the top corresponds to a cell at theborder of the physical barrier, but that can look freely ahead. In this case, R M = R MS(cid:48) = R max S .Similarly for the bottom full line. However, looking instead at the top dotted line (line 1), thecell encounters a physical barrier again at λ = R , so that R M = R and R MS(cid:48) = R + ∆. For theother two dotted lines (numbered as 2 and 3) the cell is barely touching or not encountering yetthe barrier, so that R MS(cid:48) = R max S(cid:48) .The dependence of M on time is due to the fact that the mechanical cue may change in time,for instance because of ECM degradation due to metalloproteinases or redistribution of cell mass inthe volume filling case. However, in the following, to simplify a little the notation the dependenceon time is usually dropped.In particular, we remark that 6igure 2: Example of sensing ranges limited by a physical limit of migration corresponding to alevel M th of M . Gray circles denote R M and black squares R MS(cid:48) . • when the cell center is in a point where M < M th , it means that it has not reached thephysical barrier yet and R M > • when the cell center is in a point where M > M th , then it means that it is stuck in anovercrowded region and R M = 0; • when the cell center is in a point where M = M th , unless for very peculiar cases, it meansthat it is at the border of the barrier. In this case R M vanishes if the cell is polarized towardsthe barrier and strictly positive if it wants to move away from it.We also observe that when the presence of a physical barrier limiting R max S(cid:48) to R MS(cid:48) ( t, x , ˆ v )might also limit R max S to R MS ( t, x , ˆ v ).Summarizing, in order to polarize along the most likely direction ˆ v , the cell with sensing radius R max S averages the signal S over a region that is described by a weight function γ S up to R max S ,or possibly up to a lower value R MS ( t, x , ˆ v ) if an obstacle is encountered along its sensing activity.Then, the cell polarized along ˆ v determines its speed averaging another signal S (cid:48) through a weightfunction γ S(cid:48) over the sensing radius R max S(cid:48) or, as above, up to R M(cid:48)S(cid:48) ( t, x , ˆ v ).Coming back to Eq.(11), as in the work by Loy and Preziosi (2019), it is useful to denote thetwo factors in it as B [ S ]( x , ˆ v ) = (cid:90) R MS ( x , ˆ v )0 γ S ( λ ) T ˆ v λ [ S ]( x ) dλ , (15)and Ψ[ S (cid:48) ]( x , v | ˆ v ) = (cid:90) R M(cid:48)S(cid:48) ( x , ˆ v )0 γ S(cid:48) ( λ (cid:48) ) ψ ( x , v |S (cid:48) ( x + λ (cid:48) ˆ v )) dλ (cid:48) , (16)in order to highlight the directional and the speed sensing and how they affect the choice of thedirection and of the speed independently. The turning probability in (10) then reads T [ S , S (cid:48) ]( x , v, v ) = c ( x ) B [ S ]( x , ˆ v )Ψ[ S (cid:48) ]( x , v | ˆ v ) , (17)where the normalizing constant is c ( x ) = 1 (cid:90) S d − B [ S ]( x , ˆ v )Γ S(cid:48) ( x , ˆ v ) d ˆ v , (18)with Γ S(cid:48) ( x , ˆ v ) = (cid:90) R M(cid:48)S(cid:48) ( x , ˆ v )0 γ S(cid:48) ( λ (cid:48) ) dλ (cid:48) . (19)7f R M(cid:48)S(cid:48) is independent of ˆ v , then also Γ S(cid:48) is independent of ˆ v and c ( x ) can be factorized as c ( x ) = c ( x ) / Γ S(cid:48) ( x ) with c ( x ) = 1 (cid:90) S d − (cid:90) R MS ( x , ˆ v )0 γ S ( λ ) T ˆ v λ [ S ]( x ) dλ d ˆ v , that only depends on the directional sensing.However, in general, with the definitions (15) and (16), the turning operator (10) writes J [ p ]( x , v, v ) = µ ( x ) (cid:16) ρ ( t, x ) c ( x ) B [ S ]( x , ˆ v )Ψ[ S (cid:48) ]( x , v | ˆ v ) − p ( t, x , v, v ) (cid:17) , (20)and the distribution function which nullifies it is p ( t, x , v, v ) = ρ ( t, x ) c ( x ) B [ S ]( x , ˆ v )Ψ[ S (cid:48) ]( x , v | ˆ v ) . (21)As stated by Loy and Preziosi (2019), provided that the probability distribution has initially afinite mass and energy and non-absorbing boundary conditions hold, the function (21) is a localstable asymptotic equilibrium state (Bisi et al., 2008; Pettersson, 2004). The related mean velocityis defined by U S , S(cid:48) ( x ) = c ( x ) (cid:90) V p B [ S ]( x , ˆ v )Ψ[ S (cid:48) ]( x , v | ˆ v ) v dv d ˆ v , (22)and may be rewritten as U S , S(cid:48) ( x ) = c ( x ) (cid:90) S d − Γ S(cid:48) ( x , ˆ v ) B [ S ]( x , ˆ v ) ¯ U S(cid:48) ( x | ˆ v )ˆ v d ˆ v , (23)where ¯ U S(cid:48) ( x | ˆ v ) = 1Γ S(cid:48) ( x , ˆ v ) (cid:90) R M(cid:48)S(cid:48) ( x , ˆ v )0 γ S(cid:48) ( λ (cid:48) )¯ v ( x |S (cid:48) ( x + λ (cid:48) ˆ v )) dλ (cid:48) . (24)Similarly, the variance-covariance tensor of the transition probability D S , S(cid:48) ( x ) = c ( x ) (cid:90) V p B [ S ]( x , ˆ v ) Ψ[ S (cid:48) ]( x , v | ˆ v )( v − U S , S(cid:48) ( x )) ⊗ ( v − U S , S(cid:48) ( x )) dv d ˆ v , (25)can be written as D S , S(cid:48) ( x ) = c ( x ) (cid:90) S d − Γ S(cid:48) ( x , ˆ v ) B [ S ]( x , ˆ v ) ¯ D S(cid:48) ( x | ˆ v )ˆ v ⊗ ˆ v d ˆ v − U S , S(cid:48) ( x ) ⊗ U S , S(cid:48) ( x ) , (26)where ¯ D S(cid:48) ( x | ˆ v ) = 1Γ S(cid:48) ( x , ˆ v ) (cid:90) R M(cid:48)S(cid:48) ( x , ˆ v )0 γ S(cid:48) ( λ (cid:48) ) D ( x |S (cid:48) ( x + λ (cid:48) ˆ v )) dλ (cid:48) . (27)In order to understand the meaning of ¯ U S(cid:48) in the presence of physical limits of migration, letus consider the following case: assume that ahead of a cell there is a space free of obstacles up toa distance R , e . g . , the ECM density has a density yielding a mean speed ¯ v , while after R theECM density is so large to represent a physical limit of migration with a much lower mean speed,say ¯ ε or even vanishing with ¯ v = 0 and ψ = δ ( v ). For sake of simplicity, let us assume that γ S(cid:48) isa Heaviside function. The following interesting cases are possible • If R > R max S(cid:48) , then the critical value for M is located beyond the sensing radius (as for thedashed line 1 in Fig. 2), so that R MS(cid:48) = R M = R max S(cid:48) , and finally ¯ U S(cid:48) = ¯ v ; • If R max S(cid:48) − ∆ < R < R max S(cid:48) (as for the dashed line 2 in Fig. 2), then the critical value for M is barely sensed with the cell poking a bit in the region with M > M th . In this case R M = R but R MS(cid:48) = R max S(cid:48) , so that ¯ U S(cid:48) = ¯ v R + ¯ ε ( R max S(cid:48) − R ) R max S(cid:48) ;8 If R < R max S(cid:48) − ∆ (as for the dashed line 3 in Fig. 2), then the cell pokes for a depth∆ into the region with M > M th . In this case R M = R and R MS(cid:48) = R + ∆, so that¯ U S(cid:48) = ¯ v R + ¯ ε ∆ R + ∆ ; • In particular, if R = 0, the cell has arrived at the barrier, then R M = 0 and R MS(cid:48) = ∆, sothat ¯ U S(cid:48) = ¯ ε . Hence, in the limit ∆ → ε , which can also be taken to be zero, as in the simulations to follow; • Conversely, if the cell is so close to an exit, say at a distance R exit ∈ [0 , ∆], that it can senseoutside the physical barrier, then R = 0 and R MS(cid:48) = ∆, so that ¯ U S(cid:48) = ¯ v − (¯ v − ¯ ε ) R exit ∆Hence, in the limit ∆ → v .Figure 3: Speed of a cell polarized to move to the right according to its position in presenceof a speed limit region in the gray zone ( a, b ), for instance due to too dense ECM to allow cellmigration.Figure 3 gives an example of the measured ¯ U S(cid:48) as a function of the cell position (with the cellpolarized to move to the right) when the region with M > M th is the interval ( a, b ). The cellmoves at a velocity ¯ v till it reaches a distance R max S(cid:48) from the barrier. Then it decreases its speedat first only sligthly and then faster when getting closer to the barrier to reach a very small orvanishing value. If on the other hand the cell is at the right border of the barrier it can movefreely with velocity ¯ v .In the simulations to follow, in order to test the physical limit of migration, we shall alwaysconsider ∆ = ¯ ε = 0, i . e . , the case in which the cell can not poke in the limited area. We leave toa future work the investigation of the effect of considering a non-vanishing ∆. One of the main issues about transport equations consists in recovering the appropriate macro-scopic limit allowing to highlight the driving macroscopic phenomenon. Typically these are ob-tained by identifying a small parameter (cid:15) (cid:28) τ = (cid:15) t, ξ = (cid:15) x ,τ = (cid:15)t, ξ = (cid:15) x . The proper choice of the time scale comes from a nondimensionalization of the turning operator(10). Referring to the functions introduced in (15) and (16), we suppose that, up to a nondimen-sionalization, B [ S ] and Ψ[ S (cid:48) ] may be written as B [ S ] = B [ S ] + (cid:15)B [ S ] + O ( (cid:15) ) , Ψ[ S (cid:48) ] = Ψ[ S (cid:48) ] + (cid:15) Ψ[ S (cid:48) ] + O ( (cid:15) ) , which describe different orders of bias. Then, the turning kernel T [ S , S (cid:48) ] can be written as T [ S , S (cid:48) ]( ξ , v, ˆ v ) = T [ S , S (cid:48) ] ( ξ , v, ˆ v ) + (cid:15)T [ S , S (cid:48) ] ( ξ , v, ˆ v ) + O ( (cid:15) ) , (28)where T [ S , S (cid:48) ] ( v, ˆ v ) = c ( x ) B [ S ] (ˆ v )Ψ[ S (cid:48) ] ( v | ˆ v ) , (29)and T [ S , S (cid:48) ] ( v, ˆ v ) = c ( x ) B [ S ] (ˆ v )Ψ[ S (cid:48) ] ( v | ˆ v ) + c ( x ) B [ S ] (ˆ v )Ψ[ S (cid:48) ] ( v | ˆ v ) . (30)Coherently, 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) . Similarly, the distribution function p is expanded as p = p + (cid:15)p + O ( (cid:15) ) . (31)As there is conservation of mass we have that (Othmer and Hillen, 2000) all the mass is in p , i . e . , ρ = ρ, ρ i = 0 ∀ i ≥ , (32)where ρ i = (cid:90) V p p i dv d ˆ v . Furthermore, for performing the diffusive limit we shall suppose that (cid:90) V p p i v ˆ v dv d ˆ v = 0 ∀ i ≥ (cid:90) U Ψ[ S (cid:48) ] ( v | ˆ v ) dv = (cid:90) U Ψ[ S (cid:48) ]( v | ˆ v ) dv , (cid:90) U Ψ[ S (cid:48) ] i ( v | ˆ v ) dv = 0 ∀ i ≥ , and (cid:90) S d − B [ S ] (ˆ v ) d ˆ v = (cid:90) S d − B [ S ](ˆ v ) d ˆ v , (cid:90) S d − B [ S ] i (ˆ v ) d ˆ v = 0 ∀ i ≥ , then (cid:90) V p T [ S , S (cid:48) ] ( ξ , v, ˆ v ) dv d ˆ v = 1 , (33a)and (cid:90) V p T [ S , S (cid:48) ] i ( ξ , v, ˆ v ) dv d ˆ v = 0 ∀ i ≥ , (33b)that are needed for performing a macroscopic limit (Othmer and Hillen, 2000).Denoting by J i the turning operators defined by T [ S , S (cid:48) ] i , it turns out for instance that J [ p ]( v, ˆ v ) = µ ( ρ T [ S , S (cid:48) ] ( v, ˆ v ) − p ( v, ˆ v )) = 0 , (34)and therefore the function p is the equilibrium function given by p ( v p ) = ρ T [ S , S (cid:48) ] ( v, ˆ v ) . (35)10he velocity and tensor of the second moments of the transition probability then naturallywrite as U S , S(cid:48) = U S , S(cid:48) + (cid:15) U S , S(cid:48) + O ( (cid:15) ) , (36)and D S , S(cid:48) = D S , S(cid:48) + (cid:15) D S , S(cid:48) + O ( (cid:15) ) , (37)where U i S , S(cid:48) = (cid:90) V p T [ S , S (cid:48) ] i v dv d ˆ v , 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) ) dv d ˆ v , for i ≥
0. Then, using (29) and (30), the macroscopic velocity of order 0 is U S , S(cid:48) ( ξ ) = c ( ξ ) (cid:90) S d − Γ S(cid:48) ( ξ , ˆ v ) B [ S ] (ˆ v ) ¯ U S(cid:48) ( ξ | ˆ v )ˆ v d ˆ v , (38)that of order 1 is U S , S(cid:48) ( ξ ) = c ( ξ ) (cid:90) S d − Γ S(cid:48) ( ξ , ˆ v ) (cid:16) B [ S ] (ˆ v ) ¯ U S(cid:48) ( ξ | ˆ v ) + B [ S ] (ˆ v ) ¯ U S(cid:48) ( ξ | ˆ v ) (cid:17) ˆ v d ˆ v (39)and the equilibrium diffusion tensor is D S , S(cid:48) ( x ) = c ( ξ ) (cid:90) S d − Γ S(cid:48) ( ξ , ˆ v ) B [ S ] ( ξ , ˆ v ) ¯ D S(cid:48) ( ξ | ˆ v )ˆ v ⊗ ˆ v d ˆ v − U S , S(cid:48) ( ξ ) ⊗ U S , S(cid:48) ( ξ ) , (40)where ¯ D S(cid:48) ( ξ | ˆ v ) = 1Γ S(cid:48) ( ξ , ˆ v ) (cid:90) R M(cid:48)S(cid:48) ( ξ , ˆ v )0 γ S(cid:48) ( λ (cid:48) ) D ( ξ |S (cid:48) ( ξ + λ (cid:48) ˆ v )) dλ (cid:48) . The diffusion tensor is in general anisotropic, as ¯ D S(cid:48) may be non isotropic in ˆ v .The discriminating factor on which scaling can be performed is whether T is such that U S , S(cid:48) = (41) i . e . , the leading order macroscopic velocity vanishes, or not. This is the necessary functionalcondition for performing a diffusive limit, and it is actually what we expect in a diffusive regime.Referring to Eq. (38), we observe that, even in the case in which B [ S ] is constant, correspondingfor instance to no directional bias and R M(cid:48)S(cid:48) , and therefore Γ
S(cid:48) not depending on ˆ v , Eq. (41) issatisfied if the mean velocity ¯ U S(cid:48) ( ξ | ˆ v ) is even as a function of ˆ v for a.e. ξ ∈ Ω. Moreover, if R M(cid:48)S(cid:48) depends on ˆ v , condition (41) can not be satisfied, unless for trivial cases that do not includephysical limits of migration, as we will illustrate in the next sections through some examples.Therefore, in general one can only perform a hyperbolic limit which reads ∂ρ∂τ + ∇ · (cid:0) ρ U S , S(cid:48) (cid:1) = 0 . (42)In order to see which regimes are diffusive or hyperbolic, it is instructive to see what happensif R S(cid:48) ( x , ˆ v ) is much smaller then the characteristic length of variation of S (cid:48) as done by Loy andPreziosi (2019) for the directional bias. In this case, we can expand S (cid:48) as S (cid:48) ( x + λ (cid:48) ˆ v ) = S (cid:48) ( x ) + λ (cid:48) ˆ v · ∇S (cid:48) ( x ) + O ( λ (cid:48) ) ∀ λ (cid:48) ≤ R S(cid:48) ( x , ˆ v ) (43)and we have that that the quantity S (cid:48) ( x ) + λ (cid:48) ˆ v · ∇S (cid:48) ( x ) stays positive. The density function ψ ( v |S (cid:48) ( y )), y ∈ Ω may be seen as a function of two variables˜ ψ : ( v, y ) ∈ [0 , U ] × Ω (cid:55)−→ ˜ ψ ( v, y ) = ψ ( v |S (cid:48) ( y ))11omposed with the function S (cid:48) defined on Ω. If we suppose that ˜ ψ ∈ L ([0 , U ]) × C (Ω) with (cid:90) U ˜ ψ ( v, y ) dv = (cid:90) U ψ ( v |S (cid:48) ( y )) dv = 1 , (44)then, in virtue of (43) we may write˜ ψ ( v, S (cid:48) ( x + λ (cid:48) ˆ v )) = ˜ ψ ( v, S (cid:48) ( x ) + λ (cid:48) ˆ v · ∇S (cid:48) ( x ) + O ( λ (cid:48) ))= ˜ ψ ( v, S (cid:48) ( x )) + ∂∂ S (cid:48) [ ˜ ψ ( v, S (cid:48) ( x ))] λ (cid:48) ˆ v · ∇S (cid:48) ( x ) + O ( λ (cid:48) ) ∀ λ (cid:48) ≤ R S(cid:48) ( x , ˆ v ) . Therefore, recalling Eqs.(17) and (18) the transition probability writes T ( x , v, ˆ v ) = B [ S ](ˆ v ) (cid:90) S d − B [ S ](ˆ v )Γ S(cid:48) ( x , ˆ v ) d ˆ v (cid:90) R S(cid:48) ( x , ˆ v )0 γ S(cid:48) ( λ (cid:48) ) ψ ( v |S (cid:48) ( x + λ (cid:48) ˆ v )) dλ (cid:48) ≈ B [ S ](ˆ v )Γ S(cid:48) ( x , ˆ v ) (cid:90) S d − B [ S ](ˆ v )Γ S(cid:48) ( x , ˆ v ) d ˆ v (cid:20) ψ ( v |S (cid:48) ( x )) + Λ (cid:48) ( x , ˆ v ) ∂∂ S (cid:48) [ ψ ( v |S (cid:48) ( x ))] ∇S (cid:48) ( x ) · ˆ v (cid:21) where Λ (cid:48) ( x , ˆ v ) = 1Γ S(cid:48) ( x , ˆ v ) (cid:90) R S(cid:48) ( x , ˆ v )0 γ S(cid:48) ( λ (cid:48) ) λ (cid:48) dλ (cid:48) . The quantity ∂∂ S (cid:48) ¯ v ( x |S (cid:48) ( x ))is the derivative of the master curve relating ¯ v to the signal S (cid:48) .Considering the rescaling ξ = (cid:15) x , one has that T [ S , S (cid:48) ] = B [ S ](ˆ v )Γ S(cid:48) ( ξ , ˆ v ) (cid:90) S d − B [ S ](ˆ v )Γ S(cid:48) ( ξ , ˆ v ) d ˆ v ψ ( v |S (cid:48) ( ξ ))and T [ S , S (cid:48) ] = B [ S ](ˆ v )Γ S(cid:48) ( ξ , ˆ v )Λ (cid:48) ( ξ , ˆ v ) (cid:90) S d − B [ S ](ˆ v )Γ S(cid:48) ( ξ , ˆ v ) d ˆ v ∂∂ S (cid:48) [ ˜ ψ ( v, S (cid:48) ( ξ ))] ∇S (cid:48) ( ξ ) · ˆ v . Integrating over the velocity space V p it can be readily checked that T satisfies condition (33a)and T satisfies condition (33b) observing that (cid:90) U ∂∂ S (cid:48) ψ ( v |S (cid:48) ( ξ )) dv = ∂∂ S (cid:48) (cid:90) U ψ ( v |S (cid:48) ( ξ )) dv = 0 . One then has that U S(cid:48) = ¯ v ( ξ |S (cid:48) ( ξ )) N S(cid:48) , where N S(cid:48) = (cid:90) S d − B [ S ](ˆ v )Γ S(cid:48) (ˆ v )ˆ v d ˆ v (cid:90) S d − B [ S ](ˆ v )Γ S(cid:48) (ˆ v ) d ˆ v , (45)and U S(cid:48) = ∂∂ S (cid:48) [¯ v ( ξ |S (cid:48) ( ξ ))] T S(cid:48) ∇S (cid:48) , where T S(cid:48) = (cid:90) S d − B [ S ](ˆ v )Γ S(cid:48) (ˆ v )ˆ v ⊗ ˆ v d ˆ v (cid:90) S d − B [ S ](ˆ v )Γ S(cid:48) (ˆ v ) d ˆ v .
12e then observe that the mean direction N S(cid:48) defined in Eq. (45) of the cell population may hardlybe zero, except in the case in which B and Γ S(cid:48) are even in every direction, i . e .B [ S ]( ξ , ˆ v ) = B [ S ]( ξ , − ˆ v ) and Γ S(cid:48) ( ξ , ˆ v ) = Γ S(cid:48) ( ξ , − ˆ v ) a.e. in Ω . (46)This may happen only if R S(cid:48) is even as a function of the direction ˆ v , e . g . , if it is constant whichmeans that there is no physical limit of migration. As soon as a physical barrier appears, in pointsclose to it an asymmetry appears in the evaluation of the sensing radius which make (46) notvalid.In the simulations to follow, we shall illustrate that because of these reasons, even if R max S issmall as well or in the case of isotropic polarization B = const , the mean velocity hardly vanishesclose to physical barriers and, then, Eq. (41) is not satisfied, leading to the necessity of performinga hyperbolic scaling with a parabolic scaling only possible away from the barriers.For sake of completeness, we recall that in these cases the diffusive limit would lead to ∂ρ∂τ + ∇ · (cid:16) ρ U S , S(cid:48) (cid:17) = ∇ · (cid:18) µ ∇ · (cid:0) D S , S(cid:48) ρ (cid:1)(cid:19) . (47) We shall consider the case in which there is conservation of mass in the domain of integration.Therefore, we will apply biological no-flux condition (Plaza, 2019) (cid:90) V p p ( τ, ξ , v, ˆ v )ˆ v · n ( ξ ) dv d ˆ v = 0 , ∀ ξ ∈ ∂ Ω , τ > n ( ξ ) the outer normal to the boundary ∂ Ω in the point ξ . This class of boundary conditionsis part of the wider class of non-absorbing boundary conditions. At the macroscopic level (48)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 ∂ Ω . There are two important classes of kinetic boundary conditions which satisfy (48): the regularreflection boundary operators and the non-local (in velocity) boundary operators of diffusive type.We address the reader to the works by Palcewski (1992) and by Lods (2005) for the definitionof these boundary operators. In the present work, we shall consider Maxwell-type boundaryconditions which are prescribed in the form p ( τ, ξ , v (cid:48) , ˆ v (cid:48) ) = α ( ξ ) p ( τ, ξ , v, V (ˆ v )) + (1 − α ( ξ )) M ( ξ , v, ˆ v ) (cid:90) ˆ v ∗ · n ≥ p ( τ, ξ , v ∗ , ˆ v ∗ ) | ˆ v ∗ · n | dv ∗ d ˆ v ∗ , (49)where V (ˆ v ) = − ˆ v for the bounce back reflection condition and V (ˆ v ) = ˆ v − v · n ) n for thespecular reflection. M ( ξ , v, ˆ v ) is the Maxwellian function at the wall of Ω. The numerical scheme is the same used by Loy and Preziosi (2019). We consider a computationaldomain in the form Ω × V p where in the one dimensional case Ω = [ x min , x max ] and V p symmetrizedconsidering the two only possible directions along x and − x in ˜ V p = [ − U, U ] and in the twodimensional case Ω = [ x min , x max ] × [ y min , y max ] and V p = [0 , U ] × [0 , π ]. The computational13omain is discretized with a Cartesian mesh X h × V p h , where X h and V p h are defined by (in twodimensions) X h = { x ij = ( x i , y j ) = ( x min + i ∆ x, y min + j ∆ y ) , i = 0 , . . . , n x , j = 0 , . . . , n y } V p h = { v l,k = v k (cos θ l , sin θ l ) , θ l = ( j + 1 / θ, l = 0 , . . . , n ang − , v k = v + k ∆ v, k = 0 , . . . , 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 introducethe 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,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) j =0 p ni,j,l,k . Concerning boundary conditions, in the one-dimensional case we consider regular reflectingconditions. In one dimension, the bounce-back and the specular reflection boundary conditionscoincide, that is p ( t, x = x min , v ) = p ( t, x = x min , − v ) and p ( t, x = x max , − v ) = p ( t, x = x max , 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, ˆ v ) = T [ S , S (cid:48) ]( x , v, ˆ v ) , being T the asymptotic equilibrium of the system with this class (no-flux) of boundary conditions.Concerning the relaxation step, we remark that T [ S , S (cid:48) ] is defined as in Eq. (17) where Ψ isa probability density with a minimal variance, as numerically we can not represent a Dirac delta,even because this would require a weak formulation of all the equations. We then usedΨ[ S (cid:48) ]( x , v | ˆ v ) = C exp (cid:20) − k cos (cid:18) π ( v − ¯ U S(cid:48) ) U (cid:19)(cid:21) , v ∈ [0 , U ]that is a Von Mises distribution translated on the speed interval [0 , U ] with C a normalizationconstant. Its average is exactly ¯ U S(cid:48) that we compute using the definition (24). In order toapproximate a Dirac delta we shall consider small variances of Ψ and, then, large values of k .Numerically, though, the variance of Ψ will be larger then ∆ v being ∆ v the minimum distancebetween two possible values of the speed, because if it is smaller it causes numerical spuriouserrors. 14 Random polarization
The environmental cues which may represent physical limits of migration in general affect cellspeed. Therefore, we here first focus on the particular case in which there is no bias in the decisionof the direction of motion, corresponding to a random polarization, given by B [ S ]( x , ˆ v ) = const .As said by Loy and Preziosi (2019), this does not imply isotropy, because the speed can havedifferent density distributions on every direction ˆ v , as the distribution of S (cid:48) sensed ahead alongthe direction ˆ v may be different. Furthermore, in this article, such differences of S (cid:48) may also leadto different sensing radii in different directions ˆ v .In this section we will specify the model in the random polarization case and, eventually, we willintroduce some practical examples characterized by anisotropy on the sensing radius dependingon the direction, as well as on space and time.If in Eq.(15) T ˆ v λ [ S ] is independent of ˆ v , then the transition probability (17) simplifies to T [ S , S (cid:48) ]( x , v, ˆ v ) = 1 (cid:90) S d − Γ S(cid:48) ( x , ˆ v ) d ˆ v (cid:90) R M(cid:48)S(cid:48) ( x , ˆ v )0 γ S(cid:48) ( λ (cid:48) ) ψ ( x , v |S (cid:48) ( x + λ (cid:48) ˆ v )) dλ (cid:48) , (50)where we recall that Γ S(cid:48) ( x , ˆ v ) = (cid:90) R M(cid:48)S(cid:48) ( x , ˆ v )0 γ S (cid:48) ( λ (cid:48) ) dλ (cid:48) . The macroscopic velocity of the transition probability simplifies to U S(cid:48) ( x ) = 1 (cid:90) S d − Γ S(cid:48) ( x , ˆ v ) d ˆ v (cid:90) S d − Γ S(cid:48) ( x , ˆ v ) ¯ U S(cid:48) ( x | ˆ v )ˆ v d ˆ v , (51)and the variance-covariance tensor to D S(cid:48) ( x ) = 1 (cid:90) S d − Γ S(cid:48) ( x , ˆ v ) d ˆ v (cid:90) S d − Γ S(cid:48) ( x , ˆ v ) ¯ D S(cid:48) ( x | ˆ v )ˆ v ⊗ ˆ v d ˆ v − U S(cid:48) ( x ) ⊗ U S(cid:48) ( x ) , where ¯ U S(cid:48) and ¯ D S(cid:48) are respectively given by Eqs.(24) and (27).The fact that R M(cid:48)S(cid:48) and therefore Γ
S(cid:48) may depend on ˆ v , in general leads to different speedsin different directions. This means that, unless for very special cases, the 0 th -order of Eq.(51)does not vanish identically and therefore the proper scaling is hyperbolic. So, even if there is nodirectional sensing, in presence of a barrier there may be anisotropy due to a non-homogeneoussensing possibility of the cell in the different directions. In classical volume filling models the speed substituted in the mass balance equation or in theadvection-diffusion equation is a decreasing function of the density, eventually vanishing for den-sities above a critical value ρ th . To handle a similar case in the kinetic framework, we consider S (cid:48) = M (cid:48) = ρ and M th = ρ th . For example, volume filling effects can be described by a ψ withmean ¯ v ( x | ρ ) = ¯ v M (cid:18) − ρρ th (cid:19) + , (52)where ( · ) + = max( · ,
0) represents the positive part operator and ¯ v M is the maximal speed. Forsake of simplicity, we will call R ρ the quantity R M(cid:48)S(cid:48) defined in Eq.(14) and in the following equationwe explicitly stress that in presence of volume filling effects the sensing radius R ρ depends on timebecause the density of cells depends on the evolution of the distribution function p itself.In order to understand the application of the model, let us consider a density distribution likethe one in the top of Fig. 4 and for sake of simplicity with respect to the general discussion on the15igure 4: Qualitative behaviour of cell speed (bottom) in a landscape with a density that mightovercome the volume filling threshold (top). The length of the arrow close to the cells indicatesthe magnitude of mean speed, while the blockhead arrow the sensing radius R ρ depending on thepresence of the overcrowded area.maximum sensing radius done at the end of Section 3, let us take ¯ ε = ∆ = 0. In the discussion weinitially take γ ρ ( λ (cid:48) ) = H ( R ρ − λ (cid:48) ), where H is the Heaviside function, meaning that the speed isdetermined by uniformly averaging the response to the signal ρ in the direction ˆ v looking aheadup to a distance R ρ defined by (14). We will finally assume that all cells want to move to theright.Referring to Figure 4, we can qualitatively observe that cell A has all the available space tomove and will do it with the maximum allowed speed, up to when it reaches the location x because then it is feeling an overcrowded environment ahead. Because of that, it will then slowdown stopping upon reaching the point x . We observe that cell C has a sensing radius limitedby the presence of the overcrowded region ahead x . Cells D and E can not move because theyare in the middle of the jammed area ( R ρ = 0). Actually, cell E can sense that there would beavailable space beyond the point x , but even is it is close to the border of the jammed area it isstill. On the other hand, cell F can move to the right because it perceives available space ahead.Conversely, if it were polarized to go to the left, then it would not move because it is at the borderof an overcrowded region ahead. This clearly points out that the speed of cell F depends on itspolarization, an information that in the model is introduced starting from the dependence of R ρ from ˆ v . From a qualitative point of view, the resulting speed of a cell moving to the right in thelandscape given on the top of Fig. 4 is given at the bottom of the same figure.16f, instead, γ (cid:48) ρ ( λ (cid:48) ) = δ ( λ (cid:48) − R ρ ) cells will only look at a distance R ρ ahead without consideringthe information within or beyond that distance. So, in the landscape on the top of Figure 4 cell Awill move at the maximum speed, cell B at a lower speed (even lower than in the case of a Heavisidefunction, because it is not averaging over the sensing region) and will stop upon reaching the point x − R max ρ . Cells C, D and E will not move, while cell F will start moving. The behaviour of cellsB and C could be in principle justified by a will of stopping at a distance from the overcrowdedregion. As we shall see, because of this effect choosing a Dirac delta as sensing weight will giverise to patterns often characterized by a typical wavelength of the order of R max ρ .In this case the transition probability (50) writes T [ ρ ]( t, x , v, ˆ v ) = 1 (cid:90) S d − Γ ρ ( t, x , ˆ v ) d ˆ v (cid:90) R ρ ( t, x , ˆ v )0 γ ρ ( λ (cid:48) ) ψ ( x , v | ρ ( t, x + λ (cid:48) ˆ v )) dλ (cid:48) . (53)In particular, we recall that because of the definition (14), if ∆ = 0, R ρ ( t, x , ˆ v ) is such that[0 , R ρ ( t, x , ˆ v )] = (cid:26) λ (cid:48) ∈ [0 , R maxρ ] (cid:12)(cid:12)(cid:12)(cid:12) ρ ( t, x + λ (cid:48) ˆ v ) ≤ ρ th (cid:27) . (54)Therefore, in the integration interval [0 , R ρ ( t, x , ˆ v )]¯ v ( t, x | ρ ( t, x + λ (cid:48) ˆ v )) = ¯ v M (cid:18) − ρ ( t, x + λ (cid:48) ˆ v ) ρ th (cid:19) . The mean speed will then be ¯ U ρ ( t, x | ˆ v ) = ¯ v M (cid:18) − ¯ ρ ( t, x , ˆ v ) ρ th (cid:19) , (55)where ¯ ρ ( t, x , ˆ v ) = (cid:90) R ρ ( t, x , ˆ v )0 γ ρ ( λ (cid:48) ) ρ ( t, x + λ (cid:48) ˆ v ) dλ (cid:48) Γ ρ ( t, x , ˆ v ) , (56)measures the average density in the direction ˆ v until the threshold value is reached.Looking at the macroscopic speed, we have that U ρ ( t, x ) = ¯ v M (cid:90) S d − Γ ρ ( t, x , ˆ v ) d ˆ v (cid:18)(cid:90) S d − Γ ρ ( t, x , ˆ v )ˆ v d ˆ v − ρ th (cid:90) S d − Γ ρ ( t, x , ˆ v )¯ ρ ( t, x , ˆ v ) ˆ v d ˆ v (cid:19) , (57)and it vanishes in points where the density is above the threshold value.It is instructive to examine the limit case of small R maxρ , which implies that λ (cid:48) ≤ R maxρ issmall as well. Like in Section 4 we can perform a Taylor expansion of ρ ( x + λ (cid:48) ˆ v ) in a smallneighbourhood of x with size R ρ ρ ( t, x + λ (cid:48) ˆ v ) ≈ ρ ( t, x ) + λ (cid:48) ∇ ρ ( t, x ) · ˆ v . (58)Therefore, up to re-scaling,¯ U ρ ( τ, ξ | ˆ v ) = ¯ v M (cid:18) − ρ ( τ, ξ ) ρ th − (cid:15) Λ (cid:48) ρ th ˆ v · ∇ ρ ( τ, ξ ) (cid:19) , where Λ (cid:48) ( τ, ξ , ˆ v ) = 1Γ ρ ( τ, ξ , ˆ v ) (cid:90) R ρ ( τ, ξ , ˆ v )0 γ ρ ( λ (cid:48) ) λ (cid:48) dλ (cid:48) U ρ = ¯ v M (cid:18) − ρρ th (cid:19) , ¯ U ρ = Λ (cid:48) ρ th ˆ v · ∇ ρ. (59)Therefore the macroscopic velocities of 0 th - and 1 st -order away from overcrowded areas are, re-spectively, U ρ = ¯ v M (cid:90) S d − Γ ρ ( τ, ξ , ˆ v ) ˆ v d ˆ v (cid:90) S d − Γ ρ ( τ, ξ , ˆ v ) d ˆ v (cid:18) − ρρ th (cid:19) , (60)and U ρ = − ¯ v M ρ th T ∇ ρ , with T ( τ, ξ ) = 1 (cid:90) S d − Γ ρ ( τ, ξ , ˆ v ) d ˆ v (cid:90) S d − Γ ρ ( τ, ξ , ˆ v )Λ (cid:48) ( τ, ξ , ˆ v )ˆ v ⊗ ˆ v d ˆ v , (61)which is anisotropic if Γ ρ is not isotropic.In this limit, we can perform a parabolic scaling only where R ρ is even or it does not dependon ˆ v , so that the numerator in (60) vanish. If so, then the parabolic limit reads ∂ρ∂τ − ∇ · ( ρ T ∇ ρ ) = ∇ · (cid:20) µ ∇ · (cid:0) D ρ ρ (cid:1)(cid:21) , with D ρ = (cid:90) S d − ¯ D ρ ( ξ | ˆ v )ˆ v ⊗ ˆ v d ˆ v . Of course, this is barely the case in presence of physical limits of migration (see for example Fig.7), and we should perform a hyperbolic limit leading to Eq. (42).We can also discuss the proper choice of scaling by considering a nondimensionalization andshow that because of the dependence of the sensing radius on the direction, a diffusive time scalecan be hardly chosen uniformly in Ω. Therefore, a hydrodynamic limit will be the appropriateone. Let us now introduce l ρ = max ρ |∇ ρ | as the characteristic length of variation of ρ , and theparameter η = ¯ R ρ l ρ , (62)where ¯ R ρ is a reference sensing radius. We have seen that if η (cid:28)
1, then we may write Eqs. (58),(59), (60), and Eq. (61). This is not possible if η (cid:29)
1. We shall rescale the variables as ξ = x l ρ , ˜ v = vU ref , τ = tσ t , The time scale σ t can be chosen as a diffusion time T Diff scale or a drift time scale T Drift . Ingeneral we may write (Othmer and Hillen, 2000) T Diff = l ρ D , D = U ref ¯ µ , T Drift = l ρ 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 ρ = O (cid:18) U ref (cid:15) (cid:19) (63)18hich implies the hierarchy T Drift (cid:28) T Diff . (64)In the present case this is equivalent to η = ¯ R ρ l ρ < . 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 ρ U ref = η ¯ µ as the hierarchy (5.1) does not hold anymore. Hence, the following relation holds η ∼ (cid:15) (cid:29) . We observe that the choice U ref = ¯ U ρ would make the nondimensionalization depend on thedirection. Therefore, we shall choose as in (Loy and Preziosi, 2019) U ref = ¯ R ρ ¯ µ being ¯ µ a reference frequency. The same holds true for ¯ R ρ that can not be chosen equal to R ρ as thelatter depends on the direction. On the other hand, in the present case, we can not even consider¯ R ρ = R max ρ in the choice of U ref like in (Loy and Preziosi, 2019), because it varies considerably intime and space and it may also be different at a fixed point in space as it depends on the direction.In conclusion, the relation (63) cannot be considered everywhere in Ω and, therefore, the choiceof a diffusive time scale is not the proper one.Referring to Section 4.2 for details concerning the numerical integration of the kinetic model,we will now present some simulations focusing on how the model deals with the volume fillingeffect. In particular, aiming at checking the effect of limiting the sensing radius R ρ because ofovercrowding, we will perform simulation with R ρ = R maxρ regardless of the presence of thresholdsand limiting it because of the threshold. For sake of simplicity, we will refer to the former case asthe unlimited R model . Of course, as expected, the effect becomes visible when cell density getscloser to the threshold value ρ th . A second aim is to put in evidence the difference between usinga Dirac delta and a Heaviside function as weight function to evaluate the cell density.In Fig. 5 the initial macroscopic density is a small perturbation of the constant distribution, sothat it is always below the threshold value ρ th . This implies that initially R ρ = R max ρ everywhere.However, the non-homogeneous distribution of speed leads to the formation of overcrowded areasand therefore to the limitation of R maxρ in Fig. 5(a). We observe that the model hampers thecell density from going above ρ th . On the contrary, the unlimited R model does not succeed inimposing such a limit (Fig. 5(b)) and cells go over the threshold value (Fig. 5(c)).In Fig. 6(a), the sensing function is a Heaviside function. Averaging the cell density over theinterval [ x, x + R ρ ], it leads to smoother solutions compared to using a delta function that impliesusing only the information in x + R ρ in order to determine the new speed (see Fig. 6(b),(c)).The initial condition barely touches the threshold value ρ th = 0 .
5. So, the cell in the center areslower than those away from the center, forming a wave of crowded cells. A second peak forms ata distance R max ρ because cells there sense the advancing front ahead and slow down. This leads tothe formation of a pattern of characteristic size comparable with R max ρ . The fact that cell densityin Fig. 6(a) never reaches the threshold value implies that R ρ is always equal to the the maximumpossible value R max ρ . 19 a) (b) R ρ = R maxρ (c) R ρ = R maxρ Figure 5: Pattern formation when γ ρ is a Dirac delta starting from a density distribution ρ ( x ) =0 . .
02 sin π x below the threshold value ρ th = 0 .
25 and µ = 10. In (b) the sensing radius isalways R max ρ = 2 while in (a) it is limited by ρ th . Figure (c) is the same as (b), but the densityover the threshold value ρ th = 0 .
25 is highlighted in white. (a) γ ρ = H (b) γ ρ = δ (c) γ ρ = δ , R ρ = R maxρ Figure 6: Pattern formation starting from a density distribution ρ ( x ) = 0 . (cid:104) − ( x − . . (cid:105) thatjust reaches the threshold value ρ th = 0 .
5. In addition µ = 200. In (a) the sensing kernel is aHeaviside function, while in (b) it is a Dirac delta. In addition, in (c) the sensing radius is notlimited by ρ th and is always R max ρ = 0 .
6. 20n the other hand, when the sensing function is a δ , the pattern is stronger and is characterizedby maxima that reach ρ th = 0 .
5. We also observe a zig-zag behaviour which is a characteristic ofother alignment-repulsion-attraction models like in the works by Carrillo et al. (2015) and Eftimie(2012). This is due to the following dynamics: referring to Fig. 6(b) at t ≈
80, in a certainpoint cells start clustering, reaching their maximal local density ρ th . So, cell sensing this excessivecrowding prefer to move in the opposite direction of motion clustering in turn to values close to ρ th . As in the previous figure using an unlimited R model leads to higher peaks beyond ρ th andsharper fronts.In Fig. 7 the initial condition is not symmetric with respect to the midpoint of Ω, it is smallerthan ρ th for x < . x > . R ρ is initially zero in the overcrowded region x > . x < . R ρ is close to R max ρ = 0 . U ρ does not vanish implies that (41) is not satisfied and a diffusive limit can not be performed.In Fig. 8 the initial distribution is a Gaussian with maximum above ρ th = 0 .
5, so that there isovercrowding in the central region, specifically nearly between 1 . .
2. In the region where thedensity is below the threshold, the solution diffuses fast. In fact, cells initially located in x > . R ρ = R max ρ = 0 . R max ρ = 0 . R max ρ = 0 . R max ρ ) allowing motion. TheSupplementary movie VF.mp4 shows the time evolution of the density distribution p that startsfrom a distribution that is uniformly distributed in the velocity space and reaches the spatialuniform distribution. In order to move in a three-dimensional environment, cells interact with the ECM. This is networkof fibres that on one hand are used by cells to adhere and exert traction forces and on the otherhand can constitute an obstacle when the characteristic pore size drops below a threshold value.The combination of these interactions leads to a bimodal dependence of cell speed on the stiffnessand density of the ECM (Harley et al., 2008; Peyton and Putnam, 2005; Zaman et al., 2006). Inparticular, there is a threshold value M th above which cell can not move in the ECM (Wolf et al.,2007; Friedl et al., 2011). In some migration experiments, mainly on artificial scaffolds, it is alsoclear that when there is little or no possibility of building focal adhesion with the substratum thenagain cells can not exert active traction forces and are unable to move (Goodman et al., 1989;Nam et al., 2016). In terms of ECM we can then infer that there might also be a minimal density M of ECM necessary to crawl in it. There is however a difference between the two thresholdsbecause in this last case cells can extend their protrusions beyond the region lacking of adhesionpoints to possibly reach a farther region where they can adhere and exert traction.Denoting by R M the quantity R M(cid:48)S(cid:48) defined in Eq.(14) for the ECM case and referring to theECM landscape in Fig. 9, we can identify a cell response similar to overcrowding where there is21 a) (a) (b) (b)(c) (c) (d) (d)(e) (e) (f) (f)
Figure 7: Evolution from the initial condition ρ ( x ) = . π tanh(30( x − . .
4, so that ρ ( x ) <ρ th = 0 . x < .
5. The sensing kernel is a Heaviside function. The other parametersare R max ρ = 0 . µ = 200. In the top row the spatio-temporal evolution of the macroscopicvariables ρ and U ρ are plotted. In the second row that of the sensing radii in both directions. Inthe third row the mean speed in the two directions and the macroscopic velocity. Finally, in thebottom row, the distribution densities p at time 5 and 10. Positive and negative v ’s stand for thespeed distributions of cells oriented along x and − x .22 a) γ ρ = H (b) γ ρ = H (c) γ ρ = δ (d) γ ρ = H , R ρ = R max ρ (e) γ ρ = δ , R ρ = R max ρ Figure 8: Spatio-temporal evolution from the initial condition ρ ( x ) = 1 . (cid:104) − ( x − . . (cid:105) thathas a central region above ρ th = 0 .
5. Other parameters are R max ρ = 0 . µ = 2. In (a,b,d)density for a Heaviside sensing kernel and in (c,e) for a Dirac delta. In particular, in (d,e) R ρ isno limited and is always given by R max ρ . 23n ECM density above M th . Namely, for cell A the sensing radius R M = R maxM , while for cell Dit is R M < R maxM because cell speed is set according to the ECM density up to x . Cell E cannot move because it is in the middle of the dense ECM region, while cell F can move to the rightbecause it will encounter a microstructure allowing cell motion, but not to the left because in thatdirection the pore size is too small.However, as already stated, here a new phenomena occurs that is related to the existence ofregions with scarse presence of ECM. In Fig. 9 it occurs for x < x < x and x > x . In thesesituations a cell has no limitation of the sensing radius R maxM , but when it extends its protrusionsit is not able to build focal adhesion and exert traction forces in these interval, so the relatedcontribution to cell speed vanishes. As a consequence, typically the speed decreases, as for cell B.Now two situations may occur when a cell reaches this problematic area. If it is able to extendits protrusions to find a place to anchor and adhere, like cell C, it can exert traction forces andjump beyond the interval [ x , x ]. On the other hand, if it is not able to do so, like the red cellC’ that is characterized by a smaller R maxM , then it is stuck in x . Finally, a red cell located at adistance R maxM from x , like cell D’, barely touches the border in x . So, the interval [ x , x ] ischaracterized by a vanishing speed as [ x , x ] but for different dynamics, that will be shown in thesimulations to follow. For a similar reason the region beyond x can not be reached by both theyellow cell H and a red cell in the same location.Figure 9: Qualitative behaviour of cell speed (bottom) in a landscape with an ECM densityshown above. The grey regions denote problematic regions either because of too dense ECM (thecentral one) or because of lack of ECM (the lateral ones). The length of the arrow close to the cellsindicates the magnitude of mean speed, while the blockhead arrow indicates the sensing radius R M depending on the presence of an area with prohibitively small pore sizes.. The dashed lines closeto the cell denote intervals that do not contribute to cell speed. In the bottom graph the blackdotted line refers to the mean speed to the yellow (lighter) cell that can extend longer protrusions,while the full red line the one of the red (darker) cell, that has vanishing speed in the interval[ x , x ]. Both cell types have vanishing speed for x > x .In order to take into account of all the effects mentioned above, in the simulations to follow24e will use the following specific form of¯ v = ¯ v + 4(¯ v M − ¯ v )( M th − M ) [( M − M )( M th − M )] + , representing (if ¯ v = 0) the positive part of a parabola with zeros in M and M th with maximumspeed ¯ v M achieved for M = ( M th + M ) denoted by M max in the following, for sake of simplicity.The discussion of the macroscopic limit closely follows what presented for the volume filling case.For this reason it will not be repeated here.With the aim of showing the importance of considering a sensing radius defined as in (13), inFig. 10 we present a simulation where the value M th = 0 . R maxM from the border of the dense area (that is in x = 3). So, if R M is not limited asdefined in (13), they sense beyond the physical barrier and manage to come out of the dense zoneof ECM because they have a positive speed in the direction going out of the ECM. Conversely, if(13) is used cells are blocked in the dense area on the left.Figure 10: Behaviour of a cell distribution initially within a dense ECM area when R M is notlimited as defined by (13), but is always equal to R maxM = 1 .
5. The sensing kernel is a Dirac delta.Other parameters are M = 0, M th = 0 .
4, ¯ v M = 0 . M = 0 .
1, but the ECM profile (represented by the greenline) presents a central region [
P, Q ] with density smaller then 0 .
1. In Fig. 11(a), R maxM = 0 . P because, like the red cell in Fig. 9, they can not reachwith their protrusions the point Q where the ECM assumes again values larger then M . On theother hand, in Fig. 11(b), R maxM = 2 and, so, cells manage to reach the point Q and go over thehole. In fact, referring to Fig. 11(c),(d) for the cells that have short protrusions there is a pointwhere the speed (24) vanishes, whilst cells with longer protrusions move with a slower velocitywhen approaching the ECM depression but always keep a positive speed (Fig. 11(d)), so thatthey manage to overcome the problematic area. In Fig. 11(e) the density distribution p shows thedistribution of the microscopic velocities in the physical space while in Fig. 11 (f) pρ representsthe equilibrium probability measure.In Fig. 12 we present a set of two-dimensional simulations mimicking the experimental set-upin which a Petri dish is coated with thick stripes of different ECM components. For instance, in(Goodman et al., 1989) they show the different locomotion on laminin (or E8) and on fibronectin.25 a) R maxM = 0 . R maxM = 2(c) R maxM = 0 . R maxM = 2(e) R maxM = 0 . , p (f) R maxM = 0 . , pρ Figure 11: Spatio-temporal evolution of the cell density (a,b) in presence of a region with extremelylow density of ECM. The sensing kernel is always a Heaviside function and in (c,d) the mean speedis given with M = 0 . M th = 0 .
3, ¯ v M = 0 . µ = 20. In (a), (c), (d) and (e) R maxM = 0 .
2, whilein (b) and (d) R maxM = 2. In (e)-(f) we present the density distribution p and the probabilitymeasure pρ . 26he stripes made of fibronectin do not encourage locomotion, irrespective of the level of coatingconcentration used, whilst laminin and E8 encourage locomotion. In particular the locomotoryresponse is peaked around a certain range of values of laminin (20 − f mol · cm − ) whilst itdecreases for larger and smaller values of coating concentrations. Similar experiments are alsoperformed by Nam et al. (2016). Specifically, as shown in Fig. 12(a), the density of the laminincoating is 1 where it is present (red stripes) and zero elsewhere (blue stripes). In this case M = 0 .
01 and so the cells are not able to adhere in the blue stripes, whilst M max = 1, andso the cells have a maximal speed on the stripes of laminin. Cells start from an initial uniformdistribution in y < .
5. Cells on the ECM laminin stripes rapidly move along them and avoidgoing on the stripes lacking of laminin. Cells in the region where there is no laminin cannot movealong the blue stripes. However, the sensing radius is sufficiently high to allow all of them to grabsome adhesion sites and pull themselves onto the laminin stripes. Initially, (see Fig. 12(b)) cellscloser to the laminin stripes are faster and soon jump onto them to then steer and move alongthem. This is the reason why there is a minimum of cell density on the blue stripes close to theinterface and a faster progression on the red stripes. On the other hand, cells in the center areslower and take longer to move more or less perpendicularly to the stripes. In fact, cells not in themiddle of the blue stripe can reach the laminin stripes also if not oriented perpendicularly to them.Eventually, the entire population evolves along the stripes (Fig. 12(d)). (See Supplementary movieStripes.mp4). (a) Matrix density (b) t = 1(c) t = 10 (d) t = 35 Figure 12: Spatial evolution of cell density in an environment with ECM stripes distributed as in(a). Cells are initially uniformly distributed below the dashed line. As M = 0 .
01 and M max = 1,cells prefer to move along the stripes and as R max M = 0 .
25 those located in the middle of the bluestripe initially have some difficulties in reaching them. The sensing kernel is a Heaviside functionand ¯ v M = 1. Here the boundary conditions are those defined in (49) with α = 0 . Double bias
We shall now consider the case of cell migration under the action of both a field S affecting thechoice of the direction of motion (specifically either an external chemoattractant or an internaleffect due to cell-cell adhesion) and a field S (cid:48) affecting the speed (specifically related either tovolume filling effects or to the presence of ECM). In this case both sensing radii may depend ontime, position and direction. We then consider the operator J [ p ]( t, x , v, ˆ v ) = µ ( x ) (cid:32) ρ ( t, x ) c ( x ) (cid:90) R S ( t, x , ˆ v )0 γ S ( λ ) b ( S ( x + λ ˆ v )) dλ Ψ[ S (cid:48) ] − p ( t, x , v, ˆ v ) (cid:33) , (65)where Ψ is given by (16). Let us consider the case in which cell-cell adhesion represents a mechanism of cell polarizationbiasing the otherwise random motion of cells but taking into account that cells can not come tooclose because of volume filling, i . e . , S = S (cid:48) = ρ . We shall denote R S = R ρ,adh and R S(cid:48) = R ρ,vf .In particular, as the two cues are the same it is natural to take R ρ,adh = R ρ,vf and we will denoteit by R ρ . Figure 13 highlights the differences between when R ρ is defined as in Section 5.1 (see(a) and (c)) and when it is equal to the maximum possible extension of protrusions R max ρ . Inparticular, in Fig. 13(a,b), we observe that if the sensing function for the volume filling is a Diracdelta the difference is not so remarkable, with the formation of patterns of size close to R maxρ .We also highlight the formation of a cell-free zone close to the boundary because of the actionof adhesion forces that pull the aggregate together. On the other hand, if the sensing kernel forvolume filling is a Heaviside function, the two different choices of sensing radius determine theformation of a plateau as in Fig. 13 (c), or two aggregates as in Fig. 13 (d). Moreover in the lattercase, the cell density goes over the threshold ρ th . In all the other cases the cell density remainsunder the threshold value ρ th . In this case we assume that cells are sensitive to a chemoattractant with b ( S ) = S and takeinto account of volume filling effects with R ρ ( t, x , ˆ v ) given by (54) and mean speed ¯ U ρ ( ξ | ˆ v ) by(55). In general the appropriate macroscopic limit will be a hyperbolic one, which, dropping thedependence from space and time for sake of implicity, reads ∂ρ∂τ + ∇ · ¯ v M ρ (cid:90) S d − Γ S (ˆ v )Γ ρ (ˆ v ) ¯ S (ˆ v ) (cid:18) − ¯ ρ (ˆ v ) ρ M (cid:19) ˆ v d ˆ v (cid:90) S d − Γ S (ˆ v )Γ ρ (ˆ v ) ¯ S (ˆ v ) d ˆ v = 0 (66)where ¯ ρ is given by (56) and, analogously,¯ S (ˆ v ) = 1Γ S (ˆ v ) (cid:90) R ρ (ˆ v )0 γ S ( λ ) S ( ξ + λ ˆ v ) dλ is the weighted average of the signal S .In Fig. 14, we consider the same volume filling effect as in Fig. 8 under the action of a normallydistributed chemoattractant centered in x = 2 .
5. The sensing radius for the chemoattractant isthe same as the one for the volume filling effect defined as in (54) which is affected by the thresholdvalue ρ th . Due to the presence of the chemoattractant, the cell density does not converge to theconstant solution, but it remains above the threshold value as the chemoattractant and the volumefilling effect are balanced. 28 a) γ ρ,vf = δ, R ρ,adh = R ρ,vf (b) γ ρ,vf = δ, R ρ,adh = R max ρ,adh (c) γ ρ,vf = H, R ρ,adh = R ρ,vf (d) γ ρ,vf = H, R ρ,adh = R max ρ,adh Figure 13: Evolution in presence of cell-cell adhesion and volume filling with ρ th = 0 .
25 startingfrom the initial distribution ρ ( x ) = 0 . .
02 sin π x with µ = 2. The sensing kernel for adhesion isalways a Heaviside function, while that for volume filling is a Dirac delta (top row) or a Heavisidefunction (bottom row). In the left column a limited R ρ is used while in the right column themaximum protrusion length R max ρ = 0 .
2. In all figure, densities above the threshold value ρ th =0 .
25 are coloured in white. 29igure 14: Spatio-temporal evolution from the initial condition ρ ( x ) = 1 . (cid:104) − ( x − . . (cid:105) as inFig. 8 but in presence of a chemoattractant distributed as S ( x ) = 0 . (cid:104) ( x − . . (cid:105) . The sensingkernel is a Heaviside function. Other parameters are R max S = R max ρ = 0 . ρ th = 0 . µ = 2. In this section we shall consider the motion of a cell cluster toward a chemoattractant in a stronglyheterogeneous environment. Specifically, the initial distribution of cells is ρ ( x ) = 0 . (cid:20) − | x − x | . (cid:21) with x = (1 . , . , as in Fig. 15(a), that of the chemoattractant is C ( x ) = 50 exp (cid:20) − | x − x c | . (cid:21) with x c = (3 , , as in Fig. 15(b), and that of the ECM is M ( x ) = M b + ( M m − M b ) exp (cid:20) − | x − x M | . (cid:21) with x M = (2 , , as in Fig. 15(c). So, in a homogeneous environment cells would tend to move more or less alongthe diagonal of the square domain toward the maximum concentration of chemoattractant, whilein the heterogeneous case they would tend to avoid the region with too dense ECM.In the first simulation reported in the second row of Fig. 15 the ECM density is distributed ina way that it does not represent physical limits of migration. In fact, its maximum concentration M m = 1 . M th = 1 .
39 and its minimum concentration M b = 0 . M = 0 .
01 for crawling. So, cells are not blocked,but they are however slowed down as the maximum speed is achieved where M = M max = 0 . ECM .mp M m = 1 . M th = 0 .
8. So, cells can not enter the region30 a) Inital condition (b) Chemoattractant (c) Matrix density(d) t = 5 (e) t = 10 (f) t = 50(g) t = 5 (h) t = 10 (i) t = 50(j) t = 5 (k) t = 10 (l) t = 50 Figure 15: Motion of a cell cluster initially distributed as in (a) toward a chemoattractant withdistribution as in (b) passing through a dense ECM region as in (c) (In the plotted case M b = 0 . M m = 1 . R max S = R max M = 0 .
3. The turning frequency is µ = 10. (d-f) M m = 1 . M th = 1 . M b = 0 . M = 0 .
01, and M max = 0 .
7. (g-i) M m = 1 . M th = 0 . M b = 0 . M = 0 .
2, and M max = 0 .
5. The region within the white dashed circle is then too dense because
M > M th = 0 .
8. (j-l) M m = 0 . M th = 1 . M b = 0 is below M = 0 .
1, and M max = 0 . M < M = 0 . ECM .mp M m = 0 . M th = 0 .
8. On the other hand, the value M = 0 . R max M = 0 .
3, while motion in other directions is hamperedif not completely inhibited. Cells then cluster in the region with a comfortable density of ECMtowards the north-east side of the circle, because they are attracted towards the maximum chemicalconcentration. However, ahead they sense a region lacking of ECM for adhesion and remain stuck.(see also Supplementary Movie
ECM .mp The kinetic model developed in this article is based on the observation thati) cells sense their environment collecting chemical and mechanical cues by extending protru-sions that can be much longer than the cell diameter;ii) the information acquired determine cell polarization and speed;iii) the mechanisms governing cell polarization and speed depend on different intracellular mech-anisms;iv) the valuable information can be limited by the presence of physical limits of migration, suchas cell overcrowding, cell layers, like mesothelial or endothelial linings, basal membranes, orin general ECM with pores too small to be penetrated by the cell nucleus or even by itsprotrusions.From the modelling point of view, these points respectively imply that the kinetic model is char-acterized by (i) non-local turning operators ( i . e . , the integrals over λ and λ (cid:48) ), (ii) a probabilitydistribution that depends on a speed v and an orientation unit vector ˆ v , (iii) with a turning op-erator split in a part influencing v and a part influencing ˆ v depending on different sensing kernels( γ S and γ S(cid:48) ) and cues ( S and S (cid:48) ) (iv) operating on domains (identified by R MS and R M(cid:48)S(cid:48) ) that candepend on the presence of physical limits of migration, such as those mentioned above.The article shows in particular how important it is to handle and model properly situationsthat might look extreme but actually characterize many physiological and pathological situationsleading to cell aggregation and collective migration, to cell compartmentalization by basal mem-branes, to cell invasion when the membranes rupture or cells acquire a phenotype that allows themto pass through their narrow pores, leading to intravasation and extravasation of metastasis.The model is very flexible and was applied to situations taking into account of volume fillingeffects and cell-cell adhesion. Cell-matrix interaction was also considered both in the case of thickECM and when a lack of ECM might hamper the formation of focal adhesions that are essentialfor cell crawling. In the latter case it was shown that if the cell is able to extend its protrusionbeyond the problematic area to reach a region where focal adhesions can be formed again, then,due to the non-local character of the model, it is able to cross over the region with poor ECM.Otherwise, it is segregated close to the border of the area lacking of ECM. A virtual experimentof motion along stripes of laminin is also performed.Simulations show how pattern may form spontaneously from nearly homogeneous configura-tions, especially when the sensing kernel is a Dirac delta. The characteristic size of the pattern isrelated to the sensing radius. The presence of this effect calls for a stability analysis that will beperformed in a future article. 32nother important topic to be addressed that is not included in the present model is how tohandle different cues governing either sensing or speed. The former aspect refers, for instance, tocases in which cells adhere to each other while under the action of chemotaxis and/or haptotaxis.Being able to deal con temporarily with the two phenomena, together with volume filling, is funda-mental to deal with collective chemotaxis. The latter aspect refers, for instance, to overcrowdingdue to strongly heterogeneous distributions of ECM and, in particular, volume filling in presenceof basal membranes.
Acknowledgements
This work was partially supported by Istituto Nazionale di Alta Matematica, Ministry of Educa-tion, Universities and Research, through the MIUR grant Dipartimenti di Eccellenza 2018-2022and Compagnia di San Paolo that finances NL’s Ph.D. scholarship.
References
Abercrombie, M. and Heaysman, J. E. (1953). Observations on the social behaviour of cellsin tissue culture: I. Speed of movement of chick heart fibroblasts in relation to their mutualcontacts.
Experimental Cell Research , 5(1):111 – 131.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.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.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.Bitsouni, V. and Eftimie, R. (2018). Non-local parabolic and hyperbolic models for cell polarisationin heterogeneous cancer cell populations.
Bulletin of Mathematical Biology , 80(10):2600–2632.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.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. 33olombi, 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.Davidson, P. M., Denais, C., Bakshi, M. C., and Lammerding, J. (2014). Nuclear deformabilityconstitutes a rate-limiting step during cell migration in 3-d environments.
Cellular and MolecularBioengineering , 7(3):293–306.Devreotes, P. and Janetopoulos, C. (2003). Eukaryotic chemotaxis: Distinctions between direc-tional sensing and polarization.
Journal of Biological Chemistry , 278(23):20445–20448.Eftimie, R. (2012). Hyperbolic and kinetic models for self-organized biological aggregations andmovement: a brief review.
Journal of Mathematical Biology , 65(1):35–75.Eftimie, R., Perez, M., and Buono, P.-L. (2017). Pattern formation in a nonlocal mathematicalmodel for the multiple roles of the tgf- β pathway in tumour dynamics. Mathematical Biosciences ,289:96 – 115.Friedl, P., Wolf, K., and Lammerding, J. (2011). Nuclear mechanics during cell migration.
CurrentOpinion in Cell Biology , 1:55–64.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.Giverso, C., Grillo, A., and Preziosi, L. (2014). Influence of nucleus deformability on cell entryinto cylindrical structures.
Biomechanics and Modeling in Mechanobiology , 13:481–502.Goodman, S. L., Risse, G., and Mark, K. (1989). The E8 subfragment of laminin promoteslocomotion of myoblasts over extracellular matrix.
The Journal of Cell Biology , 109:799–809.Harley, B., Kim, H., Zaman, M., Yannas, I., Lauffenburger, D., and Gibson, L. J. (2008). Mi-croarchitecture of three-dimensional scaffolds influences cell migration behavior via junctioninteractions.
Biophysical journal , 95(8):40134024.Hillen, T. (2006). M5 mesoscopic and macroscopic models for mesenchymal motion.
Journal ofMathematical Biology , 53:585–616.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.Lods, B. (2005). Semigroup generation propertiesof streaming operators with noncontractiveboundary conditions.
Mathematical and Computer Modelling , 42:1441–1462.Loy, N. and Preziosi, L. (2019). Kinetic models with non-local sensing determining cell polar-ization and speed according to independent cues.
Journal of Mathematical Biology-accepted(arXiv:1906.11039v3) .Nam, K.-H., Kim, P., K. Wood, D., Kwon, S., P. Provenzano, P., and Kim, D.-H. (2016). Multiscalecues drive collective cell migration.
Scientific Reports , 6:29749.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:311–338.34thmer, H. G., Dunbar, S. R., and Alt, W. (1988). Models of dispersal in biological systems.
Journal of Mathematical Biology , 26(3):263–298.Othmer, H. G. and Hillen, T. (2000). The diffusion limit of transport equations derived fromvelocity-jump processes.
SIAM Journal of Applied Mathematics , 61:751–775.Painter, J. K. and Hillen, T. (2002). Volume-filling and quorum-sensing in models for chemosen-sitive movement.
Canadian Applied Mathematics Quarterly , 10:501–543.Painter, K. J., Armstrong, N. J., and Sherratt, J. A. (2010). The impact of adhesion on cellularinvasion processes in cancer and development.
Journal of Theoretical Biology , 264(3):1057–1067.Painter, K. J., Bloomfield, M. J., Sherratt, J. A., and Gerisch, A. (2015). A nonlocal model forcontact attraction and repulsion in heterogeneous cell populations.
Bulletin of MathematicalBiology , 77:1132–1165.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.Pettersson, R. (2004). On solutions to the linear Boltzmann equation for granular gases.
TransportTheory and Statistical Physics , 33(5-7):527–543.Peyton, S. R. and Putnam, A. J. (2005). Extracellular matrix rigidity governs smooth muscle cellmotility in a biphasic fashion.
Journal of Cellular Physiology , 204(1):198–209.Plaza, R. G. (2019). Derivation of a bacterial nutrient-taxis system with doubly degenerate cross-diffusion as the parabolic limit of a velocity-jump process.
Journal of Mathematical Biology ,78:1681–1711.Schmeiser, C. and Nouri, A. (2017). Aggregated steady states of a kinetic model for chemotaxis.
Kinetic and Related Models , 10(1):313 – 327.Schoumacher, M., Goldman, R. D., Louvard, D., and Vignjevic, D. M. (2010). Actin, microtubules,and vimentin intermediate filaments cooperate for elongation of invadopodia.
The Journal ofCell Biology , 189(3):541–556.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.Scianna, M., Preziosi, L., and Wolf, K. (2013). A cellular Potts model simulating cell migrationon and in matrix environments.
Mathematical Biosciences and Engineering , 10:235–261.Shankar, J., Messenberg, A., Chan, J., Underhill, T. M., Foster, L. J., and Nabi, I. R. (2010).Pseudopodial actin dynamics control epithelial-mesenchymal transition in metastatic cancercells.
Cancer Research , 70(9):3780–3790.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.35osin, 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.
The Jour-nal of Cell Biology , 201:1069–1084.Wolf, K., Wu, Y. I., Liu, Y., Geiger, J., Tam, E., Overall, C., Stack, M. S., and Friedl, P. (2007).Multi-step pericellular proteolysis controls the transition from individual to collective cancercell invasion.
Nature Cell Biology , 9(8):893–904.Zaman, M. H., Trapani, L. M., Sieminski, A. L., Mackellar, D., Gong, H., Kamm, R. D., Wells,A., Lauffenburger, D. A., and Matsudaira, P. (2006). Migration of tumor cells in 3D matricesis governed by matrix stiffness along with cell-matrix adhesion and proteolysis.