1 Introduction
Background
Mathematical models formulated as reaction–diffusion equations with non-local reaction terms have been increasingly used to achieve a more in-depth theoretical understanding of the mechanisms underlying the spatial spread and the phenotypic evolution of populations with heterogeneous motility [Reference Arnold, Desvillettes and Prévost4, Reference Berestycki, Mouhot and Raoul9–Reference Bouin, Calvez, Meunier, Mirrahimi, Perthame, Raoul and Voituriez11, Reference Bouin, Henderson and Ryzhik13, Reference Turanova46].
In these models, the phenotypic state of each individual is described by a continuous structuring variable, and the model itself consists of a balance equation for the local population density function (i.e. the phenotypic distribution of the individuals at each spatial position). As is the case for the classical Fisher-KPP model [Reference Fisher20, Reference Kolmogorov31], individuals are assumed to undergo undirected, random movement, which translates into a linear diffusion term. Additionally, intrapopulation variability of individual motility is taken into account by letting the diffusion coefficient be a function of the structuring variable. Moreover, possible changes in individual motility are conceptualised as transitions between phenotypic states, which are modelled through an integral or a differential operator. Finally, in analogy with the non-local version of the Fisher-KPP model [Reference Berestycki, Nadin, Perthame and Ryzhik8, Reference Hamel and Ryzhik27], most of these models rely on the assumption that the population undergoes logistic growth at a rate that depends on the local number density of individuals (i.e. the integral of the solution with respect to the structuring variable), which is described via a non-local reaction term.
Amongst these models, the model for the cane toad invasion presented in [Reference Bénichou, Calvez, Meunier and Voituriez7] has received considerable attention from the mathematical community over the last few years. Analysis of this simple yet effective model has made it possible to find a robust mechanistic explanation for the empirical observation that highly motile individuals are, as such, more likely to be found at the edge of the invasion front, and has helped elucidate the way this form of spatial sorting can promote acceleration of the invasion front [Reference Phillips, Brown, Webb and Shine41–Reference Shine, Brown and Phillips43, Reference Urban, Phillips, Skelly and Shine47]. In particular, the existence of travelling-front solutions and the occurrence of spatial sorting in the case of bounded motility has been studied in [Reference Bouin and Calvez10–Reference Bouin, Henderson and Ryzhik12, Reference Turanova46], while front acceleration in the case of unbounded motility has been investigated in [Reference Berestycki, Mouhot and Raoul9, Reference Bouin, Calvez, Meunier, Mirrahimi, Perthame, Raoul and Voituriez11, Reference Bouin, Henderson and Ryzhik13]. Furthermore, an evolution equation for the dynamic of the maximum point of the local population density function along the phenotypic dimension (i.e. the dominant phenotypic trait) at the edge of the front has been formally derived in [Reference Bouin, Calvez, Meunier, Mirrahimi, Perthame, Raoul and Voituriez11].
Content of the paper
We consider a model for the dynamics of growing cell populations with heterogeneous mobility and proliferation rate. In analogy with the models considered in the aforementioned studies, intra-population heterogeneity is here captured by a continuous structuring variable representing the cell phenotypic state and the model consists of a balance equation for the local cell population density function. However, in contrast to the aforementioned studies, such a balance equation takes the form of a non-local advection–reaction–diffusion equation whereby the velocity field and the reaction term are both functions of the structuring variable and of the local cell density. This leads to the emergence of invasion fronts with compact support and brings about richer spatio-temporal dynamics of the dominant phenotypic trait throughout the front.
Outline of the paper
The remainder of the paper is organised as follows. In Section 2, we describe the model and the main underlying assumptions. In Section 3, we present the results of numerical simulations, which were obtained using the numerical methods detailed in Appendix A. In Section 4, we carry out formal asymptotic analysis of the model in order to provide an explanation for such numerical results. In Section 5, we discuss the main results of numerical simulations and formal analysis. Moreover, we briefly explain how these mathematical results may shed light on the interplay between spatial sorting and natural selection that underpins tumour growth and the emergence of phenotypic heterogeneity in glioma. Finally, we provide a brief overview of possible research perspectives.
2 Statement of the problem
A model for the dynamics of heterogeneous growing cell populations
We consider a mathematical model for the dynamics of a growing population of cells structured by a variable 
$y \in [0,Y] \subset \mathbb{R}_+$, which represents the phenotypic state of each cell and takes into account intra-population heterogeneity in cell proliferation rate and cell mobility (e.g. the variable y could represent the level of expression of a gene that controls cell proliferation and cell mobility). The population density at position 
$x \in \mathbb{R}$ and time 
$t \in [0,\infty)$ is modelled by the function n(t,x,y), the evolution of which is governed by the following non-local partial differential equation (PDE)
\begin{equation} \begin{cases}\displaystyle{\partial_t n - \alpha \, \mu(y) \, \partial_x \left(n \, \partial_x \rho(t,x) \right) = R(y,\rho(t,x)) \, n + \beta \, \partial^2_{yy} n,}\\\\\displaystyle{\rho(t,x) := \int_{0}^Y n(t,x,y) \, \textrm{d}y,}\end{cases}\quad(x,y) \in \mathbb{R} \times (0,Y),\end{equation}
subject to zero Neumann boundary conditions at 
$y=0$ and 
$y=Y$.
The second term on the left-hand side of the non-local PDE (2.1) represents the rate of change of the population density due to the tendency of cells to move towards less crowded regions (i.e. to move down the gradient of the cell density 
$\rho(t,x)$) [Reference Ambrosi and Preziosi3, Reference Byrne and Drasdo14]. The function 
$\alpha \, \mu(y)$, with 
$\alpha > 0$, models the mobility of cells in the phenotypic state y. Without loss of generality, we consider the case where higher values of y correlate with higher cell mobility and, therefore, we let 
$\mu(y)$ be a smooth function that satisfies the following assumptions:
\begin{equation} \mu(0) > 0, \quad \dfrac{\textrm{d} \mu(y)}{\textrm{d} y} > 0 \; \text{ for } y \in (0,Y].\end{equation}
Moreover, the first term on the right-hand side of the non-local PDE (2.1) represents the rate of change of the population density due to cell proliferation and death. The function 
$R(y,\rho(t,x))$ models the fitness (i.e. the net proliferation rate) of cells in the phenotypic state y at time t and position x under the local environmental conditions given by the cell density 
$\rho(t,x)$. We let 
$R(y,\rho)$ be a smooth and bounded function that satisfies the following assumptions:
\begin{equation} R(Y,0) = 0, \;\; R(0,\rho_M) = 0, \;\; \partial_\rho R(\cdot,\rho) < 0, \;\; \partial_y R(y,\cdot) < 0 \; \text{ for } y \in (0,Y],\end{equation}
with 
$0<\rho_M < \infty$ being the local carrying capacity of the cell population. Here, the assumption on 
$\partial_\rho R$ corresponds to saturating growth, while the assumption on 
$\partial_y R$ models the fact that more mobile cells may be characterised by a lower proliferation rate due to the energetic cost of migration [Reference Aktipis, Boddy, Gatenby, Brown and Maley1, Reference Alfonso, Talkenberger, Seifert, Klink, Hawkins-Daarud, Swanson, Hatzikirou and Deutsch2, Reference Gallaher, Brown and Anderson22–Reference Giese, Loo, Tran, Haskett, Coons and Berens26, Reference Hatzikirou, Basanta, Simon, Schaller and Deutsch28, Reference Orlando, Gatenby and Brown34, Reference Pham, Chauviere, Hatzikirou, Li, Byrne, Cristini and Lowengrub39]. In particular, we will focus on the case where
\begin{equation}R(y,\rho) := r(y) - \rho \quad \text{with} \quad r(Y) = 0, \quad r(0) = \rho_M, \quad \dfrac{\textrm{d} r(y)}{\textrm{d} y} < 0 \; \text{ for } y \in (0,Y],\end{equation}
with r(y) being a smooth and bounded function that models the proliferation rate of cells in the phenotypic state y.
Finally, the second term on the right-hand side of the non-local PDE (2.1) models the effects of spontaneous, heritable phenotypic changes [Reference Huang29], which occur at rate 
$\beta >0$.
Object of study
Focussing on a biological scenario whereby cell movement occurs on a slower timescale compared to cell proliferation and death, while spontaneous, heritable phenotypic changes occur on a slower timescale compared to cell movement [Reference Doerfler and Böhm18, Reference Smith, Tomfohr, Wells, Beebe, Kepler and Reichert44, Reference Wang, Hinow, Bryce, Weaver, Estrada, Arteaga and Webb49], we introduce a small parameter 
$\varepsilon >0$ and let
\begin{equation*}\alpha := \varepsilon \quad \text{and} \quad \beta := \varepsilon^2.\end{equation*}
Furthermore, in order to explore the long-time behaviour of the cell population (i.e. the behaviour of the population over many cell generations), we use the time scaling 
$t \to t / \varepsilon$ in (2.1), which gives the following non-local PDE for the population density function 
$n(\frac{t}{\varepsilon},x,y) \equiv n_{\varepsilon}(t,x,y)$:
\begin{equation} \begin{cases}\displaystyle{\varepsilon \, \partial_t n_{\varepsilon} - \varepsilon \, \mu(y) \, \partial_x \left(n_{\varepsilon} \, \partial_x \rho_{\varepsilon}(t,x) \right) = R(y,\rho_{\varepsilon}(t,x)) \, n_{\varepsilon} + \varepsilon^2 \, \partial^2_{yy} n_{\varepsilon},}\\\\\displaystyle{\rho_{\varepsilon}(t,x) := \int_{0}^Y n_{\varepsilon}(t,x,y) \, \textrm{d}y,}\end{cases}\quad(x,y) \in \mathbb{R} \times (0,Y).\end{equation}
3 Numerical simulations
In this section, we report on numerical solutions of the non-local PDE (2.5) in the case where 
$R(y,\rho_{\varepsilon})$ is defined via (2.4). We choose the following initial condition:
\begin{equation}n_{\varepsilon}(0,x,y) = C \, e^{-x^2} \, e^{-\frac{(y-a)^2}{\varepsilon}} \quad \text{with } \; C \; \text{ s.t. } \; C \, \int_0^Y e^{-\frac{(y-a)^2}{\varepsilon}} \, \textrm{d}y = 1 \; \text{ and } \; a \in (0,Y),\end{equation}
which satisfies 
$\displaystyle{n_{\varepsilon}(0,x,y) \begin{matrix}{\scriptstyle\ast}\\[-8pt]-\!\!\!-\!\!\!-\!\!\!\!\!\rightharpoonup\\[-8pt]{\scriptstyle\varepsilon \rightarrow 0}\end{matrix} \rho(0,x) \, \delta_{\overline{y}^0(x)}(y)}$, with 
$\rho(0,x) = e^{-x^2}$ and 
$\overline{y}^0(x) \equiv a$. Such an initial condition models a biological scenario whereby 
$y=a$ is the locally dominant phenotypic trait at every position x at time 
$t=0$. We use uniform discretisations of steps 
$\Delta t$, 
$\Delta x$ and 
$\Delta y$ of the intervals (0,T], (0,X) and (0,Y), respectively, as computational domains of the independent variables t, x and y. The implicit finite volume scheme employed to solve numerically (2.5) complemented with (3.1) and subject to zero-flux/Neumann boundary conditions at 
$x=0$ (we expect a constant step), 
$y=0$ and 
$y=Y$ is described in Appendix A. All numerical computations are performed in Matlab.
Travelling fronts
The plots in Figure 1 summarise the numerical results obtained in the case where
\begin{equation}Y:=1, \quad \mu(y) := y^2 + 0.01, \quad r(y):= 1 - y^2, \quad \rho_M := 1.\end{equation}
The above definitions of 
$\mu(y)$ and r(y) are such that assumptions (2.2) and (2.4) are satisfied.
Figure 1. Travelling fronts. Plots of the normalised cell population density function 
$n_{\varepsilon}(t,x,y)/\rho_{\varepsilon}(t,x)$ (left panel) and the cell density 
$\rho_{\varepsilon}(t,x)$ (right panel, solid blue lines) at three successive time instants (i.e. 
$t=4$, 
$t=6$ and 
$t=8$). The dashed cyan lines in the right panel highlight the corresponding values of 
$r(\overline{y}_{{\varepsilon}}(t,x))$, with 
$\overline{y}_{\varepsilon}(t,x)$ being the maximum point of 
$n_{\varepsilon}(t,x,y)$ at 
$x \in \textrm{supp}(\rho_{\varepsilon})$, while the inset of the right panel displays the plots of 
$x_{1\varepsilon}(t)$ (blue squares), 
$x_{2\varepsilon}(t)$ (red diamonds) and 
$x_{3\varepsilon}(t)$ (black stars) such that 
$\rho_{\varepsilon}(t,x_{1\varepsilon}(t))=0.2$, 
$\rho_{\varepsilon}(t,x_{2\varepsilon}(t))=0.6$ and 
$\rho_{\varepsilon}(t,x_{3\varepsilon}(t))=0.8$. These results were obtained solving numerically (2.5) with 
$\varepsilon := 0.01$ under assumptions (2.4), (3.1) with 
$a=0.2$, and (3.2). Moreover, 
$T=8$, 
$X=25$, 
$\Delta {t} = 0.01$, 
$\Delta {x} = 0.01$ and 
$\Delta {y} = 0.02$.
The left panel of Figure 1 displays the plots of the normalised cell population density function 
$n_{\varepsilon}(t,x,y)/\rho_{\varepsilon}(t,x)$ at three successive time instants (i.e. 
$t=4$, 
$t=6$ and 
$t=8$). These plots indicate that for all 
$x \in \textrm{supp}(\rho_{\varepsilon})$ the normalised population density function 
$n_{\varepsilon}(t,x,y)/\rho_{\varepsilon}(t,x)$ is concentrated as a sharp Gaussian with maximum at a point 
$\overline{y}_{\varepsilon}(t,x)$ [i.e. 
$n_{\varepsilon}(t,x,y)/\rho_{\varepsilon}(t,x) \approx \delta_{\overline{y}_{\varepsilon}(t,x)}(y)$ for all 
$x \in \textrm{supp}(\rho_{\varepsilon})$], and the maximum point 
$\overline{y}_{\varepsilon}(t,x)$ behaves like a compactly supported and monotonically increasing travelling front that connects 
$y=0$ to 
$y=Y$.
The right panel of Figure 1 displays the plots of the cell density 
$\rho_{\varepsilon}(t,x)$ (solid blue lines) and the function 
$r(\overline{y}_{{\varepsilon}}(t,x))$ (dashed cyan lines) at three successive time instants (i.e. 
$t=4$, 
$t=6$ and 
$t=8$). These plots indicate that 
$\rho_{\varepsilon}(t,x)$ behaves like a one-sided compactly supported and monotonically decreasing travelling front that connects 
$\rho_M$ to 0. Moreover, there is an excellent quantitative match between 
$\rho_{\varepsilon}(t,x)$ and 
$r(\overline{y}_{{\varepsilon}}(t,x))$, which means that if 
$\rho_{\varepsilon}(t,x)>0$ then the relation 
$R(\overline{y}_{\varepsilon}(t,x),\rho_{\varepsilon}(t,x))=0$ holds.
The inset of the right panel of Figure 1 displays the plots of 
$x_{1\varepsilon}(t)$ (blue squares), 
$x_{2\varepsilon}(t)$ (red diamonds) and 
$x_{3\varepsilon}(t)$ (black stars) such that 
$\rho_{\varepsilon}(t,x_{1\varepsilon}(t))=0.2$, 
$\rho_{\varepsilon}(t,x_{2\varepsilon}(t))=0.6$ and 
$\rho_{\varepsilon}(t,x_{3\varepsilon}(t))=0.8$. These plots show that 
$x_{1\varepsilon}(t)$, 
$x_{2\varepsilon}(t)$ and 
$x_{3\varepsilon}(t)$ are straight lines of slope 
$\approx 2.5$, which supports the idea that 
$\rho_{\varepsilon}$ behaves like a travelling front of speed 
$c \approx 2.5$. Such a value of the speed is coherent with the condition on the minimal wave speed 
$c^*$ given by (4.20). In fact, inserting into (4.20) the numerical values of 
$\overline{y}_{\varepsilon}(8,x)$ in place of 
$\overline{y}(z)$ and the numerical values of 
$\partial^2_{yy} u_{\varepsilon}(8,x,\overline{y}_{\varepsilon}(8,x))$ with 
$u_{\varepsilon} = \varepsilon \log(n_{\varepsilon})$ in place of 
$\partial^2_{yy} u(z,\overline{y}(z))$ gives 
$c^* \gtrapprox 2.5$.
Front edge acceleration and stretching fronts
The plots in Figure 2 summarise the numerical results obtained in the case where
\begin{equation}Y:=20, \quad \mu(y) := 0.01 + y^4, \quad r(y):= 1 - \dfrac{y}{1+y}, \quad \rho_M := 1.\end{equation}
The above definitions of 
$\mu(y)$ and r(y) are chosen so that assumptions (2.2) and (2.4) are satisfied for 
$Y \to \infty$, and condition (4.21) is met (see details below).
Figure 2. Front edge acceleration and stretching fronts. Plots of the normalised cell population density function 
$n_{\varepsilon}(t,x,y)/\rho_{\varepsilon}(t,x)$ (left panel) and the cell density 
$\rho_{\varepsilon}(t,x)$ (right panel, solid blue lines) at three successive time instants (i.e. 
$t=4$, 
$t=6$ and 
$t=8$). The dashed cyan lines in the right panel highlight the corresponding values of 
$r(\overline{y}_{{\varepsilon}}(t,x))$, with 
$\overline{y}_{\varepsilon}(t,x)$ being the maximum point of 
$n_{\varepsilon}(t,x,y)$ at 
$x \in \textrm{supp}(\rho_{\varepsilon})$, while the inset of the right panel displays the plots of 
$x_{1\varepsilon}(t)$ (blue circles), 
$x_{2\varepsilon}(t)$ (red squares), 
$x_{3\varepsilon}(t)$ (black diamonds) and 
$x_{4\varepsilon}(t)$ (pink stars) such that 
$\rho_{\varepsilon}(t,x_{1\varepsilon}(t))=0.1$, 
$\rho_{\varepsilon}(t,x_{2\varepsilon}(t))=0.25$, 
$\rho_{\varepsilon}(t,x_{3\varepsilon}(t))=0.45$ and 
$\rho_{\varepsilon}(t,x_{4\varepsilon}(t))=0.8$. These results were obtained solving numerically (2.5) with 
$\varepsilon := 0.01$ under assumptions (2.4), (3.1) with 
$a=0.2$, and (3.3). Moreover, 
$T=8$, 
$X=200$, 
$\Delta {t} = 0.002$, 
$\Delta {x} = 0.1$ and 
$\Delta {y} = 0.05$.
The left panel of Figure 2 displays the plots of the normalised cell population density function 
$n_{\varepsilon}(t,x,y)/\rho_{\varepsilon}(t,x)$ at three successive time instants (i.e. 
$t=4$, 
$t=6$ and 
$t=8$). Similarly to the case of Figure 1, these plots show that for all 
$x \in \textrm{supp}(\rho_{\varepsilon})$ the normalised population density function 
$n_{\varepsilon}(t,x,y)/\rho_{\varepsilon}(t,x)$ is concentrated as a sharp Gaussian with maximum at a point 
$\overline{y}_{\varepsilon}(t,x)$ [i.e. 
$n_{\varepsilon}(t,x,y)/\rho_{\varepsilon}(t,x) \approx \delta_{\overline{y}_{\varepsilon}(t,x)}(y)$ for all 
$x \in \textrm{supp}(\rho_{\varepsilon})$], and the maximum point 
$\overline{y}_{\varepsilon}(t,x)$ is a monotonically increasing function of x with minimal value 0 for all 
$t \in [0,8]$. However, in contrast to the case of Figure 1, here 
$\overline{y}_{\varepsilon}(t,x)$ has a jump discontinuity and its maximal value increases as t increases.
The right panel of Figure 2 displays the plots of the cell density 
$\rho_{\varepsilon}(t,x)$ (solid blue lines) and the function 
$r(\overline{y}_{{\varepsilon}}(t,x))$ (dashed cyan lines) at three successive time instants (i.e. 
$t=4$, 
$t=6$ and 
$t=8$). Similarly to the case of Figure 1, these plots indicate that 
$\rho_{\varepsilon}(t,x)$ is a monotonically decreasing function of x with maximal value 
$\rho_M$ and minimal value 0 for all 
$t \in [0,8]$. Furthermore, there is an excellent quantitative match between 
$\rho_{\varepsilon}(t,x)$ and 
$r(\overline{y}_{{\varepsilon}}(t,x))$, which means that if 
$\rho_{\varepsilon}(t,x)>0$ then the relation 
$R(\overline{y}_{\varepsilon}(t,x),\rho_{\varepsilon}(t,x))=0$ holds. However, in contrast to the case of Figure 1, we have that 
$\rho_{\varepsilon}(t,x)$ behaves like a stretching front, which suggests that the speed of the front edge increases with t.
Coherently with this, the plot of 
$x_{1\varepsilon}(t)$ (blue circles) such that 
$\rho_{\varepsilon}(t,x_{1\varepsilon}(t))=0.1$ displayed in the inset of Figure 2 shows that the value of 
$x_{1\varepsilon}$ undergoes super linear growth, which supports the idea that front edge acceleration occurs. This is also coherent with the fact that, in the case where 
$\mu(y)$ and r(y) are defined via (3.3), we have that condition (4.21) is met and, therefore, the minimal wave speed 
$c^*$ tends to 
$\infty$ as 
$Y \to \infty$.
4 Formal asymptotic analysis
In this section, we undertake formal asymptotic analysis of the non-local PDE (2.5) in order to provide an explanation for the numerical results presented in Section 3.
Building on the Hamilton–Jacobi approach presented in [Reference Barles, Mirrahimi and Perthame6, Reference Diekmann, Jabin, Mischler and Perthame17, Reference Lorz, Mirrahimi and Perthame33, Reference Perthame36, Reference Perthame and Barles37], we make the real phase WKB ansatz [Reference Barles, Evans and Souganidis5, Reference Evans and Souganidis19, Reference Fleming and Souganidis21]
\begin{equation} n_{\varepsilon}(t,x,y) = e^{\frac{u_{\varepsilon}(t,x,y)}{\varepsilon}},\end{equation}
which gives
\begin{equation*}\partial_t n_{\varepsilon} = \frac{\partial_t u_{\varepsilon}}{\varepsilon} n_{\varepsilon}, \quad \partial_x n_{\varepsilon} = \frac{\partial_x u_{\varepsilon}}{\varepsilon} n_{\varepsilon}, \quad \partial^2_{yy} n_{\varepsilon} = \left(\frac{1}{\varepsilon^2} \left(\partial_y u_{\varepsilon} \right)^2 + \frac{1}{\varepsilon} \partial^2_{yy} u_{\varepsilon} \right) n_{\varepsilon}.\end{equation*}
Substituting the above expressions into the non-local PDE (2.5) gives the following Hamilton–Jacobi equation for 
$u_{\varepsilon}(t,x,y)$:
\begin{equation}\partial_t u_{\varepsilon} - \mu(y) \left(\partial_{x} u_{\varepsilon} \, \partial_{x} \rho_{\varepsilon} + \varepsilon \, \partial^2_{xx} \rho_{\varepsilon} \right) = R(y,\rho_{\varepsilon}) + \left(\partial_y u_{\varepsilon} \right)^2 + \varepsilon \, \partial^2_{yy} u_{\varepsilon}, \quad (x,y) \in \mathbb{R} \times(0,Y).\end{equation}
Letting 
$\varepsilon \to 0$ in (4.2) we formally obtain the following equation for the leading-order term u(t,x,y) of the asymptotic expansion for 
$u_{\varepsilon}(t,x,y)$:
\begin{equation}\partial_t u - \mu(y) \, \partial_{x} \rho \, \partial_{x} u = R(y,\rho) + \left(\partial_y u \right)^2, \quad (x,y) \in \mathbb{R} \times(0,Y),\end{equation}
where 
$\rho(t,x)$ is the leading-order term of the asymptotic expansion for 
$\rho_{\varepsilon}(t,x)$.
Constraint onu
Consider 
$x \in \mathbb{R}$ such that 
$\rho(t,x) >0$, that is, 
$x \in \textrm{supp}(\rho)$, and let 
$\overline{y}(t,x)$ be a nondegenerate maximum point of u(t,x,y), that is, 
${\overline{y}(t,x) \in {\textrm{arg}}\max\limits_{y \in [0,Y]} u(t,x,y)}$ with 
$\partial^2_{yy} u(t,x,\overline{y})<0$. Since 
$R(y,\rho_{\varepsilon})$ satisfies assumptions (2.3), we have that 
$\rho_{\varepsilon}(t,x) < \infty$ for all 
$\varepsilon > 0$. Hence, letting 
$\varepsilon \to 0$ in (4.1) formally gives the following constraint for all 
$t>0$:
\begin{equation}u(t,x,\overline{y}(t,x)) = \max_{y \in [0,Y]} u(t,x,y) = 0, \quad x \in \textrm{supp}(\rho),\end{equation}
which also implies that
\begin{equation}\partial_y u(t,x,\overline{y}(t,x)) = 0 \quad \text{and} \quad \partial_x u(t,x,\overline{y}(t,x)) = 0, \quad x \in \textrm{supp}(\rho).\end{equation}
Remark 4.1 The system defined by (4.3) and (4.4) is a constrained Hamilton–Jacobi equation and 
$\rho(t,x)>0$ can be regarded as a Lagrange multiplier associated with constraint (4.4).
Relation between 
$\overline{y}(t,x)$ and 
\[\rho (t,x)\]
Evaluating (4.3) at 
$y=\overline{y}(t,x)$ and using (4.4) and (4.5) we find
\begin{equation}R(\overline{y}(t,x),\rho(t,x)) = 0, \quad x \in \textrm{supp}(\rho).\end{equation}
The monotonicity assumptions (2.3) ensure that 
$\rho \mapsto R(\cdot,\rho)$ and 
$\overline{y} \mapsto R(\overline{y},\cdot)$ are both invertible. Therefore, relation (4.6) gives a one-to-one correspondence between 
$\overline{y}(t,x)$ and 
$\rho(t,x)$.
Transport equation for 
$\overline{y}$
Differentiating (4.3) with respect to y, evaluating the resulting equation at 
$y=\overline{y}(t,x)$ and using (4.4) and (4.5) yields
\begin{equation}\partial^2_{yt} u(t,x,\overline{y}) - \mu(\overline{y}) \, \partial_{x} \rho \, \partial^2_{yx} u(t,x,\overline{y}) = \partial_{y} R(\overline{y},\rho), \quad x \in \textrm{supp}(\rho).\end{equation}
Moreover, differentiating (4.5) with respect to t and x we find, respectively,
\begin{equation*}\partial^2_{ty} u(t,x,\overline{y}) + \partial^2_{yy} u(t,x,\overline{y}) \, \partial_{t} \overline{y}(t,x) = 0 \; \Rightarrow \; \partial^2_{yt} u(t,x,\overline{y}) = - \partial^2_{yy} u(t,x,\overline{y}) \, \partial_{t} \overline{y}(t,x)\end{equation*}
and
\begin{equation}\partial^2_{xy} u(t,x,\overline{y}) + \partial^2_{yy} u(t,x,\overline{y}) \, \partial_{x} \overline{y}(t,x) = 0 \; \Rightarrow \; \partial^2_{yx} u(t,x,\overline{y}) = - \partial^2_{yy} u(t,x,\overline{y}) \, \partial_{x} \overline{y}(t,x).\end{equation}
Substituting the above expressions of 
$\partial^2_{yt} u(t,x,\overline{y})$ and 
$\partial^2_{yx} u(t,x,\overline{y})$ into (4.7) and using the fact that 
$\partial^2_{yy} u(t,x,\overline{y}) < 0$ gives the following transport equation for 
$\overline{y}(t,x)$:
\begin{equation}\partial_{t} \overline{y} - \mu(\overline{y}) \, \partial_{x} \rho \, \partial_{x} \overline{y} = \frac{1}{-\partial^2_{yy} u(t,x,\overline{y})} \partial_{y} R(\overline{y},\rho), \quad x \in \textrm{supp}(\rho),\end{equation}
which is a generalised Burgers’ equation with source term since 
$\overline{y}(t,x)$ and 
$ \rho(t,x)$ are related through (4.6).
Travelling-wave problem
Substituting the travelling-wave ansatz
\begin{equation*}\rho(t,x) = \rho(z), \quad u(t,x,y) = u(z,y) \quad \text{and} \quad \overline{y}(t,x) = \overline{y}(z) \quad \text{with} \quad z = x - c \, t, \quad c>0\end{equation*}
into (4.3)–(4.6) and (4.9) gives
\begin{equation}- \left(c + \mu(y) \rho' \right) \partial_z u = R(y,\rho) + (\partial_y u)^2, \quad (z,y) \in \mathbb{R} \times (0,Y),\end{equation}
\begin{equation}u(z,\overline{y}(z)) = \max_{y \in [0,Y]} u(z,y) = 0, \quad \partial_y u(z,\overline{y}(z)) = 0, \quad \partial_z u(z,\overline{y}(z)) = 0, \quad z \in \textrm{supp}{\left(\rho \right)},\end{equation}
\begin{equation}R(\overline{y}(z),\rho(z)) = 0, \quad z \in \textrm{supp}{\left(\rho \right)},\end{equation}
\begin{equation}- \left(c + \mu(\overline{y}) \rho' \right) \overline{y}' = \frac{1}{-\partial^2_{yy} u(z,\overline{y})} \partial_{y} R(\overline{y},\rho), \quad z \in \textrm{supp}{\left(\rho \right)}.\end{equation}
We consider travelling-front solutions 
$\overline{y}(z)$ that satisfy (4.13) subject to the following asymptotic condition
\begin{equation}\lim_{z \to - \infty} \overline{y}(z) =0,\end{equation}
so that, since 
$R(0,\rho_M)=0$ [cf. assumptions (2.3)], relation (4.12) gives 
$\displaystyle{\lim_{z \to - \infty} \rho(z) = \rho_M}$.
Monotonicity of travelling-front solutions
Differentiating (4.12) with respect to z gives
\begin{equation}\partial_y R(\overline{y}(z),\rho(z)) \overline{y}'(z) + \partial_{\rho} R(\overline{y}(z),\rho(z)) \rho'(z) = 0, \quad z \in \textrm{supp}{\left(\rho \right)}.\end{equation}
Substituting the expression of 
$\rho'$ given by (4.15) into (4.13) yields
\begin{equation*}-c \ \overline{y}' + \mu(\overline{y}) \ \dfrac{\partial_y R(\overline{y},\rho)}{\partial_{\rho} R(\overline{y},\rho)} \ \left(\overline{y}' \right)^2 = \frac{1}{-\partial^2_{yy} u(z,\overline{y})} \, \partial_{y} R(\overline{y},\rho),\end{equation*}
that is,
\begin{equation}\overline{y}' = \frac{-\partial_y R(\overline{y},\rho)}{c} \left(\frac{1}{-\partial^2_{yy} u(z,\overline{y})} + \dfrac{\mu(\overline{y}) \left(\overline{y}' \right)^2}{-\partial_{\rho} R(\overline{y},\rho)} \right), \quad z \in \textrm{supp}{\left(\rho \right)}.\end{equation}
Since 
$\partial^2_{yy} u(z,\overline{y})<0$ and 
$\partial_y R(y,\cdot)<0$ for 
$y \in (0,Y]$ [cf. assumptions (2.3)], using (4.16) and the expression of 
$\rho'$ given by (4.15) we find
\begin{equation}\overline{y}'(z) > 0 \quad \text{and} \quad \rho'(z) < 0, \quad z \in \textrm{supp}{\left(\rho \right)}.\end{equation}
Position of the front edge
Relation (4.12) and monotonicity results (4.17) along with the fact that 
$R(Y,0)=0$ [cf. assumptions (2.3)] imply that the position of the edge of a travelling-front solution 
$\overline{y}(z)$ that satisfies (4.13) subject to asymptotic condition (4.14) coincides with the unique point 
$\ell \in \mathbb{R}$ such that 
$\overline{y}(\ell)=Y$.
Minimal wave speed
Differentiating both sides of (4.10) with respect to y gives
\begin{equation*}- \left(c + \mu(y) \rho' \right) \partial^2_{yz} u(z,y) - \dfrac{\textrm{d} \mu(y)}{\textrm{d} y} \ \rho' \ \partial_{z} u(z,y) = \partial_y R(y,\rho) + 2 \ \partial_y u(z,y) \ \partial^2_{yy} u(z,y).\end{equation*}
Evaluating the above equation at 
$y=\overline{y}(z)$ using (4.11) yields
\begin{equation}\left(c + \mu(\overline{y}) \rho'(z) \right) \, \partial^2_{yz} u(z,\overline{y}) + \partial_y R(\overline{y},\rho) = 0.\end{equation}
Moreover, (4.8) implies that
\begin{equation*}\partial^2_{yz} u(z,\overline{y}) = - \partial^2_{yy} u(z,\overline{y}) \, \overline{y}'\end{equation*}
and substituting into the latter equation the expression of 
$\overline{y}'$ given by (4.15) we find
\begin{equation*}\partial^2_{yz} u(z,\overline{y}) = \partial^2_{yy} u(z,\overline{y}) \ \dfrac{\partial_{\rho} R(\overline{y},\rho)}{\partial_y R(\overline{y},\rho)} \ \rho'(z).\end{equation*}
Inserting the above expression of 
$\partial^2_{yz} u(z,\overline{y})$ into (4.18) gives
\begin{equation*}\mu(\overline{y}) \ \partial^2_{yy} u(z,\overline{y}) \ \partial_{\rho} R(\overline{y},\rho) \left(\rho'\right)^2 + c \, \partial^2_{yy} u(z,\overline{y}) \ \partial_{\rho} R(\overline{y},\rho) \ \rho' + \left(\partial_y R(\overline{y},\rho)\right)^2 = 0.\end{equation*}
In the case where 
$R(y,\rho)$ is defined via (2.4), we have that
\begin{equation*}\partial_{\rho} R(\cdot,\rho) = -1 \quad \text{and} \quad \displaystyle{\partial_y R(\overline{y},\cdot) = \dfrac{\textrm{d} r(\overline{y})}{\textrm{d} y}}.\end{equation*}
Hence, the latter equation becomes
\begin{equation}\mu(\overline{y}) \ \partial^2_{yy} u(z,\overline{y}) \ \left(\rho'\right)^2 + c \, \partial^2_{yy} u(z,\overline{y}) \ \rho' - \left(\dfrac{\textrm{d} r(\overline{y})}{\textrm{d} y}\right)^2 = 0.\end{equation}
Coherently with (4.17), the real roots of (4.19) seen as a quadratic equation for 
$\rho'$ are negative. Furthermore, the following condition has to hold for the roots to be real:
\begin{equation*}2 \ \left|\dfrac{\textrm{d} r(\overline{y})}{\textrm{d} y}\right| \sqrt{\dfrac{\mu(\overline{y})}{\left|\partial^2_{yy} u(z,\overline{y})\right|}} \leqslant c.\end{equation*}
This indicates that there is a minimal wave speed 
$c^*$, which satisfies the following condition:
\begin{equation}c^* \geqslant \sup_{z \in \textrm{supp}(r(\overline{y}))} 2 \ \left|\dfrac{\textrm{d} r(\overline{y}(z))}{\textrm{d} y}\right| \sqrt{\dfrac{\mu(\overline{y}(z))}{\left|\partial^2_{yy} u(z,\overline{y}(z))\right|}},\end{equation}
where we have used the fact that, when 
$R(y,\rho)$ is defined via (2.4), relation (4.12) gives
\begin{equation*}\rho(z) \equiv r(\overline{y}(z)), \quad z \in \textrm{supp}{\left(\rho \right)}.\end{equation*}
Condition (4.20) implies that if
\begin{equation}\left|\dfrac{\textrm{d} r(Y)}{\textrm{d} y}\right| \sqrt{\mu(Y)} \longrightarrow \infty \; \text{ as } \; Y \to \infty\end{equation}
then 
$c^* \to \infty$ as 
$Y \to \infty$.
5 Discussion, biological implications and research perspectives
Discussion of the main results
In this paper, we have reported on the results of numerical simulations of the non-local PDE (2.5) complemented with (2.2) and (2.4), and subject to zero Neumann boundary conditions at 
$y=0$ and 
$y=Y$. These numerical results indicate that
\begin{equation}\text{if } n_{\varepsilon}(0,x,y) \begin{matrix}{\scriptstyle\ast}\\[-8pt]-\!\!\!-\!\!\!-\!\!\!\!\!\rightharpoonup\\[-8pt]{\scriptstyle\varepsilon \rightarrow 0}\end{matrix} \rho(0,x) \, \delta_{\overline{y}^0(x)}(y) \; \text{ then } \; n_{\varepsilon}(t,x,y) \begin{matrix}{\scriptstyle\ast}\\[-8pt]-\!\!\!-\!\!\!-\!\!\!\!\!\rightharpoonup\\[-5pt]{\scriptstyle\varepsilon \rightarrow 0}\end{matrix} \rho(t,x) \, \delta_{\overline{y}(t,x)}(y),\end{equation}
with 
$\rho(t,x)$ and 
$\overline{y}(t,x)$ such that if 
$\rho(t,x)>0$ then the relation 
$R(\overline{y}(t,x), \rho(t,x)) =0$ holds. These numerical results also indicate that in the case where 
$Y \in \mathbb{R}^*_+$ (i.e. when 
$\mu(Y)<\infty$ and, therefore, the cell mobility is bounded), 
$\rho(t,x)$ in (5.1) behaves like a one-sided compactly supported and monotonically decreasing travelling front 
$\rho(z) \equiv \rho(x-ct)$ that connects 
$\rho_M$ to 0, while 
$\overline{y}(t,x)$ in (5.1) behaves like a compactly supported and monotonically increasing travelling front 
$\overline{y}(z) \equiv \overline{y}(x-ct)$ that connects 0 to Y. Furthermore, we have provided numerical evidence for the fact that front edge acceleration and formation of stretching fronts may occur in the case where 
$Y \to \infty$ (i.e. when 
$\mu(Y) \to \infty$ and, therefore, the cell mobility is unbounded).
In order to explain such numerical results, we have undertaken formal asymptotic analysis of the non-local PDE (2.5) complemented with (2.2) and (2.3) in the asymptotic regime 
$\varepsilon \to 0$ using a Hamilton–Jacobi approach. In particular, we have shown that 
$\overline{y}(t,x)$ satisfies a generalised Burgers’ equation with source term [see transport equation (4.9)] and 
$\rho(t,x)=R(\overline{y}(t,x),\rho(t,x))^{-1}(0)$ [see relation (4.6)]. Moreover, we have shown that travelling-front solutions 
$\overline{y}(z)$ of such transport equation, which connect 0 to Y are monotonically increasing, while the corresponding 
$\rho(z) = R(\overline{y}(z),\rho(z))^{-1}(0)$ is monotonically decreasing and connect 
$\rho_M$ to 0 [see the monotonicity results given by (4.17)]. Finally, in the case where 
$R(y,\rho)$ is defined via (2.4), we have characterised the minimal speed 
$c^*$ of such travelling-front solutions [see the result given by (4.20)] and derived sufficient conditions under which 
$c^* \to \infty$ as 
$Y \to \infty$ [see condition (4.21)].
Biological implications of the main results
From a biological point of view, 
$\overline{y}(t,x)$ represents the dominant phenotypic trait at position x and time t and the transport equation for 
$\overline{y}(t,x)$ can be seen as a generalised canonical equation of adaptive dynamics [Reference Dieckmann and Law16, Reference Diekmann, Jabin, Mischler and Perthame17], which describes the spatio-temporal evolution of the dominant phenotypic trait. Furthermore, the fact that 
$\rho(t,x)$ behaves like a monotonically decreasing travelling front 
$\rho(z)$ that connects 
$\rho_M$ to 0 represents the formation of an invasion front of cells that expands into the surrounding environment [Reference Pérez-Garca, Calvo, Belmonte-Beitia, Diego and Pérez-Romasanta35]. Hence, the fact that 
$\overline{y}(t,x)$ behaves like a monotonically increasing travelling front 
$\overline{y}(z)$ that connects 0 to Y has the following biological implications. First, the fact that the front 
$\overline{y}(z)$ is monotonic indicates that cells with different phenotypic characteristics populate different parts of the invasion front – i.e. phenotypic heterogeneity is dynamically maintained throughout the front. Second, since larger values of y correlate with a lower proliferation rate and a higher mobility, the fact that the front 
$\overline{y}(z)$ is increasing indicates that more mobile/less proliferative phenotypic variants occupy the front edge, whereas less mobile/more proliferative phenotypic variants are selected at the back of the front. This recapitulates previous theoretical and experimental results on glioma growth, which indicate that the interior of the tumour consists mainly of proliferative cells while the tumour border comprises mainly cells that are more mobile and less proliferative – see, for instance, [Reference Alfonso, Talkenberger, Seifert, Klink, Hawkins-Daarud, Swanson, Hatzikirou and Deutsch2, Reference Dhruv, McDonough Winslow, Armstrong, Tuncali, Eschbacher, Kislin, Loftus, Tran and Berens15, Reference Giese, Bjerkvig, Berens and Westphal25, Reference Giese, Loo, Tran, Haskett, Coons and Berens26, Reference Wang, Rath, Lal, Richard, Li, Goodwin, Laterra and Xia48, Reference Xie, Mittal and Berens50] and references therein.
Research perspectives
Building upon the results presented in this paper, a number of generalisations of the mathematical model given by the non-local PDE (2.1) could be considered in order to investigate the role of the concerted action between evolutionary and mechanical processes in tissue development and tumour growth. For example, a natural generalisation is the one given by the following non-local PDE
\begin{equation} \begin{cases}\displaystyle{\partial_t n - \mu(y) \, \nabla_{\textbf{x}} \cdot \left(n \, \nabla_{\textbf{x}} P(t,\textbf{x}) \right) = R(y,P(t,\textbf{x})) \, n + \beta \, \partial^2_{yy} n,}\\\\\displaystyle{P \equiv \Pi(\rho), \quad \rho(t,\textbf{x}) := \int_{0}^Y n(t,\textbf{x},y) \, \textrm{d}y,}\end{cases}\;(\textbf{x},y) \in \mathbb{R}^{d} \times (0,Y)\end{equation}
subject to zero Neumann boundary conditions at 
$y=0$ and 
$y=Y$. Here, 
$d=1,2,3$ depending on the biological problem considered, and the function 
$P(t,\textbf{x})$ is the pressure exerted by cells at position x and time t, which is defined via the barotropic relation 
$\Pi(\rho)$ that satisfies suitable assumptions.
On the basis of the knowledge we have here acquired on the behaviour of the solutions to the non-local PDE (2.1), under asymptotic scenarios relevant to applications we may expect 
$n(t,\textbf{x}, y)$ to converge to a singular measure of the form 
$\rho(t,\textbf{x}) \delta_{\overline{y}(t,\textbf{x})}(y)$. Moreover, depending on the choices of Y, 
$\mu(y)$, 
$R(y,P(t,\textbf{x}))$ and 
$\Pi(\rho)$, the cell density 
$\rho(t,\textbf{x})$ may develop into an invading front or it may exhibit interface instabilities [Reference Kim and Tong30, Reference Lorenzi, Lorz and Perthame32, Reference Pham, Turian, Liu, Li and Lowengrub40, Reference Tang, Vauchelet, Cheddadi, Vignon-Clementel, Drasdo and Perthame45]. Finally, when the following definition of 
$\Pi$ is considered
\begin{equation*}\Pi(\rho) := K_{\gamma} \, \rho^{\gamma}, \qquad K_{\gamma} >0, \; \gamma > 1,\end{equation*}
which was proposed in [Reference Perthame, Quirós and Vázquez38] in order to capture key aspects of tumour and tissue growth while ensuring analytical tractability of the model equation, one finds that 
$P(t,\textbf{x})$ satisfies a porous medium-type equation. Hence, free boundary problems may emerge in the asymptotic regime 
$\gamma \to \infty$ (i.e. the asymptotic regime whereby cells are regarded as an incompressible fluid). These are lines of research that we will be pursuing in the near future.
Acknowledgements
B.P. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 740623). T.L. gratefully acknowledges support of the project PICS-CNRS no. 07688 and the MIUR grant “Dipartimenti di Eccellenza 2018-2022”, and would like to thank Alexander Lorz for insightful discussions during the early stages of the project.
Appendix A. Numerical methods
Since 
$\rho_{\varepsilon}(t,x)$ might develop into a stiff travelling front, solving the non-local PDE (2.5) via an explicit finite volume scheme would result in a severe CFL constraint on t. In order to overcome such a limitation, we carried out numerical simulations using the implicit finite volume scheme presented here. For simplicity of notation, throughout this appendix we drop the subscript 
$\varepsilon$.
Time splitting
Adopting a time-splitting approach, which is based on the idea of decomposing the original problem into simpler subproblems that are then sequentially solved at each time step, we decompose the non-local PDE (2.5) posed on 
$\Omega := (0,T]\times(0, X)\times(0,Y)$ into two parts – viz. the diffusion–advection part corresponding to the following non-local PDE
\begin{equation}\begin{cases}\partial_t n -{\partial}_x(\mu(y) \, n \, {\partial}_x\rho)=\varepsilon \, {\partial}_{yy}^2 n,\\ \\[-7pt] \displaystyle{\rho(t,x) := \int_0^Y n(t,x,y) \, \textrm{d}y}\end{cases}\end{equation}
and the reaction part corresponding to the following integro-differential equation
\begin{equation}\begin{cases}\varepsilon \, \partial_t n = R(y,\rho) \, n,\\ \\[-7pt] \displaystyle{\rho(t,x) := \int_0^Y n(t,x,y) \, \textrm{d}y}.\end{cases}\end{equation}
We complement (A1) with zero-flux/Neumann boundary conditions at 
$x=0$ (we expect a constant step), 
$y=0$ and 
$y=Y$. Note that making the ansatz 
$n(t,x,y) = e^{\frac{u(t,x,y)}{\varepsilon}}$, as similarly done in Section 4, the integro-differential equation (A2) can be rewritten in the following alternative form:
\begin{equation}\begin{cases}{\partial}_t u = R(y, \rho),\\ \\[-7pt] \displaystyle{\rho(t,x) := \int_0^Y e^{\frac{u(t,x,y)}{\varepsilon}} \, \textrm{d}y}.\end{cases}\end{equation}
Preliminaries and notation
We denote by 
$[\![ k_1, k_2]\!]$ the set of integers between 
$k_1$ and 
$k_2$. We discretise 
$\Omega$ via a uniform structured grid of steps 
$\Delta {t}$, 
$\Delta {x}$, 
$\Delta {y}$ whereby 
$t_h=h\Delta {t}$ and the (j, k)-th cell is
\begin{equation}K_{j-\frac12,k-\frac12} = (x_{j-1},x_{j})\times (y_{k-1},y_{k})\quad \text{with}\quad x_j = j\Delta {x}, \quad y_k = k \Delta {y},\end{equation}
where 
$j \in [\![1, m_x]\!]$ and 
$k\in [\![1, m_y]\!]$, 
$\Delta {x} = \frac{X}{m_x}$, 
$\Delta {y} = \frac{Y}{m_y}$ and 
$m_x,m_y\in \mathbb{N}$. Moreover, we let 
$N_{j-\frac12,k-\frac12}^h$ be the numerical approximation of the average of 
$n(t_h, x, y)$ over the cell 
$K_{j-\frac12,k-\frac12}$ and we consider the following first-order approximation of the average of 
$\rho(t_h, x)$ over the interval 
$(x_{j-1},x_{j})$:
\begin{equation*}\rho_{j-\frac12}^{h} = \Delta {y} \sum_{k=1}^{m_y} N_{j-\frac12,k-\frac12}^{h}.\end{equation*}
Finally, we introduce the notation
\begin{equation}\boldsymbol{n}^h = \left(N_{j-\frac12,k-\frac12}^h\right)^\textrm{T}\in \mathbb{R}^{(m_x+1), (m_y+1)}, \quad \boldsymbol{\rho}^{h} = \left(\rho_{j-\frac12}^h\right)^\textrm{T} \in \mathbb{R}^{m_x}\end{equation}
with 
$j \in[\![ 1, m_x]\!]$ and 
$k \in[\![ 1, m_y]\!]$.
Numerical scheme
Step 1 We first solve numerically (A1) by using the following implicit scheme:
\begin{equation}\frac{N_{j-\frac12,k-\frac12}^{*} - N_{j-\frac12,k-\frac12}^h}{\Delta {t}} - \frac{1}{\Delta {x}} \left[F_{j,k-\frac12}^{*}-F_{j-1,k-\frac12}^{*}\right]=\varepsilon \frac{N_{j-\frac12,k+\frac12}^{*} - 2 N_{j-\frac12,k-\frac12}^{*} + N_{j-\frac12,k-\frac32}^{*}}{\Delta {y}^2},\end{equation}
where 
$F_{j, k-\frac12}^{*}$ represents the numerical flux at the boundary 
${\partial} K_{j-\frac12,k-\frac12} \cap \{x=x_{j}\}$, which is given by the following upwind approximation:
\begin{equation}F_{j, k-\frac12}^{*} = \mu_{k-\frac12}\left[-(\delta_x\rho_{j}^{*})_{-} \, N_{j-\frac12,k-\frac12}^{*}+(\delta_x\rho_{j}^{*})_{+} \, N_{j+\frac12,k-\frac12}^{*}\right].\end{equation}
Here, 
$\mu_{k-\frac12} = \mu(y_{k-\frac12})$,
\begin{equation*}\delta_x\rho_j^{*}=\dfrac{\left(\rho_{j+\frac12}^{*} - \rho_{j-\frac12}^{*}\right)}{\Delta {x}} \quad \text{with} \quad \rho_{j-\frac12}^{*} = \Delta {y} \sum_{k=1}^{m_y} N_{j-\frac12,k-\frac12}^{*},\end{equation*}
and 
$(\cdot)_-$ and 
$(\cdot)_+$ are, respectively, the negative and positive parts of 
$(\cdot)$. Analogous considerations hold for 
$F_{j-1,k-\frac12}^{*}$. We complement (A6) with boundary conditions corresponding to zero-flux/Neumann boundary conditions at 
$x=0$ (we expect a constant step), 
$y=0$ and 
$y=Y$.
Step 2 We solve numerically (A3) using the following implicit scheme:
\begin{equation}\frac{U_{j-\frac12,k-\frac12}^{h+1} - U_{j-\frac12,k-\frac12}^*}{\Delta {t}} = R\left(y_{k-\frac12}, \rho_{j-\frac12}^{h+1}\right),\end{equation}
where 
$U_{j-\frac12,k-\frac12}^* = \varepsilon \ln \left(N_{j-\frac12,k-\frac12}^{*} \right)$ and 
$N_{j-\frac12,k-\frac12}^{*}$ is obtained via (A6). Since
\begin{eqnarray}\rho_{j-\frac12}^{h+1} &=& \Delta {y} \sum_{k=1}^{m_y} \exp\left( \dfrac{U_{j-\frac12,k-\frac12}^{h+1}}{\varepsilon}\right) \nonumber\\&=& \Delta {y} \sum_{k=1}^{m_y} \exp\left( \dfrac{U_{j-\frac12,k-\frac12}^{*}+\Delta {t} R\left(y_{k-\frac12}, \rho_{j-\frac12}^{h+1}\right)}{\varepsilon}\right),\end{eqnarray}
in the case where the function R is defined via (2.4) the value of 
$\rho_{j-\frac12}^{h+1}$ can be found by solving (A9). The value of 
$\rho_{j-\frac12}^{h+1}$ so obtained is substituted into (A8), which is then solved in order to find 
$U_{j-\frac12,k-\frac12}^{h+1}$, whose value is finally used to compute 
$N_{j-\frac12,k-\frac12}^{h+1}$ via the formula
\begin{equation*}N_{j-\frac12,k-\frac12}^{h+1} = \exp\left(\frac{U_{j-\frac12,k-\frac12}^{h+1}}{\varepsilon} \right).\end{equation*}
Properties of the numerical scheme (A6)
Due to the the strong coupling between n(t,x,y) and 
$\rho(t,x)$ in the non-local PDE (A1), it remains an open problem to prove existence and uniqueness of the solution to the corresponding initial boundary value problem. Similarly, proving unique solvability of the non-linear, non-local, implicit scheme (A6) remains an open problem.
Here, assuming solvability of (A6), we prove that such a numerical scheme preserves non-negativity of n, maximum principle on 
$\rho$ and monotonicity of 
$\rho$ (cf. Proposition A.1).
Proposition A.1 Consider the scheme (A6) only. If the numerical scheme (A6) is uniquely solvable, then the following properties hold:
-
(i) [nonnegativity]
if
$\boldsymbol{n}^h\geqslant 0$ then $\boldsymbol{n}^*\geqslant 0$; -
(ii) [maximum principle on
$\rho$]if
$0\leqslant \boldsymbol{\rho}^h \leqslant \rho_M$ then $0\leqslant \boldsymbol{\rho}^* \leqslant \rho_M$; -
(iii) [monotonicity of
$\rho$]if
$\boldsymbol{\rho}^h$ is monotonically decreasing then $\boldsymbol{\rho}^*$ is monotonically decreasing.
Proof.
-
(i) The implicit scheme (A6) can be rewritten as
(A10)\begin{equation}-a_{j-1,k}^{*} N_{j-\frac32,k-\frac12}^{*}+b_{j,k}^{*} N_{j-\frac12,k-\frac12}^{*}-c_{j+1,k}^{*} N_{j+\frac12,k-\frac12}^{*} +\end{equation}
(A11)where\begin{equation}\varepsilon\frac{\Delta {t}}{(\Delta {y})^2} \big(-N_{j-\frac12,k-\frac32}^{*} + 2N_{j-\frac12,k-\frac12}^{*} - N_{j-\frac12,k+\frac12}^{*} \big)=N_{j-\frac12,k-\frac12}^h,\end{equation}
The system of equations (A10) can be written in matrix form as\begin{align*}&a_{j,k}^{*} = \frac{\Delta {t}}{\Delta {x}}\mu_{k-\frac12}(\delta_{x}\rho_{j}^{*})_{-} \geqslant 0,\quad c_{j,k}^{*} = \frac{\Delta {t}}{\Delta {x}}\mu_{k-\frac12}(\delta_x\rho_{j-1}^{*})_{+} \geqslant 0,\\[3pt] &b_{j,k}^{*} = 1+\frac{\Delta {t}}{\Delta {x}}\mu_{k-\frac12}\left[(\delta_x\rho_{j-1}^{*})_{+}+(\delta_x\rho_{j}^{*})_{-}\right] = 1 + a_{j,k}^{*} + c_{j,k}^{*}.\end{align*}
where\begin{equation*}\boldsymbol{M}^{*} \boldsymbol{n}^{*} = \boldsymbol{n}^h, \nonumber\end{equation*}
$\boldsymbol{M}^{*}$ is a matrix containing the terms $a_{j,k}^{*}$’s, $b_{j,k}^{*}$’s and $c_{j,k}^{*}$’s with $j \in [\![1, m_x]\!]$, $k\in [\![1, m_y]\!]$. Since the matrix $\boldsymbol{M}^{*}$ is strictly diagonally dominant by columns, it is invertible and all elements of $\left(\boldsymbol{M}^{*}\right)^{-1}$ are positive. This ensures that $\boldsymbol{n}^{*}$ is non-negative if $\boldsymbol{n}^h$ is non-negative.
-
(ii) Summing (A6) over all
$k\in [\![ 1, m_y ]\!]$, we find
(A12)where\begin{equation}\frac{\rho_{j-\frac12}^{*} - \rho_{j-\frac12}^h}{\Delta {t}} - \frac{1}{\Delta {x}} \left[<\mu_{k-\frac12} N_{j,k-\frac12}^ {*,\rm{upwind}}>\delta_x \rho_{j}^ {*}-<\mu_{k-\frac12} N_{j-1,k-\frac12}^ {*,\rm{upwind}}>\delta_x \rho_{j-1}^ {*}\right]=0,\end{equation}
$\delta_x \rho_{j}^ {*} = \big(\rho_{j+\frac12}^ {*}-\rho_{j-\frac12}^ {*}\big) / \Delta {x}$, $<\mu_{k-\frac12} N_{j,k-\frac12}^{*,\rm{upwind}} >= \Delta {y} \sum_{k=1}^{m_y}\mu_{k-\frac12} N_{j,k-\frac12}^{*,\rm{upwind}}$ and
(A13)For simplicity of notation, we define\begin{equation}N_{j,k-\frac12}^{*,\rm{upwind}} =\begin{cases}N_{j-\frac12,k-\frac12}^{*} \quad&\text{ if } \delta_x\rho_j^{*} < 0, \\\\[-7pt] N_{j+\frac12,k-\frac12}^{*} \quad&\text{ if } \delta_x\rho_{j}^{*} \geqslant 0.\end{cases}\end{equation}
$\displaystyle{d_j^{*} = \frac{\Delta {t}}{\Delta {x}^2}<\mu_{k-\frac12} N_{j,k-\frac12}^ {*,\rm{upwind}}>}.$ Notice that $d_j^{*}\geqslant 0$. Then, the system of equations (A12) can be rewritten as
(A14)Assume that\begin{equation}(1 + d_{j-1}^{*} + d_{j}^ {*}) \rho_{j-\frac12}^ {*} - d_{j-1}^{*} \rho_{j-\frac32}^ {*}-d_{j}^{*} \rho_{j+\frac12}^ {*}=\rho_{j-\frac12}^ h.\end{equation}
$\displaystyle{\rho_{j_0-\frac12}^ {*} = \min_j \{\rho_{j-\frac12}^{*}\}}$, we claim that $\rho_{j_0-\frac12}^{*} \geqslant 0$. In fact, if $\rho_{j_0-\frac12}^{*} < 0$, we have
(A15)which is a contradiction. Hence,\begin{equation}\rho_{j_0-\frac12}^h = \rho_{j_0-\frac12}^{*} + d_{j_0-1}^{*}\left(\rho_{j_0-\frac12}^{*} - \rho_{j_0-\frac32}^{*}\right) + d_{j_0}^{*}\left(\rho_{j_0-\frac12}^{*} - \rho_{j_0+\frac12}^{*}\right) \leqslant \rho_{j_0-\frac12}^{*}<0,\end{equation}
$\boldsymbol{\rho}^{*}\geqslant 0$. Similarly, one can prove that $\boldsymbol{\rho}^{*} \leqslant \rho_M$.
-
(iii) Introducing the notation
$w_j^{*} = \rho_{j+\frac12}^{*} - \rho_{j-\frac12}^{*}$, we rewrite (A14) as
(A16)Changing all subscripts in (A16) from j to\begin{equation}\rho_{j-\frac12}^{*} + d_{j-1}^{*} w_{j-1}^{*} - d_{j}^{*} w_{j}^{*} = \rho_{j-\frac12}^{h}.\end{equation}
$j+1$, after a little algebra we find
(A17)Writing the system of the equations (A17) in matrix form and using arguments similar to those used in part (i), it is possible to prove that\begin{equation}(1+2d_j^{*})w_j^{*} - d_{j-1}^{*} w_{j-1}^{*} - d_{j+1}^{*} w_{j+1}^{*} = w_j^h.\end{equation}
$\boldsymbol{w}^{*}\leqslant 0$ if $\boldsymbol{w}^{h}\leqslant 0$.
Properties of the numerical scheme (A8)
The numerical scheme (A8) satisfies the properties established by Proposition A.2.
Proposition A.2 Consider the scheme (A8) only. If 
$R(y,\rho)$ satisfies assumptions (2.3) and 
$\boldsymbol{n}^*\geqslant 0$, then the following properties hold:
-
(i) [existence and uniqueness and non-negativity]
the scheme (A8) admits a unique solution such that
$\boldsymbol{n}^{h+1}\geqslant 0$; -
(ii) [maximum principle on
$\rho$]if
$0\leqslant \boldsymbol{\rho}^* \leqslant \rho_M$, then $0 \leqslant \boldsymbol{\rho}^{h+1} \leqslant \rho_M$.
Proof.
-
(i) It is sufficient to prove existence and uniqueness of
$\rho_{j-\frac12}^{h+1}$. Let
Since\begin{equation*}f(\rho) = \rho - \Delta {y} \sum_{k=1}^{m_y} \exp\left(\frac{U_{j-\frac12,k-\frac12}^{*}+\Delta {t} R(y_{k-\frac12}, \rho)}{\varepsilon}\right).\end{equation*}
$f'(\rho) >0$, $f(0)<0$ and $\displaystyle{\lim_{\rho\to\infty} f(\rho) =\infty}$, equation (A9) has a unique positive root, which is $\rho_{j-\frac12}^{h+1}$. From this, existence, uniqueness and non-negativity of $N_{j-\frac12,k-\frac12}^{h+1}$ immediately follow.
-
(ii) Noticing that
$f'(\rho)>0$, $f(0) <0$ and
(A18)we conclude that equation (A9) has a unique solution in the interval\begin{equation}f(\rho_M) \geqslant \rho_M - \Delta {y} \sum_{k=1}^{m_y} \exp\left(\frac{U_{j-\frac12,k-\frac12}^{*}}{\varepsilon}\right) = \rho_M -\rho_{j-\frac12}^* \geqslant 0,\end{equation}
$[0, \rho_M]$. This implies that $0\leqslant \boldsymbol{\rho}^{h+1} \leqslant \rho_M$.
