Nonlocal and local models for taxis in cell migration: a rigorous limit procedure

A rigorous limit procedure is presented which links nonlocal models involving adhesion or nonlocal chemotaxis to their local counterparts featuring haptotaxis and classical chemotaxis, respectively. It relies on a novel reformulation of the involved nonlocalities in terms of integral operators applied directly to the gradients of signal-dependent quantities. The proposed approach handles both model types in a unified way and extends the previous mathematical framework to settings that allow for general solution-dependent coefficient functions. The previous forms of nonlocal operators are compared with the new ones introduced in this paper and the advantages of the latter are highlighted by concrete examples. Numerical simulations in 1D provide an illustration of some of the theoretical findings.


Introduction
Macroscopic equations and systems describing the evolution of populations in response to soluble and insoluble environmental cues have been intensively studied and the palette of such reaction-diffusion-taxis models is continuously expanding. Models of such form are motivated by problems arising in various contexts, a large part related to cell migration and proliferation connected to tumor invasion, embryonal ME and CS were supported in part by the research initiative Mathematics Applied to Real World Challenges (MathApp) of the TU Kaiserslautern. CS also acknowledges funding by the Federal Ministry of Education and Research (BMBF) in the project GlioMaTh 05M2016.
Extended author information available on the last page of the article development, wound healing, biofilm formation, insect behavior in response to chemical cues, etc. We refer, e.g. to Bellomo et al. (2015) for a recent review also containing some deduction methods for taxis equations based on kinetic transport equations.
Apart from such purely local PDE systems with taxis, several spatially nonlocal models have been introduced over the last two decades and are attracting ever increasing interest. They involve integro-differential operators in one or several terms of the featured reaction-diffusion-drift equations. Their aim is to characterize interactions between individuals or signal perception happening not only at a specific location, but over a whole set (usually a ball) containing (centered at) that location. In the context of cell populations, for instance, this seems to be a more realistic modeling assumption, as cells are able to extend various protrusions (such as lamellipodia, filopodia, cytonemes, etc.) into their surroundings, which can reach across long distances compared against cell size, see González-Méndez et al. (2017) and Sáenz-de Santa-María et al. (2017) and references therein. Moreover, the cells are able to relay signals they perceive and thus transmit them to cells with which they are not in direct contact, thereby influencing their motility, see e.g., Eom and Parichy (2017) and Garcia and Parent (2008). Cellcell and cell-tissue adhesion are essential for mutual communication, homeostasis, migration, proliferation, sorting, and many other biological processes. A large variety of models for adhesive behavior at the cellular level have been developed to account for the dynamics of focal contacts, e.g. Bell (1978), Bell et al. (1984) and Ward and Hammer (1993) and to assess their influence on cytoskeleton restructuring and cell migration, e.g. DiMilla et al. (1991), Dickinson and Tranquillo (1993), Kuusela and Alt (2008), Uatay (2019). Continuous, spatially nonlocal models involving adhesion were introduced more recently ( Armstrong et al. 2006) and are attracting increasing interest from the modeling (Bitsouni et al. 2018;Buttenschön et al. 2018;Carrillo et al. 2019;Domschke et al. 2014;Gerisch and Chaplain 2008;Gerisch and Painter 2010;Murakawa and Togashi 2015;Painter et al. 2010), analytical (Chaplain et al. 2011;Dyson et al. 2013Dyson et al. , 2010Sherratt et al. 2009;Hillen et al. 2018), and numerical (Gerisch 2010) viewpoints. Yet more recent models (Domschke et al. 2017;Engwer et al. 2017) also take into account subcellular level dynamics, thus involving further nonlocalities (besides adhesion), with respect to some structure variable referring to individual cell state. Thereby, multiscale mathematical settings are obtained, which lead to challenging problems for analysis and numerics. Another essential aspect of cell migration is the directional bias in response to a diffusing signal, commonly termed chemotaxis. A model of cell migration with finite sensing radius, thus featuring nonlocal chemotaxis has been introduced in Othmer and Hillen (2002) and readdressed in Hillen et al. (2007) from the perspective of well-posedness, long time behaviour, and patterning. We also refer to Loy and Preziosi (2019) for further spatially nonlocal models and their formal deduction.
For adhesion and nonlocal chemotaxis models, a gradient of some nondiffusing or diffusing signal is replaced by a nonlocal integral term. Here we are only interested in this type of model, and refer to Chen et al. (2020); Eftimie (2018), Kavallaris and Suzuki (2018) for reviews on settings involving other types of nonlocality. Specifically, following Armstrong et al. (2006), Gerisch and Chaplain (2008), Hillen et al. (2007) and Othmer and Hillen (2002), we consider the subsequent systems, whose precise mathematical formulations will be specified further below: 1. a prototypical nonlocal model for adhesion ∂ t c r = ∇ · (D c (c r , v r )∇c r − c r χ(c r , v r )A r (g(c r , v r ))) + f c (c r , v r ), (1.1a) is referred to as the adhesion velocity, and the function F r describes how the magnitude of the interaction force depends on the interaction range |ξ | within the sensing radius r . We require this function to satisfy Assumptions 1.1 (Assumptions on F r ) (i) (r , ρ) → F r (ρ) is continuous and positive in [0, r 0 ] 2 for some r 0 > 0; (ii) F 0 (0) = n + 1. 1 The quantity is often referred to as the total adhesion flux, possibly scaled by some constant involving the typical cell size or the sensing radius, see e.g., Armstrong et al. (2006) and Buttenschön et al. (2018). Here we also include a coefficient χ(c r , v r ) that depends on cell and tissue (extracellular matrix, ECM) densities, which can be seen as characterizing the sensitivity of cells towards their neighbours and the surrounding tissue. It will, moreover, help provide in a rather general framework a unified presentation of this and the subsequent local and nonlocal model classes for adhesion, haptotactic, and chemotactic behavior of moving cells. System (1.1) is a simplification of the integro-differential system (4) in Gerisch and Chaplain (2008). The main difference between the two settings is that in our case we ignore the so-called matrix-degrading enzymes (MDEs). Instead, we assume the cells to degrade the tissue directly: this fairly standard simplification (e.g., Painter et al. 2010) effectively assumes that proteolytic enzymes remain localised to the cells, and helps simplify the analysis. On the other hand, (1.1) can also be viewed as a nonlocal version of the haptotaxis model with nonlinear diffusion: (c, v)∇c − cχ(c, v)∇g (c, v)) + f c (c, v), (1.3b) 2. a prototypical nonlocal chemotaxis-growth model with the nonlocal gradient System (1.4) can be seen as a nonlocal version of the chemotaxis-growth model (1.5b) where χ (c, v) is the chemotactic sensitivity function. As mentioned above, in order to have a unified description of our systems (1.3) and (1.5) and of their respective nonlocal counterparts (1.1) and (1.4), we later introduce a more general version of the nonlocal chemotaxis flux, similar to the above adhesion velocity A r .
Here and below B r and S r denote the open r -ball and the r -sphere in R n , both centred at the origin, and are the usual mean values of a function u over B r and S r , respectively. The nonlocal systems (1.3) and (1.5) are stated for Unless the spatial domain Ω is the whole R n , suitable boundary conditions are required. In the latter case, usually periodicity is assumed, which is not biologically realistic in general. Still, this offers the easiest way to properly define the output of the nonlocal operator in the boundary layer where the sensing region is not fully contained in Ω. Very recently various other boundary conditions have been derived and compared in the context of a single equation modeling cell-cell adhesion in 1D (Hillen and Buttenschön 2020). Few previous works focus on solvability for models with nonlocality in a taxis term. Some of them deal with single equations that only involve cell-cell adhesion (Dyson et al. 2010(Dyson et al. , 2013Hillen and Buttenschön 2020), others study nonlocal systems of the sort considered here for two  or more components (Engwer et al. 2017). The global solvability and boundedness study in Hillen et al. (2018) is obtained for the case of a nonlocal operator with integration over a set of sampling directions being an open, not necessarily strict subset of R n . The systems studied there include settings with a third equation for the dynamics of diffusing MDEs. Conditions which secure uniform boundedness of solutions to such cell-cell and cell-tissue adhesion models in 1D were elaborated in Sherratt et al. (2009).
Some heuristic analysis via local Taylor expansions was performed in Gerisch and Chaplain (2008) and Hillen (2007) showing that as r → 0 the outputs A r u and∇ r u, respectively, converge pointwise to ∇u for a fixed and sufficiently smooth u. In Hillen et al. (2007) it was observed that it would be interesting to study rigorously the limiting behaviour of solutions of the nonlocal problems involving∇ r u. The authors ask in which sense, if at all, do these solutions converge to solutions of the corresponding local problem as r → 0. Numerical results appeared to confirm that, in certain cases, the answer is positive. Still, to the best of our knowledge, no rigorous analytical study of this issue has as yet been performed. Clearly, any approach based on representations using Taylor polynomials requires a rather high order regularity of solution components and a suitable control on the approximation errors, and that uniformly in r . This is difficult or even impossible to obtain in most cases, particularly when dealing with weak solutions. In this work we propose a different approach based on the representation of the input u in terms of an integral of ∇u over line segments. This leads to a new description of the nonlocal operators A r and∇ r in terms of nonlocal operators applied to gradients (see Sect. 3 below). Moreover, it turns out that redefining their outputs inside the vanishing boundary layer in a suitable way allows one to perform a rigorous proof of convergence: Under suitable assumptions on the system coefficients and other parameters, appropriately defined sequences of solutions to nonlocal problems involving the mentioned modified nonlocal operators converge for r → 0 to those of the corresponding local models (1.3) and (1.5), respectively. Our convergence proof is based on estimates on c r and v r which are uniform in r and on a compactness argument. The two models (1.1) and (1.4) are chosen as illustrations, however our idea can be further applied to other integro-differential systems with similar properties.
The rest of the paper is organised as follows. Section 2 introduces some basic notations to be used throughout this paper. In Sect. 3 we introduce the aforementioned adaptations of the nonlocal operators A r and∇ r and study their limiting properties as r becomes infinitesimally small. This turns out to be useful for our convergence proof later. We also establish in Sect. 4 the well-posedness for a certain class of equations including such operators. In the subsequent Sect. 5 we introduce a couple of nonlocal models that involve the previously considered averaging operators, prove the global existence of solutions of the respective systems, and investigate their limit behaviour as r → 0. Section 6 provides some numerical simulations comparing various nonlocal and local models considered in this work in the 1D case. Finally, Sect. 7 contains a discussion of the results and a short outlook on open issues.

Basic notations and function spaces
We denote the Lebesgue measure of a set A by |A|. Let Ω ⊂ R n be a bounded domain with smooth enough boundary.
For a function w : Ω → R n we assume, by convention, that For r > 0 we introduce the following subdomain of Ω Partial derivatives, in both classical and distributional sense, with respect to variables t and x i , will be denoted respectively by ∂ t and ∂ x i . Further, ∇, ∇· and Δ stand for the spatial gradient, divergence and Laplace operators, respectively. ∂ ν is the derivative with respect to the outward unit normal of ∂Ω.
We assume the reader to be familiar with the definitions and the usual properties of such spaces as: the standard Lebesgue and Sobolev spaces, spaces of functions with values in these spaces, and with anisotropic Sobolev spaces. In particular, we denote by C w ([0, T ]; L 2 (Ω)) the space of functions u : [0, T ] → L 2 (Ω) which are continuous w.r.t. the weak topology of L 2 (Ω).
Throughout the paper ·, · X * ,X denotes a duality paring between a space X and its dual X * .
Finally, we make the following useful convention: For all indices i, the quantity C i denotes a positive constant or, alternatively, a positive function of its arguments. Moreover, unless explicitly stated, these constants do not depend upon r .

Operators A r and∇ r and averages of ∇
In this section we study the applications of the non-local operators A r and∇ r to fixed, i.e. independent of r , functions u. Our focus is on the limiting behaviour as r → 0. Formal Taylor expansions performed in Gerisch and Chaplain (2008), Hillen et al. (2007) anticipate that the limit is the gradient operator in both cases. This we prove here rigorously under rather mild regularity assumptions on u. To be more precise, we replace A r and∇ r by certain integral operators T r and S r (see (3.2) and (3.7) below) applied to ∇u and show that these operators are pointwise approximations of the identity operator in the L p spaces.
We start with the operator A r . For r ∈ (0, r 0 ], u ∈ C 1 (Ω), and x ∈ Ω r we compute that Formula (3.1) extends to arbitrary u ∈ W 1,1 (Ω) by means of a density argument. Motivated by (3.1) we introduce the averaging operator In Sect. 3.1 we check that T r w(x) is well-defined for all w ∈ (L 1 (Ω)) n and a.a.
x ∈ Ω. In this notation, for all r ∈ (0, r 0 ] and u ∈ W 1,1 (Ω) identity (3.1) takes the form In the limiting case r = 0 we have that In the final step we used Assumptions 1.1(ii) which says that F 0 (0) = n + 1 (this explains our choice) and the trivial identity -B 1 |y| dy = n n + 1 .
(3.4) Thus, we have just proved the following lemma: Lemma 3.1 (Adhesion velocity vs. T r ) Let u ∈ W 1,1 (Ω). Then it holds that In a very similar manner one can establish a representation for∇ r . For this purpose we define the averaging operator (3.7) The corresponding result then reads: Lemma 3.2 (Non-local gradient vs. S r ) Let u ∈ W 1,1 (Ω). Then it holds that (3.9) The proof of Lemma 3.2 is very similar to that of Lemma 3.1 and we omit it here. Next, we observe that identity (3.5) was established for Ω r . In the boundary layer Ω\Ω r the definition (1.2) of the adhesion velocity allows various extensions. For example, one could keep (1.2) by assuming (as done, e.g., in Engwer et al. (2017)) that u := 0 in R n \Ω. An alternative would be to average over the part of the r -ball that lies inside the domain. Let us have a closer look at the first option (the second can be handled similarly). Consider the following example: Example 3.3 (A r vs. T r (∇·) in 1D) Let Ω = (−1, 1), r 0 = 1, F r ≡ 2, and u ≡ 1. In this case, u ≡ 0, hence For A r one readily computes by assuming u = 0 in R\(−1, 1) that for x ∈ (−1, 1) Thus, Example 3.3 supports our idea to average ∇u instead of u itself. The same applies to∇ r u vs. S r (∇u).
Averaging w.r.t. y ∈ B 1 and then also w.r.t. s ∈ (0, 1) might appear superfluous in the definition of the operator T r . The following example compares the effect of T r with that of an operator which averages w.r.t. to y only.

Example 3.4
Let Ω = R n , n ≥ 2, and r > 0, F r ≡ n + 1. In this case Consider also the operator It is easy to see that both operators are well-defined, linear, continuous, and self-adjoint in the space L 2 (R n ). Moreover, they map the dense subspace C 0 (R n ; R n ) into itself. This suggests the following natural extension to (C 0 (R n ; R n )) * : T r μ, ϕ (C 0 (R n ;R n )) * ,C 0 (R n ;R n ) := μ, T r ϕ (C 0 (R n ;R n )) * ,C 0 (R n ;R n ) , T r μ, ϕ (C 0 (R n ;R n )) * ,C 0 (R n ;R n ) := μ, T r ϕ (C 0 (R n ;R n )) * ,C 0 (R n ;R n ) .
Let, for instance, w := δ 0 e 1 , δ 0 and e 1 mean the usual Dirac delta and the vector (1, 0, . . . , 0), respectively. One readily computes that For n ≥ 2, the operator T r retains the singularity at the origin, however making it less concentrated, while T r eliminates that singularity entirely and produces instead jump discontinuities all over S r .

Properties of the averaging operators T r and S r
In this section we collect some properties of the averaging operators T r and S r .
Lemma 3.5 (Properties of T r ) Let F r satisfy Assumptions 1.1 and let r ∈ (0, r 0 ]. Then: The corresponding operator norm satisfies (ii) Let p, p * ∈ [1, ∞] be such that p * = p p−1 . For all w 1 ∈ (L p (Ω)) n and w 2 ∈ L p * (Ω) n it holds: (3.11) (iii) Let p ∈ [1, ∞). For all w ∈ (L p (Ω)) n it holds that (3.12) (iv) For p = 2 it holds that Remark 3.6 Due to the assumptions on F r we have in the limit that (3.14) Proof of Lemma 3.5 (i) Since w is measurable and ρ → F r (ρ), (x, s, y) → x +rsy, (y, z) → (z · y) y |y| are continuous, we have that is well-defined a.e. in Ω × B 1 × (0, 1) and is measurable. Let p ∈ (1, ∞) and p * = p p−1 . Using Hölder's inequality, Fubini's theorem, and our convention that w vanishes outside Ω, we deduce for all w ∈ (L p (Ω)) n that This implies that for all p ∈ (1, ∞) operator T r is well-defined in (L p (Ω)) n and satisfies (3.10). It is also clearly linear. Taken together we then have that T r ∈ L((L p (Ω)) n ) and (3.10) holds. The cases p = 1 and p = ∞ can be treated similarly.
(ii) Let w 1 ∈ (L p (Ω)) n and w 2 ∈ L p * (Ω) n . We compute by using Fubini's theorem, the symmetry of B 1 , and simple variable transformations that Thereby we used our convention that each function defined in Ω is assumed to be prolonged by zero outside Ω. Comparing (3.15) and (3.16) we obtain (3.11). (iii) We apply the Banach-Steinhaus theorem. Due to (i) and (3.14), {T r } r ∈(0,r 0 ] is a family of uniformly bounded linear operators in the Banach space (L p (Ω)) n . Thus, as C c (Ω; R n ) is dense in (L p (Ω)) n for p < ∞, we only need to check (3.12) for w ∈ C c (Ω; R n ). But for such w we can directly pass to the limit under the integral and thus obtain using (3.3) and the dominated convergence theorem that (iv) Here we make use of the Fourier transform, which we denote by the hat symbol.
A straightforward calculation shows that (3.17) Combining (3.17) with the Plancherel theorem and using our convention that w vanishes outside Ω, we can estimate as follows: Here |M| 2 denotes the spectral norm of a matrix M ∈ R n×n . Further, observe that Consequently, denoting by e 1 the first canonical vector of R n and appropriately constructing an orthogonal matrix O in order for Oξ = |ξ |e 1 to hold, we obtain that is a diagonal matrix, its spectral norm is given by the spectral radius. Estimating the right-hand side of (3.21) we then conclude that Finally, the pointwise convergence (3.12) and the Banach-Steinhaus theorem imply that concluding the proof.
A similar result holds for S r : Lemma 3.7 (Operator S r ) Let r ∈ [0, r 0 ]. Then: The corresponding operator norm satisfies Proof The proof almost repeats that of Lemma 3.5. Therefore, we only check (3.24) and omit further details. Let p ∈ [1, ∞) and p * = p p−1 . Using Hölder's inequality, Fubini's theorem, and our convention that w vanishes outside Ω we deduce for all w ∈ (L p (Ω)) n that which means that (3.25) The proof in the case p = ∞ follows the same steps, or, alternatively, one passes to the limit as p → ∞ in (3.25).

Remark 3.8
The constants in (3.10) for any n ≥ 1 and in (3.24) for n ≥ 2 are not necessarily optimal. For p = 2 it remains open whether or not The answer may depend upon Ω and p.

Well-posedness for a class of evolution equations involving T r or S r
In this section we establish the existence and uniqueness of solutions to a certain class of single evolution equations involving T r or S r . This result is an important ingredient for our analysis of nonlocal systems in Sect. 5. Thus, we consider the following initial boundary value problem: Here and for ε ≥ 0 we set A standard calculation shows that G ε is globally Lipschitz with a Lipschitz constant 1.

Remark 4.1
Observe that for ε = 0 equation (4.1a) is linear, whereas for ε > 0 the nonlocal part of the flux is a priori bounded. The latter helps us to construct nonnegative solutions in Sect. 5.
We make the following assumptions: (4.7) To shorten the notation, we introduce a pair of constants Due to assumptions (4.3)-(4.5) it is clear that We introduce a family of operators Lemma 4.2 Let (4.3)-(4.5) be satisfied. Then: is well-defined, monotone, hemicontinuous, and satisfies the bounds is well-defined, monotone, hemicontinuous, and satisfies the bounds Proof The assumptions on the coefficients a i together with the Lipschitz continuity of G ε readily imply that for a.a. t ∈ [0, T ] the operator M(t, ·) is well-defined and satisfies (4.9). Moreover, due to (4.3) and G ε Lipschitz, it is also clear that M(·, u) : is continuous on R, the latter meaning that M(t, ·) is hemicontinuous. Using Hölder's inequality, the fact that G ε is Lipschitz with Lipschitz constant 1, the assumptions on the a i 's, and the properties of R r , we compute that for u, v ∈ H 1 (Ω), which proves monotonicity. Further, taking v = 0 in (4.10) and using M(t, 0) = 0 yields (4.8). Part (i) is thus proved. A proof of (ii) can be done similarly; we omit further details.
Using the properties of the averaging operators proved in Sect. 3.1 we can define weak solutions to (4.1) in a manner very similar to that for the classical, purely local case (i.e., when a 2 ≡ 0): (ii) c r satisfies (4.1a)-(4.1b) in the following sense: for all ϕ ∈ H 1 (Ω) and a.a.
Using standard theory one readily proves the following existence result: Lemma 4.4 Let (4.3)-(4.7) hold. Then there exists a unique weak solution to (4.1) in terms of Definition 4.3. The solution satisfies the following estimates: Proof The existence of a unique weak solution to (4.1) is a direct consequence of Lemma 4.2(i) and the standard theory of evolution equations with monotone operators, see, e.g. Showalter (1997), Chapter III Proposition 4.1). It remains to check the bounds (4.12), and (4.13). Taking ϕ := c r in the weak formulation (4.11) and using (Temam 2001, Chapter III Lemma 1.2), (4.8), and the Young inequality, we obtain that which yields (4.12) due to the Gronwall lemma. Finally, using (4.9), we obtain from the weak formulation (4.11) that Together with (4.12) this implies (4.13).

