2.1 Description of cell motion

Velocity-jump process

We describe the motion of a cell labelled by an index i as a run-and-tumble process with run time $\tau_i \in \mathbb{R}^*_+$ and running probability $\psi(\mathbf{x}_i,\tau_i)$ , where $0< \psi(\cdot,\cdot) \leqslant 1$ and $\partial_{\tau_i}\psi(\cdot,\cdot) \leqslant 0$ . The running probability $\psi(\mathbf{x}_i,\tau_i)$ correlates with the stopping rate $\beta(\mathbf{x}_i,\tau_i)$ through the relations given by the following definition [Reference Estrada-Rodriguez and Gimperlein11]:

(2.1) \begin{equation}\psi(\mathbf{x}_i,\tau_i)\;:=\; \exp\left(\int_0^{\tau_i}\beta(\mathbf{x}_i,s)\mathop{}\!\mathrm{d} s \right) \ , \quad \beta = \frac{\varphi}{\psi} \quad \text{with} \quad \varphi \;:=\; -\partial_{\tau_i}\psi \ .\end{equation}

Hence, starting at position $\mathbf{x}_i$ at time t, the $i{\rm th}$ cell will continue moving along a straight path in the direction given by the vector $\mathbf{v}_i$ with constant speed c for a period of time $\tau_i$ , after which it may stop with rate $\beta(\mathbf{x}_i,\tau_i)$ . The cell will then instantaneously resume moving in a new randomly selected direction given by a vector $\bar{\mathbf{v}}_i$ , which is prescribed by a turning kernel $\ell(\mathbf{x}_i,t,\mathbf{v}_i;\bar{\mathbf{v}}_i)$ , that is, cells undergo a velocity-jump process [Reference Othmer, Dunbar and Alt22].

Running probability

The running probability $\psi(\mathbf{x}_i,\tau_i)$ determines the distribution of the running time $\tau_i$ and depends on the way in which the $i{\rm th}$ cell moves. Note that the running probability is here assumed to be independent from the cell velocity $\mathbf{v}_i$ . On the basis of experimental evidence reported in [Reference Boissonnas, Fetler, Zeelenberg, Hugues and Amigorena5,Reference Engelhardt, Boldajipour, Beemiller, Pandurangi, Sorensen, Werb, Egeblad and Krummel10], we assume that inactive CTLs move in a non-local search pattern corresponding to trajectories that are characterised by a strong presence of long runs, which enable them to cover larger areas. On the other hand, activated CTLs and DCs ^{Footnote 1} move in a more localised search pattern. In particular, building upon the modelling approach presented in [Reference Macfarlane, Lorenzi and Chaplain20], we describe the motion of activated CTLs and DCs as a Brownian motion, whereas we let inactive CTLs undergo superdiffusive motion consistent with a Lévy walk, whereby the mean square displacement grows nonlinearly with time. In particular, the mean square displacement at time t is proportional to $t^{{2}/{\alpha}}$ , where $\alpha \in (1,2)$ is the Lévy exponent. We recall that $\alpha=1$ and $\alpha=2$ would correspond to ballistic motion (i.e. a form of motion whereby the mean square displacement at time t is proportional to $t^2$ ) and classical diffusion, respectively.

Under these assumptions, if the $i{\rm th}$ cell belongs to population A or population D, we let the value of the running time $\tau_i$ follow a Poisson distribution [Reference Othmer and Hillen23]. Hence, under the additional simplifying assumption that cells in populations A and D are characterised by the same stopping rate, which is assumed to be constant and thus modelled by a parameter $b \in \mathbb{R}^*_+$ , we use the following definition of the running probability:

(2.2) \begin{equation}\psi(\mathbf{x}_i,\tau_i) \equiv \psi(\tau_i) \;:=\; \exp\left(-b \, \tau_i \right), \quad \varphi(\mathbf{x}_i,\tau_i) \equiv \varphi(\tau_i) \;:=\; b \, \exp\left(-b \, \tau_i \right) \ .\end{equation}

On the other hand, if the $i{\rm th}$ cell belongs to population I, we let the value of the running time $\tau_i$ follow a long-tailed distribution, and we define the running probability along the lines of [Reference Estrada-Rodriguez and Gimperlein11] as:

(2.3) \begin{equation}\psi(\mathbf{x}_i,\tau_i)\;:=\;\Bigl(\frac{\tau_0(\mathbf{x}_i)}{\tau_0(\mathbf{x}_i)+\tau_i} \Bigr)^\alpha\ , \quad \varphi(\mathbf{x}_i,\tau_i)\;:=\;\frac{\alpha \, \tau_0(\mathbf{x}_i)^\alpha}{(\tau_0(\mathbf{x}_i)+\tau_i)^{\alpha+1}}, \quad \alpha\in(1,2) \ .\end{equation}

Here, the function $\tau_0(\mathbf{x}_i) \geqslant 0$ captures possible spatial inhomogeneities in the running time distribution.

Turning kernel and turning operator

We consider the case where the new direction of cell motion given by $\bar{\mathbf{v}}_i$ is symmetrically distributed with respect to the original direction given by $\mathbf{v}_i$ and, therefore, we let the turning kernel $\ell(\mathbf{x}_i,t,\mathbf{v}_i;\bar{\mathbf{v}}_i)$ satisfy the following assumptions [Reference Alt2]:

(2.4) \begin{equation}\ell(\mathbf{x}_i,t,\mathbf{v}_i;\bar{\mathbf{v}}_i) \equiv \ell(\mathbf{x}_i,t,|\bar{\mathbf{v}}_i-\mathbf{v}_i|) \ , \quad \int_{\rm V} \ell(\cdot,\cdot,|\mathbf{v}_i-\mathbf{e}_1|)\mathop{}\!\mathrm{d} \mathbf{v}_i=1 \ ,\end{equation}

where $\mathbf{e}_1 = (1, 0, \ldots, 0) \in \mathbb{R}^n$ is a unit vector.

Moreover, we let the integral operator $\mathcal{T}$ be a turning operator such that for all test functions $\phi(\mathbf{v}_i)$ :

(2.5) \begin{equation} \mathcal{T}[\phi](\cdot, \cdot, \bar{\mathbf{v}}_i) =\int_{\rm V} \ell(\cdot,\cdot,\mathbf{v}_i;\bar{\mathbf{v}}_i) \phi(\mathbf{v}_i) \mathop{}\!\mathrm{d} \mathbf{v}_i \ , \end{equation}

where $\ell$ is the turning kernel defined via (2.4). Since $\displaystyle{\int_{\rm V} \ell(\cdot,\cdot,\cdot;\bar{\mathbf{v}}_i) \mathop{}\!\mathrm{d} \bar{\mathbf{v}}_i =1}$ , we have

(2.6) \begin{equation} \int_{\rm V}(\unicode{x1D7D9}-\mathcal{T})[\phi](\cdot, \cdot,\bar{\mathbf{v}}_i) \mathop{}\!\mathrm{d}\bar{\mathbf{v}}_i=0\ ,\end{equation}

where $\unicode{x1D7D9}$ is the identity operator.

Finally, we recall that in n-dimensions the surface area of the unit sphere ${\rm V} $ is

(2.7) \begin{equation}|{\rm V}|=\begin{cases}\dfrac{2\pi^{{n}/{2}}}{\Gamma\left(\frac{n}{2}\right)} \ , & \text{for } {n} \text{ even},\\[9pt]\dfrac{\pi^{{n}/{2}}}{\Gamma\left(\frac{n}{2}+1\right)} \ , & \text{for } {n} \text{ odd}\end{cases}\end{equation}

where $\Gamma(\cdot)$ is the Gamma function, and we also recall some useful properties of the spectrum of the turning operator $\mathcal{T}$ [Reference Alt2]:

Lemma 1 If the turning kernel $\ell(\cdot,\cdot,|\bar{\mathbf{v}}_i-\mathbf{v}_i|)$ is continuous, then $\mathcal{T}$ is a symmetric compact operator. In particular, there exists an orthonormal basis of $L^2({\rm V})$ consisting of eigenfunctions $\{ \phi_k,\ k\geqslant 0\}$ of $\mathcal{T}$ . Using the notation $\mathbf{v}_i=({v}_0^i,{v}^i_1,\dots,{v}^i_{n-1}) \in {\rm V}$ , we have

\begin{align*}\begin{aligned}\phi_{0}(\mathbf{v}_i) & =\frac{1}{|{\rm V}|} & & \text{is an eigenfunction associated with the eigenvalue} & & \iota_{0}=1,\\[5pt]\phi_{1}^j(\mathbf{v}_i) & =\frac{n{v}^i_j}{|{\rm V}|} & & \text{are eigenfunctions associated with the eigenvalue:} & &\end{aligned}\end{align*}

(2.8) \begin{equation}\iota_{1}=\int_{\rm V}\ell(\cdot,\cdot,|\bar{\mathbf{v}}_i-\mathbf{e}|)\bar{v}^{i}_1\mathop{}\!\mathrm{d}\bar{\mathbf{v}}_i<1,\end{equation}

where $\mathbf{e} = (1, 1, \ldots, 1) \in \mathbb{R}^n$ is the vector with all components equal to 1. Moreover, any function $p_i\in L^2(\mathbb{R}^n\times \mathbb{R}_+\times {\rm V})$ admits a unique decomposition of the form:

(2.9) \begin{equation}p_i=\frac{1}{|{\rm V}|}\left(\rho_i+n\mathbf{v}_i\cdot {w}_i \right)+\hat{z},\end{equation}

where $\hat{z}$ is orthogonal to all linear polynomials in $\mathbf{v}_i$ ,

\begin{equation*}\rho_i(\mathbf{x}_i,t)=\int_{\rm V}p_i(\mathbf{x}_i,t,\mathbf{v}_i)\phi_0(\mathbf{v}_i) \mathop{}\!\mathrm{d}\mathbf{v}_i,\quad {w}_i^j(\mathbf{x}_i,t)=\int_{\rm V} p_i(\mathbf{x}_i,t,\mathbf{v}_i)\phi_1^j(\mathbf{v}_i) \mathop{}\!\mathrm{d}\mathbf{v}_i,\end{equation*}

and ${w}_i = ({w}^i_0, \dots, {w}^i_{n-1})$ .

2.2 Description of the interactions between cells

Building on previous work on individual-based models of interaction dynamics between DCs and CTLs [Reference Macfarlane, Chaplain and Lorenzi19,Reference Macfarlane, Lorenzi and Chaplain20], we consider only the effects of binary cell–cell interactions, thus neglecting interactions that involve more than two cells.

