Kinetic models with non-local sensing determining cell polarization and   speed according to independent cues

Nadia Loy; Luigi Preziosi

arXiv:1906.11039·q-bio.CB·July 2, 2019

Kinetic models with non-local sensing determining cell polarization and speed according to independent cues

Nadia Loy, Luigi Preziosi

PDF

TL;DR

This paper develops a kinetic model for cell movement that incorporates non-local sensing cues influencing both direction and speed, analyzing how environmental factors affect cell polarization and motility.

Contribution

It introduces a novel velocity-jump process model with double bias based on non-local environmental cues, linking microscopic sensing to macroscopic cell behavior.

Findings

01

The model captures how external cues influence cell polarization.

02

Analysis of cell size and sensing mechanisms on population dynamics.

03

Comparison between kinetic and macroscopic models validates the approach.

Abstract

Cells move by run and tumble, a kind of dynamics in which the cell alternates runs over straight lines and re-orientations. This erratic motion may be influenced by external factors, like chemicals, nutrients, the extra-cellular matrix, in the sense that the cell measures the external field and elaborates the signal eventually adapting its dynamics. We propose a kinetic transport equation implementing a velocity-jump process in which the transition probability takes into account a double bias, which acts, respectively, on the choice of the direction of motion and of the speed. The double bias depends on two different non-local sensing cues coming from the external environment. We analyze how the size of the cell and the way of sensing the environment with respect to the variation of the external fields affect the cell population dynamics by recovering an appropriate macroscopic limit…

Tables2

Table 1. Table 1: Summary of models

[-10pt] Biological aspect	Section	$𝒥$	$\frac{Macro}{limit}$	Application	Figure
Orientation
Local sensing	5.2.1	(80)	(83)
$\frac{Encounter rate depending}{on the incoming velocity}$	5.2.2	(85)	(88)
Non-local sensing	5.2.3	(89)		Chemotaxis	1, 2, 4
				Cell adhesion	5
				Durotaxis	7
(small sensing limit)		(96)	(98)–(99)
Comparative sensing	5.2.4	(107)		Chemotaxis	3
(small sensing limit)		(112)	(98)–(113)
Speed
Random polarization	6.1.1	(130)	(131)	ECM steryc hindrance	9
(Non-local speed sensing)
Orientation and Speed
Non-local orientation	6.1.3	(133)	(134)	Chemotaxis and	11
and speed sensing				ECM steryc hindrance

Table 2. Table 2: Chemotactic motion.

[-10pt] $l_{𝒮}$	$R_{𝒮}$	$γ_{𝒮}$	$Λ$	$η$	relative $L_{\infty}$ error
17.964	0.02	$δ$	0.02	0.001	0.03
17.964	0.02	H	0.01	0.001	0.03
17.964	0.2	$δ$	0.2	0.01	0.03
17.964	0.2	H	0.1	0.01	0.05
17.964	2	$δ$	2	0.1	0.13
17.964	2	H	1	0.1	0.53

Equations373

\frac{\partial p}{\partial t} (t, x, v_{p}) + v \cdot \nabla p (t, x, v_{p}) = J [p] (t, x, v_{p})

\frac{\partial p}{\partial t} (t, x, v_{p}) + v \cdot \nabla p (t, x, v_{p}) = J [p] (t, x, v_{p})

ρ (t, x) = \int_{V_{p}} p (t, x, v_{p}) d v_{p};

ρ (t, x) = \int_{V_{p}} p (t, x, v_{p}) d v_{p};

U (t, x) = \frac{1}{ρ ( t , x )} \int_{V_{p}} p (t, x, v_{p}) v d v_{p};

U (t, x) = \frac{1}{ρ ( t , x )} \int_{V_{p}} p (t, x, v_{p}) v d v_{p};

P (t, x) = \int_{V_{p}} (v - U) \otimes (v - U) p (t, x, v_{p}) d v_{p};

P (t, x) = \int_{V_{p}} (v - U) \otimes (v - U) p (t, x, v_{p}) d v_{p};

E (t, x) = \int_{V_{p}} \frac{∣ v - U ∣ ^{2}}{2} p (t, x, v_{p}) d v_{p} .

E (t, x) = \int_{V_{p}} \frac{∣ v - U ∣ ^{2}}{2} p (t, x, v_{p}) d v_{p} .

\int_{V_{p}} f (v, \hat{v}) d v_{p} = \int_{0}^{U} \int_{S^{d - 1}} f (v, \hat{v}) d \hat{v} d v .

\int_{V_{p}} f (v, \hat{v}) d v_{p} = \int_{0}^{U} \int_{S^{d - 1}} f (v, \hat{v}) d \hat{v} d v .

\int_{V_{p}} f (v, \hat{v}) d v_{p} = \int_{0}^{U} f^{1} (v) d v \int_{S^{d - 1}} f^{2} (\hat{v}) d \hat{v} .

\int_{V_{p}} f (v, \hat{v}) d v_{p} = \int_{0}^{U} f^{1} (v) d v \int_{S^{d - 1}} f^{2} (\hat{v}) d \hat{v} .

\int_{S^{d - 1}} f^{2} (\hat{v}) d \hat{v} = \int_{0}^{2 π} f^{2} (cos (θ), sin (θ)) d θ

\int_{S^{d - 1}} f^{2} (\hat{v}) d \hat{v} = \int_{0}^{2 π} f^{2} (cos (θ), sin (θ)) d θ

J [p] (x, v_{p}) = G [p] (x, v_{p}) - L [p] (x, v_{p}),

J [p] (x, v_{p}) = G [p] (x, v_{p}) - L [p] (x, v_{p}),