Nonlocal models involving averaging operators T r and S r
In this section we study the following model IBVP: Here, as in the previous section, R r stands for any of the two averaging operators: We assume that the diffusion coefficient D v is either a positive number, or it is zero. Equations (5.1a)-(5.1b) are closely related to (1.1) and (1.4) in Sect. 1, the difference being that the terms involving the adhesion velocity/non-local gradient are now replaced by those including the averaging operators T r /S r from Sect. 3. Our motivation for introducing this change is twofold. First of all, due to (3.5) and (3.8) it affects the points in the boundary layer Ω\Ω r , at the most. On the other hand, Example 3.3 indicates that including, e.g., A r can lead to limits with unexpected blow-ups on the boundary of Ω.
System (5.1) is a non-local version of the hapto-/chemotaxis system In this case, the actual diffusion and haptotactic sensitivity coefficients are so that in the classical formulation (5.2a) takes the form The main goal of this section is to establish, under suitable assumptions on the system coefficients which are introduced in Sect. 5.1, a rigorous convergence as r → 0 of solutions of the nonlocal model family (5.1) to those of the local model (5.2), see Theorem 5.8. This is accomplished in the final Sect. 5.4. Since we are dealing here with a new type of nonlocal system, we establish for (5.1) the existence of nonnegative solutions in Sects. 5.2, and 5.3.