Moreover, given that the focus of this work is on modelling the switch in T cell migration modes mediated by interactions between inactive CTLs and DCs, we explicitly model the effects of interactions between cells of population I and cells of population D, while for simplicity, we neglect the effects of intrapopulation cell–cell interactions and interactions between cells of population I and cells of population A.

Furthermore, the spatial dynamics of DCs are primarily affected by interactions with inactive CTLs [Reference Bousso6,Reference Gardner and Ruffell14,Reference Rothoeft, Balkow, Krummen, Beissert, Varga, Loser, Oberbanscheidt, van den Boom and Grabbe26]. Hence, for simplicity, we explicitly model the effect of interactions between cells of population D and cells of population A on the motion of A cells, while we neglect the effect of these interactions on the motion of D cells, since we take it to be negligible compared to that of interactions with cells of population I.

On the basis of these considerations, we incorporate into the model only the effects of interactions between pairs of cells that are summarised by the schematics in Figure 1, which correspond to the following definitions and assumptions.

Figure 1.

Schematics of cell–cell interactions corresponding to Assumptions 1 and 2. Prime symbols indicate a change in cell velocity upon interaction.

We allow interactions between a cell i in the phase-space state $(\mathbf{x}_i, \mathbf{v}_i)$ and a cell j in the phase-space state $(\mathbf{x}_j, \mathbf{v}_j)$ to occur when the cell j is in the domain of interaction of the cell i, which is defined as the set:

(2.10) \begin{equation}\Omega_j(\mathbf{x}_i) \;:=\; \{ \mathbf{x}_j\in\mathbb{R}^n: |\mathbf{x}_i-\mathbf{x}_j|\geqslant\varrho\} \equiv \mathbb{R}^n \setminus {\rm B}_\varrho(\mathbf{x}_i) \ ,\end{equation}

where ${\rm B}_\varrho(\mathbf{x}_i)$ denotes the ball of radius $\varrho$ centred at $\mathbf{x}_i$ (see the schematics in Figure 2). If a cell i acquires a new velocity upon interaction with a cell j, the new velocity is defined, for simplicity, as the following post-interaction velocity:

(2.11) \begin{equation}\mathbf{v}'_i=\mathbf{v}_i - 2 \, (\mathbf{v}_i\cdot \nu) \, \nu \quad \text{with} \quad \nu \;:=\; \frac{\mathbf{x}_i-\mathbf{x}_j}{|\mathbf{x}_i-\mathbf{x}_j|} \ ,\end{equation}

where $\nu$ is the normal vector at the point of interaction (i.e. $\nu$ is the unit normal that points outward from $\Omega_j(\mathbf{x}_i)$ and inward to ${\rm B}_\varrho(\mathbf{x}_i)$ ) [Reference Cercignani, Illner and Pulvirenti7].

Figure 2.

Schematics of the interaction domain defined in (2.10).

3.1 Preliminaries, assumptions and notation

The state of the system at time t is described by the N-particle distribution function $f^N\left(\mathbf{x}_1, \ldots, \mathbf{x}_N, t, \mathbf{v}_1, \ldots, \mathbf{v}_N, \mathbf{\tau}_1, \ldots, \mathbf{\tau}_N\right)$ [Reference Cercignani, Illner and Pulvirenti7,Reference Villani27]. In the case where cell dynamics at the microscopic scale obey the rules presented in Section 2, the evolution of $f^N$ is governed by the following transport equation [Reference Kennard16]:

(3.1) \begin{equation} \partial_t f^N+\sum_{i=1}^{N}\Bigl(\partial_{\tau_i}\,f^N+c \, \mathbf{v}_i\cdot\nabla_{\mathbf{x}_i}\, f^N \Bigr)=-\sum_{i=1}^{N}\beta \, f^N \end{equation}

posed on $\Omega^{N} \times \mathbb{R}^*_+ \times {\rm V}^{N} \times \mathbb{R}^{* N}_+$ , with

\begin{align*}\Omega^{N} \;:=\; \{(\mathbf{x}_1, ..., \mathbf{x}_{N}) \in \mathbb{R}^{n\times {N}}:\ |\mathbf{x}_i-\mathbf{x}_j|\geqslant\varrho \ \forall i,j \} \ .\end{align*}

We consider the transport equation (3.1) subject to smooth, compactly supported initial conditions at $t=0$ , boundary conditions corresponding to elastic interactions on $\partial\Omega^{N}$ , and suitable Dirichlet boundary conditions at $\tau_i=0$ linked to the running probability $\psi$ for $i=1,\ldots,N$ . In the mathematical framework given by (3.1), the probability of finding at position $\mathbf{x}_1$ and at time t the cell labelled by the index 1 that is moving in direction $\mathbf{v}_1$ for a period of time $\tau_1$ is related to the one-particle marginal:

\begin{align*}{f}(\mathbf{x}_1,t,\mathbf{v}_1,\tau_1)&= \frac{1}{|{\rm V}|^{{ N-1}}}\int_{[0,t]^{{ N-1}}}\int_{\Omega_{N-1}(\mathbf{x}_1)}\int_{{\rm V}^{{ N-1}}}f^{{ N}}(\mathbf{x}_1, \ldots, \mathbf{x}_N, t, \mathbf{v}_1, \ldots, \mathbf{v}_N, \mathbf{\tau}_1, \ldots, \mathbf{\tau}_N)\nonumber\\[5pt] &\ \times\mathop{}\!\mathrm{d}\mathbf{v}_{2} \mathop{}\!\mathrm{d}\mathbf{x}_{2} \mathop{}\!\mathrm{d}\tau_{2} \ldots \mathop{}\!\mathrm{d} \mathbf{v}_{N} \mathop{}\!\mathrm{d} \mathbf{x}_{N} \mathop{}\!\mathrm{d} \tau_{N}\ .\end{align*}

Here, $|{\rm V}|$ denotes the surface area of the unit sphere V and $\Omega_{N-1}(\mathbf{x}_1)\;:=\; \{(\mathbf{x}_{2}, \ldots,\mathbf{x}_{N}) \in \mathbb{R}^{n \times N-1} : (\mathbf{x}_1,\mathbf{x}_{2},\ldots,\mathbf{x}_{N})\in\Omega^{N} \}$ .

A comprehensive description of cell dynamics would in principle require considering possible interactions between all cells. However, as mentioned earlier, building on previous work on the mathematical modelling of the interaction dynamics between DCs and CTLs [Reference Macfarlane, Chaplain and Lorenzi19,Reference Macfarlane, Lorenzi and Chaplain20], we consider only the effect of binary cell–cell interactions, thus neglecting interactions that involve more than two cells. Therefore, as per the scaling and assumptions introduced in Section 4.1, which are similar to those typically considered in low-density regimes [Reference Cercignani, Illner and Pulvirenti7,Reference Othmer, Dunbar and Alt22,Reference Villani27], we truncate the hierarchy of equations corresponding to (3.1) at the second order by integrating out cells $3, \ldots ,N$ from the N-particle distribution function $f^N\left(\mathbf{x}_1, \ldots, \mathbf{x}_N, t, \mathbf{v}_1, \ldots, \mathbf{v}_N, \mathbf{\tau}_1, \ldots, \mathbf{\tau}_N\right)$ .

Two-particle distribution functions

Let $f_{hk}(\mathbf{x}_h,\mathbf{x}_k,t,\mathbf{v}_h,\mathbf{v}_k,\tau_h,{\tau_k})$ with $h, k \in \{A,D,I\}$ and $k \neq h$ denote the two-particle distribution function associated with:

• – a cell of population h in the generic phase-space state $(\mathbf{x}_h,\mathbf{v}_h) \in \mathbb{R}^n \times {\rm V}$ , with generic run time $\tau_h \in [0,t]$ and stopping rate $ \beta_h(\mathbf{x}_h,\tau_h)$ defined via (2.1);

• – a cell of population k in the generic phase-space state $(\mathbf{x}_k,\mathbf{v}_k) \in \mathbb{R}^n \times {\rm V}$ , with generic run time $\tau_k \in [0,t]$ and stopping rate $ \beta_k(\mathbf{x}_k,\tau_k)$ defined via (2.1).

Truncating the hierarchy of equations corresponding to (3.1) at the second order, we obtain the following transport equation for $f_{hk}(\mathbf{x}_h,\mathbf{x}_k,t,\mathbf{v}_h,\mathbf{v}_k,\tau_h,{\tau_k})$ :

(3.2) \begin{equation}\begin{aligned} (\partial_t+\partial_{\tau_h}+\partial_{\tau_k} + c \, \mathbf{v}_h\cdot\nabla_{\mathbf{x}_h}+c \, \mathbf{v}_k\cdot\nabla_{\mathbf{x}_k})f_{hk}=-(\beta_h+\beta_k)f_{hk}\end{aligned}\end{equation}

posed on $\Omega^2 \times \mathbb{R}_+ \times {\rm V}^{2} \times \mathbb{R}^{* 2}_+$ , with

(3.3) \begin{equation}\Omega^2 \;:=\; \{(\mathbf{x}_h, \mathbf{x}_{k}) \in \mathbb{R}^{n\times {2}}:\ |\mathbf{x}_h-\mathbf{x}_k|\geqslant\varrho \ \ \forall h,k \} \ .\end{equation}

This equation is subject to a smooth, compactly supported initial condition at $t=0$ , specular reflective boundary conditions corresponding to elastic interactions on $\partial \Omega^2$ , and with boundary conditions at $\tau_k=0$ and $\tau_h=0$ given by:

(3.4) \begin{equation}\begin{aligned} f_{hk}(\mathbf{x}_h,\mathbf{x}_k,t,\mathbf{v}_h,\mathbf{v}_k,\tau_h,\tau_k=0)&=\mathcal{T}\int_0^t\beta_kf_{hk}(\mathbf{x}_h,\mathbf{x}_k,t,\mathbf{v}_h,\mathbf{v}_k,\tau_h,{\tau_k})\mathop{}\!\mathrm{d} \tau_k\ ,\\[5pt] f_{hk}(\mathbf{x}_h,\mathbf{x}_k,t,\mathbf{v}_h,\mathbf{v}_k,\tau_h=0,\tau_k)&=\mathcal{T}\int_0^t\beta_hf_{hk}(\mathbf{x}_h,\mathbf{x}_k,t,\mathbf{v}_h,\mathbf{v}_k,\tau_h,{\tau_k})\mathop{}\!\mathrm{d} \tau_h\ .\end{aligned}\end{equation}

One-particle distribution functions

Given the two-particle distribution function:

(3.5) \begin{equation}\tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_k,t,\mathbf{v}_h,\mathbf{v}_k)\;:=\;\int_0^t\int_0^tf_{hk}\mathop{}\!\mathrm{d}\tau_h \mathop{}\!\mathrm{d}\tau_k \ , \end{equation}