G [p] (x, v_{p}) = \int_{V_{p}} μ (x, v_{p}^{'}) T [S, S^{'}] (x, v_{p} ∣ v_{p}^{'}) p (t, x, v_{p}^{'}) d v_{p}^{'},

G [p] (x, v_{p}) = \int_{V_{p}} μ (x, v_{p}^{'}) T [S, S^{'}] (x, v_{p} ∣ v_{p}^{'}) p (t, x, v_{p}^{'}) d v_{p}^{'},

L [p] (x, v_{p}) = \int_{V_{p}} μ (x, v_{p}) T [S, S^{'}] (x, v_{p}^{'}^{'} ∣ v_{p}) p (t, x, v_{p}) d v_{p}^{'}^{'},

L [p] (x, v_{p}) = \int_{V_{p}} μ (x, v_{p}) T [S, S^{'}] (x, v_{p}^{'}^{'} ∣ v_{p}) p (t, x, v_{p}) d v_{p}^{'}^{'},

\int_{V_{p}} T [S, S^{'}] (x, v_{p}^{'}^{'} ∣ v_{p}) d v_{p}^{'}^{'} = 1, \forall x \in Ω, \forall v_{p} \in V_{p} .

\int_{V_{p}} T [S, S^{'}] (x, v_{p}^{'}^{'} ∣ v_{p}) d v_{p}^{'}^{'} = 1, \forall x \in Ω, \forall v_{p} \in V_{p} .

J [p] (x, v_{p}) = \int_{V_{p}} μ (x, v_{p}^{'}) p (t, x, v_{p}^{'}) d v_{p}^{'} T [S, S^{'}] (x, v_{p}) - μ (x, v_{p}) p (t, x, v_{p}),

J [p] (x, v_{p}) = \int_{V_{p}} μ (x, v_{p}^{'}) p (t, x, v_{p}^{'}) d v_{p}^{'} T [S, S^{'}] (x, v_{p}) - μ (x, v_{p}) p (t, x, v_{p}),

\mathcal{J}[p](t,{\bf x},{\bf v}_{p})=\mu({\bf x})\,\Big{(}\rho(t,{\bf x})T[\mathcal{S},\mathcal{S}^{\prime}]({\bf x},{\bf v}_{p})-p(t,{\bf x},{\bf v}_{p})\Big{)}\,.

\mathcal{J}[p](t,{\bf x},{\bf v}_{p})=\mu({\bf x})\,\Big{(}\rho(t,{\bf x})T[\mathcal{S},\mathcal{S}^{\prime}]({\bf x},{\bf v}_{p})-p(t,{\bf x},{\bf v}_{p})\Big{)}\,.

p (t, x, v_{p}) = ρ (t, x) T [S, S^{'}] (x, v_{p}),

p (t, x, v_{p}) = ρ (t, x) T [S, S^{'}] (x, v_{p}),

\int_{V_{p}} J [p] (t, x, v_{p}) d v_{p} = 0,

\int_{V_{p}} J [p] (t, x, v_{p}) d v_{p} = 0,

\int_{V_{p}}\mathcal{J}[p](t,{\bf x},{\bf v}_{p})\,{\bf v}\,d{\bf v}_{p}=\mu({\bf x})\rho(t,{\bf x})\big{(}{\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}({\bf x})-{\bf U}(t,{\bf x})\big{)}\,

\int_{V_{p}}\mathcal{J}[p](t,{\bf x},{\bf v}_{p})\,{\bf v}\,d{\bf v}_{p}=\mu({\bf x})\rho(t,{\bf x})\big{(}{\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}({\bf x})-{\bf U}(t,{\bf x})\big{)}\,

U_{S, S^{'}} (x) = \int_{V_{p}} T [S, S^{'}] (x, v_{p}) v d v_{p} .

U_{S, S^{'}} (x) = \int_{V_{p}} T [S, S^{'}] (x, v_{p}) v d v_{p} .

D_{S, S^{'}} (t, x) = \int_{V_{p}} (v - U_{S, S^{'}}) \otimes (v - U_{S, S^{'}}) T [S, S^{'}] (x, v_{p}) d v_{p} .

D_{S, S^{'}} (t, x) = \int_{V_{p}} (v - U_{S, S^{'}}) \otimes (v - U_{S, S^{'}}) T [S, S^{'}] (x, v_{p}) d v_{p} .

γ_{S} : R_{+} ⟼ R_{+} and γ_{S^{'}}^{'} : R_{+} ⟼ R_{+},

γ_{S} : R_{+} ⟼ R_{+} and γ_{S^{'}}^{'} : R_{+} ⟼ R_{+},

T [S, S^{'}] (x, v_{p}) = c (x) \int_{R_{+}} γ_{S} (λ) T_{λ}^{\hat{v}} [S] (x) d λ \int_{R_{+}} γ_{S^{'}} (λ^{'}) ψ (x, v ∣ S^{'} (x + λ^{'} \hat{v})) d λ^{'} .

T [S, S^{'}] (x, v_{p}) = c (x) \int_{R_{+}} γ_{S} (λ) T_{λ}^{\hat{v}} [S] (x) d λ \int_{R_{+}} γ_{S^{'}} (λ^{'}) ψ (x, v ∣ S^{'} (x + λ^{'} \hat{v})) d λ^{'} .

\mathcal{T}_{\lambda}^{\hat{{\bf v}}}[\mathcal{S}]({\bf x})=b\big{(}\mathcal{S}({\bf x}+\lambda\hat{{\bf v}})\big{)}\,,

\mathcal{T}_{\lambda}^{\hat{{\bf v}}}[\mathcal{S}]({\bf x})=b\big{(}\mathcal{S}({\bf x}+\lambda\hat{{\bf v}})\big{)}\,,

c (x) = \frac{1}{\int _{S^{d - 1}} Γ _{0}^{'} \int _{R_{+}} γ _{S} ( λ ) T _{λ}^{\hat{v}} [ S ] ( x ) d λ d v ^}

c (x) = \frac{1}{\int _{S^{d - 1}} Γ _{0}^{'} \int _{R_{+}} γ _{S} ( λ ) T _{λ}^{\hat{v}} [ S ] ( x ) d λ d v ^}

\begin{array}[]{lr}\tau=\epsilon^{2}t,\quad{\bf\boldsymbol{\xi}}=\epsilon{\bf x},\\[8.0pt] \tau=\epsilon t,\quad{\bf\boldsymbol{\xi}}=\epsilon{\bf x}.\end{array}

\begin{array}[]{lr}\tau=\epsilon^{2}t,\quad{\bf\boldsymbol{\xi}}=\epsilon{\bf x},\\[8.0pt] \tau=\epsilon t,\quad{\bf\boldsymbol{\xi}}=\epsilon{\bf x}.\end{array}

T [S, S^{'}] (ξ, v_{p}) = T [S, S^{'}]_{0} (ξ, v_{p}) + ϵ T [S, S^{'}]_{1} (ξ, v_{p}) + O (ϵ^{2})

T [S, S^{'}] (ξ, v_{p}) = T [S, S^{'}]_{0} (ξ, v_{p}) + ϵ T [S, S^{'}]_{1} (ξ, v_{p}) + O (ϵ^{2})

\int_{V_{p}} T [S, S^{'}]_{0} (ξ, v_{p}) d v_{p} = 1,

\int_{V_{p}} T [S, S^{'}]_{0} (ξ, v_{p}) d v_{p} = 1,

\int_{V_{p}} T [S, S^{'}]_{i} (ξ, v_{p}) d v_{p} = 0 \forall i \geq 1.

U_{S, S^{'}} = U_{S, S^{'}}^{0} + ϵ U_{S, S^{'}}^{1} + O (ϵ^{2})

U_{S, S^{'}} = U_{S, S^{'}}^{0} + ϵ U_{S, S^{'}}^{1} + O (ϵ^{2})

D_{S, S^{'}} = D_{S, S^{'}}^{0} + ϵ D_{S, S^{'}}^{1} + O (ϵ^{2})

D_{S, S^{'}} = D_{S, S^{'}}^{0} + ϵ D_{S, S^{'}}^{1} + O (ϵ^{2})

U_{S, S^{'}}^{i} = \int_{V_{p}} T [S, S^{'}]_{i} v d v_{p}

U_{S, S^{'}}^{i} = \int_{V_{p}} T [S, S^{'}]_{i} v d v_{p}

D_{S, S^{'}}^{i} = \int_{V_{p}} T [S, S^{'}]_{i} (v - U_{S, S^{'}}^{i}) \otimes (v - U_{S, S^{'}}^{i}) d v_{p}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Kinetic models with non-local sensing determining cell polarization and speed according to independent cues

Nadia Loy

Luigi Preziosi

Nadia Loy Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy, and Department of Mathematics “G. Peano”, Via Carlo Alberto 10 ,10123 Torino, Italy ([email protected])

Luigi Preziosi Department of Mathematical Sciences “G. L. Lagrange”, Dipartimento di Eccellenza 2018-2022, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy ([email protected])

Abstract

Cells move by run and tumble, a kind of dynamics in which the cell alternates runs over straight lines and re-orientations. This erratic motion may be influenced by external factors, like chemicals, nutrients, the extra-cellular matrix, in the sense that the cell measures the external field and elaborates the signal eventually adapting its dynamics.

We propose a kinetic transport equation implementing a velocity-jump process in which the transition probability takes into account a double bias, which acts, respectively, on the choice of the direction of motion and of the speed. The double bias depends on two different non-local sensing cues coming from the external environment. We analyze how the size of the cell and the way of sensing the environment with respect to the variation of the external fields affect the cell population dynamics by recovering an appropriate macroscopic limit and directly integrating the kinetic transport equation. A comparison between the solutions of the transport equation and of the proper macroscopic limit is also performed.

1 Introduction

It is well known that a classical migration mode of bacteria consists in alternating runs or swims over straigth lines and tumbles (Berg,, 1983). Also eukariotic cells alternate persistent crawlings along a polarization axis with reorientation phases in which they re-polarize as a result of several external stimuli. These can be represented by nutrients, chemical factors, or noxious substances (Block et al.,, 1983), and also by environmental cues, such as density, stiffness, and structure of the extra-cellular matrix (ECM). In addition, the presence of other cells along the way can act in a two fold-way, either attractive due to the interaction of adhesion molecules (e.g, cadherin complexes) expressed on the cellular membrane, or repulsive when the region is becoming too overcrowded.

Cells measure all these signals by transmembrane receptors located on their protrusions that can extend up to several cell diameters. The captured chemical or physical signals activate in turn transduction pathways downstream that lead to the cell response. This includes in particular the polarization of the cell with the formation of a “head” and a “tail” and the activation of adhesion molecules and traction forces leading to motion (Adler,, 1966). In the framework of kinetic models, in the present paper we will focus on how the environmental sensing over a finite radius can be translated in the choice of the direction of motion and of the speed of the cell. In fact, in order to do that we develop a kinetic model in which the distribution function depends on cell speed and orientation, respectively a scalar and a unit vector, in addition to the usual dependence on space and time. In the tumbling phase, then, the new orientation and subsequent speed are chosen as a result of a sensing over a finite neighbourhood of the cell, giving the kinetic model a non-local character.

Focusing on mesoscopic models of cell migration and referring to Hillen and Painter, (2008) for an extensive and more general review on PDE models of chemotaxis, we can observe that the run and tumble dynamics can be modeled by a stochastic process called velocity-jump process that was introduced by Stroock, (1974). In this work, the author derived a linear transport equation from the velocity-jump process. Such an equation describes the movement of a single-particle distribution function like in the Boltzmann equation (Cercignani,, 1987). The main elements of such a process are the tumbling frequency, the mean speed, and the transition probability that describes the probability of choosing a certain velocity after re-orientation. The mean speed, mean runtime (which is the inverse of the frequency), and the tumbling probability may be measured from individual patterns of members of the population. Another advantage of such models is that they admit finite propagation speed. Alt, (1980) and Othmer et al., (1988) introduced the bias induced by an external stimuli in the linear transport equation. Dickinson, (2000) generalized the model to a transport equation in an anisotropic environment and applied to an environment with a gradient of a stimulus.

Hillen, (2006) proposed a Boltzmann-like model taking into account of the interaction between cells and ECM, also including the degradation of ECM fibers as a function of the angle between the ECM fiber and the cell velocity. In particular, the author recovered macroscopic limits (diffusive or hyperbolic) according to the structure of the matrix. Painter, (2008) used a similar model to describe ECM remodelling resulting from the migration of cells in a heterogeneous environment by integrating numerically the kinetic equation.

An extension of this Boltzmann-like model taking also into account of the interaction between cells and between cells and ECM is proposed by Chauviere et al., 2007a and Chauviere et al., 2007b . A further extension is suggested by Chauviere and Preziosi, (2010) where the operators describing interactions between cells and between cells and ECM include the bias due to an external cue. Surulescu and coworkers (Engwer et al.,, 2017; Engwer et al., 2015a, ; Engwer et al., 2015b, ; Engwer et al.,, 2016; Stinner et al.,, 2015) instead discussed in more details cell-ECM adhesion phenomena, and, focusing on tumour growth, mainly gliomas, also included proliferation and therapies.

Many of the papers above determine for the proposed mesoscopic models the related macroscopic, usually diffusive, limits and implement them numerically. The diffusive limit of transport equations is deeply described by Hillen and Othmer, (2000). More in detail, Othmer and Hillen, (2002) investigated various forms and orders of the bias which may be included in the transition probability. However, in some regimes a hyperbolic limit is more suitable (for example in presence of formation of networks, as studied by Chalub et al., (2004), Hwang et al., (2005), Filbet et al., (2005)).

Coming to the sensing strategies and then to the determination of cell re-orientation and speed, they have been considered in several models. In position jump processes, the transition probabilities model the strategy of sensing and may typically be (i) local, if the information at the present position is considered, (ii) neighbour based, if the information at the target jump site is taken into account, (iii) gradient based, if the local difference in information between the target and the local site is considered, (iv) local average, if the average of the information between the present and target site as done by Othmer and Stevens, (2001) and by Painter and Sherratt, (2003).

In general, a way of including cell sensing is to consider a non-local average of the external fields. For what concerns bacteria and cells dynamics, Othmer and Hillen, (2002) and Hillen et al., (2007) introduced a finite sampling radius and they defined a non-local gradient as the average of the external field on a surface which represents the membrane of the cell. In particular, the Authors derived a macroscopic diffusive model with a non-local gradient both from a position jump process and from a velocity jump process by postulating that the non-local sensing is a bias of higher order.

The notion of non-local sensing was also used for cell adhesion and haptotaxis by Armstrong et al., (2006) yielding a macroscopic integro-differential equation. This model was recently derived from a space jump process by Buttenschön et al., (2018). Similar equations were proposed in 2D set-ups by Colombi et al., (2017), Colombi et al., (2015), and applied also to model crowd dynamics and traffic flow by Tosin and Frasca, (2011).

Eftimie, (2012) and Carrillo et al., (2015) proposed non-local kinetic models including repulsion, alignement and attraction. They distinguish cells going along the two directions in the one-dimensional set-up with the possibility of switching between the two directions with a constant speed. The one-dimensional kinetic model takes then a discrete structure that is then integrated directly.

As already stated, in this paper we want to include the sensing strategy in a velocity-jump process and to model the tactic and kinetic response of a cell as a consequence of the non-local sensing of the environment. In particular, we focus on the fact that different fields can influence respectively cell direction and speed. Because of these characteristics we name the model as a double-bias non-local kinetic model. To be specific, a chemical factor may influence the orientation of a cell (taxis), while the cell population density or the matrix density may influence the speed of the cell once the direction is chosen. For instance, a volume filling effect can hamper the motion in a certain direction due to cell overcrowding. Similarly, fiber networks with too narrow pores can hamper or even completely inhibit migration, a phenomenon called physical limit of migration by te Boekhorst et al., (2016) and Wolf et al., (2013).

Therefore, cell response is related to the choice of the transition probability evaluated in a non-local way over the neighborhood of the cell. Several examples are given, some in order to show the link with existing models (Buttenschön et al.,, 2018; Hillen et al.,, 2007; Painter et al.,, 2010, 2015; Painter and Hillen,, 2002), others to point out the novelties of the proposed model. In particular, from the numerical point of view, we simulate directly the kinetic equation taking into account of different phenomena, ${\it e.g.}$ , chemotaxis, ECM steryc hindrance, adhesion induced aggregation, durotaxis.

We also discuss how the choices made on the transition probability and on the size of the sampling volume determine the type of macroscopic limit that can be performed. Then, we compare the outcome of the numerical integration of the transport equation with the proper macroscopic limits and the related analytical solution that can be obtained for the stationary case. This shows a satisfactory agreement when the hypothesis necessary to perform the limits are satisfied.

The plan of the article is as follows. After a section introducing the general and preliminary concepts, Section 3 introduces the general structure of the double-bias non-local kinetic model. Then, in Section 4 the macroscopic limits are discussed, distinguishing when it is possible to perform a parabolic scaling or a hyperbolic scaling according to the modelling assumptions. In Section 5 we start with a simpler single-bias model in which there is no cue affecting the speed, that is determined through a given distribution function. Several sensing strategies determining cell orientation are discussed, starting from local responses to non-local averages and strategies comparing the level of chemicals in different points of the neighbourhood, also including the possible dependence of the encounter rate from taxis factors. In Section 6 the double-bias model is discussed. Again, for clarity we first focus on the simpler case in which there is no cue influencing cell orientation, so that it is randomly chosen. We finally give an example for the complete model, already presented in its general form in Section 3.

2 Preliminaries

The statistical description of the cell population is performed through the distribution density $p=p(t,{\bf x},{\bf v}_{p})$ which is parametrized by the time $t>0$ , the position ${\bf x}\in\Omega\subseteq\mathbb{R}^{d}$ and the velocity pair ${\bf v}_{p}=(v,\hat{{\bf v}})\in V_{p}=[0,U]\times\mathbb{S}^{d-1}$ , being $\mathbb{S}^{d-1}$ the unit sphere boundary in $\mathbb{R}^{d}$ and $U$ the maximal cell speed, so that $\hat{{\bf v}}$ is the unit vector and $v$ is the speed, ${\it i.e.}$ , the modulus of the velocity vector ${\bf v}=v\hat{{\bf v}}$ . The choice of representing the distribution function $p$ depending on the direction and on the speed of the velocity, instead of the velocity vector ${\bf v}$ , lies in the need of separating the mechanisms governing cell polarization and motility, for instance in response of chemotaxis and in presence of other influencing factors (Alt,, 1980). The mesoscopic model consists in the transport equation for the cell distribution:

[TABLE]

where the operator $\nabla$ denotes the spatial gradient. The term $\mathcal{J}[p](t,{\bf x},{\bf v}_{p})$ , named turning operator, is an integral operator that describes the change in velocity which is not due to the free-particle transport. It may describe the classical run and tumble behaviors, contact guidance phenomena, or cell-cell interactions. For the moment we will consider the classical run and tumble, ${\it e.g.}$ , random re-orientations, which, however, may be biased by external cues. Therefore, our turning operator will be the implementation of a velocity-jump process in a kinetic transport equation as introduced by Alt, (1980).

A macroscopic description for the cell population can be classically recovered through the definition of moments of the distribution function $p$ :

the cell number density

[TABLE]

the cell mean velocity

[TABLE]

the cell variance-covariance matrix

[TABLE]

the cell speed variance

[TABLE]

We remark that, because of the definition of $V_{p}$ , the integrals over the velocity space are

[TABLE]

In particular, if the dependence of $f$ on $v$ and $\hat{{\bf v}}$ can be factorized, ${\it i.e.}\ f(v,\hat{{\bf v}})=f^{1}(v)f^{2}(\hat{{\bf v}})$ , we have that

[TABLE]

The integral over $\mathbb{S}^{d-1}$ , which we remind to be the boundary of the unit sphere, is not a surface integral. It has to be interpreted as

[TABLE]

in 2D and similarly in 3D.

Computing the moments of the transport equation (1) allows to obtain evolution equations for the moments above. Macroscopic limits can then be achieved by classical procedures. In particular, the diffusive limit of transport equations with velocity jump processes is deeply treated by Hillen and Othmer, (2000), Othmer and Hillen, (2002), and Hillen, (2006), where the hyperbolic scaling is treated as well. The interest here is focused on when, according to the hypothesis done on the transition probability and on the sensing radius, it is possible to perform either a parabolic or a hyperbolic scaling, leading to the related diffusive and hyperbolic limits.

3 The turning operator

In taxis processes cells are capable of detecting and measuring external signals through membrane receptors located along cell protrusions that can extend over a finite radius. Information of both mechanical and chemical origin are then transduced and act as control factors for the dynamics of cells. Therefore, in the turning operator of the kinetic model we will include the evaluation of mean fields in a neighborhood of the re-orientation position. In particular, we will consider two different control factors $\mathcal{S}$ and $\mathcal{S}^{\prime}$ which are positive quantities defined on $\Omega$ . They respectively bias cell polarization (and therefore orientation) and speed, once cell orientation is set. Hence, we will have a transition probability and a frequency of tumbling which depend on both $\mathcal{S}$ and $\mathcal{S}^{\prime}$ .

The general form of the turning operator which implements a velocity jump processes is

[TABLE]

where the gain term is

[TABLE]

the loss term is

[TABLE]

and ${\bf v}_{p}^{\prime}$ is the pre-turning velocity of the gain term and ${\bf v}_{p}\!\!\!^{\prime}{}^{\prime}$ is the post-turning velocity of the loss term. The so-called turning kernel $T[\mathcal{S},\mathcal{S}^{\prime}]({\bf x},{\bf v}_{p}|{\bf v}_{p}^{\prime})$ is the probability for a cell in ${\bf x}$ of choosing the velocity ${\bf v}_{p}$ after a re-orientation biased by the external fields $\mathcal{S}$ and $\mathcal{S}^{\prime}$ given the pre-turning velocity ${\bf v}_{p}^{\prime}$ . Being a transition probability, it satisfies

[TABLE]

Also the turning frequency $\mu$ may depend on the external signal and/or on its gradient, and on the microscopic orientation, ${\it e.g.}$ through its polarization, typically through ${\bf v}\cdot\nabla\mathcal{S}$ (see for example the work by Othmer and Hillen, (2002)). In fact, it is proved that the speed and the turning rates of individuals depend not only on the magnitude of an external signal $\mathcal{S}$ but also on its temporal and spatial variations as highlighted by Berg, (1983), Fisher et al., (1989), Koshland, (1981), Soll and Wessels, (1998).

As done by Chauviere et al., 2007a and Chauviere et al., 2007b , in the following, we will assume that cells retain no memory of their velocity prior to the re-orientation, ${\it i.e.}$ , $T=T[\mathcal{S},\mathcal{S}^{\prime}]({\bf x},{\bf v}_{p})$ . The independence from the pre-tumbling velocity lies in the fact that the choice of the new velocity is linked to the slow interaction process also related to cell ruffling and sensing which is responsible for the biased re-orientation. However, the assumption might be restrictive in some cases, as it does not include, for instance, the case in which the sensing region depends on the incoming velocity through a polarization-dependent expression of transmembrane proteins. It also excludes persistence effects in which the re-orientation direction depends on the pre-tumbling polarization of the cell. Under such an assumption, the operator (8) simplifies to

[TABLE]

where $\int_{V_{p}}\mu({\bf x},{\bf v}_{p}^{\prime})p(t,{\bf x},{\bf v}_{p}^{\prime})\,d{\bf v}_{p}^{\prime}$ represents the fraction of re-orientating cells per time unit in ${\bf x}$ . If also the frequency $\mu({\bf x})$ does not depend on the microscopic velocity, then we have that

[TABLE]

However, in the following, we will also consider a ’weak’ dependence of the frequency on the incoming velocity and in particular on cell polarization.

We observe that the distribution function nullifying the operator in (13) is simply

[TABLE]

which, in general, is time dependent, but it may also be a stationary equilibrium state.

As a consequence of (11), for any function $T[\mathcal{S},\mathcal{S}^{\prime}]$ the operator $\mathcal{J}[p]$ satisfies the property

[TABLE]

that expresses mass conservation. On the other hand, the evaluation of the first moment of $\mathcal{J}[p]$ with respect to the velocity variable yields

[TABLE]

that indicates no preservation of momentum in general. Similarly, also energy is not preserved. In (16) the vector ${\bf U}$ is the mean velocity of the cells defined by (3), while ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}$ is a macroscopic velocity due to the variations of $\mathcal{S},\,\mathcal{S}^{\prime}$ , defined by

[TABLE]

It is the mean outgoing velocity after a re-orientation biased by $\mathcal{S}$ and $\mathcal{S}^{\prime}$ . Another important quantity is the variance-covariance matrix of $T$ , which is defined as

[TABLE]

Bisi et al., (2008) and Pettersson, (2004) deal with an operator which has a structure similar to (13). They prove that, provided that the probability distribution has initially a finite mass and energy and non-absorbing boundary conditions hold, the function (14), which makes (13) vanish, is a stable asymptotic equilibrium state. It is stationary, which means that $\rho$ is stationary. We shall consider no-flux boundary conditions ( ${\it i.e.}$ (59)) which are non-absorbing boundary conditions. Moreover, ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}$ and $\mathbb{D}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}$ are the mean velocity and variance-covariance tensor of the cells at the equilibrium.

3.1 Structure of the turning probability $T$

From the physiological point of view, the decision process determining the new cell velocity can be split in two parts. First the cell decides where to go and polarizes, forming the so-called cell “head” and “tail”. Then, it starts to move in that direction with a certain speed. These two processes are mainly independent and might be influenced by different chemical and mechanical cues. For instance, environment microstructure, availability of intracellular space and of adhesion sites, expression of integrins, activity of motor proteins, calcium influx, availability of ATP all influence cell speed, while gradients of free and bound chemical factors, of ECM components and ECM stiffness all influence cell polarization. A cell can even be polarized but unable to move. This, for instance, occurs if the cell is treated with latrunculin as shown by Devreotes and Janetopoulos, (2003), or because of overcrowding due to the presence of too many cells in the direction where it would like to move, giving rise to the so-called volume filling effect (see the work by Painter and Hillen, (2002) and references therein), or because of a too dense ECM, giving rise to the so-called physical limit of migration (Wolf et al.,, 2013; Arduino and Preziosi,, 2015; Giverso et al.,, 2017; Scianna and Preziosi,, 2013, 2014; Scianna et al.,, 2013).

As already stated, the choice of the new velocity is a result of a sensing activity of the cell in its neighborhood performed to monitor the quantity of two generally different fields $\mathcal{S}$ and $\mathcal{S}^{\prime}$ , respectively determining cell polarization and speed. This information is then weighted by two functions

[TABLE]

both with compact support that is related to the finite size of the influencing neighborhood related respectively to the fields $\mathcal{S}$ and $\mathcal{S}^{\prime}$ around the position ${\bf x}$ and weighting the information according to the distance from the cell center. They may be

•

Dirac deltas, ${\it e.g.}$ , $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=\delta(\lambda-R_{\scriptscriptstyle\mathcal{S}})$ , if the cells only measure the information perceived on a surface of given radius $R_{\scriptscriptstyle\mathcal{S}}$ ,

•

Heaviside functions, ${\it e.g.}$ , $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=H(R_{\scriptscriptstyle\mathcal{S}}-\lambda)$ if the cells explore the whole volume of the sphere centered in ${\bf x}$ with radius $R_{\scriptscriptstyle\mathcal{S}}$ and weight the information uniformly, or

•

decreasing functions of the distance $\lambda$ from the cell center ${\bf x}$ taking for instance into account that the probability of making longer protrusion decreases with the distance, so the sensing of closer regions is more accurate.

The distance $R_{\scriptscriptstyle\mathcal{S}}$ is, therefore, a measure of the capability of sensing and detecting of the cell. It may be the measure of cell protrusions, or it may be bounded by $\frac{U}{\mu}$ which is the mean linear tract travelled between two consecutive re-orientations. Berg and Purcell, (1977) have shown that the sampling volume in which the signal is significantly detected depends on how rapidly the receptors on the cell’s membrane process the signal. The same considerations are valid for $\gamma_{\scriptscriptstyle\mathcal{S}^{\prime}}$ with a radius $R_{\scriptscriptstyle\mathcal{S}^{\prime}}$ which in general, might be different from $R_{\scriptscriptstyle\mathcal{S}}$ .

The above considerations justify the following splitting

[TABLE]

The quantity $\mathcal{T}_{\lambda}^{\hat{{\bf v}}}[\mathcal{S}]({\bf x})$ is a functional which acts on $\mathcal{S}$ and describes the way the cell measures the quantity $\mathcal{S}$ around ${\bf x}$ along the direction $\hat{{\bf v}}$ and, therefore, the bias intensity in the direction $\hat{{\bf v}}$ . In some cases, we shall write for instance

[TABLE]

where $b$ is a quantity depending on the value of $\mathcal{S}$ in a point at a distance $\lambda$ from ${\bf x}$ along the direction $\hat{\bf v}$ . So, if $b$ is an increasing function of $\mathcal{S}$ and the signal is stronger in the direction $\hat{\bf v}$ , then there will be a higher probability for the cell to move along $\hat{\bf v}$ than along $-\hat{\bf v}$ . We shall also consider other forms for $\mathcal{T}_{\lambda}^{\hat{{\bf v}}}[\mathcal{S}]({\bf x})$ as well.

The quantity $\psi({\bf x},v|\mathcal{S}^{\prime}({\bf x}+\lambda^{\prime}\hat{{\bf v}}))$ is a conditional probability distribution function of the cell speed given the value of $\mathcal{S}^{\prime}$ in a neighborhood of ${\bf x}$ in the direction $\hat{\bf v}$ that is then weighted by $\gamma_{\scriptscriptstyle\mathcal{S}^{\prime}}$ . Again, we shall discuss this further in Section 6. Finally, the factor $c({\bf x})$ is a normalizing factor which is given by

[TABLE]

where $\Gamma_{0}^{\prime}=\displaystyle\int_{\mathbb{R}_{+}}\gamma_{\scriptscriptstyle\mathcal{S}^{\prime}}(\lambda^{\prime})\,d\lambda^{\prime}$ , so that (11) is satisfied.

4 Macroscopic limits

In order to highlight the driving macroscopic phenomenon which may be diffusion or convection, we will discuss the appropriate macroscopic limits of the kinetic equations. A way to do that is to perform, respectively, a parabolic scaling and then a diffusive limit, or a hyperbolic scaling and then a hydrodynamic limit. In these approaches it is assumed that there is a small parameter $\epsilon\ll 1$ that allows to suitably rescale the transport equation. The two scalings are, respectively,

[TABLE]

Following the work by Othmer and Hillen, (2002), we also assume that $T$ can be expanded as

[TABLE]

This means that there are different orders of bias as will be better illustrated in the examples to follow. If we assume that $\mu=\mathcal{O}(1)$ , we may denote $\mathcal{J}^{0}$ and $\mathcal{J}^{1}$ the related operators defined by $T_{0}$ and $T_{1}$ respectively, and we assume that

[TABLE]

Coherently, the velocity and tensor of the second moments of the transition probability may be written as

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

for $i\geq 1$ . In both limits we consider an expansion of the distribution function $p$ in the form

[TABLE]

where the function $p_{0}$ is the equilibrium function given by the solution to the leading order part of the turning operator

[TABLE]

As $\mathcal{J}^{0}[p_{0}]({\bf v}_{p})=\mu(\rho_{0}T[\mathcal{S},\mathcal{S}^{\prime}]_{0}({\bf v}_{p})-p_{0}({\bf v}_{p}))$ , Eq.(27) is satisfied if and only if

[TABLE]

4.1 The diffusive limit

The diffusive limit of velocity jump processes is deeply treated by Hillen and Othmer, (2000) and Othmer and Hillen, (2002). In particular, in (Othmer and Hillen,, 2002) macroscopic models for chemosensitive movement are systematically derived for different orders of $\epsilon$ of the turning probability. In the same spirit as in the works above, because of the conservation of mass, we have that

•

all the mass is in $p_{0}$ , ${\it i.e.}$ ,

[TABLE]

where $\displaystyle{\rho_{i}=\int_{V_{p}}p_{i}d{\bf v}}$ ;

•

[TABLE]

.

In the works by Hillen and Othmer, (2000) and Othmer and Hillen, (2002) the turning probability depends on the pre-tumbling velocity, ${\it i.e.},T=T({\bf v}_{p}|{\bf v}_{p}^{\prime})$ and it is required that

[TABLE]

Thanks to this property and to other regularity assumptions on $T$ Hillen and Othmer, (2000) show that the function identically equal to $1$ is an eigenfunction of the turning operator with eigenvalue 0 that is a simple eigenvalue. In the present work there is no dependence on the pre-tumbling velocity and (31) does not hold true. Furthermore, because of the structure of the turning operator, the density $p_{0}$ making the leading part of the turning operator vanish (28) is an eigenfunction of the velocity jump operator and it depends on the velocity. Moreover, as all the functions nullifying the turning operators are proportional to $T$ up to a multiplicative constant, we have that the eigenvalue 0 is simple and the subspace of $L^{2}(V_{p})$ defined by

[TABLE]

is the kernel of the turning operator and it has dimension equal to one. Therefore, the appropriate scalar product allowing to determine uniquely the solvability condition is the one proposed by Hillen, (2006)

[TABLE]

Whence the equilibrium function $p_{0}$ belongs to the subspace generated by the eigenfunction of $\mathcal{J}^{0}$ with eigenvalue [math], ${\it i.e.}$ , $\langle T[\mathcal{S},\mathcal{S}^{\prime}]_{0}(\boldsymbol{\xi},{\bf v}_{p})\rangle$ . This means that the equilibrium function has mean and the variance-covariance tensor equal to ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S^{\prime}}}^{0}$ and $\mathbb{D}_{\scriptscriptstyle\mathcal{S},\mathcal{S^{\prime}}}^{0}$ respectively. For $i\geq 1$ , the functions $p_{i}\in\langle T[\mathcal{S},\mathcal{S}^{\prime}]_{0}(\boldsymbol{\xi},{\bf v}_{p})\rangle^{\bot}\subseteq L^{2}(V_{p})$ for $a.e.\,\boldsymbol{\xi}$ , where $\langle T[\mathcal{S},\mathcal{S}^{\prime}]_{0}(\boldsymbol{\xi},{\bf v}_{p})\rangle^{\bot}$ is the orthogonal subspace to the subspace generated by $T[\mathcal{S},\mathcal{S^{\prime}}]_{0}(\boldsymbol{\xi},{\bf v}_{p})$ . Thanks to the scalar product (32) we have that $g\in\langle T[\mathcal{S},\mathcal{S}^{\prime}]_{0}(\boldsymbol{\xi},{\bf v}_{p})\rangle^{\bot}$ if and only if $\int_{V_{p}}g({\bf v}_{p})d{\bf v}_{p}=0$ . This justifies the assumptions (29).

The dependence of $p_{0}$ on ${\bf v}_{p}$ has also consequences on the conditions for the isotropy of $T$ and of the diffusion tensor are not the same as the ones considered by Hillen and Othmer, (2000) and Othmer and Hillen, (2002). The case in which $T$ does not depend on the pre-tumbling velocity is treated by Hillen, (2006) where it is assumed that the turning probability $T$ is of order zero in $\epsilon$ . Therefore, we illustrate the diffusive limit procedure in the case in which the turning probability does not depend on the pre-tumbling velocity and it can be expressed as in (22), and we also illustrate the necessary conditions which are needed in order to perform it.

The discriminating factor is whether $T$ is such that

[TABLE]

${\it i.e.}$ , the leading order macroscopic velocity vanishes, or not. In this case we are led to perform a parabolic scaling, ${\it i.e.}$ , $\tau=\epsilon^{2}t$ , ${\bf\boldsymbol{\xi}}=\epsilon{\bf x}$ , $\epsilon\ll 1$ .

The rescaled equation is

[TABLE]

We suppose that $\mu\sim\mathcal{O}(1)$ . From now on, to simplify the notation, we will drop the dependencies on $\tau$ and $\boldsymbol{\xi}$ . Comparing terms of equal order in $\epsilon$ , we then obtain the following system of equations:

In $\epsilon^{0}$ :

[TABLE]

In $\epsilon^{1}$ :

[TABLE]

In $\epsilon^{2}$ :

[TABLE]

As already discussed, Eq. (35) implies (28), that is the equilibrium state of order zero. From Eq. (36), by inverting $\mathcal{J}^{0}$ , we get

[TABLE]

The solvability condition for inverting $\mathcal{J}^{0}$ is that $p_{1}\in\langle T[\mathcal{S},\mathcal{S^{\prime}}]_{0}({\bf v}_{p})\rangle^{\bot}$ . Because of (32) this means that

[TABLE]

Equation (39) is satisfied because

[TABLE]

that is (33), and

[TABLE]

because of (23b).

The evolution equation for $\rho\,(=\rho_{0})$ is obtained from the solvability condition at $\mathcal{O}(\epsilon^{2})$ , ${\it i.e.}$ , the integral on $V_{p}$ of Eq. (37) which is given by

[TABLE]

where we used (30). The term ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S^{\prime}}}^{1}$ is the chemotactic velocity, which is present if there is a bias of order $\epsilon$ as in the work by Othmer and Hillen, (2002). In the following sections, we will show that under appropriate conditions of sensing of the cells, microscopic rules may allow to recover the chemotaxis equation in which the chemotactic velocity is proportional to the gradient of $\mathcal{S}$ and the chemotactic sensitivity depends on the sensitivity function, which captures the cell capability of sensing and detecting the signal.

4.2 Diffusive limit for encounter rates weakly depeding on the incoming velocity

We now generalize the above procedure in the case it is possible to write $\mu({\bf x},{\bf v}_{p})=\mu_{0}({\bf x})+\epsilon\mu_{1}({\bf x},{\bf v}_{p})+\mathcal{O}(\epsilon^{2})$ . In this case the rescaled transport equation reads

[TABLE]

where also $p$ has to be intended as expanded in terms of $\epsilon$ . Comparing terms of equal order in $\epsilon$ , we obtain the same equation as before in $\epsilon^{0}$ , while at the first order

[TABLE]

and at second order

[TABLE]

From equation (44), we get

[TABLE]

where

[TABLE]

In this case the solvability condition for inverting $\mathcal{J}^{0}$ becomes

[TABLE]

that is again satisfied because

[TABLE]

and

[TABLE]

The first two equations vanish as proved in the preceding section, while recalling the definition (47) and (23a), we have that (51) is trivially satisfied.

As in the previous section, the solvability condition at $\mathcal{O}(\epsilon^{2})$ , ${\it i.e.}$ , the integral on $V_{p}$ of (45) gives the evolution equation for $\rho$

[TABLE]

where ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{1}=\displaystyle\int_{V_{p}}T[\mathcal{S},\mathcal{S^{\prime}}]_{1}{\bf v}\,d{\bf v}_{p}$ and

[TABLE]

where the first term vanished because of (33). So, the chemotactic velocity presents a correction to the term ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{1}$ given by ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{1,\mu_{1}}$ . In the particular case in which $\mu_{1}$ is independent of ${\bf v}_{p}$ we obviously recover (42) as ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{1,\mu_{1}}={\bf 0}$ .

4.3 The hyperbolic limit

The discriminating factor on whether a diffusive or a hyperbolic limit can be performed lies in the properties of the chosen turning probability $T$ and specifically on the satisfaction of (33). So, at variance with the previous section, we here assume that $T$ does not satisfy (33). In this case, a hyperbolic scaling, ${\it i.e.}$ $\tau=\epsilon,{\bf\boldsymbol{\xi}}=\epsilon{\bf x}$ , $\epsilon\ll 1$ can be performed. The rescaled transport equation is

[TABLE]

As the equilibrium state is the same as before, we consider a Chapman-Enskog expansion of $p$ in the form

[TABLE]

where $g\in\langle T[\mathcal{S},\mathcal{S^{\prime}}]_{0}({\bf v}_{p})\rangle^{\bot}\subseteq L^{2}(V_{p})$ with the same scalar product as before. Substituting (55) in (54) and integrating on $V_{p}$ the equation at the order $\epsilon^{1}$ , we obtain

[TABLE]

thanks to (29), where (with the dependencies)

[TABLE]

In this case we have then a drift-driven phenomenon.

If we consider a frequency weakly depending on the velocity, the rescaled transport equation modifies into

[TABLE]

and, thanks to (29) and (47), we recover, again, (56) with ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{0}$ given by (57) .

Isotropy and anisotropy

As we said at the beginning of the section, the necessary and sufficient conditions for the isotropy of the diffusion tensor are not the same as the ones considered by Hillen and Othmer, (2000), as the transition probability does not depend on the incoming velocity. This implies that the mean outgoing velocity for a given incoming velocity is always zero and, therefore, the adjoint persistence is zero. In our case, we have that $\mathbb{D}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{0}$ is isotropic if and only if $T_{0}[\mathcal{S},\mathcal{S^{\prime}}]$ is isotropic as a function of $\hat{{\bf v}}$ , ${\it i.e.}$ , if $T_{0}[\mathcal{S},\mathcal{S^{\prime}}]$ is invariant with respect to rotations of $\frac{\pi}{2}$ as a function of $\hat{{\bf v}}$ . This also implies that ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{0}$ is zero and that the directional persistence is zero. Viceversa, if ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{0}$ is not zero, $\mathbb{D}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}^{0}$ is anisotropic.

4.4 Boundary conditions

We shall consider the following biological no-flux condition (Plaza,, 2019)

[TABLE]

being ${\bf n}(\boldsymbol{\xi})$ the outer normal to the boundary $\partial\Omega$ in the point $\boldsymbol{\xi}$ , and we remind that ${\bf v}_{p}=(v,\hat{{\bf v}})$ . This choice is motivated by the fact that in experiments in vitro there is no exchange of biological material between the system and its environment. Equation (59) implies that there is no normal mass flux across the boundary (Lemou and Mieussens,, 2008). At the macroscopic level (59) gives (Plaza,, 2019)

[TABLE]

for the diffusive limit, whilst for the hyperbolic limit the corresponding boundary condition is

[TABLE]

Two important classes of kinetic boundary conditions which satisfy (59) are the regular reflection boundary operators and the non-local (in velocity) boundary operators of diffusive types defined by Palcewski, (1992) (see also the work by Lods, (2005)). Let us denote the boundary operator

[TABLE]

Two main examples of the regular reflection boundary conditions are

•

bounce back reflection condition

[TABLE]

and $p(\tau,\boldsymbol{\xi},v^{\prime},\hat{{\bf v}}^{\prime})=p(\tau,\boldsymbol{\xi},v,\hat{{\bf v}})$ elsewhere.

•

specular reflection boundary condition

[TABLE]

and $p(\tau,\boldsymbol{\xi},v^{\prime},\hat{{\bf v}}^{\prime})=p(\tau,\boldsymbol{\xi},v,\hat{{\bf v}})$ elsewhere.

Diffusive boundary conditions are prescribed in the form

[TABLE]

where $h=h(\boldsymbol{\xi},{\bf v}_{p},{\bf v}_{p}^{*})$ is a non-negative measurable function called the Gaussian equilibrium function satisfying

[TABLE]

A linear combination of a regular reflection with a diffusive boundary operator is called a Maxwell-type boundary operator and reads

[TABLE]

where $M(\boldsymbol{\xi},v,\hat{{\bf v}})$ is the Maxwellian of the wall and $\mathcal{V}(\hat{{\bf v}})=-\hat{{\bf v}}$ for the bounce back reflection condition and $\mathcal{V}(\hat{{\bf v}})=\hat{{\bf v}}-(\hat{{\bf v}}\cdot{\bf n}){\bf n}$ for the specular reflection. We shall consider regular reflection boundary conditions for the one-dimensional simulations and Maxwell type boundary conditions for the two-dimensional simulations (see section 5.3.1).

5 Directional bias

In this section we will specialize the above kinetic models and macrosopic limits to several situations of biological interest, that are listed in Table 1 in order to guide the reader through them. For sake of clarity, for the moment we will forget about the $\mathcal{S}^{\prime}$ field, assuming that once the cell is polarized its speed is determined through a distribution function over the speed $\psi({\bf x},v)$ . The case in which the distribution function depends on a signaling cue $\mathcal{S}^{\prime}$ will be denoted as double bias and will be treated in Section 6. Therefore, $\psi({\bf x},v)$ describes the random unbiased probability for a cell at position ${\bf x}$ of having modulus $v$ after a random re-orientation. In this case we expect a factorized turning probability, as the distributions for the (biased) direction and for the speed are independent.

We will apply the model to several non-local sensing dynamics, comparing the results with other existing models.

5.1 Directional turning probability

Let us introduce the quantities

[TABLE]

and

[TABLE]

where $c_{2}({\bf x})=\dfrac{1}{\displaystyle\int_{\mathbb{S}^{d-1}}\int_{\mathbb{R}_{+}}\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)\mathcal{T}_{\lambda}^{\hat{{\bf v}}}[\mathcal{S}]({\bf x})d\lambda\,d\hat{{\bf v}}}$ is the normalization constant.

In this case, the turning probability can then be factorized as

[TABLE]

We also introduce

[TABLE]

that is the mean speed and

[TABLE]

that is the variance of the distribution function of the speed. The final expression of the operator $\mathcal{J}[p]$ then becomes

[TABLE]

The distribution which makes (68) vanish corresponds to a macroscopic velocity

[TABLE]

and to a diffusion tensor

[TABLE]

The expansion (22) of $T$ reflects now on the expansion of $\bar{B}$

[TABLE]

with

[TABLE]

and

[TABLE]

We remark that if (33) holds true ( ${\it i.e.}$ , if ${\bf U}_{\scriptscriptstyle\mathcal{S}}^{0}$ vanishes)

[TABLE]

The diffusive limit in the case in which the turning frequency does not depend on the microscopic velocity, ${\it i.e.}$ Eq. (42), now reads

[TABLE]

while the diffusive limit in the case in which $\mu_{1}$ depends on the microscopic velocity, ${\it i.e.}$ Eq. (52), now reads

[TABLE]

where

[TABLE]

while the hyperbolic limit is

[TABLE]

5.2 Examples

The function $\mathcal{T}^{\lambda}_{\hat{{\bf v}}}$ , introduced to account for a bias in the choice of the direction, has to be chosen according to the specific taxis process to be considered. In particular, as the cell scouts and detects the signal around itself, the functional $\mathcal{T}^{\lambda}_{\hat{{\bf v}}}$ will depend on $\mathcal{S}$ through the quantity $\mathcal{S}({\bf x}+\lambda\hat{{\bf v}})$ as introduced by Othmer and Hillen, (2002) and later treated by Hillen et al., (2007), with $\lambda=R_{\scriptscriptstyle\mathcal{S}}$ . For us, $\lambda\in[0,R_{\scriptscriptstyle\mathcal{S}}]$ , where $R_{\scriptscriptstyle\mathcal{S}}$ accounts for the size of the sensing neighborhood and the information is weighted through $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)$ . This quantity will be considered in order to define for every ${\bf x}$ the probability of going in direction $\hat{{\bf v}}$ , which is $\bar{B}[\mathcal{S}({\bf x})](\hat{{\bf v}})$ . According to the characteristics of this probability, we will recover a macroscopic model, depending on the value of

[TABLE]

which discriminates between a diffusion driven or a drift driven phenomena, according to whether (33) is satisfied or not.

5.2.1 Local sensing

In the first example we consider a transition rate that depends only on the information given by the control factor at that site, ${\it i.e.},\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=\delta(\lambda-0)$ . In this case, the term $\bar{B}[\mathcal{S}]$ in the transition probability is simply

[TABLE]

where $|\cdot|$ denotes the measure of the set. Then the turning operator reads

[TABLE]

This is what we expect, as the measure of $\mathcal{S}$ in ${\bf x}$ does not affect the choice of the direction. In this model, $\mathcal{S}({\bf x})$ may only affect the frequency of turning. As (33) is satisfied and $B[\mathcal{S}]_{1}=0$ , we perform a diffusive limit. In addition, since

[TABLE]

$\mathbb{D}_{\scriptscriptstyle\mathcal{S}}({\bf\boldsymbol{\xi}})$ is isotropic

[TABLE]

and, if $\mu\sim\mathcal{O}(1)$ , we have that

[TABLE]

Therefore, there is no directional bias which comes from the transition probability, coherently with the fact that no non-local information is taken into account.

5.2.2 Encounter rate dependent taxis

In order to consider situations in which organisms are too small to perform a non-local measurement, like E.Coli, following the ideas by Block et al., (1983) and Othmer and Hillen, (2002), we allow here the encounter rate to weakly depend on a signal $S$ , still working with local models. For sake of simplicity, we assume a time-independent signal, or however changes over a time scale much larger than the free-fly time. In fact, E.Coli can measure a pathwise gradient in time, ${\it i.e.}$ , $\dfrac{D\mathcal{S}}{Dt}=\dfrac{\partial}{\partial t}\mathcal{S}+{\bf v}\cdot\nabla\mathcal{S}$ . Block et al., (1983) show that the movement of E.Coli may be modelled with a Poisson process with turning frequency depending on the temporal gradient of $\mathcal{S}$

[TABLE]

being $\mu_{r}$ the turning rate in absence of the external signal. The turning operator reads

[TABLE]

For instance, one can take the function proposed by Lapidus and Schiller, (1976) ${\it i.e.}$ $f(\mathcal{S})=\dfrac{C_{1}K_{D}}{(K_{D}+\mathcal{S})^{2}}$ where $K_{D}$ is the dissociation constant for the attractant. Hence, in the spirit of the work by Othmer and Hillen, (2002) by a parabolic scaling, we can expand (84) as

[TABLE]

and set

[TABLE]

and

[TABLE]

We then obtain ${\bf U}_{\mathcal{S}}^{1}={\bf 0}$ because $B[\mathcal{S}]_{1}=0$ and, recalling (76),

[TABLE]

and the macroscopic equation reads

[TABLE]

that is a chemotaxis model also in the case of local sensing, when an organism is too small for measuring a gradient, but it is able to measure the temporal gradient of $\mathcal{S}$ along its path.

5.2.3 Non-local sensing average

We now consider the case in which the choice of the new velocity depends on a suitable average of the value of the signal $\mathcal{S}$ perceived through transmembrane receptors by a cell that extending its protrusions can scout its neighborhood.

We shall consider in this case

[TABLE]

that models the fact that in order to decide its new direction, the cell measures $\mathcal{S}$ in ${\bf x}+\lambda\hat{{\bf v}}$ and evaluates an average of a quantity which is linked to this measure ( $b$ ). Therefore, recalling (68) the operator for biased random motion including taxis is

[TABLE]

if the turning frequency does not depend on the microscopic velocity. In order to understand what to expect from the integration of (89) that will be performed in Section 5.3, we can look at the asymptotic limit.

For a general $b$ and $R_{\scriptscriptstyle\mathcal{S}}$ , the macroscopic velocity of order zero is

[TABLE]

As, in this case, (33) is generally not satisfied, we have a non-zero drift term and a hyperbolic limit needs be performed, leading to the integro-differential equation

[TABLE]

In particular, if $b(\mathcal{S})=\mathcal{S}$ , this drift term measures the dominant direction in the extracellular factor in the neighbourhood of the cell as

[TABLE]

However, it is possible to simplify the turning operator (89), if $R_{\scriptscriptstyle\mathcal{S}}$ is much smaller than the characteristic length $l_{\mathcal{S}}$ of variation of $\mathcal{S}$ , ${\it e.g.}$ , $l_{\mathcal{S}}=1/\max\frac{|\nabla\mathcal{S}|}{\mathcal{S}}$ . In fact, we can expand $b$ in a Taylor series

[TABLE]

Therefore, if $b(\mathcal{S}({\bf x}))$ does not vanish, the operator may be approximated by

[TABLE]

where

[TABLE]

For instance, if only the signals at a membrane having a distance $R_{\scriptscriptstyle\mathcal{S}}$ from ${\bf x}$ is taken into account, ${\it i.e.}$ $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=\delta(\lambda-R_{\scriptscriptstyle\mathcal{S}})$ , then $\Lambda=R_{\scriptscriptstyle\mathcal{S}}$ . If instead all the signals between the cell center and its membrane are uniformly mediated, ${\it i.e.}$ $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=H(R_{\scriptscriptstyle\mathcal{S}}-\lambda)$ , then $\Lambda=\dfrac{R_{\scriptscriptstyle\mathcal{S}}}{2}$ .

In this case, the turning operator is

[TABLE]

We remark that the smallness of $R_{\scriptscriptstyle\mathcal{S}}$ and the approximation (93) localizes the non-local integral model (89) into (96). Under these limit assumptions it is possible to perform a parabolic scaling that gives

[TABLE]

Hence, $\bar{B}[\mathcal{S}]_{0}(\boldsymbol{\xi},\hat{{\bf v}})=\dfrac{1}{|\mathbb{S}^{d-1}|}$ and $\bar{B}[\mathcal{S}]_{1}(\boldsymbol{\xi},\hat{{\bf v}})=\dfrac{1}{|\mathbb{S}^{d-1}|}\Lambda\dfrac{b^{\prime}(\mathcal{S}(\boldsymbol{\xi}))}{b(\mathcal{S}(\boldsymbol{\xi}))}\nabla\mathcal{S}(\boldsymbol{\xi})\cdot\hat{{\bf v}}$ . In this case recalling (74) we obtain

[TABLE]

where

[TABLE]

So, the chemotactic sensitivity $\chi$ takes into account of the kinetic response through $\bar{U}$ and of the sensing capability through $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)$ contained in $\Lambda$ . More importantly, different signal-dependent sensitivity models can be obtained according to the choice of $b$ . Viceversa, any relation on the chemotactic sensitivity can be obtained by a proper $b$ given by

[TABLE]

Trivially, if $b$ is independent of $\mathcal{S}$ , one has no chemotaxis. If $b$ is proportional to $\mathcal{S}$ (and $\mathcal{S}\neq 0$ , always), then

[TABLE]

Let us now introduce the parameter

[TABLE]

We have seen that if the sensing radius is smaller then the characteristic length of variation of the chemotactic field, ${\it i.e.}\,\eta\ll 1$ , then (89) can be simplified to (96), whereas this is not possible if the sensing radius is larger then $l_{\scriptscriptstyle\mathcal{S}},{\it i.e.}\,\eta\gg 1$ . This leads to different macroscopic limits (98) and (91) respectively. This different macroscopic behavior and the proper choice of the scaling may be justified thanks to a nondimensionalization argument. We shall rescale the variables as

[TABLE]

where we shall choose

[TABLE]

being $\bar{\mu}$ a reference frequency. The time scale $\sigma_{t}$ can be chosen as a diffusion time $T_{Diff}$ scale or a drift time scale $T_{Drift}$ . In general we may write (Hillen and Othmer,, 2000)

[TABLE]

The regime is diffusive - and we will choose $\sigma_{t}=T_{Diff}$ - if the frequency $\bar{\mu}$ is very large with respect to the reference time scale $\sigma_{t}$ , ${\it i.e.}$ if we can find a small parameter $\epsilon$ such that

[TABLE]

The latter is equivalent to

[TABLE]

which implies the hierarchy

[TABLE]

In the present case this is equivalent to

[TABLE]

On the other hand, the macroscopic regime is hyperbolic, and we choose a faster time scale if

[TABLE]

In this case the appropriate choice will be

[TABLE]

as the hierarchy (102) does not hold anymore. Hence, the following relation holds

[TABLE]

With respect to (100), different chemotactic coefficients can be obtained by different choices of the bias function $b(\mathcal{S})$ : for instance, if $b(\mathcal{S})=k+\mathcal{S}^{n}$ with $k>0$ , then we have a saturating chemotactic sensitivity

[TABLE]

Another example is the receptor binding process. In this case we may consider a saturating dependence on the signal concentration

[TABLE]

yielding

[TABLE]

We notice that in this case, if we choose $n=1$ and $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=\delta(\lambda-R_{\scriptscriptstyle\mathcal{S}})$ , ${\it i.e.}$ , a membrane sensing, we recover the model proposed by Hillen et al., (2007).

Generally speaking, $b$ increasing with $\mathcal{S}$ will give rise to chemoattraction, while $b$ decreasing with $\mathcal{S}$ will lead to a repulsive effect. For instance, if

[TABLE]

we have (103) but with the opposite sign, ${\it i.e.}$ , alignment in the opposite direction with respect to the gradient of $\mathcal{S}$ .

At this stage we did not specify whether $\mathcal{S}$ represents a diffusing or a matrix bound chemical. Even more, one can deal with durotaxis in a completely analogous way. In this case the signal $\mathcal{S}$ is the perceived stiffness of the ECM and the sensing of the mechanical properties of the ECM around the cell will induce motion toward stiffer region of the ECM if $b$ is an increasing function of $\mathcal{S}$ , as we shall see as last application in Section 5.3.

5.2.4 Comparative sensing

In some cases the signals perceived by the cell in differently localized receptors on its membrane is amplified by a polarization dynamics involving PTEN and the phosforillation of PIP2 into PIP3 (Ambrosi et al.,, 2005). This causes the formation of a “head” and a “tail” in the cell that chooses the direction accordingly. In order to mimick this phenomenon, we assume here that the turning rate depends on what is sensed in ${\bf x}+\lambda\hat{{\bf v}}$ and in ${\bf x}-\lambda\hat{{\bf v}}$ , ${\it i.e.}$ ,

[TABLE]

where $b\in(-1,1)$ is a quantity comparing the values $\mathcal{S}({\bf x}-\lambda\hat{{\bf v}})$ and $\mathcal{S}({\bf x}+\lambda\hat{{\bf v}})$ in a way that keeps $\mathcal{T}_{\hat{{\bf v}}}^{\lambda}$ always positive. Hence,

[TABLE]

being $\Lambda_{0}=\displaystyle\int_{\mathbb{R}^{+}}\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)d\lambda$ . In the attractive case, the cell is most likely to go where there is a larger concentration of chemical, and therefore we might take

[TABLE]

On the other hand, in the repulsive case, the cell tends to go where there is a smaller concentration of chemical, and therefore we might take

[TABLE]

As done in the previous section, if we can assume that the size of the neighborhood providing information through signaling is small, one can perform a Taylor expansion of the function $\mathcal{T}_{\lambda}^{\hat{{\bf v}}}$ in $\lambda=0$ and write

[TABLE]

in the attractive case (the repulsive case is similar with a minus sign, or equivalently a negative $\beta$ ). The biased transition probability becomes

[TABLE]

that for a logarithmic dependence of $b$ from $\mathcal{S}$ has the same structure as (97). In this limit the operator for biased random motion including taxis becomes

[TABLE]

Equation (111) satisfies (33) and, therefore, a diffusive limit can be performed. As, $\bar{B}[\mathcal{S}]_{0}(\boldsymbol{\xi},\hat{{\bf v}})=\psi({\bf\boldsymbol{\xi}};v)\dfrac{1}{|\mathbb{S}^{d-1}|}$ and $\bar{B}[\mathcal{S}]_{1}(\boldsymbol{\xi},\hat{{\bf v}})=\dfrac{\beta}{\alpha}\Lambda\dfrac{\nabla\mathcal{S}({\bf x})\cdot\hat{{\bf v}}}{k+\mathcal{S}({\bf x})}$ , we get a drift-diffusion equation (98) with

[TABLE]

Even if in the limit of small sensing radii the comparative sensing is almost the same as the non-local sensing introduced in the previous section (in the sense that the diffusive limit are similar), it allows to add some details. In fact, the ratio of the coefficients $\beta$ and $\alpha$ measures the different weight of the diffusive and the advective (chemotactic) contributions. Further, we may include in $\beta$ a dependence on other substances like $\mathcal{S}$ itself or auto-inducers to model quorum sensing. For example in the case of the cellular slime mold Dictyostelium discoideum, the response to the auto-inducer Netrin-1 may be attractive in high concentrations of cAMP, whilst it may be decreasing in case of low levels of cAMP (Deery and Gomer,, 1999). Painter and Hillen, (2002) consider a response in the form $\beta(w)=1-\frac{w}{w^{*}}$ , being $w^{*}$ a saturation level. In addition, when the sensing radius is not small the kinetic models are different, reflecting the fact that the mechanisms are fundamentally different, with the comparative sensing giving rise to a stronger response, especially when considering that the reception of the signal can be amplified by feedback mechanisms arising from the activation of proper protein cascades inside the cells, giving rise, for instance to the segregation of PIP2 and PIP3 in different parts of the cell.

5.3 Numerical simulations

We simulate the kinetic model with turning operator given by (89) with $b(\mathcal{S})=\mathcal{S}$ both in one and two dimensions. In order to do that, we use the numerical scheme proposed by Filbet and Vauchelet, (2010) in which a kinetic model for chemotaxis is simulated in two-dimensions using a van Leer scheme for the space transport.

5.3.1 Numerical resolution of the kinetic model

We consider a computational domain that will be in the form $[x_{min},x_{max}]\times[0,U]$ in the one dimensional case and $[x_{min},x_{max}]\times[y_{min},y_{max}]\times V_{p}$ in the two dimensional case. The computational domain is discretized with a Cartesian mesh ${\bf X}_{h}\times{\bf V}_{p_{h}}$ , where ${\bf X}_{h}$ and ${\bf V}_{p_{h}}$ are defined by (in two dimensions)

[TABLE]

where $\Delta x=\dfrac{x_{max}-x_{min}}{n_{y}},\Delta y=\dfrac{y_{max}-y_{min}}{n_{y}},\Delta v=\dfrac{U}{n_{v}},\Delta\theta=\dfrac{2\pi}{n_{ang}}$ . Denoting by $p^{n}_{i,j,l,k}$ an approximation of the distribution function $p(t^{n},{\bf x}_{i,j},{\bf v}_{p_{l,k}})$ , where ${\bf v}_{p_{l,k}}=(v_{k},\hat{{\bf v}}_{l})$ . We introduce the first order splitting

[TABLE]

where $h=(\Delta t,\Delta x,\Delta y)$ , ${\bf v}\cdot\nabla_{{\bf x},h}p_{i,j,l,k}^{n}$ is an approximation of the transport operator ${\bf v}\cdot\nabla p$ computed with a Van Leer scheme. It is a high resolution monotone, conservative scheme which is second order if the solution is smooth and first order near the shocks. $T[\mathcal{S},\mathcal{S^{\prime}}]_{i,j,l,k}$ is the discretization of the transition probability $T$ . We observe that as the turning operator preserves mass and the turning probability is known and does not depend on $p$ , the relaxation step may be implicit and we may consider the density at time $n+1$ . In particular the density is computed by using a trapezoidal rule

[TABLE]

Concerning boundary conditions, in the one dimensional case we consider regular reflecting conditions. In one dimension, in a domain $[0,L]$ , the bounce-back and the specular reflection boundary conditions coincide, that is $p(t,x=0,v)=p(t,x=0,-v)$ and $p(t,x=L,-v)=p(t,x=L,v)$ . We do not consider Maxwell type conditions as only the outgoing speed would be affected. In the two dimensional case, the regular reflection is biologically unrealistic, as cells do not bounce back nor they collide with the wall as hard spheres. Therefore, Maxwell type boundary conditions are more realistic, and we shall consider for the Maxwellian to the wall

[TABLE]

being $T$ the asymptotic equilibrium of the system with this class (no-flux) of boundary conditions.

Chemotactic motion

In the first simulation we consider a fixed gaussian distribution of chemoattractant

[TABLE]

and a constant initial condition for the cell population, as shown as in Figure 1 (a). From Figures 1(b), 2, we can observe that cells tend to assume a profile similar to that of the chemoattractant (see also Figure 2). This is not surprising since for a turning rate $\mu$ and distribution function $\psi(v)$ independent of $x$ (and therefore constant $\bar{U}$ and $D$ ) in one dimension the stationary solution of (98) is given by

[TABLE]

and the constant $C$ given by the initial mass, ${\it i.e.}$ ,

[TABLE]

The maximum of the cell stationary state in response to the chemoattractant depends, as shown in Figure 2, both on the size of the sensing neighborhood and on the kind of sensing. A larger sensing radius leads to a stronger motility of cells and, therefore, to a higher peak in the steady state. This is because a bigger sensing radius allows to scout a bigger neighborhood. Furthermore, the sensing of a volume ( $\gamma_{\scriptscriptstyle\mathcal{S}}=H$ ) leads to an average of $\mathcal{S}$ over a bigger region, with respect to $\gamma_{\scriptscriptstyle\mathcal{S}}=\delta$ . This trend is also confirmed by the diffusive limit. In fact, the maximum density depends on the power $m$ defined in (114). Having set $\mu=50$ and chosen a uniform distribution in $[0,2]$ , so that $\bar{U}=1$ and $D=\dfrac{1}{3}$ , then $m=150\Lambda$ with $\Lambda=0.02$ for $\gamma_{\scriptscriptstyle\mathcal{S}}=\delta(\lambda-R_{\scriptscriptstyle\mathcal{S}})$ and $\Lambda=0.01$ for $\gamma_{\scriptscriptstyle\mathcal{S}}=H(\lambda-R_{\scriptscriptstyle\mathcal{S}})$ .

Referring to Table 2, when $R_{\scriptscriptstyle\mathcal{S}}=0.02\ll l_{\mathcal{S}}$ , the collision operator (89) is well approximated by (96) and we can compare the solution between the kinetic model and Equation (98) with $\chi(S)=\Lambda\frac{\bar{U}}{\mathcal{S}}$ . As shown in Figures 2(a),(b), the macroscopic density of the stationary solution of the kinetic model with turning operator given by (89) is very close to the analytic solution (114) of the diffusive limit. In addition, one can observe that a larger sensing radius gives rise to a stronger sensitivity as reflected by a larger $\Lambda$ and therefore a larger $\chi(\mathcal{S})$ in the advection-diffusive equation (98) and a larger $m$ in (114). The comparison between the two solutions remains quite good for $R_{\scriptscriptstyle\mathcal{S}}=0.2$ , which corresponds to $\eta=0.01$ . Instead, as shown in Figure 2 and Table 2, an $R_{\scriptscriptstyle\mathcal{S}}=2$ is not much smaller than $l_{\scriptscriptstyle\mathcal{S}}$ , yielding to larger discrepancies. The trends obtained for $\gamma_{\scriptscriptstyle\mathcal{S}}=\delta$ and $\gamma_{\scriptscriptstyle\mathcal{S}}=H$ are similar. However, the transient time is larger if the sensing function is the Heavyside function.

In the second simulation reported in Figure 3, we use a comparative sensing kernel (112). Here $R_{\scriptscriptstyle\mathcal{S}}=0.02$ and we consider different values of $\alpha$ and $\beta$ . We may observe that if $\alpha$ is much larger than $\beta$ , then the diffusive dynamics is dominant.

Finally, in the set of simulations shown in Figure 4, we start from a macroscopic gaussian distribution for cells, that moves due to the perception of the chemoattractant which is distributed linearly. Cells gradually move to the right until they reach the right boundary where the specular reflection boundary condition, corresponding to a no-flux condition for the macroscopic density, keeps them in the domain.

Cell-cell adhesion

Using the same framework, we may also consider cellular adhesion by considering $\mathcal{S}=\rho$ as a biasing signal. In this case, taking for instance $b(\rho)=\rho$ and $\gamma_{\rho}$ a Dirac delta, the turning operator reads

[TABLE]

As initial condition we take $p(0,x,{\bf v}_{p})=\dfrac{0.2}{U}(1+0.09\sin x)$ moving both to the right and to the left, that corresponds to a small perturbation of the constant initial condition with unitary mass that would stay constant. However, cell-cell adhesion triggers an instability so that the small perturbation amplifies, reaching a peaked distribution in the center of the domain as shown in Figure 5.

In two dimensions we show cell-cell adhesion for a four peaks distribution. Boundary conditions are of Maxwell-type, with $\alpha=0.5$ in (62). We may observe that at the beginning the distributions peak up. Then the sensing radius leads to aggregation of the single peaks in a nearly circular crown structure. Peaks get closer and eventually aggregate.

Durotaxis

We now consider the case in which cells are guided by the rigidity of the extracellular matrix moving towards stiffer regions. We consider a stiffness profile as in Fig. 7 (red line) and $b(\mathcal{S})=\mathcal{S}$ . The sensing function is again a Dirac delta and the sensing radius in the simulation shown in the top row of Figure 7 (red line) is $R_{\scriptscriptstyle\mathcal{S}}=0.02$ . So, the cells within that radius start moving to the stiffer region. Eventually, a stationary state is reached that depends on the sensing radius. In particular, Figure 7 (bottom row) clearly shows that for a small radius $R_{\scriptscriptstyle\mathcal{S}}=0.02$ (which corresponds to $\eta=0.02$ ) the analytic solution of the advection-diffusion equation (115) is similar to the macroscopic density of the solution to the kinetic model. If, instead, $R_{\scriptscriptstyle\mathcal{S}}=2.5$ (which corresponds to $\eta=3$ ), we have that (115) fails in approximating the solution to the kinetic equation, while the solution to the hyperbolic macroscopic equation well approximates the profile of the cell density.

In Figure 8 we present the corresponding simulations in two dimensions. We may observe that the qualitative behavior is the same. Here we show the results with Maxwell type boundary conditions ( $\alpha=0$ in (62)). We checked that the macroscopic stationary state is the same with regular reflection boundary conditions ( $\alpha=1$ in (62)).

6 Double bias

In the previous section we have assumed that the external signal only affects the choice of the turning direction keeping the distribution function $\psi$ determining the speed always the same and in particular unaffected by environmental cues or external chemical signals. Instead, in cell motility one should distinguish polarization mechanisms from the ability to move as they can depend on different factors. In fact, as stated in the introduction, a cell can even be polarized but unable to move. For this reason, in this section we consider a double bias, due to a second factor $\mathcal{S}^{\prime}$ evaluated in a non-local way that influences the speed, once the direction $\hat{{\bf v}}$ is chosen according to the bias ruled by $\mathcal{S}$ . These can be represented not only by different external free and bound chemical factors, or subcellular cues involved in cell motility, but also by cell and ECM density.

In order to do that, we allow the distribution function $\psi({\bf x},v)$ to depend on the level of the signal $\mathcal{S}^{\prime}$ . In order to stress this dependence, we use the following notations $\psi=\psi({\bf x},v|\mathcal{S}^{\prime}({\bf x}+\lambda^{\prime}\hat{{\bf v}}))$ for the distribution function, $\bar{v}({\bf x}|\mathcal{S}^{\prime}({\bf x}+\lambda^{\prime}\hat{\bf v}))$ for its mean, and $D({\bf x}|\mathcal{S}^{\prime}({\bf x}+\lambda^{\prime}\hat{\bf v}))$ for its variance, where $|\mathcal{S}^{\prime}$ reads as ”given $\mathcal{S}^{\prime}$ in a point ${\bf x}+\lambda^{\prime}\hat{{\bf v}}$ in the neighborhood of the cell”. Then, the cell polarized along $\hat{\bf v}$ determines its speed averaging the signal over the sensing radius $R_{\scriptscriptstyle\mathcal{S}^{\prime}}$ through a weight function $\gamma_{\scriptscriptstyle\mathcal{S}^{\prime}}$ .

In the same spirit as Eq. (63), introducing also

[TABLE]

the turning probability in (12) or (13) reads

[TABLE]

The normalization term $c$ can be factorized as

[TABLE]

with

[TABLE]

and

[TABLE]

In this way,

[TABLE]

and

[TABLE]

are both distribution functions and the turning probability reads $T=\bar{B}\bar{\Psi}$ .

Going back to the general case, the turning operator reduces to

[TABLE]

and the distribution function which nullifies the turning operator is

[TABLE]

Now, ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}$ in Eq. (17) reads

[TABLE]

where

[TABLE]

and the diffusion tensor is

[TABLE]

Let us now assume that, like $B$ , $\Psi[\mathcal{S}^{\prime}]$ may be written, up to re-scaling, as

[TABLE]

The means and variances of $\Psi[\mathcal{S}^{\prime}]_{0}$ and $\Psi[\mathcal{S}^{\prime}]_{1}$ will be, respectively, denoted by $\bar{U}_{\scriptscriptstyle\mathcal{S}^{\prime}}^{0}$ , $\bar{U}_{\scriptscriptstyle\mathcal{S}^{\prime}}^{1}$ , $\bar{D}_{\scriptscriptstyle\mathcal{S}^{\prime}}^{0}$ , $\bar{D}_{\scriptscriptstyle\mathcal{S}^{\prime}}^{1}$ .

Up to re-scaling, the re-scaled turning probability becomes $T[\mathcal{S},\mathcal{S}^{\prime}]({\bf v}_{p})=T[\mathcal{S},\mathcal{S}^{\prime}]_{0}({\bf v}_{p})+\epsilon T[\mathcal{S},\mathcal{S}^{\prime}]_{1}({\bf v}_{p})$ , where, now

[TABLE]

and

[TABLE]

In this case, Eqs.(23a) and (23b) are satisfied if

[TABLE]

Therefore, the macroscopic velocity of order 0 is

[TABLE]

while the macroscopic velocity of order 1 is

[TABLE]

and the equilibrium diffusion tensor

[TABLE]

The diffusion tensor is in general anisotropic, as $\bar{D}^{0}$ may be not isotropic in $\hat{{\bf v}}$ .

If (33) is satisfied, the equilibrium tensor is

[TABLE]

where

[TABLE]

In this case, we get as a diffusive limit Eq. (42), or (52) if the frequency depends weakly on the polarization. Otherwise, if (33) is not satified, we get (56) as a hyperbolic limit.

6.1 Examples

First of all, we observe that if $\gamma_{\scriptscriptstyle\mathcal{S}^{\prime}}(\lambda)=\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=\delta(\lambda-0)$ are Dirac deltas, then the local model with $T[\mathcal{S},\mathcal{S}^{\prime}]=\frac{1}{|\mathbb{S}^{d-1}|}\psi({\bf x},v)$ is recovered. Hence ${\bf U}_{\scriptscriptstyle\mathcal{S},\mathcal{S}^{\prime}}$ is zero and $\bar{U}_{\scriptscriptstyle\mathcal{S}^{\prime}}({\bf x}|\hat{\bf v})$ and $\bar{D}_{\scriptscriptstyle\mathcal{S}^{\prime}}({\bf x}|\hat{\bf v})$ are the mean speed and variance of $\psi({\bf x},v)$ as defined in Section 5.1.

In this section we shall at first consider $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=\delta(\lambda-0)$ , ${\it i.e.}$ , when there is no bias in the decision of the direction of motion. We shall do this in order to analyze the influence of the signal $\mathcal{S}^{\prime}$ . In particular, we will see that if, for instance, $\mathcal{S^{\prime}}$ is the ECM density then in the limit we recover a model describing the steryc hindrance of motion. Finally, we will analyse models where both $\mathcal{S}$ and $\mathcal{S}^{\prime}$ have a non-local influence in determining cell motion.

6.1.1 Random polarization

In order to understand the meaning of the transport model we here focus on the particular case in which there is no bias in the decision of the direction of motion. This does not mean that the situation is isotropic, because the speed can have different distribution functions in the randomly chosen direction of motion because the distribution of $\mathcal{S}^{\prime}$ sensed ahead may be different.

In this case the turning operator (117) simplifies to

[TABLE]

where $\bar{\Psi}$ is defined in (119), the macroscopic velocity of the transition probability is

[TABLE]

and the variance-covariance matrix is

[TABLE]

Up to re-scaling, we obviously have that

[TABLE]

and we can consider the following cases:

$\bullet$

If $\bar{\Psi}[\mathcal{S}^{\prime}]_{0}(\boldsymbol{\xi},v|\hat{{\bf v}})$ is even as a function of $\hat{{\bf v}}$ , ${\it i.e.}$ if $\bar{U}$ is even as a function of $\hat{{\bf v}}$ , we have that ${\bf U}_{\mathcal{S}^{\prime}}^{0}(\boldsymbol{\xi})=\boldsymbol{0}$ . We can therefore perform a diffusive limit to get

[TABLE]

if $\bar{\Psi}[\mathcal{S}^{\prime}]_{1}$ is zero and

[TABLE]

otherwise. The tensor $\mathbb{D}_{\scriptscriptstyle\mathcal{S}^{\prime}}^{0}$ may be anisotropic, as $\bar{U}$ and $\bar{D}$ may be anisotropic as functions of $\hat{{\bf v}}$ .

$\bullet$

If $\bar{\Psi}[\mathcal{S}^{\prime}]_{0}(\boldsymbol{\xi},v|\hat{{\bf v}})$ is not even as a function of $\hat{{\bf v}}$ , ${\it i.e.}$ if $\bar{U}$ is not even as a function of $\hat{{\bf v}}$ , we can perform a hyperbolic limit to get

[TABLE]

Therefore, we see that, even if there is no directional sensing, as the speed sensing is conditioned by the direction, there may be a drift driven phenomenon or an anisotropic diffusive one.

6.1.2 ECM steryc hindrance

We shall now focus on the effect of the density $M({\bf x})$ of extracellular matrix. Following the experimental results by Palecek et al., (1997), Paul et al., (1993), Wolf et al., (2013), we take into account that there is an optimal ECM density or pore cross section for cell migration. In fact, for lower densities cell motion results hampered for the lack of adhesion sites that the cell can use to exert traction. On the other hand, for higher densities cell motion is hampered as well by excessive adhesion and lack of space among the fibers. Actually, it is known that there is a maximum value of ECM density, or better a minimum value of pore area cross section, above which the fiber network is too tight to be penetrated (the so-called physical limit of migration (Wolf et al.,, 2013; Scianna and Preziosi,, 2013; Giverso et al.,, 2017). It is more controversial whether there is a minimum value of ECM density denoted in the following by $M_{0}$ below which cells are unable to migrate because of the absence of ECM fibers on which to crawl. The existence of these physical limits of migration will be examined in more detail in a future work.

One can then take such experimentally determined functions as an indication of the dependence of the mean $\bar{v}$ of the distribution function $\psi$ from the density $M$ of extracellular matrix. In the following simulation, to mimick such a behaviour we assume that

[TABLE]

and choose $\psi$ to be

[TABLE]

We consider an initial condition which is uniformly distributed in the velocity space and with macroscopic density a gaussian centered in $x=1.5$ and the heterogeneous distribution of ECM density ranging in $[M_{0},M_{1}]=[0.1,1]$ shown in Fig. 9 (a) and (b)(green line). We take here $\gamma_{\scriptscriptstyle M}(\lambda^{\prime})=H_{M}(R_{\scriptscriptstyle M}-\lambda^{\prime})$ as sensing function, meaning that the speed is determined by uniformly averaging the density of ECM in the direction $\hat{\bf v}$ looking ahead up to a distance $R_{\scriptscriptstyle M}$ that in the simulations is $0.04$ .

If a matrix density distributed as in Fig. 9 (a) and (b) is present, cells change their speed accordingly. In particular, if $M_{max}<M_{0}=0.1$ (as in the simulations reported in Fig. 9 (a)) cell mean speed is on the decreasing branch of (132), i.e., always smaller than $v_{max}\frac{M_{0}}{M_{max}}\exp\left[-(M_{0}-M_{max})/M_{max}\right]$ . Hence, they are strongly slowed down by the ECM. On the other hand, if $M_{max}>M_{1}=1$ (as in the simulations reported in Fig. 9 (b)) cell mean speed is on the increasing branch of (132), i.e., always larger than the value given above. So, they progressively invade the whole spatial domain and crawl across of the ECM.

In two dimensions, if a matrix density is distributed as in Fig. 10 (a), we observe that cells tend to go where the value of the matrix density is $0.5*M_{max}$ , ${\it i.e.}$ where the speed is the largest. Upon reaching this zone with maximum mean speed (see Fig.10(e)) diffusion looks anisotropic because cells tend to follow the region with optimal ECM density staying away from the regions where matrix density gets closer to $M_{max}$ where motility is hampered.

6.1.3 ECM steryc hindrance and non-local chemotaxis

We shall now consider a real double bias, which means that there is a non-local sensing of the field $\mathcal{S}$ for the directional bias, ${\it i.e.}$ , $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=\delta(\lambda-R_{\scriptscriptstyle\mathcal{S}})$ or $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=H(R_{\scriptscriptstyle\mathcal{S}}-\lambda)$ . In addition we take $\mathcal{T}_{\lambda}^{\hat{{\bf v}}}[\mathcal{S}]({\bf x})=b\big{(}\mathcal{S}({\bf x}+\lambda\hat{{\bf v}})\big{)}$ . Then the turning operator reads

[TABLE]

As (33) is hardly satisfied, the macroscopic limit is advective

[TABLE]

On the other hand, if $R_{\scriptscriptstyle\mathcal{S}}$ can be considered small, the turning operator reads

[TABLE]

In this case

[TABLE]

and

[TABLE]

Condition (33) is satisfied if $\Psi[\mathcal{S}^{\prime}]_{0}$ is even as a function of $\hat{{\bf v}}$ for all ${\bf x}$ and the diffusive limit is (42) with advective velocity

[TABLE]

being

[TABLE]

It may be anisotropic as, in general, the mean speed is anisotropic.

If $\Psi[\mathcal{S}^{\prime}]_{0}$ is not even as a function of $\hat{{\bf v}}$ , we must perform a hyperbolic limit that gives

[TABLE]

This shows that even if the directional sensing gives rise to an even distribution function of order zero in $\hat{{\bf v}}$ and, then, to a diffusive behavior, the speed sensing may change this tendency and we may have a drift driven dynamics.

As an application we shall now consider the case of cell migration in the extra-cellular matrix under the action of a chemoattractant. We then integrate the transport equation (1) with the operator (133) with $b(\mathcal{S})=\mathcal{S}$ , the mean $\bar{v}$ given by (132). We take here the same sensing function as before, ${\it i.e.},$ $\gamma_{\scriptscriptstyle M}(\lambda^{\prime})=H_{M}(R_{\scriptscriptstyle M}-\lambda^{\prime})$ . Also $\gamma_{\scriptscriptstyle\mathcal{S}}(\lambda)=H(R_{\scriptscriptstyle\mathcal{S}}-\lambda)$ .

In Figure 11, the initial gaussian distribution of cells, then, starts travelling to the right under the action of the linearly distibuted chemoattractant. However, at variance with what happend in Figure 4, when cells encounter the denser area of ECM they slow down. If $M_{max}$ is large ( $M_{max}=1/0.3$ ), the mean speed is large and the cells are driven by the chemoattractant, even if the radius of the direction sensing is small ( $R_{\scriptscriptstyle\mathcal{S}}=0.04$ ). So, eventually the cells succeed in going through the ECM as shown in the third row of Figure 8. The evolution of the distribution function is also given in the bottom row. If $M_{max}$ is small ( $M_{max}=1/10$ ), even if the influence of the chemoattractant is very strong ( $R_{\scriptscriptstyle\mathcal{S}}=1$ ), cells are strongly hampered by the ECM as shown in the second row of Figure 8.

7 Conclusion

The kinetic model developed in this paper is based on several observations, all related to a particular feature of the model. Specifically,

in order to move cells react to external stimuli by polarizing and deciding their direction of motion, forming a “head” and “tail”. Then they try to move in that direction forming adhesion sites and polimerizing actin at the head and removing adhesion sites at the tail. It is therefore convenient to write the kinetic model with the orientation and speed, separately, as independent variables; 2. 2.

Both orientation and speed depend on cell sensing of their environment over a finite radius related, for instance, to the length of their protrusion. This gives the kinetic model a non-local character; 3. 3.

Orientation and speed can be related to different chemical and mechanical cues, with the former influencing mainly orientation and the latter mainly speed (though this is not always the case). This induces the inclusion of different cues determining the changes of orientation and speed in the turning operator; 4. 4.

While chemical cues are measured by the cells extending their protrusions, many mechanical cues are related to the fact that the nucleus, which represents the stiffest organelle of the cell, has difficulty in passing through narrow pores in the ECM or through densely packed cells. So, the sensing radii involved in the identification of the orientation and of the speed might be different with the former being larger than the latter.

A particular attention was devoted to the identification of the proper macroscopic limit according to the properties of the turning operator. However, all simulations integrate the kinetic model. In some cases, the solution is compared with the solution of the related macroscopic equations. The model was in fact applied to several sensing strategies, including when the signal is perceived locally, averaged over a region, or perceived at the cell “head” and “tail” and then compared. Several applications were also considered as chemotaxis (or equivalently haptotaxis), durotaxis, cell-cell adhesion. The effect of the presence of an heterogeneous ECM was also considered in absence of any cue affecting cell orientation, then yielding a random polarization, and in presence of a chemical cue.

We have left to a coming paper the case in which the density of cell or of ECM is so large that they represent physical limits for migration. This introduces a dependence of the sensing radius determining the speed that depends on the density of cells or of ECM. In fact, for instance, if the characteristic pore size of the ECM is very small, then cell can protrude its cytoskeleton across the dense ECM and scout for chemical signals and polarize, but the cell cannot penetrate the ECM because the nucleus is stuck. This holds true even in the case in which the layer of dense ECM is very thin, like in basal membrane or in cell lining, e.g., endothelial or mesothelial lining. The cell would have no difficulties in moving after passing through the membrane but cannot cross it. One should then address the question of deducing a kinetic model capable of dealing with such situations that are of great interest because of their relationship with the spread of metastasis.

Aknowledgements

This work was partially supported by Istituto Nazionale di Alta Matematica, Ministry of Education, Universities and Research, through the “MIUR grant Dipartimenti di Eccellenza 2018-2022” and Compagnia San Paolo that finances NL PhD scholarship.

Bibliography61

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Adler, (1966) Adler, J. (1966). Chemotaxis in bacteria. Science , 153(3737):708–116.
2Alt, (1980) Alt, W. (1980). Biased random walk models for chemotaxis and related diffusion approximations. Journal of Mathematical Biology , 9(2):147–177.
3Ambrosi et al., (2005) Ambrosi, D., Gamba, A., and Serini, G. (2005). Cell directional persistence and chemotaxis in vascular morphogenesis. Bulletin of Mathematical Biology , 67(1):195–195.
4Arduino and Preziosi, (2015) Arduino, A. and Preziosi, L. (2015). A multiphase model of tumour segregation in situ by a heterogeneous extracellular matrix. International Journal of Non-Linear Mechanics , 75:22–30.
5Armstrong et al., (2006) Armstrong, N. J., Painter, K. J., and Sherratt, J. A. (2006). A continuum approach to modelling cell-cell adhesion. Journal of Theoretical Biology , 243(1):98–113.
6Berg and Purcell, (1977) Berg, H. and Purcell, E. (1977). Physics of chemoreception. Biophysical Journal , 20(2):193 – 219.
7Berg, (1983) Berg, H. C. (1983). Random Walks in Biology . Princeton University Press, rev - revised edition.
8Bisi et al., (2008) Bisi, M., Carrillo, J. A., and Lods, B. (2008). Equilibrium solution to the inelastic boltzmann equation driven by a particle bath. Journal of Statistical Physics , 133(5):841–870.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Kinetic models with non-local sensing determining cell polarization and speed according to independent cues

Abstract

1 Introduction

2 Preliminaries

3 The turning operator

3.1 Structure of the turning probability TTT

4 Macroscopic limits

4.1 The diffusive limit

4.2 Diffusive limit for encounter rates weakly depeding on the incoming velocity

4.3 The hyperbolic limit

Isotropy and anisotropy

4.4 Boundary conditions

5 Directional bias

5.1 Directional turning probability

5.2 Examples

5.2.1 Local sensing

5.2.2 Encounter rate dependent taxis

5.2.3 Non-local sensing average

5.2.4 Comparative sensing

5.3 Numerical simulations

5.3.1 Numerical resolution of the kinetic model

Chemotactic motion

Cell-cell adhesion

Durotaxis

6 Double bias

6.1 Examples

6.1.1 Random polarization

6.1.2 ECM steryc hindrance

6.1.3 ECM steryc hindrance and non-local chemotaxis

7 Conclusion

Aknowledgements

3.1 Structure of the turning probability $T$