Problem setting and main result of the section
We begin with several general assumptions about the coefficients of system (5.1).
Assume that the coefficients satisfy the following bounds: Further, we assume that the initial values satisfy then assumption (5.5) can be replaced by a weaker one, such as v 0 ∈ L 2 (Ω).
We keep (5.5) in order to simplify the exposition.
In addition, we will later choose one of the following assumptions on f c and the nonlocal operator: Assumptions 5.3 (Further assumptions on f c ) One of the following conditions holds: Assumptions 5.4 (Assumptions on R r ) One of the following holds: (a) for a given fixed r ∈ (0, r 0 ] and assume that Then, it holds a priori that for any v which solves (5.1b). Therefore it suffices to consider the coefficient functions For C 7 :=μ c (K c + 1 + η c K v ), C 8 :=μ c (K c + 1) and C 9 :=μ c we can estimate on Further, holds. Thus, Assumptions 5.1, 5.3(b) and 5.4 (b) are fulfilled if This choice of coefficient functions can be used to describe a population of cancer cells which interact among themselves and with the surrounding extracellular matrix (ECM) tissue. Both interaction types are due to adhesion, whether to other cells (cellcell adhesion) or to the matrix (cell-matrix adhesion). The interaction force F r (ρ) is taken to diminish with increasing interaction range ρ and/or of the sensing radius r : cells too far apart/out of reach hardly interact in a direct way. Function g(c, v) characterises effective interactions. Here the coefficients S cc and S cv represent cellcell and cell-matrix adhesion strengths, respectively. Our choice of g accounts for some adhesiveness limitation imposed by high local cell and tissue densities. It is motivated by the fact that overcrowding may preclude further adhesive bonds, e.g. due to saturation of receptors. The diffusion coefficient D c (c, v) is chosen to be everywhere positive and increase with a growing population density, thus enhancing diffusivity under population pressure, but, further, limited by excessive cell-tissue interaction. The latter also applies to the choice of the sensitivity function χ . Indeed, there is evidence that tight packing of cells and ECM limits diffusivity and the advective effects of haptotaxis (Lu et al. 2012). Thereby the constant b > 0 is assumed to be rather small. Finally, f c and f v describe growth of cells and tissue limited by concurrence for resources.