the one-particle distribution function of population h is given by:

(3.6) \begin{equation}p_h(\mathbf{x}_h,t,\mathbf{v}_h)\;:=\;\frac{1}{|{\rm V}|}\int_{\Omega_k(\mathbf{x}_h)}\int_{\rm V}\tilde{\tilde{f}}_{hk}\mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\mathbf{x}_k \ , \end{equation}

with the set $\Omega_k(\mathbf{x}_h)$ defined via (2.10). The function $p_h(\mathbf{x}_h,t,\mathbf{v}_h)$ describes the density of cells of population h which at position $\mathbf{x}_h$ and time t move with velocity $\mathbf{v}_h$ (i.e. the quantity $p_h(\mathbf{x}_h,t,\mathbf{v}_h) \mathop{}\!\mathrm{d}\mathbf{v}_h\mathop{}\!\mathrm{d}\mathbf{x}_h$ is the number of cells of population h in the volume element $ \mathop{}\!\mathrm{d}\mathbf{v}_h\mathop{}\!\mathrm{d}\mathbf{x}_h$ centred at the point $(\mathbf{x}_h, \mathbf{v}_h)$ of the phase space). Moreover, we will consider the weighted two-particle distribution function given by:

(3.7) \begin{equation}\tilde{\tilde{f}}^{\beta_z}_{hk}(\mathbf{x}_h,\mathbf{x}_k,t,\mathbf{v}_h,\mathbf{v}_k)\;:=\;\int_0^t\int_0^t \beta_z \ f_{hk}\mathop{}\!\mathrm{d}\tau_h \mathop{}\!\mathrm{d}\tau_k \ , \quad z \in \{h, k \} \ , \end{equation}

and the weighted one-particle distribution function given by:

(3.8) \begin{equation}p^{\beta_h}_{h}(\mathbf{x}_h,t,\mathbf{v}_h) \;:=\;\frac{1}{|{\rm V}|}\int_{\Omega_k(\mathbf{x}_h)} \ \int_{\rm V} \ \tilde{\tilde{f}}^{\beta_h}_{hk} \mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\mathbf{x}_k \ . \end{equation}

Here, $|{\rm V}|$ denotes the surface area of the unit sphere V.

3.2 Derivation of a system of transport equations

Transport equations for two-particle distribution functions

The dynamics of the two-particle distribution functions $f_{ID}$ and $f_{AD}$ are governed by the following specific forms of transport equation (3.2):

(3.9) \begin{align} (\partial_t+\partial_{\tau_I}+\partial_{\tau_D} + c \, \mathbf{v}_I\cdot\nabla_{\mathbf{x}_I}+c \, \mathbf{v}_D\cdot\nabla_{\mathbf{x}_D})\,f_{ID}&=-(\beta_I+\beta_D)\,f_{ID}\ ,\ \end{align}

(3.10) \begin{align} (\partial_t+\partial_{\tau_A}+\partial_{\tau_D}+c \, \mathbf{v}_A\cdot\nabla_{\mathbf{x}_A}+c \, \mathbf{v}_D\cdot\nabla_{\mathbf{x}_D})\,f_{AD}&= -(\beta_A+\beta_D)\,f_{AD}\ , \end{align}

which are posed on $\Omega^2 \times \mathbb{R}^*_+ \times {\rm V}^{2} \times \mathbb{R}^{* 2}_+$ . The boundary conditions at $\tau_A=0,\ \tau_D=0$ and $\tau_I=0$ are analogous to (3.4). Starting from transport equations (3.9)–(3.10) and using the method employed in [Reference Estrada-Rodriguez and Gimperlein11], it is possible to show (see 6) that the two-particle distribution functions $\tilde{\tilde{f}}_{ID}$ and $\tilde{\tilde{f}}_{AD}$ given by (3.5) satisfy the following transport equations:

(3.11) \begin{align} (\partial_t+c \ \mathbf{v}_I\cdot\nabla_{\mathbf{x}_I} +c \ \mathbf{v}_D\cdot\nabla_{\mathbf{x}_D})\,\tilde{\tilde{f}}_{ID} =&-(\unicode{x1D7D9}-\mathcal{T}_I)\left[\tilde{\tilde{f}}^{\beta_I}_{ID}\right] -(\unicode{x1D7D9}-\mathcal{T}_D)\left[\tilde{\tilde{f}}^{\beta_D}_{ID}\right]\ , \end{align}

(3.12) \begin{align} (\partial_t+c \ \mathbf{v}_A\cdot\nabla_{\mathbf{x}_A} +c \ \mathbf{v}_D\cdot\nabla_{\mathbf{x}_D})\,\tilde{\tilde{f}}_{AD} =&-(\unicode{x1D7D9}-\mathcal{T}_A)\left[\tilde{\tilde{f}}^{\beta_A}_{AD}\right] -(\unicode{x1D7D9}-\mathcal{T}_D)\left[\tilde{\tilde{f}}^{\beta_D}_{AD}\right]\ , \end{align}

posed on $\Omega^2 \times \mathbb{R}^*_+ \times {\rm V}^{2}$ . Here, $\mathcal{T}_I$ , $\mathcal{T}_D$ and $\mathcal{T}_A$ are the turning operators defined via (2.5), and $\tilde{\tilde{f}}^{\beta_h}_{hk}$ and $\tilde{\tilde{f}}^{\beta_k}_{hk}$ are the weighted two-particle distribution functions given by (3.7).

Transport equation for $p_h$

Starting from transport equation (3.2) and building upon the method presented in [Reference Estrada-Rodriguez and Gimperlein11], it is possible to show (see 7) that the one-particle distribution function $p_h(\mathbf{x}_h,t,\mathbf{v}_h)$ given by (3.6) satisfies the following transport equation:

(3.13) \begin{equation} \partial_tp_h+c \, \mathbf{v}_h\cdot\nabla_{\mathbf{x}_h}p_h=-(\unicode{x1D7D9}-\mathcal{T}_h)\left[p^{\beta_h}_{h}\right] + \mathcal{Q}_{hk}, \quad \mathbf{x}_h \in \mathbb{R}^n, t \in \mathbb{R}_+, \mathbf{v}_h \in {\rm V} \ . \end{equation}

Here, the turning operator $\mathcal{T}_h$ is defined via (2.5), the weighted one-particle distribution function $p^{\beta_h}_{h}(\mathbf{x}_h,t,\mathbf{v}_h)$ is given by (3.8) and

(3.14) \begin{equation}\mathcal{Q}_{hk}(\mathbf{x}_h,t,\mathbf{v}_h) \;:=\;\dfrac{c}{|{\rm V}|}\int_{\partial {\rm B}_\varrho(\mathbf{x}_h)}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\tilde{\tilde{f}}_{hk} \mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\sigma \ .\end{equation}

In (3.14), $\nu$ is the unit normal defined in (2.11) and $\mathop{}\!\mathrm{d} \sigma$ denotes the surface element.

The first term on the right-hand side of transport equation (3.13) represents the rate of change of the one-particle distribution function due to cell movement, while the term $\mathcal{Q}_{hk}$ is the rate of change due to interactions between cells. The specific forms of these terms depend, respectively, on the way in which cells move and the interactions they undergo, as discussed in the remainder of this section.

Expressions for $p^{\beta_h}_{h}$

The specific form of the first term on the right-hand side of transport equation (3.13) depends on the expression for $p^{\beta_h}_{h}$ which, in turn, will depend on the definition of the stopping rate $\beta_h$ .

When cells move in a local search pattern (i.e. for $h=A$ and $h=D$ ), the stopping rate $\beta_h$ is defined via (2.1) and (2.2). In this case, inserting the definition of $\beta_h$ into (3.8) yields

(3.15) \begin{equation}p_h^{\beta_h}(\mathbf{x}_h,t,\mathbf{v}_h)=b \, p_h(\mathbf{x}_h,t,\mathbf{v}_h) \ .\end{equation}

On the other hand, when cells move in a non-local search pattern (i.e. for $h=I$ ), the stopping rate $\beta_h$ is defined via (2.1) and (2.3). In this case, it is possible to show (see 8) that

(3.16) \begin{equation} p^{\beta_h}_h(\mathbf{x}_h,t,\mathbf{v}_h)=\mathcal{B}[p_h](\mathbf{x}_h,t,\mathbf{v}_h) \ ,\end{equation}

where $\mathcal{B}$ is a convolution operator such that

(3.17) \begin{equation}\mathcal{B}[p_h](\mathbf{x}_h,t, \mathbf{v}_h) = \int_0^t B(\mathbf{x}_h,t-s)p(\mathbf{x}_h-(c\mathbf{v}_h+b)(t-s),s,\mathbf{v}_h)\mathop{}\!\mathrm{d} s \ ,\end{equation}

with B being defined through its Laplace transform in time $\hat{B}$ as:

(3.18) \begin{equation} \hat{B}(\mathbf{x}_h,\lambda+b+c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})=\frac{\hat{\varphi}_h(\mathbf{x}_h,\lambda+b+c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})}{\hat{\psi}_h(\mathbf{x}_h,\lambda+b+c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})}\ .\end{equation}

Here, $\lambda$ is the Laplace variable, $\hat{\varphi}_h$ and $\hat{\psi}_h$ are the Laplace transforms in $\tau_h$ of the functions $\varphi_h$ and $\psi_h$ defined via (2.3), and the parameter b is defined via (2.2).

Expressions for $\mathcal{Q}_{hk}$

Following [Reference Estrada-Rodriguez and Gimperlein11,Reference Franz, Taylor-King, Yates and Erban13], we first note that when a cell in the phase-space state $({\bf x}_h, {\bf v}_h)$ interacts with a cell in the phase-space state $({\bf x}_k, {\bf v}_k)$ , we have $|{\bf x}_h - {\bf x}_k| = \varrho$ . Hence, the normal vector at the point of physical contact between the interacting cells, $\nu \in {\rm V}$ , defined via (2.11) can be written as $\nu=(\mathbf{x}_h-\mathbf{x}_k)/\varrho$ , that is, $\mathbf{x}_k=\mathbf{x}_h-\nu\varrho$ . As a result, using the fact that ${\rm B}_\varrho = \varrho {\rm V}$ along with the change of variable $\nu \mapsto -\nu$ , we rewrite (3.14) as:

(3.19) \begin{align}\mathcal{Q}_{hk}(\mathbf{x}_h,t,\mathbf{v}_h) \;:=\; & \frac{c}{|{{\rm V}}|}\int_{\partial {\rm B}_\varrho(\mathbf{x}_h)}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\tilde{\tilde{f}}_{hk} \mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\sigma \nonumber \\[5pt] = & -\frac{c}{|{{\rm V}}|}\varrho^{n-1}\int_{{\rm V}}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h+\nu\varrho,t,\mathbf{v}_h,\mathbf{v}_k)\mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu\ .\end{align}

Following [Reference Estrada-Rodriguez and Gimperlein11], we also note that ${\rm V} \equiv {\rm V}^+_{hk}\cup {\rm V}^-_{hk}$ with

(3.20) \begin{equation}\begin{aligned}{\rm V}^+_{hk} &\;:=\; \{\nu \in {\rm V} : \nu\cdot(\mathbf{v}_h-\mathbf{v}_k)>0 \}\ , \\[5pt]{\rm V}^-_{hk} &\;:=\; \{\nu \in {\rm V} : \nu\cdot(\mathbf{v}_h-\mathbf{v}_k)<0 \} \equiv \{-\nu \in {\rm V} : \nu\cdot(\mathbf{v}_h-\mathbf{v}_k)>0 \}\ .\end{aligned}\end{equation}

Therefore, a cell moving in direction ${\bf v}_h$ and a cell moving in direction ${\bf v}_k$ will move towards each other if $\nu \in {\rm V}^+$ and away from each other if $\nu \in {\rm V}^-$ .

Under Assumptions 1–2, denoting the post-interaction directions corresponding to $\mathbf{v}_h$ and $\mathbf{v}_k$ by $\mathbf{v}'_{\!h}$ and $\mathbf{v}'_{\!k}$ , which are defined via (2.11), we consider two different types of interactions between cells: conservative interactions (cf. Definition 2.1), between a cell of population h and a cell of population k, whereby both cells remain in their original populations upon interaction and acquire the post-interaction velocities; population-switching interactions (cf. Definition 2.2), between a cell in population h and a cell in population k, whereby the h cell switches from its original population to a different one upon interaction.

From (3.19), we define the rate of change of the one-particle distribution function $p_h(\mathbf{x}_h,t,\mathbf{v}_h)$ due to conservative interactions as:

(3.21) \begin{align} \mathcal{K}_{hk}(\mathbf{x}_h,t,\mathbf{v}_h) & \;:=\;-\frac{c}{|{{\rm V}}|}\varrho^{n-1} \, N_k(t) \, \left[\int_{{\rm V}^+_{hk}}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h+\nu\varrho,t,\mathbf{v}_h,\mathbf{v}_k)\mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu\right.\nonumber\\[5pt] & \qquad \left.+\int_{{\rm V}^-_{hk}}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h+\nu\varrho,t,\mathbf{v}_h,\mathbf{v}_k)\mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu\right]\nonumber\\[5pt] &=\frac{c}{|{{\rm V}}|}\varrho^{n-1} \, N_k(t) \, \int_{{\rm V}_{hk}^+}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\Bigl[\tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h-\nu\varrho,t,\mathbf{v}'_h,\mathbf{v}'_k)\nonumber\\[5pt] &\qquad -\tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h+\nu\varrho,t,\mathbf{v}_h,\mathbf{v}_k)\Bigr]\mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu\ , \end{align}

with $N_k(t)$ being the number of cells in population k at time t. The second equality in (3.21) is obtained by using the normal vector $-\nu$ and the post-interaction directions $\mathbf{v}_h'$ and $\mathbf{v}_k'$ in $\tilde{\tilde{f}}_{hk}$ over the set ${\rm V}_{hk}^-$ . Notice that the following property holds

(3.22) \begin{equation}\int_{\rm V} \mathcal{K}_{hk}(\cdot,\cdot,\mathbf{v}_h) \mathop{}\!\mathrm{d}\mathbf{v}_h = 0 \ ,\end{equation}

which ensures that the density of cells in population h will be preserved in the course of such interactions.

Moreover, based on~(3.19) and~(3.21), we define the rate of change of the one-particle distribution function $p_h(\mathbf{x}_h,t,\mathbf{v}_h)$ due to population-switching interactions leading the cell to leave population h as:

(3.23) \begin{equation} \mathcal{J}_{hk}(\mathbf{x}_h,t,\mathbf{v}_h) \;:=\; - \dfrac{c}{|{\rm V}|} \, \varrho^{n-1} \, N_k(t) \, \int_{{\rm V}^+_{hk}}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k) \, \tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h+\nu\varrho, t,\mathbf{v}_h,\mathbf{v}_k)\mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu \ . \end{equation}

Analogously, we define the rate of change of $p_h(\mathbf{x}_h,t,\mathbf{v}_h)$ due to population-switching interactions leading a cell to leave a generic population $l \neq h$ and enter population h as:

(3.24) \begin{align}\mathcal{J}^h_{lk}(\mathbf{x}_h,t,\mathbf{v}_h) \;:=\; &\dfrac{c}{|{\rm V}|} \, \varrho^{n-1} \, N_k(t) \, \int_{ \Omega_l(\mathbf{x}_k)}\int_{\rm V} \delta(\mathbf{x}_l - \mathbf{x}_h) \, \delta(\mathbf{v}_l - \mathbf{v}_h)\nonumber \\[5pt]&\quad\times \int_{{\rm V}^+_{lk}}\int_{\rm V}\nu\cdot(\mathbf{v}_l-\mathbf{v}_k) \tilde{\tilde{f}}_{lk}(\mathbf{x}_l,\mathbf{x}_l+\nu\varrho, t, \mathbf{v}_l,\mathbf{v}_k) \mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu \mathop{}\!\mathrm{d}\mathbf{v}_l\mathop{}\!\mathrm{d}\mathbf{x}_l \ ,\end{align}

with $\delta({\mathbf{z}} - {\mathbf{z}}^*)$ being the Dirac delta distribution centred at ${\mathbf{z}}^*$ . Definition (3.24) ensures that the density of cells that leave population l due to such interactions will appear in population h. In fact, we have

\begin{equation*} \mathcal{J}_{lk}^h(\mathbf{x}_h,t,\mathbf{v}_h)=-\int_{\Omega_l(\mathbf{x}_k)}\int_{\rm V}\delta(\mathbf{x}_l-\mathbf{x}_h)\delta(\mathbf{v}_l-\mathbf{v}_h)\mathcal{J}_{lk}(\mathbf{x}_l,t,\mathbf{v}_l)\mathop{}\!\mathrm{d}\mathbf{v}_l\mathop{}\!\mathrm{d}\mathbf{x}_l\ .\end{equation*}

In summary, the term $\mathcal{Q}_{hk}$ in transport equation (3.13) is defined in terms of (3.21)–(3.24) in different possible ways depending on the cell–cell interactions that are considered.

Under Assumptions 1–2 and definitions (3.21), (3.23) and (3.24), the rates of change of the one-particle distribution functions $p_I(\mathbf{x}_I,t,\mathbf{v}_I)$ , $p_A(\mathbf{x}_A,t,\mathbf{v}_A)$ and $p_D(\mathbf{x}_D,t,\mathbf{v}_D)$ due to cell–cell interactions will be, respectively,

(3.25) \begin{align}\mathcal{Q}_{ID}(\mathbf{x}_I,t,\mathbf{v}_I)& = (1-\zeta)\mathcal{K}_{ID}(\mathbf{x}_I,t,\mathbf{v}_I) {+} \zeta \mathcal{J}_{ID}(\mathbf{x}_I,t,\mathbf{v}_I) \ ,\end{align}

(3.26) \begin{align}\mathcal{Q}_{AD}(\mathbf{x}_A,t,\mathbf{v}_A) &= \mathcal{K}_{AD}(\mathbf{x}_A,t,\mathbf{v}_A) + \zeta \, \mathcal{J}^A_{ID}(\mathbf{x}_A,t,\mathbf{v}_A) \ , \end{align}

(3.27) \begin{align}\mathcal{Q}_{DI}(\mathbf{x}_D,t,\mathbf{v}_D) &= \mathcal{K}_{DI}(\mathbf{x}_D,t,\mathbf{v}_D) \ .\end{align}

Substituting (3.15), (3.16) and (3.25)–(3.27) into transport equation (3.13), we obtain the following transport equations for $p_I(\mathbf{x}_I,t,\mathbf{v}_I)$ , $p_A(\mathbf{x}_A,t,\mathbf{v}_A)$ and $p_D(\mathbf{x}_D,t,\mathbf{v}_D)$ :

(3.28) \begin{align} \partial_tp_I+c \, \mathbf{v}_I\cdot\nabla_{\mathbf{x}_I}p_I =&\underbrace{-(\unicode{x1D7D9}-\mathcal{T}_I)\mathcal{B}[p_I]}_{{\scriptsize{\rm cell\ motion}}} + \underbrace{(1-\zeta)\mathcal{K}_{ID}}_{{\scriptsize{\rm interactions}}} \nonumber \\[5pt] & \qquad \qquad \qquad \underbrace{{+} \, \zeta \mathcal{J}_{ID},}_{\substack{{\scriptsize{\rm outflow\ due}}\\[5pt] {\scriptsize{\rm to\ activation}}}} \quad \mathbf{x}_I \in \mathbb{R}^n, t \in \mathbb{R}^*_+,\mathbf{v}_I \in {\rm V} \ , \end{align}

(3.29) \begin{align} \partial_t{p}_A+c \, \mathbf{v}_A\cdot\nabla_{\mathbf{x}_A}{p}_A=& \underbrace{-b \, (\unicode{x1D7D9}-\mathcal{T}_A)[p_{A}]}_{{\scriptsize{\rm cell\ motion}}} + \underbrace{\mathcal{K}_{AD}}_{{\scriptsize{\rm interactions}}} \nonumber \\[5pt] & \quad \underbrace{+ \zeta \, \mathcal{J}^A_{ID},}_{\substack{{\scriptsize{\rm inflow\ due}}\\[5pt] {\scriptsize{\rm to\ activation}}}} \quad \mathbf{x}_A \in \mathbb{R}^n, t \in \mathbb{R}^*_+,\mathbf{v}_A \in {\rm V}\ , \end{align}

(3.30) \begin{align} \partial_t {p}_D +c\mathbf{v}_D\cdot\nabla_{\mathbf{x}_D}{p}_D = &\underbrace{-b \, (\unicode{x1D7D9}-\mathcal{T}_D)[p_{D}]}_{{\scriptsize{\rm cell\ motion}}} \nonumber \\[5pt] & \quad \underbrace{+\mathcal{K}_{DI},}_{{\scriptsize{\rm interactions}}} \quad \mathbf{x}_D \in \mathbb{R}^n, t \in \mathbb{R}^*_+,\mathbf{v}_D \in {\rm V} \ . \end{align}

The terms on the right-hand sides of (3.28)–(3.30) represent the rate of change of the one-particle distributions due to the biophysical phenomena specified below each term.

4.1 Preliminaries, assumptions and notation

Scaling

We assume the mean run time $\bar \tau$ to be small compared to the characteristic temporal scale for the dynamics of the macroscopic cell densities, which is represented by the parameter $T \in \mathbb{R}^*_+$ , that is, we make the assumption:

\begin{align*}\dfrac{\bar \tau}{T} \;=:\; \varepsilon \ll 1 \ .\end{align*}

Moreover, we let $X \in \mathbb{R}^*_+$ represent the characteristic spatial scale for the dynamics of the macroscopic cell densities and introduce the rescaled quantities:

\begin{align*}\hat{t} = \dfrac{t}{T}, \quad \hat{{\bf x}} = \dfrac{{\bf x}}{X}, \quad \hat{\tau} = \dfrac{\bar \tau}{T}, \quad \hat{c} = c \, \dfrac{T}{X} \ .\end{align*}

As similarly done in [Reference Alt2,Reference Estrada-Rodriguez and Gimperlein11], in order to obtain a mathematical model for the dynamics of the cells at the macroscopic scale, we consider the scaling:

(4.1) \begin{equation} (\mathbf{x},t,c,\tau)\mapsto(\hat{\mathbf{x}}/\varepsilon,\hat{t}/\varepsilon,\hat{c}/\varepsilon^\gamma,\hat{\tau}/\varepsilon^\mu) \ ,\end{equation}

with

(4.2) \begin{equation}\gamma, \mu \in \mathbb{R}^*_+ \ , \quad \gamma < 1 \quad \text{and} \quad \mu > 1 - \gamma \ .\end{equation}

Throughout the rest of the paper, we will drop the carets from (4.1) and we will study two-dimensional cell dynamics (i.e. we assume $n=2$ ).

Furthermore, noting that the diameter of the cells is small compared to the characteristic spatial scale for the dynamics of the macroscopic cell densities, and considering a biological scenario where the number of cells in the system is large and activation of CTL occurs with a small probability $\zeta$ , we assume

(4.3) \begin{equation}\varrho = \varepsilon^{\xi} \ , \quad N_I(t) \equiv \varepsilon^{-\vartheta}, \quad N_D = \varepsilon^{-\vartheta}\ , \quad \zeta = \varepsilon^\kappa \ , \quad \xi, \vartheta, \kappa \in \mathbb{R}^*_+ \ .\end{equation}

In particular, we will be focussing on a biological scenario corresponding to the following assumptions:

(4.4) \begin{equation}\gamma \;:=\; \dfrac{1}{2} \ , \quad \xi-\vartheta \;:=\; 1 - \dfrac{\gamma}{\alpha -1} \quad \text{and} \quad \kappa = - (\xi-\vartheta) + \dfrac{3}{2} > 0 \ .\end{equation}

Notice that $\xi-\vartheta<0$ when $\alpha<3/2$ . In the case where cells follow a Brownian motion (i.e. for $h=A$ and $h=D$ ) we have $\alpha=2$ and, therefore, $\xi-\vartheta=1-\gamma=1/2$ . Under scaling (4.1) definitions (2.3) become

(4.5) \begin{equation} \psi^\varepsilon(\mathbf{x}_i,\tau_i)=\Bigl(\frac{\varepsilon^\mu\tau_0({\bf x}_i)}{\varepsilon^\mu\tau_0({\bf x}_i)+\tau_i}\Bigr)^\alpha\ , \quad \varphi^\varepsilon(\mathbf{x}_i,\tau_i)\;:=\;\frac{\alpha \, \varepsilon^\mu \, \tau_0(\mathbf{x}_i)^\alpha}{(\varepsilon^\mu \, \tau_0(\mathbf{x}_i)+\tau_i)^{\alpha+1}}, \quad \alpha\in(1,2) \ . \end{equation}

Moreover, under assumption (4.3) on $\varrho$ we have

(4.6) \begin{equation} \tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h \pm \nu \rho,t,\mathbf{v}_h,\mathbf{v}_k) \equiv \tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h \pm\varepsilon^{\unicode{x03BE}}\nu,t,\mathbf{v}_h,\mathbf{v}_k) \ .\end{equation}

‘Molecular chaos’ assumption

Considering a biological scenario where cell densities are sufficiently low, we assume the velocities of any two cells which are about to interact to be uncorrelated, that is, we make the so-called ‘molecular chaos’ assumption, which holds at low densities and is commonly used in kinetic theory [Reference Cercignani, Illner and Pulvirenti7,Reference Othmer, Dunbar and Alt22,Reference Villani27]. Under this assumption, the two-particle distribution function $\tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h \pm\varepsilon^{\unicode{x03BE}}\nu,t,\mathbf{v}_h,\mathbf{v}_k)$ can be approximately expressed as the product of the corresponding one-particle distribution functions, that is,

(4.7) \begin{equation} \tilde{\tilde{f}}_{hk}(\mathbf{x}_h,\mathbf{x}_h \pm\varepsilon^{\unicode{x03BE}}\nu,t,\mathbf{v}_h,\mathbf{v}_k)=p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h) \, p_{k}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_k)+\mathcal{O}(\varepsilon^{\unicode{x03BE}})\ .\end{equation}

We draw the attention of the reader to the fact that, throughout the rest of the paper, superscripts and subscripts related to the scaling should not be confused with cell population indices.

Under scaling (4.1) and assumptions (4.3), using (4.6), (4.7) and assuming $n=2$ , the interaction terms defined via (3.21), (3.23) and (3.24) read as:

(4.8) \begin{align} \mathcal{K}_{hk}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)& = \frac{1}{|{{\rm V}}|} \varepsilon^{\xi -\vartheta - \gamma} \,c \, \int_{{\rm V}_{hk}^+}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\Bigl[p_{h}^\varepsilon(\mathbf{x}_h,t,{\bf v}'_h) p_{k}^\varepsilon(\mathbf{x}_h,t,{\bf v}'_k) \nonumber\\[5pt] &\qquad\qquad - p_{h}^\varepsilon(\mathbf{x}_h,t,{\bf v}_h) p_{k}^\varepsilon(\mathbf{x}_h,t,{\bf v}_k)\Bigr]\mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu\ ,\end{align}

(4.9) \begin{align} \mathcal{J}_{hk}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)&= - \frac{1}{|{{\rm V}}|} \varepsilon^{\xi -\vartheta - \gamma} \, c \, \int_{{\rm V}^+_{hk}}\int_{\rm V}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k) \, p_{h}^\varepsilon(\mathbf{x}_h,t,{\bf v}_h) p_{k}^\varepsilon(\mathbf{x}_h,t,{\bf v}_k)\mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu \end{align}

and

(4.10) \begin{align}{}_{\varepsilon}\mathcal{J}^h_{lk }(\mathbf{x}_h,t,\mathbf{v}_h) = &\dfrac{1}{|{\rm V}|} \, \varepsilon^{\xi -\vartheta - \gamma} \, c \, \int_{ \Omega_l(\mathbf{x}_k)}\int_{\rm V} \delta(\mathbf{x}_l - \mathbf{x}_h) \, \delta(\mathbf{v}_l - \mathbf{v}_h)\nonumber \\[5pt]&\quad\times \int_{{\rm V}^+_{lk}}\int_{\rm V}\nu\cdot(\mathbf{v}_l-\mathbf{v}_k) p_{l}^\varepsilon(\mathbf{x}_l,t,{\bf v}_l) p_{k}^\varepsilon(\mathbf{x}_l,t,{\bf v}_k) \mathop{}\!\mathrm{d}\mathbf{v}_k\mathop{}\!\mathrm{d}\nu \mathop{}\!\mathrm{d}\mathbf{v}_l\mathop{}\!\mathrm{d}\mathbf{x}_l \ .\end{align}

Expansion of $p_{h}^\varepsilon$ and macroscopic cell quantities

Exploiting the results established by Lemma 1 in the case where $n=2$ , we expand the one-particle distribution function $p_{h}^\varepsilon$ in terms of its zeroth moment $\rho_{h}^\varepsilon$ (i.e. the macroscopic cell density) and its first moment $w_{h}^\varepsilon$ (i.e. the local macroscopic direction of cell motion). This is possible because, as one can see from the right-hand side of transport equation (4.13), the interaction terms are of higher order in $\varepsilon$ (cf. the scaling used in (4.8)–(4.10)) and, therefore, we can write

(4.11) \begin{equation} p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)=\frac{1}{|{\rm V}|} \Big(\rho_{h}^\varepsilon(\mathbf{x}_h,t) + \varepsilon^{\gamma} \,{ 2} \, \mathbf{v}_h\cdot w_{h}^\varepsilon(\mathbf{x}_h,t) \Big) + o(\varepsilon^{\gamma}) \ , \quad h \in \{A, D, I \} \ ,\end{equation}

where

(4.12) \begin{equation}\rho_{h}^\varepsilon(\mathbf{x}_h,t) \;:=\; \int_{\rm V} p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)\mathop{}\!\mathrm{d} \mathbf{v}_h\ ,\ \ \ w_{h}^\varepsilon(\mathbf{x}_h,t) \;:=\; \frac{1}{\varepsilon^\gamma}\int_{\rm V}\mathbf{v}_h \, p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)\mathop{}\!\mathrm{d}\mathbf{v}_h \ .\end{equation}

We refer the reader to [Reference Othmer and Hillen23,Reference Othmer, Maini and Murray24] and the seminal work [Reference Alt2] for a complete derivation in the case of no interactions and to [Reference Estrada-Rodriguez and Gimperlein11,Reference Franz, Taylor-King, Yates and Erban13] for the case of velocity-jump models with interacting particles. The appropriate choice of scaling for the local macroscopic direction of motion is found by first inserting (4.11) into (4.13) and then integrating over V in order to obtain a suitable macroscopic equation (see transport equation (4.2)).

4.2 Derivation of a macroscopic scale system

Transport equation for $p_{h}^\varepsilon$

Under scaling (4.1) and assumptions (4.3), using (4.6) and (4.7) and assuming $n=2$ , we rewrite transport equation (3.13) for the one-particle distribution function $p_h(\mathbf{x}_h,t,\mathbf{v}_h)$ as:

(4.13) \begin{equation}\varepsilon \partial_tp_{h}^\varepsilon + \varepsilon^{1-\gamma} c \, \mathbf{v}_h\cdot\nabla_{\mathbf{x}_h}p_{h}^\varepsilon=-(\unicode{x1D7D9}-\mathcal{T}_h)[{}_{\varepsilon}p^{\beta_h}_{h}] + \mathcal{Q}_{hk}^\varepsilon \ , \end{equation}