Remark 5.7
Observe that for r = 0 we obtain a corresponding solution definition for the local system (5.2).
Our main result now reads: This Theorem is proved in Sect. 5.4.

Notation 5.9
Dependencies upon such parameters as the space dimension n, domain Ω, function c, the norms of the initial data c 0 and v 0 , norms and bounds for the coefficient functions are mostly not indicated in an explicit way.

Global existence of solutions to (5.1): the case of f c Lipschitz
In this subsection we address the existence of solutions to the nonlocal model (5.1) for the case when f c satisfies Assumptions 5.3(a). The main result of the subsection is as follows: Theorem 5.10 Let Assumptions 1.1, 5.1, and 5.3(a) hold and let r satisfy Assumption 5.4(a). Then there exists a global weak-strong solution to (5.1) in terms of Definition 5.6 with ∂ t c r ∈ L 2 (0, T ; (H 1 (Ω)) * ).
Since we aim at constructing nonnegative solutions, it turns out to be helpful to consider first the following family of approximating problems: where G ε was defined in (4.2). In order to obtain existence for the original problem, i.e., for ε = 0, we first prove existence of nonnegative solutions for the cases when ε, D c > 0. This corresponds to a chemotaxis problem with a nonlocal flux-limited drift. Weak-strong solutions to (5.9) are understood as in Definition 5.6, with the obvious modification of the weak formulation, which now reads: Lemma 5.11 Let Assumptions of Theorem 5.10 be satisfied. Assume further that Then there exists a global weak-strong solution to (5.9) with ∂ t c r ε ∈ L 2 (0, T ; (H 1 (Ω)) * ).
Proof To begin with, we extend the coefficients: for c < 0 These coefficients still satisfy Assumptions 5.1, 5.3(a), and 5.4(a) if we consider all suprema over c ∈ R instead of c ∈ R + 0 . Our approach to proving existence is based on the classical Leray-Schauder principle (Zeidler 1986, Chapter 6, §6.8, Theorem 6.A). In order to apply this theorem we first 'freeze' c r ε in the system coefficients of (5.9), replacing it byc r ε . Correspondingly, we obtain the following weak formulation in place of (5.10): For all ϕ ∈ H 1 (Ω) and a.a. t > 0 Let T > 0 and letc r ε ∈ L 2 (0, T ; L 2 (Ω)). Since f v is assumed to be Lipschitz, we can make use of the standard theory (Ladyzhenskaya et al. 1968) which implies that the semilinear parabolic initial boundary value problem (5.11c)-(5.11e) possesses a unique global strong solution 0 ≤ v r ε ∈ L 2 (0, T ; H 2 (Ω)) with ∂ t v r ε ∈ L 2 (0, T ; L 2 (Ω)), and satisfying the estimate Here and further in the proof we omit the dependence of constants upon D v . Set Due to our assumptions about D c , χ, g, and f c , these coefficients a i and f satisfy the requirements of Lemma 4.2. Consequently, there exists a unique global weak solution c r to problem (4.1) with these coefficients. We estimate for the corresponding constants α r and M r introduced in Lemma 4.2: α r ≥ C 5 C 10 (r ) =: C 15 (r ), (5.13) M r ≤ C 6 + C 12 C 13 R r L((L 2 (Ω)) n ) =: C 16 (r ), (5.14) and, due to (5.12), ( 5.15) Combining (4.12)-(4.13) and (5.13)-(5.15), we obtain the following bounds for c r : Next, we verify that Φ is closed in L 2 (0, T ; L 2 (Ω)). Consider a sequence By the definition of Φ we have thatc r εm and c r εm satisfy: for all ϕ ∈ H 1 (Ω) and a.a. t ∈ (0, T ) (5.26e) From (5.12) and (5.21) we conclude that the sequence {v r εm } is uniformly bounded in L 2 (0, T ; H 2 (Ω)) and ∂ t v r εm ∈ L 2 (0, T ; (L 2 (Ω)). Hence the Lions-Aubin lemma and the Banach-Alaoglu theorem imply that there exists v r ε s.t. (after switching to a subsequence, if necessary) and this v r ε satisfies equation (5.11c) forc r ε as well as the initial and boundary conditions in the required sense. Further, due to (5.24), and (5.25) we have in the usual way that In particular, i.e. the initial condition is satisfied.
(5.29) Next, from (5.23) and (5.27) we conclude using the boundedness and continuity of functions G ε , ∇g, ∇ f c , and (c, v) → cχ(c, v) over R × R + 0 and of operator R r in L 2 (Ω) and the dominated convergence theorem that so that due to (5.24) and the compensated compactness Observe that the weak formulation (5.26a) is equivalent to As c r ε satisfies (5.32), it follows from the last inequality that holds for all w ∈ L 2 (0, T ; H 1 (Ω)).

Remark 5.12
Observe that c r ε cannot be replaced by −(c r ε ) − inside the nonlocal operator. This is why we introduced the flux-limitation. Now we are ready to prove Theorem 5.10.

Proof of Theorem 5.10
We start with the case Lemma 5.11 gives the existence of solutions (c r ε , v r ε ) to (5.9). Setting ϕ = c r ε in (5.10), using the facts that f c is Lipschitz and |G ε (x)| ≤ |x|, we can estimate similarly to Theorem 5.13 below and obtain upper bounds of the form (5.41)-(5.47), which are independent from ε (with p = q = 2 there). Applying the Lions-Aubin lemma and the Banach-Alaoglu theorem, we conclude the existence of a pair of nonnegative functions c r and v r having the regularity stated in Definition 5.6 and such that for a sequence ε m → m→∞ 0 it holds that where the last term tends to 0 as ε m → m→∞ 0. As the term inside the integral is moreover bounded in L 2 (0, T ; L 2 (Ω)) by a constant independent from ε m , we conclude by using a result from Evans (1990, p. 6 From this and the boundedness of ∇c r ε m L 2 (0,T ;(L 2 (Ω)) n ) , (5.37)-(5.39), Lemmas 3.5 or 3.7(i) and (ii), respectively, the fact that |G ε (x)| ≤ |x|, the continuity of ∂ c g, χ, (5.3), (5.4), compensated compactness, the dominated convergence theorem, and the Hölder inequality, we obtain that for all ψ ∈ L 2 (0, T ; H 1 (Ω)) it holds that The convergence to the remaining terms in (5.8a) and the rest of (5.8) can be obtained in a way either completely analogous or very similar to the corresponding parts of the proof of Lemma 5.11. In order to prove existence for the case 1). Estimating similarly to the proof of Theorem 5.13 below and performing a standard limit procedure based on the Banach-Alaoglu theorem, the dominated convergence theorem, the Lions lemma (Lions 1969, Lemma 1.3), and the compensated compactness, one readily obtains a solution (c r 0 , v r 0 ) for D v = 0 in the sense of Definition 5.6. Observe that this time the gradient of v-component enters linearly, so that no strong convergence is required. We omit further details.

Global existence of solutions to (5.1): the case of f c dissipative
In this subsection we provide an extension of the existence Theorem 5.10 from Sect. 5.2: Then there exists a global weak-strong solution to (5.1) in terms of Definition 5.6, with ∂ t c r ∈ L q (0, T ; (W 1,q * (Ω)) * ) and satisfying the following estimates: For all T > 0 ||c r || L ∞ (0,T ;L 2 (Ω)) ≤ C 22 (T , R r L((L 2 (Ω)) n ) ), (5.41) where η k is a cut-off function: Since f ck is Lipschitz, Theorem 5.10 implies the existence of a solution (c rk , v rk ) in terms of Definition 5.6 with ∂ t c rk ∈ L 2 (0, T ; (H 1 (Ω)) * ), which corresponds to f c = f ck . Our next aim is to prove that (c rk , v rk ) satisfies the same bounds as in the statement of the Theorem with some constant C 22 (T , R r L((L 2 (Ω)) n ) ) which does not depend upon k. Set Taking ϕ := c rk in (5.8a) written for c rk and using Assumptions 5.1, 5.3(b), 5.4(a) and the Hölder and Young inequalities, we compute Next, we estimate v rk . If D v > 0, then standard theory (Ladyzhenskaya et al. 1968) yields that for all Here and further in the proof we omit the dependence of constants upon D v . If D v = 0, then we get the ODE Hence, the assumptions on f v and the solution components together with the chain rule imply that Computing the gradient on both sides of (5.52), multiplying by ∇v rk throughout, integrating over Ω, and using Assumptions 5.1 and the Young inequality, we obtain that ∇c rk (L 2 (Ω)) n ∇v rk (L 2 (Ω)) n ≤ C 28 ∇v rk 2 (L 2 (Ω)) n + C 29 ∇c rk 2 (L 2 (Ω)) n . (5.53) Applying the Gronwall inequality to (5.53) yields From (5.6) and (5.58), the embedding of Lebesgue spaces, and η k ∈ [0, 1] we conclude that so that (5.47) holds for f ck (c rk , v rk ). Combining (5.41) and (5.42) for c rk with (5.51) or (5.56) and (5.57) (depending on the sign of D v ) and using the equation for v rk yields such bounds as (5.44)-(5.46) and (5.48) for c rk and v rk . Finally, combining Assumptions 5.1 with bounds on ∇c rk , ∇v rk , and f ck (c rk , v rk ), the weak formulation (5.8a), and estimating in a standard way yields (5.43) for ∂ t c rk .
Since (c rk , v rk ) satisfy (5.41)-(5.48) uniformly in k, a standard limit procedure based on the Banach-Alaoglu theorem, the dominated convergence theorem, the Lions lemma, and the compensated compactness yields the existence of a weak-strong solution (c r , v r ) to (5.8) which satisfies (5.41)-(5.48).

Limiting behaviour of the nonlocal model (5.1) as r → 0
In this subsection we finally prove our main result concerning convergence for r → 0.
Since for each such r m the Assumptions 5.4(a) are satisfied, Theorem 5.13 is applicable and yields the existence of solutions (c r m , v r m ) which satisfy (5.41)-(5.48). Replacing R r by C 11 in C 22 (T , R r L((L 2 (Ω)) n ) ) makes the constant in (5.41)-(5.48) independent of m. Using the Lions-Aubin lemma and the Banach-Alaoglu theorem we conclude (by possibly switching to a subsequence) that c r m m→∞ c, v r m m→∞ v in L 2 (0, T ; H 1 (Ω)). (5.60) Using standard arguments based on the Banach-Alaoglu theorem, the dominated convergence theorem, the Lions lemma, and assumptions on χ and g we conclude from (5.59) and (5.60) that Observe that for any ψ ∈ L ∞ (0, T ; W 1,∞ (Ω)) the following estimate holds: Now, using (5.61) together with Lemma 3.5(i) and (iii) and (3.6) or (3.7)(i) and (iii) and (3.9), respectively, we conclude that the right hand side of (5.64) tends to zero, hence Thus, using Lemma 3.5(ii) or Lemma 3.7(ii), respectively, and compensated compactness, we obtain from (5.62) and (5.65) that The convergence in the remaining terms, equations, and conditions follows by means of a standard limit procedure based on the Banach-Alaoglu theorem, the dominated convergence theorem, the Lions lemma, and the compensated compactness. We omit these details.

Numerical simulations in 1D
We perform numerical simulations to investigate on the one hand the effect of differences between hitherto choices of nonlocal operators and our novel ones proposed in Sect. 3, and on the other hand convergence between nonlocal and local formulations.
For compactness, our current study restricts to the prototypical nonlocal model for cellular adhesion (1.1), its reformulation as (5.1), and the corresponding local model (5.2).
Thus, for (5.1) we take the operator form R r = T r , with T r as in (3.2). These models can be interpreted in the context of a population of cells invading an adhesion-laden ECM/tissue environment and, with this in mind, we initially concentrate cells at the centre of a one-dimensional domain Ω = [0, L] and impose an initially homogeneous ECM. Specifically, we set for the ECM v 0 (x) = 1, x ∈ Ω (6.1) and consider for the cell population a Gaussian-shaped aggregate where we set x c = L/2 or x c = 0. The numerical scheme follows that described in Gerisch (2010), which we refer to for details. Briefly, a Method of Lines approach is invoked whereby equations are first discretised in space (in conservative form, via a finite volume method) to yield a highdimensional system of ODEs, which are subsequently integrated in time. Discretisation of advective terms follows a third order upwinding scheme, augmented by flux limiting to preserve positivity of solutions and the resulting scheme is (approximately) secondorder accurate in space. Time integration has been performed with standard Matlab ODE solvers: our default is "ode45" with absolute and relative error tolerances set at 10 −6 , but simulations have been compared for varying space discretisation step, ODE solver, and error tolerances. To measure the difference between two distinct solutions over time we define a distance function as follows: where u 1 and u 2 denote the two solutions that are being compared.

Comparison of nonlocal operator representations
We first explore the correspondence between forms of nonlocal operator representation: we choose the prototypical nonlocal model for cell/matrix adhesion (1.1) and its reformulation (5.1), therefore taking for the latter the operator form R r = T r with T r as in (3.2). In what follows, solutions to (1.1) are denoted c A and v A and those for (5.1) denoted c T and v T . For simplicity we restrict in this section to a minimalist formulation in which D c = constant, χ = 1, f c = 0. Cell-matrix interactions are defined by g(c, v) = S cc c + S cv v and f v (c, v) = −μcv, where S cc and S cv respectively represent cell-to-cell and cell-to-matrix adhesion strengths and f v simplistically describes (direct) proteolytic degradation of matrix by cells parametrised by degradation rate μ. Figure 1 shows the computed solutions under (a-c) negligible cell-cell adhesion (S cc = 0) and (d-f) moderate cell-cell adhesion (S cc = S cv /4). The equivalence of the two formulations is revealed through the negligible difference between solutions, with the distance magnitude attributable to the subtly distinct numerical implementation.
Both simulations describe an invasion/infiltration process, in which matrix degradation by the cells generates an adhesive gradient that pulls cells into the acellular surroundings. The impact of cell-cell adhesion is manifested in the compaction of cells at the leading edge into a tight aggregate.
However, as pointed out in Sect. 3, differences in the nonlocal formulations can emerge in the vicinity of boundaries. To highlight this we consider an equivalent formulation to Fig. 1a-c, but with the cells initially placed at the left boundary [x c = 0 in (6.2)], e.g. suggesting a tumor mass which is concentrated there and whose cells are expected to detach and migrate into the considered 1D domain, travelling from left to right. As stated earlier we impose zero-flux boundary conditions at x = 0 (and x = L), and further suppose c = v = 0 and ∇c = ∇v in the extradomain region (R\Ω). Representative simulations are shown in Fig. 2. They are in agreement with our observation in Example 3.3. Indeed, for this scenario, in the prototypical nonlocal model (1.1)-(1.2) there is a very large adhesion velocity modulus at x = 0; the cells are crowded within the tumor mass and their mutual interactions are maintained during the invasion process in a sufficiently strong manner to ensure a collective shift of the still concentrated cell aggregate, with a correspondingly strong tissue degradation in its wake. In the reformulation (5.1)-(3.2), rather, the adhesion magnitude at x = 0 is for the same initial condition much lower -suggesting a tumor whose cells are readier to detach and migrate individually. This results in a more diffusive spread, with accordingly less degradation of tissue, and with cell mass remaining available at the original site over a larger time span. The latter scenario is different from the former one, but it seems nevertheless reasonable, as a tumor mass would very often not move as a whole from its original location to another in a relatively short time;  Fig. 1ac, but with the cells initially concentrated at the boundary. d, f Comparison of the two forms of nonlocal operator corresponding to the simulations represented in (a-c). The operators are practically identical sufficiently far from the boundary, but can diverge significantly for distances < r from the boundaries moreover, the active cells in a sufficiently large tumor (releasing substantial amounts of acidity) are known to preferentially adopt a migratory phenotype and perform EMT (epithelial-mesenchymal transition), see e.g., Gupta (2015), Peppicelli et al. (2014) and Prieto-García et al. (2017), which supports the idea of cells moving in a loose way rather than in compact, highly aggregated assemblies. 3 As such, our simulations suggest that, within this particular function-and parameter setting, choosing the adhesion operator in the form (1.2) instead of (3.2) might possibly overestimate the tumor invasion speed and associated healthy tissue degradation, thereby predicting a spatially concentrated tumor and neglecting regions with lower cell densities which can nevertheless trigger tumor recurrence if untreated.

Comparison between nonlocal and local formulation
Having compared together the original, (1.1), and the new, (5.1), nonlocal formulations, we next consider the extent to which their dynamics can be captured by the classical local formulation (5.2). Note that for nonlocal model simulations we will restrict to the original formulation (1.1), so that we can avail ourselves of an already well-established efficient (in terms of computational time) numerical scheme (Gerisch 2010). Here we use c L and v L to denote solutions to the local formulation and c Ar and v Ar to denote solutions to the nonlocal model with sensing radius r . We remark that a large number of related local and nonlocal models have been numerically studied to describe the invasion-type process considered here (e.g. Perumpanani et al. 1996;Anderson et al. 2000;Gerisch and Chaplain 2008;Painter et al. 2010): here the specific focus is to explore the convergence of nonlocal to local form as r → 0, which, as far as we are aware, has not been systematically investigated.
As in the first test we use the initial values (6.1) and (6.2), choosing x c = L/2, α = 10 in the latter, and consider the coefficients and functions as proposed in Example 5.5. Under these choices the resultant nonlinear diffusion coefficient for the c-equation in the classical local formulation (compare (5.2a)) becomes Notably, this potentially becomes negative under an injudicious combination of adhesive strengths S cc , S cv , and of a, b. Likewise, the actual haptotaxis sensitivity function takes the formχ (6.4) Again, depending on the relationship between S cc and S cv , this can become negative, which would lead to repellent haptotaxis: cells effectively moving away from regions with large ECM gradients, a rather unexpected behaviour. This suggests that celltissue adhesions should dominate over cell-cell adhesions, 4 as 'usual' haptotaxis, i.e. towards the increasing tissue gradient, is known to be an essential component of cell migration, this applying to several types of cells moving through the ECM (tumor cells, mesenchymal stem cells, fibroblasts, endothelial cells, etc.) see e.g. Lamalice et al. (2007), Pickup et al. (2014) and Wen et al. (2015) and references therein.
Simulations are plotted in Fig. 3 where we show cell densities for the local model (c L ) and nonlocal model under three sensing radii: In this first set of simulations we assume negligible cell-cell adhesion (S cc = 0), which automatically ensures positivity for the diffusion coefficient of the equivalent local model,D c (c, v). We note that matrix renewal is absent (μ v = 0) in the left-hand column and present (μ v > 0) in the central column. In the right-hand column we show the greater generality of the results under vastly simplified kinetics, specifically setting f c (c, v) = 0 and f v (c, v) = −cv (with the other functional forms as in Example 5.5). Simulations highlight the convergence between local and nonlocal models as r → 0: for r = 0.1, the solution differences become negligible. However, distinctions emerge for large r , where we can expect significant discrepancy between the solutions. This suggests that the local model fails to accurately predict the behaviour in cases where cells sample over relatively large regions of their local environment.
Next, we extend to include a degree of cell-cell adhesion, setting functions and parameters as in Fig. 3, except now S cc > 0. Notably this raises the possibility of a negative diffusion coefficient in the classical formulation and subsequent illposedness. Solutions under a representative set of parameters are shown in Fig. 4. For t below  ([0, 20]), the bottom row magnifies a relevant portion for clarity. Solutions to local and nonlocal models under the functional forms proposed in Example 5.5 for r = 0.01, 0.1, 0.3, 1.0 at a t = 3, b t = 3.5 and c t = 5. In (a) solutions to the local model continue to exist and we observe convergence between local and nonlocal formulations. In (b, c) the solutions to the local model are noncomputable. Nonlocal models, however, can destabilise into a pattern of aggregates. Parameters: a = 0.01, b = 1, μ c = 0.01, K c = 2, η c = 1, μ v = 0, λ v = 1 and adhesion parameters as above . Negligible cell-cell adhesion, S cc = 0, S cv = 10: solutions shown at (left) t = 2.5 and (middle) t = 5, with the distance between solutions to the nonlocal and local model shown in the right panel some critical time we observe convergence as before, with the nonlocal formulation converging to solutions of the local model as r → 0. However, continued matrix degradation further depletes v, with the result that (6.3) can become negative. At this point (in this case t ≈ 3.2 . . .) the local model becomes illposed and its solutions become incomputable (implying nonexistence of solutions). However, the nonlocal formulation appears to preserve wellposedness, consistent with previous theoretical studies where extending to a nonlocal formulation regularises a singular local model (e.g. Hillen et al. 2007). Solutions to the nonlocal model instead destabilise into a quasi-periodic pattern of cell aggregations, maintained through the cell-cell adhesion, and with a wavelength shrinking as r → 0.
Finally, we remark that convergence of solutions extends beyond the specific functional forms and, as a representative example, we consider a minimalist setting based on linear/constant forms. Specifically, we set D c = a (constant), χ = 1, f c = 0, g(c, v) = S cc c + S cv v and f v (c, v) = −μcv. In this scenario, the diffusion and haptotaxis coefficients for the classical local formulation (5.2) reduce tõ D c (c, v) = a − S cc c andχ(c, v) = S cv . (6.5) Positivity is only guaranteed under appropriate parameter selection. Such a case is illustrated in Fig. 5 where we assume negligible cell-cell adhesion (S cc = 0). Clearly, we observe convergence between the nonlocal and local formulations as r → 0. Inappropriate parameter selection, however, generates backward diffusion in the local model and solutions are consequently incomputable. In all cases considered in this test the cells do not reach the boundary region where the difference between the nonlocal formulations (1.1) and (5.1) can play a role. Thus, we expect the same solution if reformulation (5.1) is applied instead.

Discussion
In this work we provide a rigorous limit procedure which links nonlocal models involving adhesion or a nonlocal form of chemotaxis gradient to their local counterparts featuring haptotaxis, respectively chemotaxis in the usual sense. As such, our paper closes a gap in the existing literature. Moreover, it offers a unified treatment of the two types of models and extends the previous mathematical framework to settings allowing for more general, solution dependent, coefficient functions (diffusion, tactic sensitivity, adhesion velocity, nonlocal taxis gradient, etc.). Finally, we provide simulations illustrating some of our theoretical findings in 1D.
Our reformulations in terms of T r and S r reveal the tight relationship between the nonlocal operators A r and∇ r and the (local) gradient. This suggests that both nonlocal descriptions (adhesion, chemotaxis) actually encompass the dependence on the signal gradients rather than on the signal concentration/density itself, which is in line with the biological phenomenon. Indeed, through their transmembrane elements (e.g. receptors, ion channels etc.) the cells are mainly able to perceive and respond to differences in the signal at various locations or within more or less confined areas rather than measure effective signal concentrations. Along with the mentioned solution dependency of the nonlocal model coefficients, the influence of the gradient possibly reflects into contributions of the adhesion/nonlocal chemotaxis to the (nonlinear) diffusion in the local setting obtained through the limiting procedure.
The set Ω r (as introduced in Sect. 2) can be regarded as the 'domain of restricted sensing', meaning that there cells a priori sense only what happens inside Ω, the domain of interest. The measure of this subdomain is a decreasing function of the sensing radius r . When r → 0 the set Ω r tends to cover the whole domain Ω, whereas as r increases the cells can sense at increasingly larger distances; correspondingly, Ω r shrinks. For r > diam(Ω) the restricted sensing domain is empty: everywhere in Ω the cells can perceive signals not only from any point within Ω but potentially also from the outside. In this paper, however, we look at models with no-flux boundary conditions. This corresponds, e.g., to the impenetrability of the walls of a Petri dish or that of comparatively hard barriers limiting the areas populated by migrating cells, e.g. bones or cartilage material. As a result, the cells in the boundary layer Ω\Ω r have a much reduced ability to stretch their protrusions outside Ω and thus gain little information from without. To simplify matters, we assume in this work that there is no such information or it is insufficient to trigger any change in their behaviour. In the definitions of T r and S r this corresponds to the integrands being set to zero in Ω\Ω r .
It is important to note that for points x ∈ Ω\Ω r the influence of a signal p in a direction y ∈ S 1 is not taken into account by∇ r at all if x + r y / ∈ Ω. If S r is used instead, then its contribution to the average is given bỹ y := n 1 0 χ Ω ∇ p(x + rsy) ds · y y.
Thus, thanks to integration w.r.t. s, the resulting vectorỹ assembles the impact of those parts of the segment connecting x and x + r y which are contained in Ω. It is parallel to y, and it may have the same or the opposite orientation. In particular this means that although for a certain range of directions large parts of the sensing region of a cell are actually outside Ω, this may still strongly influence the speed and actual direction of the drift. The effect of integration w.r.t. s in T r is less obvious, since in this case the average w.r.t. y is computed over the ball B 1 . This already achieves the covering of the whole sensing region by allowing a cell to gather information about the signal not only in any direction y/|y|, but also at any distance less than r . The additional integration over the path x + rsy, s ∈ [0, 1], appears to mean that cells at x ∈ Ω r are able to measure the average of the signal gradient all along such line segment rather than its value directly at the ending point. Indeed, from a biological viewpoint this description seems to make more sense, as cells do not jump from one position to another, nor do they send out their protrusions in a discontinuous way bypassing certain space points along a chosen direction. Averages over cell paths are then averaged w.r.t. y, which finally determines the direction of population movement. Example 3.4 indicates that the effect of even an extremely concentrated signal gradient is mollified by averaging. This agrees with our expectations from using non-locality. In higher dimensions n ≥ 2, the two-stage averaging in T r (w.r.t. s and y) produces a direction field which is smooth away from the concentration point and also weakens but still keeps the singularity there. In contrast, averaging only w.r.t. y leads instead to jump discontinuities at a unit distance from the accumulation point. Moreover, we remark that without integrating w.r.t. s in T r (∇·) one cannot regain A r .
The effect observed in Example 3.3 further supports the conjecture that the nonlocal operators which act directly on the signal gradients might actually be a more appropriate modelling tool. While inside the subdomain Ω r there is no difference (recall Lemmas 3.1 and 3.2), inside the boundary layer Ω\Ω r the limiting behaviour as r → 0 is qualitatively distinct. Indeed, Example 3.3 shows that using, e.g., A r , leads, for r → 0, to unnatural sharp singularities at the boundary of Ω even in the absence of signal gradients, whereas this does not happen if T r is used instead. Simulations in Sect. 6.1 (see Fig. 2) confirm our theoretical findings and show a substantial difference between the solutions obtained with the two nonlocal formulations involving (1.2) and (3.2), respectively. The choice (3.2) is motivated above all from a mathematical viewpoint (as it enables a rigorous, well-justified passage to the limit for r → 0), but it also seems to make sense biologically, as our above comments and the simulations performed for the particular setting in Sect. 6.1 suggest.
In this work we have only dealt with models that include a nonlocality in the chemotaxis or cell-cell and/or cell-tissue adhesion terms and assumed the diffusion to be local. This is in line with most of the previously developed nonlocal models for cell migration, albeit they usually cover just linear diffusion. If cell-cell adhesion is present, this means that the cell flux contains the local cell gradient, as well as some averaging of it. The latter is described in our case by a suitably chosen operator T r . A possible model extension could involve a diffusion flux which is also nonlocal and has a similar form. This would mean that the cell flux is completely devoid of the local gradient. From the modelling point of view this could be seen as a population pressure acting 5 in a nonlocal manner: each cell is sensing the population mass not only at its current position, but over a whole region (of radius r ) around that location. This is actually true in vivo, where cells sample their biological environment by extending protrusions as far as several cell lengths. While cell-cell adhesions certainly play a role in this process and contribute to self-diffusion (as in the example handled in Sect. 6.2), there might be yet other ways of interaction by which the cells are able to perceive smaller or larger aggregates of their own kind. In this context one could think about replacing the local gradient by a nonlocal operator, e.g. of the form T r (∇). However, the analysis of such a model would be considerably more involved and it is to expect that existence of solutions can be established only under rather restrictive assumptions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.