where $\mathcal{Q}_{hk}^\varepsilon$ is defined in terms of $\mathcal{K}_{hk}^\varepsilon$ , $ \mathcal{J}_{hk}^\varepsilon$ and ${}_{\varepsilon}\mathcal{J}^h_{lk}$ as per (3.25)–(3.27), that is,

(4.14) \begin{align}\mathcal{Q}_{ID}^\varepsilon(\mathbf{x}_I,t,\mathbf{v}_I) &= (1-\varepsilon^\kappa)\mathcal{K}_{ID}^\varepsilon(\mathbf{x}_I,t,\mathbf{v}_I) + \varepsilon^\kappa \mathcal{J}_{ID}^\varepsilon(\mathbf{x}_I,t,\mathbf{v}_I) \ ,\end{align}

(4.15) \begin{align}\mathcal{Q}_{AD}^\varepsilon(\mathbf{x}_A,t,\mathbf{v}_A) &= \mathcal{K}_{AD}^\varepsilon(\mathbf{x}_A,t,\mathbf{v}_A) + \varepsilon^\kappa \, {}_{\varepsilon}\mathcal{J}^A_{ID}(\mathbf{x}_A,t,\mathbf{v}_A)\end{align}

and

(4.16) \begin{equation}\mathcal{Q}_{DI}^\varepsilon(\mathbf{x}_D,t,\mathbf{v}_D) = \mathcal{K}_{DI}^\varepsilon(\mathbf{x}_D,t,\mathbf{v}_D) \ .\end{equation}

We recall that in the case where cells move in a local search pattern (i.e. for $h=A$ and $h=D$ ), $\beta_h$ is defined via (2.1) and (2.2), and thus ${}_{\varepsilon}p^{\beta_h}_{h}(\mathbf{x}_h,t,\mathbf{v}_h)$ is given as in (3.15). On the other hand, in the case where cells move in a non-local search pattern (i.e. for $h=I$ ), $\beta_h$ is defined via (2.1) and (2.3), and thus ${}_{\varepsilon}p^{\beta_h}_{h}(\mathbf{x}_h,t,\mathbf{v}_h)$ is given by (3.16)

with

\begin{align*}\mathcal{B}^{\varepsilon}[p_{h}^\varepsilon](\mathbf{x}_h,t,\mathbf{v}_h) = \int_0^t B^{\varepsilon}(\mathbf{x}_h,t-s)p_{h}^\varepsilon(\mathbf{x}_h-(c\mathbf{v}_h+b)(t-s),s)\mathop{}\!\mathrm{d} s \ .\end{align*}

As before, $B^{\varepsilon}$ is defined through its Laplace transform in time $\hat{B}^{\varepsilon}$ and, in particular, under assumptions (4.2), we make the approximation:

\begin{equation*}\hat{B}^\varepsilon(\mathbf{x}_h,\varepsilon\lambda+\varepsilon^\mu b+\varepsilon^{1-\gamma}c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})\simeq \hat{B}^\varepsilon(\mathbf{x}_h,\varepsilon^{1-\gamma}c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h}) \ .\end{equation*}

Using the properties of the Laplace transform of a convolution, we write

\begin{equation*}\int_0^t B^\varepsilon(\mathbf{x}_h,t-s)p_{h}^\varepsilon(\mathbf{x}_h-(c\mathbf{v}_h+b)(t-s),s,\mathbf{v}_h)\mathop{}\!\mathrm{d} s\simeq \hat{B}^\varepsilon(\mathbf{x}_h,\varepsilon^{1-\gamma}c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})p_h^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)\ ,\end{equation*}

with

(4.17) \begin{equation} \hat{B}^\varepsilon(\mathbf{x}_h,\varepsilon^{1-\gamma}c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})=\frac{\hat{\varphi}_h^\varepsilon(\mathbf{x}_h,\varepsilon^{1-\gamma}c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})}{\hat{\psi}_h^\varepsilon(\mathbf{x}_h,\varepsilon^{1-\gamma}c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})} \ .\end{equation}

Analogous calculations are fully detailed in 8. Substituting the expressions of $\hat{\varphi}_{h}^\varepsilon$ and $\hat{\psi}_{h}^\varepsilon$ into (4.17), calculations similar to those carried out in [Reference Estrada-Rodriguez and Gimperlein11,Reference Estrada-Rodriguez, Gimperlein and Painter12] allow one to show that

(4.18) \begin{align} \hat{B}^\varepsilon&(\mathbf{x}_h,\varepsilon^{1-\gamma}c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})= \frac{\alpha-1}{d_{\varepsilon}}-\frac{\varepsilon^{1-\gamma}}{2-\alpha}c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h}\nonumber\\[5pt] &-d_{\varepsilon}^{\alpha-2}\varepsilon^{(1-\gamma)(\alpha-1)}(c\mathbf{v}_h\cdot\nabla_{\mathbf{x}_h})^{\alpha-1}(\alpha-1)^2\Gamma(-\alpha+1)+\mathcal{O}(d_{\varepsilon}^{\alpha-1}\lambda^\alpha)\ \ .\end{align}

In (4.18), $d_{\varepsilon}(\mathbf{x}_h) \;:=\; \tau_0(\mathbf{x}_h) \, \varepsilon^{\mu}$ , where $\tau_0(\mathbf{x}_h)$ is defined via (2.3).

Transport equations for $\rho_{I}^\varepsilon$ , $\rho_{A}^\varepsilon$ and $\rho_{D}^\varepsilon$

Integrating both sides of transport equation (4.13) with respect to ${\bf v}_h$ over the set V and using the fact that the turning operator $\mathcal{T}_h$ satisfies (2.6), we find that the macroscopic cell density $\rho_{h}^\varepsilon(\mathbf{x}_h,t)$ given by (4.12) satisfies the following transport equation:

(4.19) \begin{align} \partial_t \rho_{h}^\varepsilon + 2c \, \nabla_{\mathbf{x}_h} \cdot w_{h}^\varepsilon &= \varepsilon^{-1} \int_{\rm V} \mathcal{Q}_{hk}^\varepsilon \mathop{}\!\mathrm{d}\mathbf{v}_h \ , \quad \mathbf{x}_h \in \mathbb{R}^n, t \in \mathbb{R}^*_+ \ . \end{align}

Moreover, substituting the expressions for $p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)$ and $p_{k}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_k)$ given by (4.11) into the definitions of $\mathcal{K}_{hk}^\varepsilon$ , $ \mathcal{J}_{hk}^\varepsilon$ and ${}_{\varepsilon}\mathcal{J}^l_{hk}$ given by (4.8)–(4.10), we find

(4.20) \begin{equation} \int_{\rm V} \mathcal{K}_{hk}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h) \mathop{}\!\mathrm{d}\mathbf{v}_h = 0 \end{equation}

and, neglecting higher-order terms, we also obtain

(4.21) \begin{equation}\int_{\rm V} \mathcal{J}_{hk}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h) \mathop{}\!\mathrm{d}\mathbf{v}_h = -\varepsilon^{\xi -\vartheta - \gamma} \, c \, M \, \rho_{h}^\varepsilon \, \rho_{k}^\varepsilon\ , \quad \int_{\rm V} {}_{\varepsilon}\mathcal{J}^h_{lk}(\mathbf{x}_h,t,\mathbf{v}_h) \mathop{}\!\mathrm{d}\mathbf{v}_h = \varepsilon^{\xi -\vartheta - \gamma} \, c \, M \, \rho_{l}^\varepsilon \, \rho_{k}^\varepsilon \ , \end{equation}

where M is given by:

(4.22) \begin{equation}M\;:=\;\frac{1}{|{\rm V}|^3}\int_{\rm V}\int_{\rm V}\int_{{\rm V}^+_{hk}}\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\mathop{}\!\mathrm{d}\nu\mathop{}\!\mathrm{d}\mathbf{v}_h\mathop{}\!\mathrm{d}\mathbf{v}_k\ ,\quad h, k \in \{A, D, I \}, \ h \neq k \ .\end{equation}

Notice that relation (4.19) is obtained using property (3.22).

In conclusion, using (4.19) and (4.20) along with (4.14)–(4.16), from transport equation (4.2) we obtain the following equations for the macroscopic cell densities $\rho_{I}^\varepsilon(\mathbf{x}_I,t)$ , $\rho_{A}^\varepsilon(\mathbf{x}_A,t)$ and $\rho_{D}^\varepsilon(\mathbf{x}_D,t)$ :

(4.23) \begin{align} \partial_t\rho_{I}^\varepsilon+ 2c \, \nabla_{\mathbf{x}_I}\cdot w_{I}^\varepsilon&=- \, c \, M \, \rho_{I}^\varepsilon\rho_{D}^\varepsilon \ ,&& \mathbf{x}_I \in \mathbb{R}^2, t \in \mathbb{R}^*_+ \ ,\end{align}

(4.24) \begin{align} \partial_t\rho_{A}^\varepsilon + 2c \, \nabla_{\mathbf{x}_A}\cdot w_{A}^\varepsilon &= \, c \, M \, \rho_{I}^\varepsilon \rho_{D}^\varepsilon \ , && \mathbf{x}_A \in \mathbb{R}^2, t \in \mathbb{R}^*_+ \ ,\end{align}

(4.25) \begin{align} \partial_t\rho_{D}^\varepsilon + 2c \, \nabla_{\mathbf{x}_D}\cdot w_{D}^\varepsilon&=0 \ , && \mathbf{x}_D \in \mathbb{R}^2, t \in \mathbb{R}^*_+ \ .\end{align}

Here, we have used the scaling relations in (4.4) for the parameter $\kappa$ . On the right-hand side of (4.22), we have the density of cells that are leaving the state I (due to interactions with cells in the population D) and are appearing in the new state A in (4.23).

Transport equations for $w_{I}^\varepsilon$ , $w_{A}^\varepsilon$ and $w_{D}^\varepsilon$

Multiplying both sides of transport equation (4.13) by ${\bf v}_h$ and then integrating both sides of the resulting equation with respect to ${\bf v}_h$ over the set V, we find that the local macroscopic direction of cell motion $w_{h}^\varepsilon(\mathbf{x}_h,t)$ given by (4.12) satisfies the following transport equation:

(4.26) \begin{align} \varepsilon^{1+\gamma}2\partial_tw_{h}^\varepsilon +\varepsilon^{1-\gamma} c\nabla_{\mathbf{x}_h}\int_{\rm V}\mathbf{v}_h\otimes\mathbf{v}_h \, {p}_{h}^\varepsilon\mathop{}\!\mathrm{d}\mathbf{v}_h = &- \int_{\rm V} \mathbf{v}_h (\unicode{x1D7D9}-\mathcal{T}_h)\left[{}_{\varepsilon}p^{\beta_h}_{h}\right] \mathop{}\!\mathrm{d}\mathbf{v}_h \nonumber \\[5pt] & + \int_{\rm V} \mathbf{v}_h \, \mathcal{Q}_{hk}^\varepsilon \mathop{}\!\mathrm{d}\mathbf{v}_h \ . \end{align}

In the case where $\beta_h$ is defined via (2.1) and (2.2), using (4.11) and (3.15), and the properties of the turning operator $\mathcal{T}_h$ established by Lemma 1, we find that the first term on the right-hand side of (4.2) is given by:

(4.27) \begin{equation}\int_{\rm V} \mathbf{v}_h (\unicode{x1D7D9}-\mathcal{T}_h)\left[{}_{\varepsilon}p^{\beta_h}_{h}\right] \mathop{}\!\mathrm{d}\mathbf{v}_h= \frac{2\varepsilon^\gamma}{|{\rm V}|} b{(1-\iota_1)}w_h^\varepsilon \ . \end{equation}

Here, $\iota_1(\mathbf{x}_h,t)$ is the first non-zero eigenvalue of the turning operator $\mathcal{T}_h$ , which is given by (2.8). On the other hand, when $\beta_h$ is defined via (2.1) and (2.3), using (4.11) and (3.16) and the properties of the turning operator $\mathcal{T}_h$ established by Lemma 1, it was proved in [Reference Estrada-Rodriguez, Gimperlein and Painter12] that the following approximate expression of the first term on the right-hand side of (4.2) holds

(4.28) \begin{align}\int_{\rm V} &\mathbf{v}_h (\unicode{x1D7D9}-\mathcal{T}_h)[{}_{\varepsilon}p^{\beta_h}_{h}] \mathop{}\!\mathrm{d}\mathbf{v}_h = \varepsilon^{1-\frac{\gamma}{\alpha-1}} \left(g_\alpha\nabla^{\alpha-1}\rho^\varepsilon_h{-}\frac{2(\alpha-1)}{\tau_0|{\rm V}|} (\iota_1-1) w_h^\varepsilon\right) +\text{l.o.t.}\ , \end{align}

where

(4.29) \begin{equation}g_\alpha(\mathbf{x}_h,t)\;:=\;\frac{\pi\tau_0^{\alpha-2}(1-\alpha)^2c^{\alpha-1}}{\sin(\pi\alpha)\Gamma(\alpha)}\left(\frac{4\iota_1-|{\rm V}|}{|{\rm V}|}\right)\ \ \text{for}\ \ \Gamma(-\alpha+1)=\frac{\pi}{\sin(\pi\alpha)\Gamma(\alpha)} \ .\end{equation}

Notice that $g_\alpha(\cdot,\cdot)>0$ since $\sin(\pi\alpha)<0$ for $\alpha\in(1,2)$ and $4\iota_1-|{ {\rm V}}|<0$ by using (2.7) for $n=2$ and recalling that $\iota_1<1$ .

Moreover, as similarly done in [Reference Estrada-Rodriguez and Gimperlein11], using the fact that $(\cdot)':{\rm V}\mapsto {\rm V}$ is a bijection and $\mathbf{v}_h'\cdot\nu=-\mathbf{v}_h\cdot\nu$ , whence $\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)=-\nu\cdot(\mathbf{v}_h'-\mathbf{v}_k')$ , we find

(4.30) \begin{align}\int_{\rm V} & \mathbf{v}_h \, \mathcal{K}_{hk}^\varepsilon \, \mathop{}\!\mathrm{d}\mathbf{v}_h \nonumber \\[5pt]& = \frac{1}{|{{\rm V}}|} \varepsilon^{\xi -\vartheta - \gamma} \, c \, \left(\int_{\rm V} \int_{\rm V} \int_{{\rm V}^+_{hk}} {\bf v}_h \, p_{h}^\varepsilon(\mathbf{x}_h,t,{\bf v}_h') p_{k}^\varepsilon(\mathbf{x}_h,t,{\bf v}_k') \nu\cdot(\mathbf{v}_h-\mathbf{v}_k) \mathop{}\!\mathrm{d}\nu \mathop{}\!\mathrm{d}\mathbf{v}_k \mathop{}\!\mathrm{d}\mathbf{v}_h \right. \nonumber\\[5pt] &\left.- \int_{\rm V} \int_{\rm V} \int_{{\rm V}^+_{hk}} {\bf v}_h \, p_{h}^\varepsilon(\mathbf{x}_h,t,{\bf v}_h) p_{k}^\varepsilon(\mathbf{x}_h,t,{\bf v}_k) \nu\cdot(\mathbf{v}_h-\mathbf{v}_k) \, \mathop{}\!\mathrm{d}\nu\mathop{}\!\mathrm{d}\mathbf{v}_k \mathop{}\!\mathrm{d}\mathbf{v}_h \right) \nonumber\\[5pt]&=- \frac{1}{|{{\rm V}}|} \varepsilon^{\xi -\vartheta - \gamma} \, c \, \int_{\rm V}\int_{\rm V}\int_{{\rm V}^+_{hk}}(\mathbf{v}'_h)'p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h')p_{k}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_k')\nu\cdot(\mathbf{v}_h'-\mathbf{v}_k')\mathop{}\!\mathrm{d}\nu\mathop{}\!\mathrm{d}\mathbf{v}_h'\mathop{}\!\mathrm{d}\mathbf{v}_k'\nonumber\\[5pt] &\ \ \ -\int_{\rm V}\int_{\rm V}\int_{{\rm V}^+_{hk}}\mathbf{v}_h p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)p_{k}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_k)\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\mathop{}\!\mathrm{d}\nu\mathop{}\!\mathrm{d}\mathbf{v}_h\mathop{}\!\mathrm{d}\mathbf{v}_k\nonumber{}\\[5pt] &= \frac{1}{|{{\rm V}}|} \varepsilon^{\xi -\vartheta - \gamma} \, c \, \int_{\rm V}\int_{\rm V}\int_{{\rm V}^+_{hk}}(\mathbf{v}_h'-\mathbf{v}_h)p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)p_{k}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_k)\nu\cdot(\mathbf{v}_h-\mathbf{v}_k)\mathop{}\!\mathrm{d}\nu\mathop{}\!\mathrm{d}\mathbf{v}_h\mathop{}\!\mathrm{d}\mathbf{v}_k\nonumber\\[5pt] & =- \frac{1}{|{{\rm V}}|} \varepsilon^{\xi -\vartheta - \gamma} \, c \, \frac{4}{3}\int_{\rm V}\int_{\rm V}|\mathbf{v}_h-\mathbf{v}_k| \, \mathbf{v}_h \, p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h) \, p_{k}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_k)\mathop{}\!\mathrm{d}\mathbf{v}_h\mathop{}\!\mathrm{d}\mathbf{v}_k\ . \end{align}

The last equality in (4.28) is obtained using the fact that $\mathbf{v}_h'-\mathbf{v}_h=-2(\mathbf{v}_h\cdot\nu)\nu$ . Substituting the expressions for $p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)$ and $p_{k}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_k)$ given by (4.11) into (4.28) and into definitions (4.9) and (4.10) of $ \mathcal{J}_{hk}^\varepsilon$ and ${}_{\varepsilon}\mathcal{J}^h_{lk}$ , neglecting higher-order terms we find

(4.31) \begin{equation} \int_{\rm V} \mathbf{v}_h \, \mathcal{K}_{hk}^\varepsilon \mathop{}\!\mathrm{d}\mathbf{v}_h = - \varepsilon^{\xi -\vartheta} \, c \, \frac{8}{3} \, \dfrac{1}{|{\rm V}|^3} \, q_h \, \rho_k^\varepsilon w_h^\varepsilon\ , \end{equation}

and

(4.32) \begin{equation}\int_{\rm V} \mathbf{v}_h \, \mathcal{J}_{hk}^\varepsilon \mathop{}\!\mathrm{d}\mathbf{v}_h = - \frac{\varepsilon^{\xi -\vartheta}2c}{|{\rm V}|^2} \tilde{q}_h\rho_k^\varepsilon w_h^\varepsilon \ , \quad \int_{\rm V} \mathbf{v}_h \, {}_{\varepsilon}\mathcal{J}^h_{lk} \mathop{}\!\mathrm{d}\mathbf{v}_h = \frac{\varepsilon^{\xi -\vartheta}2c}{|{\rm V}|^2} \tilde{q}_h\rho_k^\varepsilon w_h^\varepsilon \ , \end{equation}

where $q_h$ and $\tilde{q}_h$ are defined as:

(4.33) \begin{equation}q_h\;:=\;\int_{{\rm V}}|\mathbf{v}_h-\mathbf{v}_k|\mathop{}\!\mathrm{d}\mathbf{v} \ , \quad { \tilde{q}_h\;:=\;\int_{{\rm V}^+_{hk}}|\nu||\mathbf{v}_h-\mathbf{v}_k|\mathop{}\!\mathrm{d} \nu}\ , \quad h \in \{A, D, I\} \ .\end{equation}

Finally, substituting the expression of $p_{h}^\varepsilon(\mathbf{x}_h,t,\mathbf{v}_h)$ into the second term on the left-hand side of (4.2) and neglecting higher-order terms yields

(4.34) \begin{equation}c \, \nabla_{\mathbf{x}_h}\int_{\rm V}\mathbf{v}_h\otimes\mathbf{v}_h \, {p}_{h}^\varepsilon \mathop{}\!\mathrm{d}\mathbf{v}_h = C_h \, \nabla_{\mathbf{x}_h} {\rho}_{h}^\varepsilon \ , \end{equation}

with $C_h$ being defined as:

(4.35) \begin{equation}C_h \;:=\; \dfrac{c}{|{\rm V}|}\int_{\rm V}\mathbf{v}_h\otimes\mathbf{v}_h\mathop{}\!\mathrm{d}\mathbf{v}_h \ , \quad h \in \{A, D, I\} \ .\end{equation}

In conclusion, using (3.15) and (3.16), (4.29), (4.30) and (4.32) along with (4.14)–(4.16), from transport equation (4.2) we obtain the following transport equations for the local macroscopic directions of cell motion $w_{I}^\varepsilon(\mathbf{x}_I,t)$ , $w_{A}^\varepsilon(\mathbf{x}_A,t)$ and $w_{D}^\varepsilon(\mathbf{x}_D,t)$ :

(4.36) \begin{align} &\varepsilon^{1+\gamma}2\partial_tw_{I}^\varepsilon+\varepsilon^{1-\gamma} C_I \, \nabla_{\mathbf{x}_I} \rho_{I}^\varepsilon ={ -} \varepsilon^{1-\frac{\gamma}{\alpha-1}}\left(g_\alpha\nabla^{\alpha-1}\rho_{I}^\varepsilon{ -}\frac{2(\alpha-1)}{\tau_0|{\rm V}|} (\iota_1-1) w_{I}^\varepsilon\right)\nonumber\\[5pt] & - \varepsilon^{\unicode{x03BE}-\vartheta} \, c \, \left((1-\varepsilon^{\kappa})\frac{8}{3}\frac{1}{|{\rm V}|^3} q_I \, \rho_{D}^\varepsilon w_{I}^\varepsilon +{ \frac{2c}{|{\rm V}|^2}\tilde{q}_I\rho_D^\varepsilon w_I^\varepsilon}\right) \ , \quad \mathbf{x}_I \in \mathbb{R}^2, t \in \mathbb{R}^*_+ \ , \end{align}

(4.37) \begin{align} \varepsilon^{1+\gamma}2\partial_tw_{A}^\varepsilon+\varepsilon^{1-\gamma} &C_A \, \nabla_{\mathbf{x}_A} \rho_{A}^\varepsilon = -\frac{2\varepsilon^\gamma}{|{\rm V}|} b{ (1-\iota_1)}w_{A}^\varepsilon \nonumber\\[5pt] & - \varepsilon^{\unicode{x03BE}-\vartheta} \, c \, \left(\frac{8}{3}\frac{1}{|{\rm V}|^3} q_A \, \rho_{D}^\varepsilon w_{A}^\varepsilon- {\frac{2c}{|{\rm V}|^2}\tilde{q}_A\rho_D^\varepsilon w_A^\varepsilon}\right) \ , \quad \mathbf{x}_A \in \mathbb{R}^2, t \in \mathbb{R}^*_+ \ , \end{align}

and

(4.38) \begin{align} \varepsilon^{\gamma+1}2\partial_t w_{D}^\varepsilon +\varepsilon^{1-\gamma} C_D \, \nabla_{\mathbf{x}_D} & \rho_{D}^\varepsilon = -\frac{2\varepsilon^{\gamma}}{|{\rm V}|} b(\iota_1-1)w_{D}^\varepsilon\nonumber\\[5pt] & - \varepsilon^{\unicode{x03BE}-\vartheta} \, c \, \frac{8}{3} \frac{1}{|{\rm V}|^3} q_D \, \rho_{I}^\varepsilon w_{D}^\varepsilon \ , \quad \mathbf{x}_D \in \mathbb{R}^2, t \in \mathbb{R}^*_+ \ .\end{align}

Macroscopic scale model

Noting that $1-\dfrac{\gamma}{\alpha-1} < 1 -\gamma$ since $\alpha \in (1,2)$ and using assumptions (4.4) for the scaling parameters, letting $\varepsilon \to 0$ in (4.34)–(4.36), we formally find the following expressions for the leading-order terms $w_I({\bf x},t)$ , $w_A({\bf x},t)$ and $w_D({\bf x},t)$ of the asymptotic expansions for the local macroscopic directions of cell motion $w_{I \varepsilon}({\bf x},t)$ , $w_{A \varepsilon}({\bf x},t)$ and $w_{D \varepsilon}({\bf x},t)$ :

• – From (4.34) and choosing the scaling parameters as $1-\frac{\gamma}{\alpha-1}=\xi-\vartheta$ , we have

  (4.39) \begin{equation} w_I=-\frac{g_{\alpha}}{H_I}\nabla_{\mathbf{x}}^{\alpha-1}\rho_I\ , \quad \text{where} \quad H_I \;:=\;\frac{2(\alpha-1)(1-\iota_1)}{\tau_0|{\rm V}|}+\frac{8c}{|3{\rm V}|^3}q_I\rho_D>0\ . \end{equation}

• – From (4.35) and choosing $\xi-\vartheta=\gamma$ in agreement with (4.4), we obtain

  (4.40) \begin{equation} w_A=-\frac{C_A}{H_A}\nabla_{\mathbf{x}}\rho_A\ , \quad \text{where}\quad H_A \;:=\;\frac{2b}{|{\rm V}|}(1-\iota_1)+\frac{8c}{3|{\rm V}|^3}q_A\rho_D>0\ . \end{equation}

• – Finally, using the same scaling rules as in the previous case, we obtain from (4.36)

  (4.41) \begin{equation} w_D=-\frac{C_D}{H_D}\nabla_{\mathbf{x}}\rho_D\ ,\quad \text{where}\quad H_D \;:=\;\frac{2b}{|{\rm V}|}(1-\iota_1)+\frac{8c}{3|{\rm V}|^3}q_D\rho_I>0\ .\end{equation}

Furthermore, under assumptions (4.4), letting $\varepsilon \to 0$ in (4.22)–(4.24) and using(4.37)–(4.39), we formally obtain the following balance equations for the leading-order terms $\rho_I({\bf x},t)$ , $\rho_A({\bf x},t)$ and $\rho_D({\bf x},t)$ of the asymptotic expansions for the macroscopic cell densities $\rho_{I}^\varepsilon({\bf x},t)$ , $\rho_{A}^\varepsilon({\bf x},t)$ and $\rho_{D}^\varepsilon({\bf x},t)$ :

(4.42) \begin{align} \partial_t\rho_I-\nabla_{\mathbf{x}}\cdot\Bigl(D_I \, \nabla^{\alpha-1}_{\mathbf{x}}\rho_I \Bigr)&=- a \, \rho_I \,\rho_D \ , \quad \alpha \in (1,2) \ , && \mathbf{x} \in \mathbb{R}^2, t \in \mathbb{R}^*_+\ ,\end{align}

(4.43) \begin{align} \partial_t\rho_A-\nabla_{\mathbf{x}}\cdot\Bigl(D_A \, \nabla_{\mathbf{x}}\rho_A \Bigr)&= a \, \rho_I \,\rho_D \ , && \mathbf{x} \in \mathbb{R}^2, t \in \mathbb{R}^*_+ \ ,\end{align}

(4.44) \begin{align} \partial_t\rho_D-\nabla_{\mathbf{x}}\cdot \Bigl(D_D \, \nabla_{\mathbf{x}}\rho_D\Bigr)&=0\ , && \mathbf{x} \in \mathbb{R}^2, t \in \mathbb{R}^*_+ \ , \end{align}

where

\begin{align*}D_I \;:=\; \frac{2 \, c \, g_\alpha}{H_I}, \quad D_A \;:=\; \frac{2 \, c \, C_A}{H_A}, \quad D_D \;:=\; \frac{2 \, c \, C_D}{H_D}, \quad a \;:=\; c \, M.\end{align*}

Considerations on the macroscopic scale model (4.40)–(4.42)

The functions $\rho_I({\bf x},t)$ , $\rho_A({\bf x},t)$ and $\rho_D({\bf x},t)$ model, respectively, the density of CTLs in the pre-activation state, activated CTLs and DCs presenting a tumour antigen on their surface at position x and time t. The spatio-temporal coevolution of CTLs and DCs is modelled through the coupled system of balance equations (4.40)–(4.42), which governs the dynamics of the cell density functions.

The mathematical model defined by (4.40)–(4.42) provides a macroscopic description of cell dynamics that takes explicitly into account the effects of cell–cell interactions and the characteristics of cell motion that are encapsulated in the parameters c (i.e. the magnitude of the cell velocity, which is assumed to be constant), $\alpha$ (i.e. the characteristic exponent of the long-tailed distribution followed by the running time of CTLs in the pre-activation state) and b (i.e. the characteristic exponent of the Poisson distribution followed by the running time of activated CTLs and DCs).

This model effectively captures the fact that interactions between DCs presenting a tumour antigen on their surface and CTLs in the pre-activation state lead to CTL activation. In particular, the term on the right-hand side of (4.40) models the decay in the density of CTLs in the pre-activation state at position x and time t due to contact interactions with DCs which result in CTL activation, while the term on the right-hand side of (4.41) models the corresponding growth in the density of activated CTLs. As one would expect, these two terms differ only in their signs and are proportional to the product between the cell density functions $\rho_{I}({\bf x},t)$ and $\rho_{D}({\bf x},t)$ . The factor of proportionality a increases with the value of the parameter c. This is coherent with the observation that higher cell motilities may increase the encounter rate of CTLs in the pre-activation state with DCs.

The model captures also the fact that CTLs in the pre-activation state move in a non-local search pattern, while the search pattern of activated CTLs is more localised. In fact, the rate of change of the density of CTLs in the pre-activation state due to cell movement (i.e. the second term on the left-hand side of (4.40)) is a fractional diffusion term, while that of the density of activated CTLs (i.e. the second term on the left-hand side of (4.41)) is a classical diffusion term. The function modelling the diffusivity of CTLs in the pre-activation state (i.e. the function $D_I$ ) and the function modelling the diffusivity of activated CTLs (i.e. the function $D_A$ ) are proportional to the parameter c. This is coherent with the observation that, ceteris paribus, a higher magnitude of the cell velocity correlates with a higher cell motility. Both $D_I$ and $D_A$ are monotonically decreasing functions of the cell density function $\rho_{D}({\bf x},t)$ , which means that, all else being equal, the higher the density of DCs at a given position, the lower the diffusivity of CTLs. This reflects the fact that higher densities of DCs will make it more likely that interactions between CTLs and DCs occur and, since these interactions force CTLs to change their direction of movement at the mesoscopic scale, this will ultimately result in a lower cell diffusivity at the macroscopic scale. Moreover, $D_A$ is an increasing function of b. This is coherent with the fact that larger values of this parameter correspond to larger mean values of the cell running times.

The fact that the right-hand side of (4.42) is zero translates in mathematical terms the idea that we are not taking into account the effects of division and death of DCs. Moreover, coherently with the fact that the motion of DCs is here described as a Brownian motion, the rate of change of the density of DCs due to cell movement (i.e. the second term on the left-hand side of (4.42)) is a classical diffusion term. Considerations analogous to those made above about the dependence of the function $D_A$ on the parameters c and b apply to the function modelling the diffusivity of DCs (i.e. the function $D_D$ ) as well. Furthermore, considerations similar to those made above about the dependence of $D_I$ on the density function $\rho_{D}({\bf x},t)$ hold for the dependence of $D_D$ on the density function $\rho_{I}({\bf x},t)$ .