As motivated by studies of cellular motility driven by  spatiotemporal chemotactic gradients in microdevices, we develop a framework for constructing approximate analytical solutions for the location, speed and cellular densities for cell chemotaxis waves in heterogeneous fields of chemoattractant from the underlying partial differential equation models. In particular, such chemotactic waves are not in general translationally invariant travelling waves, but possess a spatial variation that evolves in time, and may even may oscillate back and forth in time, according to the details of the chemotactic gradients. The analytical framework exploits the observation that unbiased cellular diffusive flux is typically small compared to chemotactic fluxes and is first developed and validated for a range of exemplar scenarios. The framework is subsequently applied to more complex models considering the full dynamics of the chemoattractant and how this may be driven and controlled within a microdevice by considering a range of boundary conditions. In particular, even though solutions cannot be constructed in all cases, a wide variety of scenarios can be considered analytically, firstly providing global insight into the important mechanisms and features of cell  motility in complex spatiotemporal fields of chemoattractant. Such analytical solutions also provide a means of rapid evaluation of model predictions, with the prospect of application in computationally demanding investigations relating theoretical models and experimental observation, such as Bayesian parameter estimation.

§ INTRODUCTION

Most biological processes integrate different cell populations, extracellular matrix (ECM) properties, chemotactic gradients and physical cues,  constituting  a  complex, dynamic and interactive microenvironment [1, 2, 3, 4, 5]. Cells are also constantly adjusting   to accommodate  their surroundings,  particularly for homeostatically maintaining the intracellular and extracellular environment within  physiological constraints [6]. In response to the different chemical and physical external stimuli, cells can modify their shape, location, internal structure and genomic expression, as well as their capacity to proliferate, migrate, differentiate, produce ECM or other biochemical substances, changing, in turn, the surrounding medium as well as sending new signals to other cells [7, 8, 9]. This two-way interaction between cells and their environment is crucial in physiological processes such as embryogenesis, organ development, homeostasis, repair, and long-term evolution of tissues and organs among others, as well as in pathological processes such as atherosclerosis or cancer [10, 11, 12, 13]. Furthermore, developing novel frameworks to investigate and elucidate  these mechanisms and interactions is key to developing novel therapeutic strategies aiming at promoting (blocking) desirable (undesirable) cellular behaviours [14].

In particular, due to the underlying complexity, in vivo research – both in humans and in animals – is impeded by the fact it is difficult to control and isolate effects. Thus, a simpler alternative is using in vitro experiments. Nevertheless, the predictive power  currently available, whether in vivo or in vitro, is still poor, as demonstrated the continuous drop in the number of new drugs appearing annually, despite billion-dollar investments [15, 16]. Indeed, structural three-dimensionality is one of the most important characteristics of biological processes  [17], but in vitro cells are mostly cultured in a traditional Petri dish (2D culture), where cell behaviour is dramatically different from   real tissues [18]. Recently, microfluidics has arisen as a powerful tool to recreate the complex microenvironment that governs tumour dynamics [19, 20]. This technique allows the  reproduction of numerous important features that are lost in 2D cultures, as well as testing drugs in a much more reliable and efficient way [21, 22, 23, 24, 25].

In addition to such in vitro models, mathematical in silico models  are a powerful tool for dealing with many problems in physics, engineering, and biology. In particular, cell population evolution models based on transport partial differential equations (PDEs) have been widely used to study many biological processes, including cancer [26, 27]. For instance, tumour development is a key example of a highly dynamic and complex biological process that originates from external signals or stimuli modulated by the particular microenvironment. Furthermore, when a given treatment is applied (surgery, chemotherapy, radiotherapy, immunotherapy, hormones or a combination thereof), the tumour and its microenvironment undergo significant alterations. This leads tumour cells to proliferate and generate microenvironments that promote the death of surrounding cell types and the survival of tumour cells that are more adaptable and resistant. That is why, when modelling the enormous variety and complexity of a tumour and its microenvironment, the resulting differential equations are highly non-linear and strongly coupled [28, 29, 30]. The numerical resolution of the equations, especially in the era of high performance computing, has been extensively utilised in the simulation of  “what if” scenarios and the study of effects and hypotheses in isolation, something that is often impossible to do with in vivo and in vitro models [31, 32]. In turn, the construction and exploitation  of in silico experiments is thus being increasingly used in the early stages of designing drugs and therapies against tumours.

A particular niche of interest is  Glioblastoma (GBM), the most common and aggressive primary brain tumour [33], with  extensive studies dedicated to    mathematical modelling  its evolution [34],   reproducing aspects of GBM  histopathology [29] and incorporating the influence of  tumour microenvironment (TME) chemical and mechanical cues   [35]. It has been demonstrated that GBM progression is extensively  controlled by the local oxygen concentrations and gradients [36], motivating many studies to incorporate  the role of oxygen gradients and hypoxia in tumour progression [37, 38, 39]. Some models of GBM have  reproduced cell culture evolution under different experimental configurations [40], using a go-or-grow transition switch, governed  by  nonlinear activation functions for the chemotaxis and growth. Such studies therefore implicate the balance between cell migratory and proliferative activity, together with  their relation to the different TME stimuli, as playing a key role in GBM evolution.

Nevertheless, the complexity of the equations to be solved often require numerical simulations that are  impractical, due to the high computational cost, especially in the resolution of inverse problems such as parameter estimation, model selection, the design of experiments, sensitivity studies, model structural analysis and Uncertainty Quantification (UQ). Although many modern techniques as Reduced Order Models (ROM) and metamodels using Artificial Intelligence (AI) have been developed in recent years [41], the existence of analytical solutions, although approximate, provides key information to test and validate numerical algorithms, inform a mechanism based understanding across parameter space and to  allow initial predictions of Quantities of Interest (QoI),  such as travelling  fronts, equilibria, the ranges of variation in the solution across parameter space and parameter sensitivities, among others. Indeed, some works in the last years have focused on the use of these techniques for analysing GBM progression [42, 43, 44].

The interaction of cells with a chemoattractant leads to a type of Keller - Segel (K-S) model [45], which generally have a rich structure as reviewed by G. Arumugam and J. Tyagi [46]. One of the main interests concerning the K-S model is the existence and characteristics of travelling waves (see for instance [47]). In this work, we explore the dynamics of cell populations under gradients of a chemotactic agent for one-dimensional problems. We move beyond travelling waves to investigate evolving wave solutions in the heterogeneous environments that are often found in microdevices and physiology. This general class of problems allows the treatment of a wide variety of situations related to the evolution of tumours, while the  analysis of the associated PDEs enables the quantification of histopathological characteristics, such as the spread of pseudopalisades and the response of the population to oscillatory stimuli. This knowledge can be used for the design of experiments, to speed up the characterisation processes of cell populations and to validate or rule out possible models.

In particular, we are interested in modelling cell motility dynamics in microfluidic devices, which are experimental platforms that have been demonstrated to accurately recreate biomimetic physiological conditions [48], with application  in bioengineering and biomedical research [19]. In many situations, the chemotactic agent concentration may be assumed as known, either because it can be directly measured, or because its concentration can be computed by solving a diffusion problem that is, to good approximation, decoupled from the cell population field.

To proceed, we first describe the structure of the mathematical problem associated with the response of cell populations to chemotactic gradients, together with  the general assumptions and hypotheses about the underlying mechanisms. We derive pertinent features about the solution field, for instance that it possesses a migratory structure with a transition zone wavefront. We are also able to estimate the wave front evolution and the shape of the solution profile. In particular, we compute an analytical solution for  specific exemplar cases associated with specific relevant experimental situations. These include a constant spatial gradient of chemoattractant, together with temporal oscillations associated with a fluctuating source, quadratic profiles of chemoattractant, offering additional nonlinear features, and an exponential profile of chemoattractants corresponding to the diffusion from a localised source. We further apply the general results to the analysis of a cell culture microfluidic experiment, representing a slightly simplified version for an  in vitro model of GBM progression, as developed in [40], showing how the methods presented here can generate analytical results for the simulation of microdevice representations for migratory tumour cell dynamics.

§ METHODS

§.§ The model

We study a broad class of problems that are related to the dynamics  of a cell culture in microfluidic devices under the influence of a chemotattractant, when the concentration of the agent can be computed or measured. A schematic view of this situation is represented in Fig. <ref>.

Scheme of the experimental configuration.

1D approximation of the cell culture.

[Typical experimental configuration for modelling cell cultures]Typical experimental configuration for modelling cell cultures. Due to the much larger length of the lateral channels relative to the width of the chamber, the domain geometry of the model is assumed one-dimensional, with length given by the length of the chamber, $L$. The cell concentrations are associated with a continuum field $u = u(x,t)$, with $t$ denoting time and $x$ the spatial coordinate along the chamber, as indicated. At the  lateral edges of the channel, that is $x=0,L$, zero flux boundary conditions are imposed, corresponding to the inability of cells to pass through these boundaries. Image created with BioRender.com.

The non-dimensional equation for the cell population concentration, $u$, represents cellular  diffusion and chemotaxis in response to a heterogeneous field of chemoattractant that  generates an advective flux  $\alpha(t,x)u$, and  also impacts cellular proliferation via $\beta(t,x)$ so as to  generate the governing equation

\begin{equation} \label{eq:eq_governing}
u_t + \left(\alpha(t,x)u\right)_x = D u_{xx} + \beta(t,x)u(1-u),
\end{equation}

with $D>0$ the non-dimensional cellular diffusion coefficient. This governing equation is also supplemented by zero-flux boundary conditions, given by

\begin{align}
\left.Du_x-\alpha(t,x)u\right|_{x=0} &= 0, \\
\left.Du_x-\alpha(t,x)u\right|_{x=L} &= 0,
\end{align}

and initial conditions

\begin{equation} \label{eq:eq_initial}
\end{equation}

§.§ Computation of the general solution for small diffusion

§.§.§ Outer solution

The main hypothesis, which is usually true for cellular motility due to weak cell-based random motility, is that the non-dimensional diffusion coefficient satisfies $D \ll 1$, as we will verify below in the case of GBM cells in microdevices. Hence, away from boundary layers,  diffusion may be neglected  compared to the influence of  growth and chemoattractant driven migration. Then,  Eq. (<ref>) may be approximated by:

\begin{equation} \label{eq:eq_governing_simp}
u_t + (\alpha(t,x)u)_x = \beta(t,x)u(1-u).
\end{equation}

This is a first-order hyperbolic PDE, amenable to  the method of characteristics. If we know the initial data  $u(0,x) = u_0(x)$, we can parameterise the initial data via  $s$ with the relation  $(t,x,u) = (0,s,u_0(s))$. Also, there is another family of characteristic curves, emerging from the $(t,x)$ points where $x=0$ and $t>0$, with the imposition of the $x=0$ boundary condition of no flux,  Eq. (<ref>). Assuming that $\alpha(0,t) \neq 0$, and that the boundary is away from the transition region of the cellular wavefront, so that to excellent approximation $u_x= 0$ since spatial gradients are small,  we conclude from the boundary contition, Eq. (<ref>), that $u(0,t)=0$ to the same level of approximation. Therefore, this boundary condition can be parameterised via $s$ with the relation $(t,x,u) = (s,0,0)$. Hence, at $t=0$ and $x=0$ there is an emerging singular characteristic that splits the domain in two regions. The geometric interpretation of the method of characteristics is shown in Fig. <ref>, which shows the projection of the characteristic curves onto the plane $(t,x)$.

[black, ->] (0,0) – (0,5);
[left] at (0,5) $t$;
[black, ->] (-0.3,0) – (10,0);
[below] at (10,0) $x$;
[left] at (-0.3,0) $t_0$;
[black, dashed] (0,2) – (10,2);
[left] at (-0.3,2) $t^*$;
[black, dashed] (4,0) – (4,4);
(0,0) .. controls (2,3) and (3,0) .. (4,2);
(4,2) .. controls (5,4) and (6,5) .. (8,4);

[red, dashed] (0,1) .. controls (2,4) and (3,1) .. (4,3);
[red, dashed] (4,3) .. controls (5,5) and (6,6) .. (8,5);

[red, dashed] (0,2) .. controls (2,5) and (3,2) .. (4,4);
[red, dashed] (4,4) .. controls (5,6) and (6,7) .. (8,6);

[blue, dashed] (1,0) .. controls (3,2) and (3.5,0.5) .. (4,1.5);
[blue, dashed] (4,1.5) .. controls (5,3) and (6,3) .. (8,3);

[blue, dashed] (2,0) .. controls (3,1) and (3.5,0) .. (4,0.5);
[blue, dashed] (4,0.5) .. controls (5,3) and (6,3) .. (8,2);

[Projection of the characteristic curves]Projection of the characteristic curves. The two families of characteristic curves are shown in blue and red.

We have, therefore, two families of characteristic curves. For the first one:

\begin{align}
\frac{\mathrm{d}t}{\mathrm{d}\tau} &= 1, \quad t(0) = 0, \label{eq:cF1a_general} \\
\frac{\mathrm{d}x}{\mathrm{d}\tau} &= \alpha(t,x), \quad x(0) = s, \label{eq:cF1b_general}\\
\frac{\mathrm{d}u}{\mathrm{d}\tau} &= u\left(\beta(t,x)(1-u)-\alpha_x(t,x)\right), \quad u(0) = u_0(s), \label{eq:cF1c_general}
\end{align}

and for the second:

\begin{align}
\frac{\mathrm{d}t}{\mathrm{d}\tau} &= 1, \quad t(0) = s, \label{eq:cF2a_general} \\
\frac{\mathrm{d}x}{\mathrm{d}\tau} &= \alpha(t,x), \quad x(0) = 0, \label{eq:cF2b_general}\\
\frac{\mathrm{d}u}{\mathrm{d}\tau} &= u\left(\beta(t,x)(1-u)-\alpha_x(t,x)\right), \quad u(0) = 0. \label{eq:cF2c_general}
\end{align}

Noting the uniqueness of the solution to Eqs. (<ref>) courtesy of Picard's theorem,  we have by inspection that these equations only generate  the trivial solution $u(x,t)=0$.

However, the solution to the first family,
given by Eqs. (<ref>), is typically more complex, though one always has
\begin{equation} \label{eq:t_sol}
t = \tau.
\end{equation}

Progress can be readily made when

* $\alpha(t,x)$ is linear in $x$, such that $\alpha(t,x)=a(t)x+b(t)$

* $\alpha(t,x)$ is separable, that is $\alpha(t,x)=f(x)g(t)$.

In particular when  $\alpha(t,x)$ is linear in $x$ we have

\begin{eqnarray} \label{sollin} x = s \mathrm{e}^{\int_0^t a(\eta)\mathrm{d} \eta} +
\int_0^t b(\eta)\mathrm{e}^{\int_{\eta}^t a(\xi)\mathrm{d} \xi}   \, \mathrm{d} \eta =: F(t;s),
\end{eqnarray}

which defines $F(t;s)$ for  $\alpha(t,x)=a(t)x+b(t)$. We generalise this definition so that $x=F(t;s)$ on the characteristic curve given by the value of $s$, and whenever this relation can be uniquely inverted for $s$, we write $s=G(t;x)$.

In contrast when $\alpha(t,x)=f(x)g(t)$ Eqs. (<ref>) and (<ref>) may be integrated to obtain

\begin{equation} \label{eq:x_sol}
\int_s^x\frac{\mathrm{d}z}{f(z)} = \int_0^\tau g(\eta)\,\mathrm{d}\eta .
\end{equation}

In turn, Eq. (<ref>) generates the relation  $x=F(t;s)$ on the characteristic curve.
Further motivation for examples of linear and separable chemotactic response functions for $\alpha(t,x)$ are given in Section  <ref> below.

In both cases, or even in the more general case where $F(t;s)$ cannot readily be determined analytically for all relevant $t,s$, the location of the transition from $u=0$, and thus the location of the transition region for the cellular  wavefront, $x^* = x^*(t)$, is given by the characteristic with $s=0$ – hence   $x^*(t)=F(t;0)$. This is particularly  informative about the general behaviour of the solution, for instance in determining the wavespeed. With respect to Eq. (<ref>), we proceed using the change of variable $r = 1/u$, we obtain the ODE in terms of the $\tau = t$ variable:

\begin{equation} \label{eq:edo}
r'(\tau)+\left(\beta(\tau,x(\tau;s))-H(\tau;s)\right)r =\beta(\tau,x(\tau;s)),
\end{equation}

where for the linear case $H(\tau; s) = a(\tau)$ and $H(\tau;s) = f'(F(\tau;s))g(\tau)$ for the separable case.

Noting the integration is along a characteristic, and thus  $s$ is fixed,  this equation is of the form $r'+p(\tau)r=q(\tau)$ for $p(\tau) =\beta(\tau,x(\tau;s))-H(\tau;s)$ and $q(\tau) = \beta(\tau,x(\tau;s))$, with $s$ fixed, so a general expression is given by

\begin{equation} \label{eq:edo_sol}
r(\tau,s) = \exp\left(-\int_0^\tau p(\eta,s) \, \mathrm{d}\eta \right)\left[\frac{1}{u_0(s)} + \int_0^\tau q(\eta,s) \exp\left(\int_0^\eta p(\xi,s) \, \mathrm{d}\xi \right)\, \mathrm{d}\eta \right],
\end{equation}

where $u_0(s)$ is indeed the value of $u_0$ at the location of the characteristic when $\tau=t=0$.

Recapping, suppose $x = F(t;s)$ may be inverted to give $~s=G(t;x)$. Then, noting
$$ x(0,s)=F(0;s)=s,$$
by the parameterisation of the initial data, we have

$$ u(x,t) = u(x = F(t;s),t) = 1/r(t;s = G(t;x)),$$
and, in particular
$$ u_0(s) = u_0(F(0,G(t;x))) = u_0(G(t;x)).$$

Combining these expressions with Eq. (<ref>), we obtain a general expression for $r$ and therefore for $u = u(x,t)$:

\begin{equation} \label{eq:u_sol_prev}
u(x,t) = \left.\left( \frac{u_0(s)\exp\left(\int_0^t p(\eta,s) \, \mathrm{d}\eta \right)}{1 +  u_0(s)\int_0^t q(\eta,s) \exp\left(\int_0^\eta p(\xi,s) \, \mathrm{d}\xi \right)\, \mathrm{d}\eta}\right)\right|_{s = G(t;x)},
\end{equation}

where $s=G(t;x) $ is fixed on each characteristic curve and

\begin{align}
p(\eta,s) &= \beta(\eta,F(\eta;s))-H(\eta;s), \nonumber \\
q(\eta,s) &= \beta(\eta,F(\eta;s)).  \nonumber
\end{align}

Eq. (<ref>) may be also written in terms of $x$ and $t$ directly, obtaining

\begin{equation} \label{eq:u_sol}
u(x,t) = \frac{u_0(G(t;x))\exp\left(\int_0^t p(\eta,G(t;x)) \, \mathrm{d}\eta \right)}{1 +  u_0(G(t;x))\int_0^t q(\eta,G(t;x)) \exp\left(\int_0^\eta p(\xi,G(t;x)) \, \mathrm{d}\xi \right)\, \mathrm{d}\eta}.
\end{equation}

A special separable case for .

Often below, when the chemoattractant flux term is independent of time, so that $\alpha(t,x)=f(x)$, we will have $u_0\ll 1$, $\beta(t,x)\equiv 1$, $g(t)=1$. In these circumstances we have the simplification

by the symmetry $(x,s,\tau)\rightarrow(s,x,-\tau)$ in Eq. (<ref>) with $g(t)=1$
and  also that

$$p(\eta,s) = 1 - f'(F(\eta;s)), $$
which gives

$$u(x,t) =  u_0(G(t;x))\mathrm{e}^t \exp\left(-\int_0^t f'(F(\eta,s)) \, \mathrm{d}\eta \right) +\mathcal{O}(u_0^2). $$

Recalling that the integration is along a characteristic, so that $s$ is fixed,  we change the integration variable via $X= F(\eta;s)$, noting from Eq.  (<ref>) and from Eq. (<ref>), with $s$ fixed and $\alpha(t,x) = f(x)$, that

$$ \frac 1 {f(X(\eta))}= \frac{\mathrm{d}\eta}{\mathrm{d}X}.$$
Hence, on further noting $s=F(0;s),~x=F(t;s)$ and $s = G(t,x)$, we have, to within $\mathcal{O}(u_0^2)$ corrections, that

\begin{equation}
u(x,t) = u_0(G(t;x)) \mathrm{e}^t \exp\left(-\int_{s}^{x} \frac{f'(X)}{f(X)} \, \mathrm{d}X\right)  = u_0(G(t;x)) \mathrm{e}^t \frac{f(G(t;x))}{f(x)}   = u_0(G(t;x)) \mathrm{e}^t \frac{f(F(-t;x))}{f(x)}. \label{u0approx}
\end{equation}

In the examples plotted below we also have $u_0(x)=u_0$ is constant, in which case $  u_0(G(t;x))$ collapses to the constant $u_0$. For clarity please note  that $u_0$ is an abbreviation for  $ u_0(G(t;x))$ in general, though this  is  constant and denoted simply by $u_0$ in all examples, separable or otherwise, plotted below.

§.§.§ Inner solution

To explore the transition layer moving with the wavefront, we introduce a scaling of coordinates such that diffusion and advection provide a leading order dominant balance in the transition layer, with $x=x^*(t)+\delta X$ for $X$ the inner variable and $\delta \ll 1 $. With the change of variables

$$ (t,x)\rightarrow(\tau, X), ~~~ ~~~  t=\delta \tau,  ~~~~~~  x=x^*(t)+\delta X, ~~~~~~  U( \tau,X)=u(t,x),$$

one has the inner solution equations

\begin{eqnarray*} \frac{1}{\delta}U_\tau &=&\frac 1 \delta \left( x_t^*-\alpha(t,x^*(t)+\delta X ) \right)U_X + \frac{D}{\delta^2}u_{XX} +\beta U(1-U) - \alpha_x(t,x^*(t)+\delta X) \\ &=&   \frac 1 \delta \left( \alpha(t,x^*(t) ) -\alpha(t,x^*(t)+\delta X ) \right)U_X + \frac{D}{\delta^2}U_{XX} +\mathcal{O}(1),
\end{eqnarray*}

where the second line uses $x_t^*=\alpha(t,x^*(t) )$. We take  $\delta =D \ll 1 $ to  bring the advective and diffusive terms into a nominal dominant balance.

Let us now assume that $\alpha$ uniformly possesses an order one derivative with respect to $x$, that is, we assume that $\alpha(t,x^*(t)) - \alpha(t,x^*(t) + \delta X) \simeq -\alpha_x(t,x^*(t))\delta X$, therefore

$$ U_{\tau} = u_{XX} +\mathcal{O}(\delta),$$

for $X\sim \mathcal{O}(1)$. The solution of this equation, when $U(X,\tau=0) = H(X)$, with $H$ the Heaviside step function is

$$U(X,\tau) =  \frac{1}{2}\left(1+\mathrm{erf}\left(\frac{X}{2\sqrt{\tau}}\right)\right).$$

§.§.§ Composite solution

We can combine the inner solution, rewritten in terms of $(x,t)$ with the outer solutions via a composite approximation to generate an approximation for the full numerical solution across the domain (away from any prospective boundary layer at $x=1$). Thus, with $u_0(x)=u_0$, constant, and $u_{\mathrm{an}}(x,t)$ the analytical characteristic solution of Eq (<ref>) we have

\begin{equation}
u(x,t) \sim
\left\{
\begin{array}{ccc}
u_{\mathrm{an}}(x^*(t),t)U\left(\frac{x-x^*(t)}{D},\frac{t}{\tau}\right) + &0, & ~~ x< x^*(t), \\
u_{\mathrm{an}}(x^*(t),t)U\left(\frac{x-x^*(t)}{D},\frac{t}{\tau}\right) + &\left(u_\mathrm{an}(x,t) -  u_\mathrm{an}(x^*(t),t)\right), & ~~ x\geq x^*(t),
\end{array} \right.
\end{equation}

so, finally:

\begin{equation}
u(x,t) \sim
\left\{
\begin{array}{cc}
\frac{1}{2}u_{\mathrm{an}}(x^*(t),t)\left(\mathrm{erf}\left(\frac{x-x^*(t)}{2\sqrt{tD}}\right)+1\right), & ~~ x< x^*(t), \\
\frac{1}{2}u_{\mathrm{an}}(x^*(t),t)\left(\mathrm{erf}\left(\frac{x-x^*(t)}{2\sqrt{tD}}\right)-1\right) + u_\mathrm{an}(x,t), & ~~ x\geq x^*(t),
\end{array} \right.
\end{equation}

For the special case discussed in the previous section,  that is, when $\alpha(t,x) = f(x)$, $\beta(t,x) \equiv 1$ and $u_0 \ll 1$, $u_\mathrm{an}(x,t) = u_0(G(t;x))\mathrm{e}^t\frac{f(G(t;x))}{f(x)}$ and in particular $u_\mathrm{an}(x^*(t),t) = u_0(0)\mathrm{e}^t\frac{f(0)}{f(x)}$ so that for this special case

\begin{equation}
u(x,t) \sim
\left\{
\begin{array}{cc}
\frac{1}{2}u_0(0)\mathrm{e}^t\frac{f(0)}{f(x)}\left(\mathrm{erf}\left(\frac{x-x^*(t)}{2\sqrt{tD}}\right)+1\right), & ~~ x< x^*(t), \\
\frac{1}{2}u_0(0)\mathrm{e}^t\frac{f(0)}{f(x)}\left(\mathrm{erf}\left(\frac{x-x^*(t)}{2\sqrt{tD}}\right)-1\right) +  u_0(G(t;x))\mathrm{e}^t\frac{f(G(t;x))}{f(x)} , & ~~ x\geq x^*(t),
\end{array} \right.
\end{equation}

with continuity assured from the constraint $G(t;x^*(t))=0$, which holds as, by construction,  both $s=0$ and $x=x^*(t)$ on the separating characteristic and $s=G(t;x)$.

§.§ Some particular cases of interest

The solution given by Eq. (<ref>) is general given that   $D \ll 1$, though even in the separable case using Eq. (<ref>) to determine the functions $F$ and $G$ can generate complicated solutions that are not readily expressed in terms of standard functions.  Further, even when $\alpha(t,x)$ is separable or linear, allowing extensive analytical progress, there is still considerable freedom in the form of $\alpha(t,x)$ and hence we firstly analyse models where $\alpha(t,x)$ is linear, quadratic, and exponential in space, before proceeding to consider an example of cellular behaviour in a  microdevice.

[Comparison of numerical and analytical solutions for $\alpha$ linear.]Comparison of numerical and analytical solutions for the linear case. The analytical and simulated profiles are compared at different times, considering $a = 2$ and $b = 1$, $\beta(t,x) = g(t) = 1$ and $u_0(x) = u_0 = 0.05$ whereas the full numerical simulation was obtained using $D=1\times 10^{-3}$.  For this plot, and for analogous plots below, there are localised boundary layer effects near the right hand edge of the domain, $x=1$, that are not captured by the presented analytical solution.

§.§.§ Linear chemotaxis

We first consider a function of the form

\begin{equation}
\alpha(t,x) = (ax+b)g(t).
\end{equation}

Then, we have $f(x)=ax+b$ and Eq. (<ref>) yields

\begin{equation}
\frac 1 a \ln\left(\frac{ax+b}{as+b}\right) = \int_0^\tau g(\eta)\,\mathrm{d}\eta =: {\cal T}(\tau).
\end{equation}

Thus defines the function ${\cal T}$, and with ${\cal T}(t) = \int_0^t g(\eta)\,\mathrm{d}\eta $ we have

\begin{align}
F(t;s) &= \frac{1}{a}\left[(as+b)\mathrm{e}^{a\intgt }-b\right], \nonumber \\
G(t;x) &= \frac{1}{a}\left[(ax+b)\mathrm{e}^{-a\intgt }-b\right]. \nonumber
\end{align}

Hence, the transition is located at

\begin{equation} \label{eq:xstar_lin}
x^*(t) = \frac{b}{a}\left(\mathrm{e}^{a\intgt  }- 1\right).
\end{equation}

Furthermore, for $\beta(t,x) = 1$ and $g(t)=1$, so that ${\cal T}(t)=t$,  the integral expression for $u(x,t)$ in Eq. (<ref>) is readily determined to reveal

$$ u(x,t) = \frac{ u_0 (a-1) \mathrm{e}^{-(a-1)t}}{a-1+u_0 \left(1- \mathrm{e}^{-(a-1)t}  \right)} =  u_0  \mathrm{e}^{-(a-1)t} +\mathcal{O}(u_0^2),$$
while $u=0$ before the transition.
Fig. <ref> shows a comparison between the numerical results with $D=1\times 10^{-3}$ and the approximate analytical solution for $\alpha(x)=2x+1$ and $\beta(t,x) = 1 = g(t)$.  For this plot, and for analogous plots below, note the accurate prediction of the evolving front, $x^*(t)$ and the general agreement between the numerical and analytical solutions for $u(x,t)$.  Furthermore, one can expect boundary layer effects near the right hand edge of the domain, $x=1$, that are not captured by the presented solution, though these are not plotted in the current Figure.

A further, and particularly relevant case, is when $g(t)=\cos (\omega t)$, whence

\begin{align}
F(t;s) &= \frac{(as+b)\exp\left(\frac{a}{\omega}\sin(\omega t)\right)-b}{a}, \\
G(t;x) &= \frac{(ax+b)\exp\left(-\frac{a}{\omega}\sin(\omega t)\right)-b}{a}.
\end{align}

Note we still have  $G(t;x) = F(-t,x)$ since $g(t)=\cos(\omega t)$ is even so that $(x,s,\tau)\rightarrow(s,x,-\tau)$ remains  a symmetry of  Eq. (<ref>).
The transition is located at

\begin{equation} \label{eq:xstar_osc}
x^*(t) = \frac{b}{a}\left[\exp\left(\frac{a}{\omega}\sin(\omega t)\right)-1\right].
\end{equation}

Fig. <ref> shows the comparison between the numerical results with $D=1\times 10^{-3}$ and the analytical solution, with $\beta(t,x) = 1$. Note that as all the domain is displayed, the boundary layer in the numerical solution near $x=1$ can be readily observed.

Comparison between analytical solutions using asymptotic theory and numerical solutions for the oscillatory gradient.  The analytical and simulated profiles at different times are compared, considering $a=1$, $b=3$, $\omega = 10$, $\beta(t,x) = 1$ and $u_0(x) = u_0 = 0.05$ whereas the full numerical simulation was again obtained using $D=1\times 10^{-3}$. One can clearly observe that the wave of cells oscillates, with the cell density spatially constant on the right of the transition, except on approaching $x=1$, where there is a boundary layer that is not accommodated in the analysis.

The oscillating gradient case is very pertinent  as it corresponds to the case where the oxygenation feed between the two channels at the microfluidic device switches, so we shall explore it in more detail. Rather than working with the general solution given by Eq. (<ref>) it can be more expedient to consider the differential equation  for $r=1/u$ given by Eq. (<ref>); noting $\beta=1$, $f(x)=(ax+b)$, $u_0(x)=u_0$, constant, for the cases considered here, this reduces to

\begin{equation} \label{eq:edo_osc_prev}
r'+\left(1-a\cos(\omega t)\right) r = 1, \quad r(0) = 1/u_0,
\end{equation}

with no $x$-dependence. Hence, on the right of the transition region the solution is constant, $u(x) = u^*$.

In  Appendix <ref> we show that  Eq. <ref> may be solved using different asymptotic methods for four different regimes:

* Slow variations of the gradients, $\omega \ll 1$.

* Fast variations of the gradients, $\omega \gg 1$.

* Dominant chemotaxis, $\beta \ll a$.

* Dominant growth, $a \ll \beta$.

For spatial locations to the right of the transition region, but away from any boundary layer at $x=1$, these approximate solutions are compared to  numerical solutions computed using standard Runge-Kutta solvers in Fig <ref>.

$\omega \ll 1$ ($a =2$, $\beta = 3$ and $\omega = 0.1$).

$\omega \gg 1$ ($a =2$, $\beta = 3$ and $\omega = 20$).

$\beta \ll a$ ($a =1$, $\beta = 0.001$ and $\omega = 1$).

$a \ll \beta$ ($a =0.01$, $\beta = 1$ and $\omega = 1$).

Comparison between analytical solutions using asymptotic theory and numerical solutions for the oscillatory gradient. The four exposed cases are analysed: slow variations of the gradiens ($\omega \ll 1$), fast variations of the gradients ($\omega \gg 1$), dominant chemotaxis ($\beta \ll a$) and dominant growth ($a \ll \beta$).

§.§.§ Quadratic chemotaxis

Now, we consider a function of the form

\begin{equation}
\alpha(t,x) = (ax^2+bx+c)g(t),
\end{equation}

and hence $f(x)=ax^2+bx+c$.

* If $\Delta = b^2-4ac = 0$, Eq. (<ref>) gives

\begin{equation}
\frac{1}{2as+b}-\frac{1}{2ax+b} = \frac{1}{2}\int_0^\tau g(\eta)\,\mathrm{d}\eta = \frac 1 2 {\cal T}(\tau),
\end{equation}

and in turn

\begin{align}
F(t;s) &= \frac{4as+2abs\intgt +b^2\intgt }{4a-2ba\intgt -4a^2s\intgt }, \\
G(t;x) &= \frac{4ax-2ab\intgt x-b^2\intgt }{4a^2 x\intgt +2ab\intgt +4a}.
\end{align}

The transition is located at

\begin{equation} \label{eq:xstar_q1}
x^*(t) = \frac{b^2\intgt }{4a-2ba\intgt }  = \frac{c\intgt}{1-b\intgt/2}.
\end{equation}

Furthermore, in the case where $g(t)=\beta(t,x)=1$ we have
to an accuracy of $\mathcal{O}(u_0^2)$ that

$$ u(x,t) = u_0 \mathrm{e}^t \frac{f(s)}{f(x)} =  u_0 \mathrm{e}^t \left. \frac{as^2+bs+c}{ax^2+bx+c}\right|_{s=G(t;x)}= \frac{4u_0\mathrm{e}^t}{(2+bt +2atx)^2},$$

via use of  Eq. (<ref>), with $u_0\ll1$ where significant, but elementary, manipulation is required to deduce the final expression.
In particular for the parameters of Fig. <ref> it is straightforward to show that at leading order the behaviour in $x$ for $t$ fixed is linear, as observed in this figure.

* If $\Delta = b^2-4ac < 0$, Eq. (<ref>) yields

\begin{equation}
\left(\arctan\left(\frac{2ax+b}{\sqrt{ -\Delta}}\right)-\arctan\left(\frac{2as+b}{\sqrt{ -\Delta}}\right)\right) = \frac{1}{2}\sqrt{ -\Delta}\int_0^\tau g(\eta)\,\mathrm{d}\eta = \frac{1}{2}\sqrt{ -\Delta}{\cal T}(\tau) ,
\end{equation}

so that

\begin{align}
F(t;s) &= \frac{1}{2a}\left[\sqrt{-\Delta}\tan\left(\frac{1}{2}\sqrt{-\Delta}\intgt +\arctan\frac{2as+b}{\sqrt{-\Delta}}\right)-b\right], \\
G(t;x) &= \frac{1}{2a}\left[\sqrt{-\Delta}\tan\left(-\frac{1}{2}\sqrt{-\Delta}\intgt +\arctan\frac{2ax+b}{\sqrt{-\Delta}}\right)-b\right],
\end{align}

with the transition is located at

\begin{equation} \label{eq:xstar_q2}
x^*(t) = \frac{1}{2a}\left[\sqrt{-\Delta}\tan\left(\frac{1}{2}\sqrt{-\Delta}\intgt +\arctan\frac{b}{\sqrt{-\Delta}}\right)-b\right].
\end{equation}

Furthermore, the analytical and full numeric solutions are plotted for  the parameters $a=b=c=2$ with $\beta(t,x)=g(t)=1$ in Fig. <ref>. For this set of parameters,  where $t\leq 0.25$ as in Fig. <ref>, we have

\begin{eqnarray} \tan \left(\frac{t\sqrt{-\Delta}}2\right) \approx \frac{t\sqrt{-\Delta}}2 \left(1+ O\left(\frac 1 3 \left( \frac{t\sqrt{-\Delta}}2  \right)^2\right)  \right) \approx  \frac{t\sqrt{-\Delta}}2,  ~~~~ \mbox{noting} ~~~~  \frac 1 3 \left( \frac{t\sqrt{-\Delta}}2  \right)^2 \lesssim 0.06 ~. \label{app1}
\end{eqnarray}

Dropping the  corrections in higher powers of  $t\sqrt{-\Delta}$, valid for sufficiently small time including the times plotted in Fig. <ref>,
we have

$$ x^*(t) =  \frac{ct}{1-bt/2} ,$$

and then also dropping terms scaling with $\mathcal{O}(u_0^2)$,  one finds

$$ u(x,t) = u_0 \mathrm{e}^t \frac{f(s)}{f(x)} =  u_0 \mathrm{e}^t \left. \frac{as^2+bs+c}{ax^2+bx+c}\right|_{s=G(t;x)}= \frac{u_0\mathrm{e}^t (4-t^2\Delta ) }{(2+bt +2atx)^2}, $$

with a relative correction of $(1+\mathcal{O}([t\sqrt{-\Delta}/2] ^2/3)).$
The latter  again gives an approximate  linear dependence in $x$ for $t$ fixed given the parameters of  Fig. <ref>, as observed.

* If $\Delta = b^2-4ac > 0$, Eq. (<ref>) gives

\begin{equation}
\ln\left(\frac{2ax+b-\sqrt{\Delta}}{2ax+b+\sqrt{\Delta}}\right)-\ln\left(\frac{2as+b-\sqrt{\Delta}}{2as+b+\sqrt{\Delta}}\right) = \sqrt{\Delta}\int_0^\tau g(\eta)\,\mathrm{d}\eta = \sqrt{\Delta}\mathcal{T}(\tau),
\end{equation}

so that

\begin{align}  F(t;s) & = \frac{1}{2a}\left[\frac{\gamma^+(2as + \gamma^-)\exp(\frac{1}{2}\sqrt{\Delta}\intgt) - \gamma^-(2as + \gamma^+)\exp(-\frac{1}{2}\sqrt{\Delta}\intgt)}
{(2as+\gamma^+)\exp(-\frac{1}{2}\sqrt{\Delta}\intgt) - (2as+\gamma^-)\exp(\frac{1}{2}\sqrt{\Delta}\intgt)}\right],  \\
G(t;x) & = \frac{1}{2a}\left[\frac{\gamma^+(2ax + \gamma^-)\exp(-\frac{1}{2}\sqrt{\Delta}\intgt) - \gamma^-(2ax + \gamma^+)\exp(\frac{1}{2}\sqrt{\Delta}\intgt)}
{(2ax+\gamma^+)\exp(\frac{1}{2}\sqrt{\Delta}\intgt) - (2ax+\gamma^-)\exp(-\frac{1}{2}\sqrt{\Delta}\intgt)}\right],
\end{align}

where we have defined $\gamma^+ = b+\sqrt{\Delta}$ and $\gamma^- = b-\sqrt{\Delta}$, with the transition location given by

\begin{equation} \label{eq:xstar_q3}
x^*(t) = 2c\left[\frac{\exp(\frac{1}{2}\sqrt{\Delta}\intgt) - \exp(-\frac{1}{2}\sqrt{\Delta}\intgt)}{\gamma^+\exp( -\frac{1}{2}\sqrt{\Delta}\intgt) - \gamma^-\exp( \frac{1}{2}\sqrt{\Delta}\intgt})\right].
\end{equation}

For  Fig <ref>, we again have $\beta(t,x)=1=g(t)$, so that ${\cal T}(t)=t$; we also have $t\sqrt{\Delta}/2 \ll 1$ throughout the simulation regime. Hence, on neglecting higher powers of $t\sqrt{-\Delta}/2 $, the transition location simplifies to

$$x^*(t) = \frac{ct}{1- bt/2 },$$

which agrees with the above as $t\sqrt{-\Delta} / 2 \rightarrow 0.$ Furthermore, under these conditions with $u_0^2\ll 1$ and  neglecting higher orders in $u_0$, one finds

$$ u(x,t) = u_0 \mathrm{e}^t \frac{f(s)}{f(x)} =  u_0 \mathrm{e}^t \left. \frac{as^2+bs+c}{ax^2+bx+c}\right|_{s=G(t;x)}=  u_0 \mathrm{e}^t(1-(b+2ax)t) , $$

with an relative error
of $(1+\mathcal{O}([t\sqrt{-\Delta}/2]^2)).$  For $t$ fixed, again we have $u(x,t)$ is approximately linearly decreasing in $x$, as observed in Fig  <ref>.

$a=1$, $b=2$, $c=1$ ($\Delta = 0$).

$a=2$, $b=2$, $c=2$ ($\Delta < 0$).

$a=1$, $b=4$, $c=1$ ($\Delta > 0$).

[Comparison of numerical and analytical solutions for $\alpha$ quadratic.]Comparison of numerical and analytical solutions for $\alpha$ quadratic. The analytical and simulated (for $D = 1 \times 10^{-3}$) profiles at different times are compared for three different $\alpha(t,x) = ax^2+bx+c$ expressions.

More generally, Fig. <ref> shows a comparison between the numerical results with $D=1\times 10^{-3}$, $\beta(t,x)=1=g(t)$ and the analytical solution for three different expressions of $\alpha(x)$ and $\beta(t,x) = 1$, with full consistency with the above solutions and approximations.

[Comparison of numerical and analytical solutions for $\alpha$ exponential.]Comparison of numerical and analytical solutions for the exponential case. The analytical and simulated profiles are compared at different times, considering $\beta(t,x)=1=g(t)$, $\alpha(t,x) = a\exp(-\lambda x)$ with $a = 2$,   $\lambda = 1 $ and the initial condition $u_0(x)=u_0=0.05 \ll 1$. In that case, we use $D = 1 \times 10^{-3}$ for computing the numerical solutions.

§.§.§ Exponential chemotaxis

We consider now a function of the type:

\begin{equation}
\alpha(t,x) = \left(a\exp(-\lambda x)+b\exp(\lambda x)\right)g(t).
\end{equation}

Then, we have $f(x)=a\exp(-\lambda x)+b\exp(\lambda x)$, with Eq. (<ref>) reducing to

\begin{equation}
\arctan\left(\frac{b\exp(\lambda x)}{\sqrt{ab}}\right) - \arctan\left(\frac{b\exp(\lambda  s)}{\sqrt{ab}}\right) = \sqrt{ab}\lambda \int_0^\tau g(\eta)\,\mathrm{d}\eta  = \sqrt{ab}\cal{T}(\tau).
\end{equation}

\begin{align}
F(t;s) &= \frac{1}{\lambda}\ln\left(\frac{\sqrt{ab} + a\tan\left(\sqrt{ab}\lambda \intgt \right)\mathrm{e}^{-\lambda s}}{\sqrt{ab}\mathrm{e}^{-\lambda s} -b\tan\left(\sqrt{ab}\lambda \intgt \right)}\right), \\
G(t;x) &= \frac{1}{\lambda}\ln\left(\frac{\sqrt{ab} - a\tan\left(\sqrt{ab}\lambda \intgt \right)\mathrm{e}^{-\lambda x}}{\sqrt{ab}\mathrm{e}^{-\lambda x}  + b\tan\left(\sqrt{ab}\lambda \intgt \right)}\right).
\end{align}

the transition is located at

\begin{equation} \label{eq:xstar_exp}
x^*(t) = \frac{1}{\lambda}\ln\left(\frac{\sqrt{ab} + a\tan\left(\sqrt{ab}\lambda \intgt \right)}{\sqrt{ab} -b\tan\left(\sqrt{ab}\lambda \intgt \right)}\right).
\end{equation}

In particular, if $b = 0$ this reduces to

\begin{align}
F(t;s) &= \frac{1}{\lambda}\ln\left(a\lambda \intgt  + \mathrm{e}^{\lambda s}\right), \\
G(t;x) &= \frac{1}{\lambda}\ln\left(-a\lambda \intgt  + \mathrm{e}^{\lambda x}\right),
\end{align}

\begin{equation} \label{eq:xstar_exp_part}
x^*(t) = \frac{1}{\lambda}\ln\left(1+\lambda a \intgt \right).
\end{equation}

Fig. <ref> shows the comparison between the numerical results and the analytical solution with $\alpha(t,x) = 2\exp(-x)$ and $\beta(t,x) = 1=g(t)$, $D=1\times 10^{-3}$. Furthermore, in this regime  on
neglecting $\mathcal{O}(u_0^2)$ corrections,  we have

$$ u(x,t) = u_0\mathrm{e}^t \left. \frac{ f(s)}{f(x)} \right|_{s=G(t;x)} = \frac{u_0 \mathrm{e}^t}{1-a\lambda t \mathrm{e}^{-\lambda x}}, $$

noting that $x\geq x_*(t)$ on this characteristic, so that

$$ 1-a \lambda t \mathrm{e}^{-\lambda x} \geq 1- a \lambda t \mathrm{e}^{-\lambda x_*(t)} = 1- \frac{a\lambda t}{1+a\lambda t}= \frac 1 {1+a\lambda t} >0,

in turn demonstrating that $u(x,t)$ does not possess a singularity with exponential decay in $x$ for $t$ fixed away from the transition to a good approximation, as may be also observed in Fig. <ref>.

§ APPLICATIONS TO MICROFLUIDIC EXPERIMENTS

§.§ The general model

We study a broad class of problems that are related to the evolution of a cell culture in microfluidic devices under a chemotactic agent, such as oxygen, when the concentration of the agent can be computed or measured, as schematically  represented in Fig. <ref>. Hence we proceed to
consider the following model of generic cell culture evolution in microfluidic devices:

\begin{align}
\frac{\partial C_n}{\partial T} &=  \frac{\partial}{\partial X}\left(D\frac{\partial C_n}{\partial X} - \chi C_n\frac{\partial B}{\partial X}\right) + \alpha_n M_g(B)C_n\left(1-\frac{C_n}{c_\mathrm{sat}}\right) - \alpha_{nd}M_d(B)C_n, \label{eq:governing_a} \\
\frac{\partial C_d}{\partial T} &=  \alpha_{nd}M_d(B)C_n, \label{eq:governing_b} \\
\frac{\partial B}{\partial T} &=  \frac{\partial}{\partial X}\left(D_B\frac{\partial B}{\partial X}\right) - \alpha_B  W(B,C_n). \label{eq:governing_c}
\end{align}

where $C_n$ and $C_d$ are respectively the alive and dead cell concentrations and $B$ is a chemotactic agent, $M_g$ and $M_d$ are nonlinear dimensionless corrections accounting for the effect of the chemo-attractant on cell growth and death. In addition $W(B,C_n)$ is the nonlinear dimensionless correction of the chemo-attractant consumption by all cells in the microdevice, noting that additional cells other than those of direct interest may be present, such as the study of metastatic tumour cells of interest migrating within a stromal cell population for instance. We also assume that dead cell concentration is  sufficiently low that it does not  compromise either live cell proliferation or migration. Note that if the chemotaxis agent is a cell nutrient, such as glucose or oxygen, $W(B,C_n) \neq 0$, whereas for unconsumed biochemical signals, $W(B,C_n)=0$.  The boundary conditions considered are

\begin{align}
\left.\frac{\partial f_n}{\partial x}\right|_{x=0} &= 0, \\
\left.\frac{\partial f_n}{\partial x}\right|_{x=L} &= 0, \\
B(x=0,t) &= B^L(t), \\
B(x=L,t) &= B^R(t),
\end{align}

with $B^L(t)$ and $B^R(t)$, the chemotactic agent concentration at the left and right channel, and $f_n = D\frac{\partial C_n}{\partial X} - \chi C_n\frac{\partial B}{\partial X}$ the alive cell flow.  As the dead cell population is slave to the other variables in this model, we neglect it henceforth  which is equivalent to taking $\alpha_{nd} \simeq 0$.

Therefore, the full model here analysed is

\begin{align}
\frac{\partial C_n}{\partial T} &=  \frac{\partial}{\partial X}\left(D\frac{\partial C_n}{\partial X} - \chi C_n\frac{\partial  B}{\partial X}\right)  +  \alpha_n M_g(B) C_n \left(1-\frac{C_n}{c_{\mathrm{sat}}}\right), \label{eq:igoverning_a} \\
\frac{\partial B}{\partial T} &=  \frac{\partial}{\partial X}\left(D_B\frac{\partial B}{\partial X}\right) - \alpha_B  W(B,C_n). \label{eq:igoverning_b}
\end{align}

In order to evaluate the relevance of the different phenomena, we define the dimensionless variables:

\begin{align}
C_n &=  c_\mathrm{sat}u, \label{eq:dimensionless_variables_u} \\
B &=  B^*v, \label{eq:dimensionless_variables_v} \\
X &=  Lx, \label{eq:dimensionless_variables_x} \\
T &=  \frac{t}{\alpha_n}, \label{eq:dimensionless_variables_t}
\end{align}

Hence Eqs. (<ref>) become

\begin{align}
\frac{\partial u}{\partial t} &=  \frac{\partial}{\partial x}\left(\Pi_1\frac{\partial u}{\partial x} - \Pi_2 u\frac{\partial v}{\partial x}\right) + m(v)u\left(1-u\right), \label{eq:governing_a_d} \\
\frac{\partial v}{\partial t} &=  \frac{\partial}{\partial x}\left(\Pi_3\frac{\partial v}{\partial x}\right) + \Pi_4 w(v,u), \label{eq:governing_b_d}
\end{align}

\begin{align}
\Pi_1 &=  \frac{D}{\alpha_n L^2}, \label{eq:Pi_1} \\
\Pi_2 &=  \frac{\chi\mathrm{O_2^*}}{\alpha_n L^2}, \label{eq:Pi_2} \\
\Pi_3 &=  \frac{D_B}{\alpha_n L^2}, \label{eq:Pi_3} \\
\Pi_4 &=  \frac{\alpha_B c_{\mathrm{sat}}}{ \alpha_n B^*}, \label{eq:Pi_4} \\
m(v)  &=  M_g(B^*v), \label{eq:r} \\
w(v,u)  &=  W(B^*v,c_{\mathrm{sat}}u). \label{eq:s} \\
\end{align}

The associated boundary conditions are

\begin{align}
\left.\frac{\partial f_n}{\partial x}\right|_{x=0} &= 0, \label{eq:boundary_d1} \\
\left.\frac{\partial f_n}{\partial x}\right|_{x=1} &= 0, \label{eq:boundary_d2} \\
v(x=0,t) &= \psi_1(t), \label{eq:boundary_d3}\\
v(x=1,t) &= \psi_2(t), \label{eq:boundary_d4}
\end{align}

where now $f_n = \Pi_1\frac{\partial u}{\partial x} - \Pi_2 u\frac{\partial v}{\partial x}$ and $\psi_1(t),~\psi_2(t)$ are prescribed functions of time corresponding to the level of nutrient or chemoattractant fixed to be at the channel edges.  In what follows, we will use $u_t$ and $u_x$ as an abbreviation for $\frac{\partial u}{\partial t}$ and $\frac{\partial u}{\partial x}$.

In particular the governing PDEs described in Eqs. (<ref>) may be reformulated as

\begin{align}
u_t &=  \Pi_1 u_{xx} -\Pi_2 \left(v_x u\right)_x + m(v)u\left(1-u\right), \label{eq:governing_dimensionless_1a} \\
v_t &=  \Pi_3 v_{xx} - \Pi_4 w(v,u). \label{eq:governing_dimensionless_1b}
\end{align}

§.§.§ The weak consumption limit

First we consider the case where the chemoattractant is not a nutrient and therefore is not consumed by cells. In that case, we can set
$$w(v,u) = 0.$$
For all the problems in which it is possible to assume $\Pi_1 \ll 1$,  whereby  random cellular motility is negligible compared to directed chemotaxis and with $\Pi_3 \gg 1$,  so that the  chemoattractant diffusion is large relative to cellular diffusion Eqs. (<ref>) reduce to

\begin{align}
u_t  + k \left(v_x u\right)_x &=  m(v)u\left(1-u\right), \label{eq:governing_dimensionless_2a} \\
v_{xx} &= 0. \label{eq:governing_dimensionless_2b}
\end{align}

where $k = \Pi_2$.

However, care is required in considering the boundary conditions for $u$ and the initial conditions for $v$ in this reduced model, due to the loss of the second spatial derivative of $u$ and the first temporal derivative of $v$. In particular, we cannot satisfy all the boundary conditions for $u$; instead we have boundary layers. We have seen in the examples these occur  at internal transitions and at the right of the domain (see Fig <ref>). Thus for the simplified system the boundary condition is enforced at $x=0$ for $u$ with the simplification of the flux to $f_n=-kv_x u$, as diffusion is treated as negligible. However,  enforcing the boundary condition at $x=1$ will require the consideration of a boundary layer that is not resolved in the simplified model as it is complicated, but does not further insight into cell migratory and chemotaxis. The two boundary conditions for $v$, as given by Eqs. <ref> and <ref>  are inherited and applied at both boundaries. We now consider the initial conditions for $v$, which cannot be satisfied. Instead there is an analogous temporal boundary layer for early time while initial transients relax, though such transients persist for such a short time that they are not of interest, and thus not resolved, here. The justification of the neglect of these fast transients is  further detailed in Appendix <ref>, where it is demonstrated that the solution of Eqn (<ref>) corresponds to the leading order composite solution in a temporal boundary layer analysis that exploits $\Pi_3\gg1$ .

Proceeding, we set $v(x=0,t) = \psi_1(t)$ and $v(x=1,t) = \psi_2(t)$, Eq. (<ref>) is immediately integrated to
\begin{eqnarray}\label{vout1}
v(t,x)=\left(\psi_2(t) - \psi_1(t)\right)x + \psi_1(t).
\end{eqnarray}

We recover therefore Eq. (<ref>) with
$$\alpha(t,x) = k(\psi_2(t) - \psi_1(t)),~~~~~~~~~\beta(t,x)  = m(\left(\psi_2(t) - \psi_1(t)\right)x + \psi_1(t)),$$
that is a special case of the linear problem $\alpha(t,x) = a(t)x+b(t)$, with $a(t) = 0$ and $b(t) = k(\psi_2(t) - \psi_1(t))$, so the different functions needed in order to compute the solution are:

\begin{align}
F(t;s) &=  s + k\int_0^t\Delta \psi(\eta)\, \mathrm{d}\eta, \label{eq:ex1_1} \\
G(t;x) &= x -  k\int_0^t\Delta \psi(\eta)\, \mathrm{d}\eta, \label{eq:ex1_2} \\
x^*(t) &=  k\int_0^t\Delta \psi(\eta)\, \mathrm{d}\eta = x-G(t;x),  \label{eq:ex1_3}
\end{align}

where we have defined $\Delta \psi (t) = \psi_2(t) - \psi_1(t)$. Also, the expression of the cell profile far from the transition is

\begin{equation}
u(x,t) =  \frac{u_0( G(t;x))\exp\left(\int_0^t K(\eta, F(\eta;G(t;x))) \, \mathrm{d}\eta\right)}{1+u_0( G(t;x))\left(\exp\left(\int_0^t K(\eta, F(\eta;G(t;x)))\, \mathrm{d}\eta\right) - 1\right)} \label{eq:ex1_4}
\end{equation}

\begin{equation}
K(\eta,X) = m\left((\psi_2(\eta) - \psi_1(\eta))X + \psi_1(\eta)\right).
\end{equation}

The evolution of the transition coordinate $x^*(t)$ and the dimensionless cell profile for different times are shown in Fig. <ref> for $k=1$, $\psi_1(t) = 0$ and different external stimuli $\psi_2(t)$. In particular, let us consider the case with $m(v) = m_0 v$, and with $\psi_1(t)=0$ and use of the change of variable $X=F(\eta;s)$, whereby on a characteristic

$$\dfrac{\mathrm{d}X}{\mathrm{d}\eta } = \alpha(\eta,X)=k\psi_2(\eta)$$

with $F(t;s)=x$ and $F(0,s)=s$. This reveals

\begin{eqnarray*} \left. \int_0^t K(\eta,F(\eta,s))\, \mathrm{d}\eta \right|_{s=G(t;x)} &=&  \left. \int_s^x \frac 1 {k\psi_2(\eta)}  m_0 \psi_2(\eta) X \, \mathrm{d} X \right|_{s=G(t;x)}= \frac{m_0}{2k} (x-G(t;x))(x+G(t;x)) \\ &=&  \frac{m_0}{2k} x^*(t)(2x-x^*(t)).
\end{eqnarray*}

Combined with Eq. (<ref>), this gives a simple expression for $u(x,t)$.  For instance

$$ u(x,t) \approx u_0 \exp\left( \frac{m_0}{k} x^*(t)(x-x^*(t)/2) \right) +\mathcal{O}(u_0^2),$$

for $u_0\ll 1$ constant;  in this case, we have an increasing function at fixed $t$ to the right of the transition given $x^*(t)>0$ and this is essentially  linear for $x^*(t)\ll 1$, as observed  in Fig. <ref>.

Value of the transition coordinate $x^* = x^*(t)$.

Cell profile at $t=2/3$.

[Solution for the case with no consumption]Approximate analytical solution for the case with no consumption. We assume $k=1$, $m(v) = v$ and the transition region is recreated with $D = 1 \times 10^{-4}$. The initial condition is set again to $u_0(x) = u_0 = 0.05$. Also, we consider three different shapes for $\psi_2$, with $\psi_1=0$. A constant oxygen level $\psi_2(t) = 1$, a trigonometric oscillatory stimulus $\psi_2(t) = 1 + \cos(\omega t)$ with $\omega = 10\pi $, and an increasing ramp stimulus $\psi_2(t) = t$.

§.§.§ Cellular consumption of chemoattractant

A more interesting case is when the chemoattractant is a nutrient and therefore it is consumed by cells. In that case, $w(v,u) \neq 0$.

With $w_{\mathrm{pc}}(v)$ denoting the non-dimensional uptake of nutrient per cell, we take $w_{\mathrm{pc}}(v)$ to be monotonic increasing with

\begin{align}
\lim_{v\rightarrow 0} w_{\mathrm{pc}}(v) = 0, \\
\lim_{v\rightarrow +\infty} w_{\mathrm{pc}}(v) = 1,
\end{align}

where the final limit is without loss of generality, with the overall scale of uptake governed by $\Pi_4$.

For instance, with Michaelis-Menten kinetics we take  [49]:

\begin{equation}
w_{\mathrm{pc}}(v) = \frac{v}{v+k_m},
\end{equation}

or, more in general, Hill-Langmuir equation for modelling the consumption kinetics [50]

\begin{equation}
w_{\mathrm{pc}}(v) = \frac{v^n}{v^n+k^n_H},
\end{equation}

In any of the aforementioned cases, there are numerous potential scenarios:

* There are no other cells at the microfluidic device besides the cell culture of our interest and we are at the low cell regime, $\Pi_4 u/\Pi_3 \ll 1$. In that case, after rapid initial transients   describing the diffusion relaxation of the nutrient, and too fast to be on the timescale of cellular motility,  Eq. (<ref>) becomes $v_{xx} = 0$, and the discussion is analogous to the case without the consumption term.

* There are other non-migrating cells within the microfluidic device in addition to the migrating cells, for instance if we are considering  a metastasis model, with the other cells at  constant concentration and in excess of the tumour cells. If we additionally have high nutrient concentrations the situation is illustrated in Fig. <ref> using oxygen as an example of nutrient. In that case, $w_{\mathrm{pc}}(v) \sim 1$ and

$$ w(u,v) = K,$$

where  $K$ is the (dimensionless) total amount of cells, essentially constant as the non-tumour cells are in excess. Then, with the definition
$$ \lambda = \Pi_4 K/\Pi_3,$$
Eq. (<ref>)

\begin{eqnarray} \label{vint1}  v_t = \Pi_3 v_{xx} - \Pi_3 \lambda . \end{eqnarray}

In addition, we assume $\Pi_3^{-1}\ll \lambda \ll \Pi_3,$ so that  $\lambda$ may be treated as order one on using asymptotic methods based on  on the leading order of approximations based on $\Pi_3\gg 1$. Then the above further  reduces to

\begin{eqnarray} \label{vfin1}  v_{xx} =   \lambda ,  \end{eqnarray}

noting, as above,  that fast initial transients are not of interest, with   further justification of Eqn (<ref>) in Appendix <ref> via a boundary layer analysis.

* There are other cells within the microfluidic device at constant concentration and in excess relative to the tumour cells, together with low nutrient concentrations. The situation is also illustrated in Fig. <ref>. Assuming the Michaelis-Menten model, $w_\mathrm{pc}(v) \sim v/k_m$ so that
$w(u,v) = Kv/k_m$ where $K$ is the effectively  constant non-dimensional total cell density and  Eq. (<ref>)
reduces to

\begin{eqnarray} \label{vint2}  v_t = \Pi_3 v_{xx} - \Pi_3 \lambda v,  \end{eqnarray}

where now $\lambda =   \Pi_4 K  /(k_m\Pi_3) $. As above,  this reduces if  $\Pi_3^{-1}\ll \lambda \ll \Pi_3,$ which we assume in order to yield
\begin{eqnarray} \label{vfin2}  v_{xx} =   \lambda v,  \end{eqnarray}

once more noting  fast initial transients are not of interest,  with additional justification of Eqn (<ref>) presented in Appendix <ref>.

High oxygen levels.

Low oxygen levels.

[Cell culture model recreating cancer cells in an oxygenated ambient.]Cell culture model recreating cancer cells in an oxygenated ambient. For illustration purposes, the oxygen is considered as the chemoattractant  nutrient of the heterogeneous growth and death. The different local and global oxygen gradients appearing in the chamber may be approximated by approximating Eq. (<ref>) when considering specific situations, such as high oxygen levels or low oxygen levels. Created with BioRender.com.

The presence of other cells in excess and high chemoattractant or oxygen levels.

With the assumptions and motivations as above in the derivation of Eq. (<ref>) we now have

\begin{align}
u_t  + k \left(v_x u\right)_x &=  m(v)u\left(1-u\right), \label{eq:governing_dimensionless_3a} \\
v_{xx} &= \lambda. \label{eq:governing_dimensionless_3b}
\end{align}

where the boundary conditions are again the ones given by Eqs. (<ref>), except for the fact that the flux is given by $f_n=-kv_x u$ as diffusion is neglected and the cell boundary condition at $x=1$, with its associated boundary layer, is no longer considered as discussed in detail previously.

If we prescribe $v(x=0,t) = \psi_1(t)$ and $v(x=1,t) = \psi_2(t)$, Eq. (<ref>)  integrates to

\begin{eqnarray} \label{vout2} v(x,t)=\frac{1}{2}\lambda x^2 + \left(\psi_2(t) - \psi_1(t) - \frac{1}{2}\lambda\right)x + \psi_1(t).
\end{eqnarray}

We have therefore
$$\alpha(t,x) = k\lambda x + k\left(\psi_2(t) - \psi_1(t)-\frac{1}{2}\lambda\right),

$$ \beta(t,x) = m\left(v(t,x)\right),$$
for cell growth.  Once more, this is a special case of $\alpha(t,x) = a(t)x+b(t)$ with $a(t) = k\lambda$ and $b(t) = k\left(\psi_2(t) - \psi_1(t)-\frac{1}{2}\lambda\right)$, so the different functions needed in order to compute the solution are:

\begin{align}
F(t;s) &= s\mathrm{e}^{k\lambda t} + \frac{1}{2} \left(1-\mathrm{e}^{k\lambda t}\right) + k\int_0^t\Delta \psi(\eta)\mathrm{e}^{-k\lambda \eta} \, \mathrm{d}\eta, \label{eq:ex2_1} \\
G(t;x) &= x\mathrm{e}^{-k\lambda t} - \frac{1}{2} \left(\mathrm{e}^{-k\lambda t}-1\right) - k\mathrm{e}^{-k\lambda t}\int_0^t\Delta \psi(\eta)\mathrm{e}^{-k\lambda \eta} \, \mathrm{d}\eta, \label{eq:ex2_2} \\
x^*(t) &= \frac{1}{2} \left(1-\mathrm{e}^{k\lambda t}\right) + k\int_0^t\Delta \psi(\eta)\mathrm{e}^{-k\lambda \eta}\, \mathrm{d}\eta. \label{eq:ex2_3}
\end{align}

The expression of the cell profile may be computed using Eq. (<ref>), which gives

\begin{equation}
u(x,t) = \frac{u_0(G(t;x))\exp(-k\lambda t)\exp\left(\int_0^t K(\eta, F(\eta;G(t;x))) \, \mathrm{d}\eta \right)}{1 + u_0(G(t;x))\int_0^t K(\eta, F(\eta;G(t;x)))\exp(-k\lambda \eta) \exp\left(\int_0^\eta K(\xi, F(\xi;G(t;x))) \, \mathrm{d}\xi \right)\, \mathrm{d}\eta} ,
\end{equation}

where now

\begin{equation}
K(\eta,X) = m\left(\frac{1}{2}X^2 + (\psi_2(\eta) - \psi_1(\eta) - \frac{1}{2}\lambda)X + \psi_1(\eta)\right).
\end{equation}

The evolution of the transition coordinate $x^*(t)$ and the dimensionless cell profile for different times are shown in Fig. <ref> for $k=1$, $\lambda=0.1$, $m(v) = m_0v$ with $m_0 = 1$, $u_0 = 0.05$ and a external stimulus given by $\psi_1(t) = 0$ and different shapes for $\psi_2(t)$.

While a fast oscillatory stimulus with $\psi_2(t)=1+\cos(10 \pi t)$  plotted in Fig. <ref> does not allow a ready approximation for cell density to the right of the transition region, we can consider the cases of $\psi_2(t)=1$ or $\psi_2(t)=t$ that are also considered in this Figure. In particular, we can use  elementary but extensive manipulation to determine and approximate

$$ \left. \int_0^t K(\eta, F(\eta,s)) \, \mathrm{d}\eta \right|_{s=G(t;x)},$$
to deduce that

$$ u_0 \approx u_0 \mathrm{e}^{-\lambda t} \exp\left( \frac {1-\mathrm{e}^{-2\lambda t}}{4\lambda} (x-x^*(t))^2 + \frac 1 3 h_1(t,\lambda)(x-x^*(t)) +  h_2(t,\lambda) \right) + \mathcal{O}(u_0^2),
in these cases.
In particular, $h_1(t,\lambda)$ is of the form
$$h_1(t,\lambda) = \left( 3t + \frac 3 2 t^2\right)-\left(\frac 3 2 + \frac 9 4 t +   t^2 \right) \lambda t +
\mathcal{O}(\lambda^2t^2),$$

for $\psi_2(t)=1$ and
$$h_1(t,\lambda) = \left(\frac{3}{2}t^2 + \frac{1}{2}t^3\right) - \left(\frac{3}{2}+ \frac{3}{4}t + \frac{1}{2}t^2+\frac{3}{8}t^3\right)\lambda t + \mathcal{O}(\lambda^2t^2),$$

for $\psi_2(t)=t$.

For  $h_2(t,\lambda)$, we have  with $\psi_2(t)=1$ that
$$ h_2(t, \lambda) =   \frac 1 2 t^2 \left( 1+\frac 1 3 t \right) - t \left(  \frac 1 2 + \frac 1 3 t + \frac 1 8 t^2  \right)\lambda t
+ \mathcal{O}(\lambda^2t^2),$$

while, in contrast, for $\psi_2(t)=t$ we note that
$$h_2(t,\lambda) =  \frac{t^4}{40} \left(5 + t\right)  - t^2\left(\frac{1}{4}+\frac{1}{16}t+\frac{1}{15}t^2+\frac{1}{36}t^3\right)\lambda t + \mathcal{O}(\lambda^2t^2).$$

For both cases note that at fixed time the cell concentration to the right of the transition is essentially the exponential of a quadratic in $x-x^*(t)$, provided that $\lambda$ is small enough.

Value of the transition coordinate $x^* = x^*(t)$.

Cell profile at $t=2/3$.

[Solution for the case with chemoattractant consumption and high chemoattractant levels]Solution for the case with chemoattractant consumption and high chemoattractant levels. Approximate analytical solution for $k=1$,   $m(v)=v$,  $\lambda = 0.1$ and the transition region is recreated with $D = 1 \times 10^{-4}$. The initial condition is set again to $u_0(x) = u_0 = 0.05$. We consider three different shapes for $\psi_2$. A constant oxygen level $\psi_2(t) = 1$, an oscillatory stimulus $ \psi_2(t) = 1 + \cos(\omega t)$ with $\omega = 10\pi$, and an increasing stimulus $\psi_2(t) = t$.

Other cells and low chemoattractant or oxygen levels. With analogous reasoning, Eqs. (<ref>) now reduce to

\begin{align}
u_t  + k \left(v_x u\right)_x &=  m(v)u\left(1-u\right), \label{eq:governing_dimensionless_4a} \\
v_{xx} &= \lambda v. \label{eq:governing_dimensionless_4b}
\end{align}

The general solution to Eq. (<ref>) generates  a chemotactic function $\alpha(t,x)$ that is neither linear in $x$ nor separable, with Eq. <ref> equivalent to a Riccati differential equation in  $\exp(x)$, for which general solutions are not known in terms of standard functions. Even though the general case is not tractable in terms of constructing solutions, two particular important configurations do allow progress, namely the gradient configuration ($v(x=0,t)=0$ and $v(x=1,t)= \psi(t)$) and the symmetric configuration (($v(x=0,t)=v(x=1,t)=\psi(t)$).

The gradient configuration is certainly the most interesting configuration. For that case, Eq. (<ref>) integrates  to

\begin{eqnarray} \label{vout3} v(t,x)=\psi(t)\sinh(\sqrt{\lambda} x),
\end{eqnarray}

$$\alpha(t,x) = k\sqrt{\lambda}\cosh(\sqrt{\lambda} x)\psi(t), ~~~~~~~~~ \beta(t,x) = m_0(v(t,x)).$$
Hence we have separability, with  $f(x) = k\sqrt{\lambda}\cosh(\sqrt{\lambda} x)$ and $g(t) = \psi(t)$. With

$$ {\cal S}(t) = \tan(k\lambda \intgt /2), ~~~~~~~ \intgt = \int_0^t g(\eta)\mathrm {d}\eta$$

we in turn have

\begin{align}
F(t;s) &= \frac 1{\sqrt{\lambda}}\ln\left(\frac{1+\mathrm{e}^{-s\sqrt{\lambda}}{\cal S}(t)}{\mathrm{e}^{-s\sqrt{\lambda}}-{\cal S}(t)}\right), \label{eq:ex3a_1} \\
G(t;x) &= \frac 1{\sqrt{\lambda}}\ln\left(\frac{1-\mathrm{e}^{- x\sqrt{\lambda}}{\cal S}(t)}{\mathrm{e}^{- x\sqrt{\lambda}}+{\cal S}(t)}\right),  \label{eq:ex3a_2} \\
x^*(t) &= \frac 1{\sqrt{\lambda}}\ln\left(\frac{1+{\cal S}(t)}{1-{\cal S}(t)}\right). \label{eq:ex3a_3}
\end{align}

As with the previous cases the expression of the cell profile may be computed using Eq. (<ref>), whereby

\begin{equation}
u(x,t) =  \frac{u_0(G(t;x))\exp\left(\int_0^t K(\eta, F(\eta;G(t;x))) \, \mathrm{d}\eta \right)}{1 +  u_0(G(t;x))\int_0^t K(\eta,F(\eta;G(t;x))) \exp\left(\int_0^\eta K(\xi, F(\xi;G(t;x))) \, \mathrm{d}\xi  \right)\, \mathrm{d}\eta }, \label{equfin}
\end{equation}

where $K(\eta,X) = m(v) - kv_{xx}$, that is

\begin{equation}
K(\eta,X) = m\left(\psi(\eta)\sinh(\sqrt{\lambda} X)\right) - k \lambda\psi(\eta)\sinh(\sqrt{\lambda} X),
\end{equation}

and, with $m(v) = m_0 v$ and $u_0(x)=u_0 \ll 1,$ constant, use of the change of variables $X=F(\eta,s)$ reduces  Eq. (<ref>) to

$$ u(x,t) = u_0 \left(\frac{\cosh(\sqrt{\lambda}G(t;x)) }{\cosh(\sqrt{\lambda} x)}   \right)^{\frac{m_0-k\lambda}{k\lambda} } +\mathcal{O}(u_0^2).

The evolution of the transition coordinate $x^*(t)$ and the dimensionless cell profile for different times are shown in Fig. <ref> for $m_0=1$, $k=1$, $\lambda = 0.1$, $ u_0(x)=u_0=0.05,$ constant. Here,  the cellular density to the right of transition further simplifies to

$$ u(x,t) =  u_0 \left( \frac{1+{\cal S}(t)^2}{(1-{\cal S}(t)^2) +2 {\cal S}(t)\sinh\left(\frac{x}{\sqrt{10}}\right) }\right)^{9} +\mathcal{O}(u_0^2), ~~~~~ {\cal S}(t) = \tan(\intgt / 20),$$
which entails for small time the spatial variation of the cellular density is that of a hyperbolic sine. Finally,  note from the expression for $x^*(t)$ and the fact $x^*(t) \leq 1$ for the transition zone to be within the domain, we have the bound ${\cal T} \leq \tanh(1/2) < 1 $ so that the denominator in the above expression is positive and there is no singularity.

Even if less interesting from the experimental point of view, we can also obtain an expression for a symmetric configuration by considering half of the domain and applying Neumann boundary conditions for $x=0$ so Eq. (<ref>) is integrated to

\begin{eqnarray} \label{vout4}v(t,x)=\psi(t)\cosh(\sqrt{\lambda} x).  \end{eqnarray}

Again we are under a separable case with $\alpha(t,x) = k\sqrt{\lambda}\sinh(\sqrt{\lambda} x)\psi(t)$ and $\beta(t,x) = 1$, so $f(x) = k\sqrt{\lambda}\sinh(\sqrt{\lambda} x)$ and $g(t) = \psi(t)$, so the different functions needed in order to compute the solution are:

\begin{align}
F(t;s) &= \frac 2 {\sqrt{\lambda}} \arctan\left(\mathrm{e}^{k\lambda\intgt}\tanh\left(\frac{s\sqrt{\lambda}}{2}\right)\right), \label{eq:ex3b_1} \\
G(t;x) &= \frac 2 {\sqrt{\lambda}}
\arctan\left(\mathrm{e}^{-k\lambda\intgt}\tanh\left(\frac{x\sqrt{\lambda}}{2}\right)\right), \label{eq:ex3b_2} \\
x^*(t) &= 0. \label{eq:ex3b_3}
\end{align}

Now in Eq. (<ref>) we shall use

\begin{equation}
K(\eta,X)= m\left(\psi(\eta)\cosh(\sqrt{\lambda} X\right) - k \lambda\psi(\eta)\cosh(\sqrt{\lambda} X).
\end{equation}

and, with $m(v) = m_0 v$ and  $u_0(x)=u_0 \ll 1,$ constant, Eq. ((<ref>)) reduces to

$$ u(x,t) = u_0 \left(\frac{\sinh(\sqrt{\lambda}x) }{\sinh(\sqrt{\lambda} G(t;x))}   \right)^{\frac{m_0-k\lambda}{k\lambda} } +\mathcal{O}(u_0^2).

Value of the transition coordinate $x^* = x^*(t)$.

Cell profile at $t=2/3$.

[Solution for the case with chemoattractant consumption and low chemoattractant levels]Solution for the case with chemoattractant consumption and low chemoattractant levels. Approximate analytical solution for $k=1$, $m(v) = v$, $\lambda = 0.1$ and the transition region is recreated with $D = 1 \times 10^{-4}$. The initial condition is set again to $u_0(x) = u_0 = 0.05$. We consider three different shapes for $\psi_1$. A constant oxygen level $\psi_1(t) = 1$, an oscillatory stimulus $\psi_1(t) = 1 + \cos(\omega t)$ with $\omega = 10\pi $, and an increasing stimulus $\psi_1(t) = t$.

§ DISCUSSION

There is an increasing use of in vitro investigations in the exploration of cellular motility, for instance with microdevice studies exploring tumour cell dynamics in response to oxygen gradients, as illustrated by GBM studies [51, 37, 40]. In turn this has motivated the main theme of this paper, namely modelling cell migration chemotaxis in heterogeneous environments, which applies for general chemoattractants, not just oxygen. In particular, the governing equations for cellular motility that have been considered are of the form

\begin{equation}
\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} - \frac{\partial }{\partial x}\left(\alpha(t,x)u_x\right) + \beta(t,x)u(1-u),
\end{equation}

and supplemented by zero flux boundary conditions and suitable, typically  constant, initial conditions.
Furthermore, the general spatiotemporal heterogeneity in the chemotactic function $\alpha(t,x)$ and the growth function $\beta(t,x)$ emerges from the chemoattractant gradients that the cells respond to, which may be manipulated extensively in microdevice experiments. However, the heterogeneity also entails that while cell migration dynamics towards high concentrations of chemoattractant is anticipated, the dynamics will not simply be that of a translationally invariant travelling wave and their associated analytical simplicity.

Hence, we have considered a  framework capable of considering spatiotemporal chemotactic gradients and the resulting wavefront and cell density dynamics for  cellular migrations in the presence of spatial and temporal heterogeneity.
In particular, the fact cell spreading in the absence of chemotaxis is generally much slower than in its presence, the non-dimensional cellular diffusion coefficient is very small, so that typically  $D\ll 1$, which we assume and exploit in this study.  With this, we have that away from boundary layers near sources, which were located at the domain edge in the examples considered, and away from cell wavefronts with sharp transitions,  an advection equation for the cellular density $u$ emerges with $D\sim 0$.
Nonetheless, this advection equation is still complex, entailing that a constraint on constructing solutions for the cell density behaviour within  the analytical framework presented here is that the chemotactic function $\alpha(t,x)$ must either be linear in $x$ or separable with respect to space and time. However, numerous cases are consistent with these constraints, as illustrated by the range of examples considered  in Section <ref>, together with  the examples that emerge from the consideration of cellular dynamics within microdevices in Section <ref>.

With these constraints, we have illustrated how the method of characteristics can be used  to construct the cell density solutions away from boundary layers and transition regions, together with the use of boundary layer methods to construct a uniformly continuous approximate solution that accommodates the transition in the wavefront of the cells. In turn this provides an analytical characterisation of the movement of cellular wavefronts and the concentrations of cells  within and either side of the wavefront. However we do not resolve boundary layers near oxygen sources, which has been the righthand  boundary in the examples considered, for the presented framework given the limited insight the boundary layer will provide for the overall cell behaviour.

Even with the restriction of linearity or separability of the chemotactic functions $\alpha(t,x)$, there is extensive freedom in the choices of $\alpha(t,x)$ and $\beta(t,x)$. Hence we first considered the predictions of the model for cell behaviour in exemplar test cases. In particular, we have explored and documented the behaviour of the cell density where the chemotactic function varies  linearly in space and non-trivially in time, as well as quadratically in space and finally  with exponential spatial behaviour.
In addition to analytical investigations, we have also verified that the analytical solutions faithfully reproduce the behaviour of direct numerical solutions of the model, as may be observed extensively in Figures <ref>-<ref>.
Furthermore the exemplar solutions are subsequently used to study models reduced from a full coupling between the cellular density and chemoattractant concentration with a microdevice setting, as documented in Section <ref>.
Once more the resulting dynamics are analytically investigated and documented, with rationally based  analytical approximation for the location of the transition region and the cellular density presented.

Such  analytical solutions are useful in numerous ways. For example, they provide oversight and insight for the system dynamics across parameter space. To illustrate, we consider the wavefront being driven by chemoattractant at the right hand boundary, so that $\psi_1(t)=0, v(t,1)=\psi_2(t)$. Then, with $k$ a measure of the strength of chemotaxis,   differentiating $x^*(t)$ from Eq. (<ref>) for the speed for the transition front gives

$$ \frac{\mathrm{d}x^*}{\mathrm{d}t}= -\frac 1 2 k\lambda \mathrm{e}^{k\lambda t} + k \psi_2(t)\mathrm{e}^{-k\lambda t}$$

for the model of section <ref>, with high chemoattractant or oxygen levels. Hence, in this case, increasing the chemoattractant/oxygen uptake, as representing by increasing $\lambda$, always slows the propagation of a rightmoving wave. In contrast, for the model with low chemoattractant or oxygen levels of  section  <ref> differentiating $x^*(t)$ from Eq. (<ref>) we obtain
$$ \frac{\mathrm{d}x^*}{\mathrm{d}t}= \frac{\sqrt{\lambda}}{\mathrm{cos}\left(k\lambda \mathcal{T}\right)}\psi_2(t)$$
therefore revealing that increasing the chemoattractant/oxygen uptake always speeds up the propagation of rightmoving wave, illustrating the general deductions that may be made from the presented analytical solutions.

As well as providing analytical characterisation of the systems behaviour the analytical approximations provide a means to very rapidly compute cell behaviour. Thus the approximate analytical solutions can support computationally intensive studies. A simple example would be a global sensitivity analyses over all parameters. In particular, rapid evaluation would be most useful for parameter estimation using experimental, often noisy, data especially if Bayesian techniques are considered as this require extensive simulation to provide posterior distributions, rather than optimisation techniques which only generate point estimates for parameters. A directly analogous example is Bayesian model selection, whereby the comparison of experimental data and model prediction is used to distinguish different model structures when these are not known a priori, such as different functional forms of $\alpha(t,x)$ or $\beta(t,x)$ representing different prospective growth and chemotactic responses for a given the tumour cell line in question. In particular both optimisation techniques and the Bayesian techniques are iterative, so that the use of a rational but rapid evaluation of  an approximation to an optimum in optimisation studies or a posterior for Bayesian techniques can in turn be used to restart the procedure with the full numerical model to further refine the results [52].

In summary, we have developed a framework for the construction of analytical approximations for front behaviours and densities for cells undergoing chemotaxis in heterogeneous environments, as characterised by the chemotactic and growth functions $\alpha$ and $\beta$ respectively.  The resulting cellular waves of migration are not simple travelling waves due to the heterogeneity induced by the chemoattractant profiles. Nonetheless, numerous  features of the wavefronts have been extracted via rational approximation, such  as the location and speed of the propagating wave, together with the cellular density profile. These have been explored and validated for exemplars scenarios as well as investigated for models fundamentally motivated by experimental microdevices for observing cellular motility under a very wide range of conditions, even if  complete generality is not feasible for progress using constructive methods. Thus, the  solutions presented here not only provide insight into the behaviour of cellular motility under the influence of spatiotemporal chemotactic heterogeneity, but also highlight general behaviours and important mechanisms and  parameters, as well as providing a means of rapid evaluation in demanding computational studies, such as Bayesian parameter estimation and model selection.

§ AKNOWLEDGEMENTS

The authors gratefully acknowledge the financial support  from the Spanish Ministry of Science and Innovation (MICINN), the State Research Agency (AEI), and FEDER, UE through the projects PID2019-106099RB-C44/AEI and PID2021-126051OB-C41, the Government of Aragon (DGA) and the “Centro de Investigación Biomedical en Red en Bioingeniería, Biomateriales y Nanomedicina (CIBER-BBN)”, financed by the Instituto de Salud Carlos III with assistance from the European Regional Development Fund (FEDER).

[1]
Christopher S Chen, Milan Mrksich, Sui Huang, George M Whitesides, and Donald E
Geometric control of cell life and death.
Science, 276(5317):1425–1428, 1997.

[2]
ASG Curtis and GM Seehar.
The control of cell division by tension or diffusion.
Nature, 274(5666):52, 1978.

[3]
Ulrich S Schwarz and Ilka B Bischofs.
Physical determinants of cell organization in soft media.
Medical Engineering & Physics, 27(9):763–772, 2005.

[4]
Stephen B Carter.
Haptotaxis and the mechanism of cell motility.
Nature, 213(5073):256, 1967.

[5]
Chun-Min Lo, Hong-Bei Wang, Micah Dembo, and Yu-li Wang.
Cell movement is guided by the rigidity of the substrate.
Biophysical Journal, 79(1):144–152, 2000.

[6]
Dennis Bray.
Cell movements: from molecules to motility.
Garland Science, 2000.

[7]
Seyed Jamaleddin Mousavi, Mohamed Hamdy Doweidar, and Manuel Doblaré.
3d computational modelling of cell migration: a
mechano-chemo-thermo-electrotaxis approach.
Journal of Theoretical Biology, 329:64–73, 2013.

[8]
Vinay Kumar, AK Abbas, and JC Aster.
Cell injury, cell death, and adaptations.
Robbins Basic Pathology, 8:1–30, 2013.

[9]
Seyed Jamaleddin Mousavi, Manuel Doblaré, and Mohamed Hamdy Doweidar.
Computational modelling of multi-cell migration in a multi-signalling
Physical Biology, 11(2):026002, 2014.

[10]
Sui Huang and Donald E Ingber.
Cell tension, matrix mechanics, and cancer development.
Cancer Cell, 8(3):175–176, 2005.

[11]
Douglas Hanahan and Robert A Weinberg.
Hallmarks of cancer: the next generation.
Cell, 144(5):646–674, 2011.

[12]
Anika Nagelkerke, Johan Bussink, Alan E Rowan, and Paul N Span.
The mechanical microenvironment in cancer: How physics affects
In Seminars in cancer biology, volume 35, pages 62–70.
Elsevier, 2015.

[13]
Daniela F Quail and Johanna A Joyce.
Microenvironmental regulation of tumor progression and metastasis.
Nature Medicine, 19(11):1423–1437, 2013.

[14]
SJ Mousavi, MH Doweidar, and M Doblaré.
Cell migration and cell-cell interaction in the presence of
Mol Cell Biomech, 10(1):1–25, 2013.

[15]
Jack W Scannell, Alex Blanckley, Helen Boldon, and Brian Warrington.
Diagnosing the decline in pharmaceutical r&d efficiency.
Nature reviews Drug Discovery, 11(3):191–200, 2012.

[16]
Garry P Nolan.
What's wrong with drug screening today.
Nature Chemical Biology, 3(4):187–191, 2007.

[17]
Rasheena Edmondson, Jessica Jenkins Broglie, Audrey F Adcock, and Liju Yang.
Three-dimensional cell culture systems and their applications in drug
discovery and cell-based biosensors.
Assay and Drug Development Technologies, 12(4):207–218, 2014.

[18]
Chengyang Wang, Zhenyu Tang, Yu Zhao, Rui Yao, Lingsong Li, and Wei Sun.
Three-dimensional in vitro cancer models: a short review.
Biofabrication, 6(2):022001, 2014.

[19]
Eric K Sackmann, Anna L Fulton, and David J Beebe.
The present and future role of microfluidics in biomedical research.
Nature, 507(7491):181–189, 2014.

[20]
Sangeeta N Bhatia and Donald E Ingber.
Microfluidic organs-on-chips.
Nature Biotechnology, 32(8):760, 2014.

[21]
Simone Bersini, Jessie S Jeon, Gabriele Dubini, Chiara Arrigoni, Seok Chung,
Joseph L Charest, Matteo Moretti, and Roger D Kamm.
A microfluidic 3d in vitro model for specificity of breast cancer
metastasis to bone.
Biomaterials, 35(8):2454–2461, 2014.

[22]
Alexandra Boussommier-Calleja, Ran Li, Michelle B Chen, Siew Cheng Wong, and
Roger D Kamm.
Microfluidics: a new tool for modeling cancer–immune interactions.
Trends in Cancer, 2(1):6–19, 2016.

[23]
Jessie S Jeon, Simone Bersini, Mara Gilardi, Gabriele Dubini, Joseph L Charest,
Matteo Moretti, and Roger D Kamm.
Human 3d vascularized organotypic microfluidic assays to study breast
cancer cell extravasation.
Proceedings of the National Academy of Sciences,
112(1):214–219, 2015.

[24]
Ioannis K Zervantonakis, Shannon K Hughes-Alford, Joseph L Charest, John S
Condeelis, Frank B Gertler, and Roger D Kamm.
Three-dimensional microfluidic model for tumor cell intravasation and
endothelial barrier function.
Proceedings of the National Academy of Sciences,
109(34):13515–13520, 2012.

[25]
Danli Wu and Patricia Yotnda.
Induction and testing of hypoxia in cell culture.
JoVE (Journal of Visualized Experiments), (54):e2899, 2011.

[26]
Helen M Byrne.
Dissecting cancer through mathematics: from the cell to the animal
Nature Reviews Cancer, 10(3):221–230, 2010.

[27]
Philipp M Altrock, Lin L Liu, and Franziska Michor.
The mathematics of cancer: integrating quantitative models.
Nature Reviews Cancer, 15(12):730–745, 2015.

[28]
Hiroaki Kitano.
Computational systems biology.
Nature, 420(6912):206, 2002.

[29]
Elaine L Bearer, John S Lowengrub, Hermann B Frieboes, Yao-Li Chuang, Fang Jin,
Steven M Wise, Mauro Ferrari, David B Agus, and Vittorio Cristini.
Multiparameter computational modeling of tumor invasion.
Cancer Research, 69(10):4493–4501, 2009.

[30]
R Eils and AE Kriete.
Computational systems biology: from molecular mechanisms to
Academic, New York, NY, 2013.

[31]
HM Byrne, T Alarcon, MR Owen, SD Webb, and PK Maini.
Modelling aspects of cancer dynamics: a review.
Philosophical Transactions of the Royal Society A: Mathematical,
Physical and Engineering Sciences, 364(1843):1563–1578, 2006.

[32]
Moriah E Katt, Amanda L Placone, Andrew D Wong, Zinnia S Xu, and Peter C
In vitro tumor models: advantages, disadvantages, variables, and
selecting the right platform.
Frontiers in Bioengineering and Biotechnology, 4:12, 2016.

[33]
Daniel J Brat.
Glioblastoma: biology, genetics, and behavior.
American Society of Clinical Oncology Educational Book,
32(1):102–107, 2012.

[34]
Haralampos Hatzikirou, Andreas Deutsch, Carlo Schaller, Matthias Simon, and
Kristin Swanson.
Mathematical modelling of glioblastoma tumour development: a review.
Mathematical Models and Methods in Applied Sciences,
15(11):1779–1794, 2005.

[35]
Yangjin Kim, Hyejin Jeon, and Hans Othmer.
The role of the tumor microenvironment in glioblastoma: A
mathematical model.
IEEE Transactions on Biomedical Engineering, 64(3):519–527,

[36]
Haralampos Hatzikirou, David Basanta, Matthias Simon, K Schaller, and Andreas
‘go or grow’: the key to the emergence of invasion in tumour
Mathematical Medicine and Biology: a journal of the IMA,
29(1):49–65, 2012.

[37]
Jose M Ayuso, Rosa Monge, Alicia Martínez-González, María
Virumbrales-Muñoz, Guillermo A Llamazares, Javier Berganzo, Aurelio
Hernández-Laín, Jorge Santolaria, Manuel Doblaré, Christopher
Hubert, et al.
Glioblastoma on a microfluidic chip: generating pseudopalisades and
enhancing aggressiveness through blood vessel obstruction events.
Neuro-Oncology, 19(4):503–513, 2017.

[38]
Alicia Martínez-González, Gabriel F Calvo, Luis A Pérez Romasanta,
and Víctor M Pérez-García.
Hypoxic cell waves around necrotic cores in glioblastoma: a
biomathematical model and its therapeutic implications.
Bulletin of Mathematical Biology, 74(12):2875–2896, 2012.

[39]
Hermann B Frieboes, Xiaoming Zheng, Chung-Ho Sun, Bruce Tromberg, Robert
Gatenby, and Vittorio Cristini.
An integrated computational/experimental model of tumor invasion.
Cancer Research, 66(3):1597–1604, 2006.

[40]
Jacobo Ayensa-Jiménez, Marina Pérez-Aliacar, Teodora Randelovic, Sara
Oliván, Luis Fernández, José Antonio Sanz-Herrera, Ignacio Ochoa,
Mohamed H Doweidar, and Manuel Doblaré.
Mathematical formulation and parametric analysis of in vitro cell
models in microfluidic devices: application to different stages of
glioblastoma evolution.
Scientific Reports, 10(1):1–21, 2020.

[41]
Marina Pérez-Aliacar, Mohamed H Doweidar, Manuel Doblaré, and Jacobo
Predicting cell behaviour parameters from glioblastoma on a chip
images. a deep learning approach.
Computers in Biology and Medicine, page 104547, 2021.

[42]
Víctor M Pérez-García, Gabriel F Calvo, Juan Belmonte-Beitia,
David Diego, and Luis Pérez-Romasanta.
Bright solitary waves in malignant gliomas.
Physical Review E, 84(2):021921, 2011.

[43]
Philip Gerlee and Sven Nelander.
Travelling wave analysis of a mathematical model of glioblastoma
Mathematical biosciences, 276:75–81, 2016.

[44]
Tracy L Stepien, Erica M Rutter, and Yang Kuang.
Traveling waves of a go-or-grow model of glioma growth.
SIAM Journal on Applied Mathematics, 78(3):1778–1801, 2018.

[45]
Evelyn F Keller and Lee A Segel.
Traveling bands of chemotactic bacteria: a theoretical analysis.
Journal of theoretical biology, 30(2):235–248, 1971.

[46]
Gurusamy Arumugam and Jagmohan Tyagi.
Keller-segel chemotaxis models: a review.
Acta Applicandae Mathematicae, 171(1):1–82, 2021.

[47]
Chuan Xue, Hyung Ju Hwang, Kevin J Painter, and Radek Erban.
Travelling waves in hyperbolic chemotaxis equations.
Bulletin of Mathematical Biology, 73(8):1695–1733, 2011.

[48]
Yoojin Shin, Sewoon Han, Jessie S Jeon, Kyoko Yamamoto, Ioannis K
Zervantonakis, Ryo Sudo, Roger D Kamm, and Seok Chung.
Microfluidic assay for simultaneous culture of multiple cell types on
surfaces or within hydrogels.
Nature Protocols, 7(7):1247–1259, 2012.

[49]
Athel Cornish-Bowden.
The origins of enzyme kinetics.
FEBS Letters, 587(17):2725–2730, 2013.

[50]
Archibald Vivian Hill.
The possible effects of the aggregation of the molecules of
haemoglobin on its dissociation curves.
J. Physiol., 40:4–7, 1910.

[51]
Jose M Ayuso, María Virumbrales-Muñoz, Alodia Lacueva, Pilar M Lanuza,
Elisa Checa-Chavarria, Pablo Botella, Eduardo Fernández, Manuel Doblare,
Simon J Allison, Roger M Phillips, et al.
Development and characterization of a microfluidic model of the
tumour microenvironment.
Scientific Reports, 6(1):1–16, 2016.

[52]
Liam V Brown, Jonathan Wagg, Rachel Darley, Andy van Hateren, Tim Elliott,
Eamonn A Gaffney, and Mark C Coles.
De-risking clinical trial failure through mechanistic simulation.
Immunotherapy Advances, 2(1):ltac017, 2022.

[53]
Frank Olver.
Asymptotics and special functions.
CRC Press, 1997.

§ OSCILLATORY SOLUTIONS FOR HOMOGENEOUS GROWTH

We consider  homogeneous growth, so that $\beta(t,x)$ is a constant independent of $x$ and $t$, together with  a temporally oscillating chemotactic response governed by $\alpha(t,x)$, so that Eq. (<ref>) reduces to

\begin{equation} \label{eq:a_ini}
r'+\left(\beta-a\cos(\omega t)\right)r =\beta, \quad r(0) = r^*.
\end{equation}

Hence the  cell concentration to the right of the transition region, but away from the boundary layer at $x=1$, is spatially constant with  a temporal oscillation, which we consider below   for four distinct parameter regimes.

§.§ Slow variations of the gradients, .

Let $T = \omega t$ and, without loss of generality, we consider the decomposition  $r(t) = y+z(\omega t)$. Then, Eq. (<ref>) becomes

\begin{equation} \label{eq:edo_truco_1}
\omega z'(\omega t) + \frac{\mathrm{d}y}{\mathrm{d}t} + \left(\beta - a\cos(\omega t)\right)z(\omega t) + \left(\beta - a\cos(\omega t)\right)y = \beta.
\end{equation}

Further, let $\omega z'(\omega t) + \left(\beta - a\cos(\omega t)\right)z(\omega t)=\beta$, so we have

\begin{equation} \label{eq:edo_f}
\omega z'(T) +\left(\beta - a\cos(T)\right)z(T) = \beta.
\end{equation}

Since $\omega \ll 1$, we can approximate the solution of Eq. (<ref>) by

\begin{equation} \label{eq:sol_f}
z(T) = \frac{\beta}{\beta -a\cos(T)} + \mathcal{O}(\omega).
\end{equation}

Then, Eq. (<ref>) yields

\begin{equation} \label{eq:edo_truco_2}
\frac{\mathrm{d}y}{\mathrm{d}t}+ \left(\beta - a\cos(\omega t)\right)y = 0,
\end{equation}

and the initial condition is

\begin{equation}
y^* = y(t=0) = r(t=0) - z(0) = r^* - \frac{\beta}{\beta - a}.
\end{equation}

As we have a slow modulation of the frequency/decay rate, we use  the  Wentzel-Kramers-Brillouin (WKB) method [53]. In terms of $T = \omega t$, Eq. (<ref>) becomes

\begin{equation} \label{eq:edo_truco_3}
\omega \frac{\mathrm{d}y}{\mathrm{d}T} + \left(\beta - a\cos(T)\right)y = 0, \quad y(0) = y^*.
\end{equation}

The WKB approximation is expressed here as

\begin{equation}\label{eq:WKBJ}
y = p\exp\left(\frac{\phi(T)}{\omega}\right)J(T), \quad J(T) = J_0 + \omega J_1 + \mathcal{O}(\omega^2).
\end{equation}

Substituting Eq. (<ref>) into Eq. (<ref>)  we obtain

\begin{equation}
p\exp\left(\frac{\phi(T)}{\omega}\right)\left[\omega\left(\frac{\dot{\phi}}{\omega}J + \dot{J}\right)+(\beta - a\cos T)J\right] = 0.
\end{equation}

Therefore, $\omega\left(\frac{\dot{\phi}}{\omega}J + \dot{J}\right)+(\beta - a\cos T)J=0$, so that

\begin{equation}
\dot{\phi}\left(J_0 + \omega J_1 + \ldots\right) + \omega\left(\dot{J}_0 + \omega \dot{J}_1 + \ldots\right) + \left(\beta - a\cos T\right)\left(J_0 + \omega J_1 + \ldots\right) = 0.
\end{equation}

The $\mathcal{O}(1)$ corresponding equation is

\begin{equation}
J_0\left(\dot{\phi} + (\beta - a \cos T)\right) = 0,
\end{equation}

and solving it for $\phi$ gives

\begin{equation}
\phi(T) = \phi^*  -\beta T + a \sin T.
\end{equation}

The $\mathcal{O}(\omega)$ corresponding equation is

\begin{equation}
J_1\left(\dot{\phi} + (\beta - a \cos T)\right) + \dot{J}_0 = 0,
\end{equation}

so, as $\dot{\phi} + (\beta - a \cos T)=0$, we obtain $J_0 = J_0^*$ (constant).

Consequently, Eq. (<ref>) becomes

\begin{equation}
y = K\exp\left(\frac{-\beta T+a\sin T}{\omega}\right),
\end{equation}

for $K = p\exp(\frac{\phi^*}{\omega})J_0$ constant. Using that $T = \omega t$ and the initial value $y(0) = y^*$, we obtain the approximation

\begin{equation}
y = \left(r^* - \frac{\beta}{\beta - a} \right)\exp\left(-\beta t +\frac{a}{\omega}\sin(\omega t)\right).
\end{equation}

Finally, as $r = y + z(\omega t)$, we have

\begin{equation}
r \sim \frac{\beta}{\beta -a\cos(\omega t)} + \left(r^* - \frac{\beta}{\beta - a} \right)\exp\left(-\beta t +\frac{a}{\omega}\sin(\omega t)\right).
\end{equation}

§.§ Fast variations of the gradients, .

We solve now the problem using the method of multiple scales. Let us assume that $r(t) = r(T_1,T_2)$, where $T_1 = t$ and $T_2 = \omega t$, so that

\begin{equation}
\frac{dr}{dt} = \frac{\partial r}{\partial T_1}\frac{\partial T_1}{\partial t} + \frac{\partial r}{\partial T_2}\frac{\partial T_2}{\partial t} = \frac{\partial r}{\partial T_1} + \omega\frac{\partial r}{\partial T_2}.
\end{equation}

If $\varepsilon = \frac{1}{\omega} \ll 1$, Eq. (<ref>) becomes

\begin{equation}
\varepsilon\left(\frac{\partial r}{\partial T_1} + \frac{\partial r}{\partial T_2}\right) + \varepsilon\left(\beta - a\cos(T_2)\right)r= \varepsilon \beta.
\end{equation}

If we use an asymptotic expansion of $r$, $r = r_0 +\varepsilon r_1$ we obtain:

\begin{equation}
\varepsilon\left(\frac{\partial r_0}{\partial T_1} + \varepsilon\frac{\partial r_1}{\partial T_1}\right) + \frac{\partial r_0}{\partial T_2} + \varepsilon \frac{\partial r_1}{\partial T_2} +  \varepsilon\left(\beta - a\cos(T_2)\right)\left(r_0 + \varepsilon r_1\right) + \mathcal{O}(\varepsilon^2) = \varepsilon \beta.
\end{equation}

Solving the equation obtained collecting the  $\mathcal{O}(1)$ terms, we find $r_0$ is a function of $T_1$ only, that is

\begin{equation}
r_0 = r_0(T_1),
\end{equation}

with the initial condition $r_0(0) = r^*.$

Now, for the equation obtained collecting the $\mathcal{O}(\varepsilon)$ terms, we have

\begin{equation} \label{eq:condition_1}
r_0'(T_1) + \frac{\partial r_1}{\partial T_2} + (\beta - a\cos(T_2))r_0(T_1) = \beta.
\end{equation}

As $r_1$ is a periodic correction, integrating Eq. (<ref>) over  $T_2\in [0,2\pi]$ we obtain

\begin{equation} \label{eq:integral}
2\pi r_0'(T_1) + 2\pi\beta r_0(T_1) = 2\pi \beta,
\end{equation}

and thus, solving for $r_0(T_1)$, we have

\begin{equation}
r_0(T_1) = r^*e^{-\beta T_1} + (1-e^{-\beta T_1}).
\end{equation}

Hence, the leading order approximation is

\begin{equation}
r(t) \sim r^*e^{-\beta t} + (1-e^{-\beta t}).
\end{equation}

Now, for the $\mathcal{O}(\varepsilon)$ equation we have:

\begin{equation} \label{eq:oeps}
\frac{\partial r_1}{\partial T_2} = a\cos(T_2)r_0(T_1).
\end{equation}

Since $r_1(T_1,T_2 = 0) = 0$,

\begin{equation} \label{eq:v1_sol}
r_1(T_1,T_2) = a\sin(T_2)r_0(T_1).
\end{equation}

and therefore the first order correction is, $r_0 + \varepsilon r_1$, that is:

\begin{equation}
r(t) \sim \left(r^*e^{-\beta t} + (1-e^{-\beta t})\right)\left(1 + \frac{a}{\omega}\sin(\omega t)\right).
\end{equation}

§.§ Dominant chemotaxis, .

We solve now the problem using the standard asymptotic expansion method. We set $\varepsilon = \beta/a \ll 1$. Then, Eq. (<ref>) reduces to

\begin{equation}
r' + a\left(\varepsilon-\cos(\omega t)\right)r = a\varepsilon.
\end{equation}

The leading order solution is obtained immediately as it is the solution to the homogeneous linear differential equation:

\begin{equation}
r_0'-a\cos(\omega t)r_0=0.
\end{equation}

Thus  the leading order approximation is

\begin{equation}
r\sim  r_0=r^*e^{\frac{1}{\omega}a\sin(\omega t)},
\end{equation}

though the next order correction generates a cumbersome expression and thus is not presented.

§.§ Dominant cell proliferation,  .

Again, we solve the problem using the standard asymptotic expansion method. Now, we set $\varepsilon = a/\beta \ll 1$. Then, Eq. (<ref>) gives

\begin{equation}
r' + \beta \left(1-\varepsilon\cos(\omega t)\right)r = \beta.
\end{equation}

The leading order solution is obtained immediately as  the solution to the inhomogeneous linear differential equation

\begin{equation}
r_0'+\beta r_0=\beta,
\end{equation}

and hence is given by

\begin{equation}
r \sim r_0=r^*e^{-\beta t} + \left(1-e^{-\beta t}\right).
\end{equation}

For the first correction, an asymptotic expansion of $r$ of the form $r=r_0 + \varepsilon r_1$ reveals that

\begin{equation}
r_1'+\beta r_1=\beta\cos(\omega t)r_0.
\end{equation}

The solution to this ODE with initial condition $r_1(0) = 0$ is given by

\begin{equation}
r_1(t) = \frac{\beta\left[(r^*-1)\gamma^2\sin(\omega t)+\omega^2e^{\beta t}\sin(\omega t) + \omega \beta e^{\beta t}\cos(\omega t)\right]}{\omega \gamma^2e^{\beta t}}-\frac{\beta^2}{\gamma^2}e^{-\beta t},
\end{equation}

where we have defined $\gamma^2 = \beta^2+\omega^2$. Hence up to the first order correction we have

\begin{equation}
r(t) \sim r_0 + \varepsilon r_1 =  1+\frac{(r^*-1)\gamma^2 \omega-A\left[\left(\gamma^2r^*-\gamma+\omega^2e^{\beta t}\right)\sin(\omega t)+ \omega \beta e^{\beta t}\cos(\omega t)-\beta\omega\right]}{\omega \gamma^2 e^{\beta t}}.
\end{equation}

§ SIMPLIFICATION OF THE TRANSPORT EQUATION FOR CHEMOTAXIS

With $v(t,x)$ denoting the concentration of the chemoattractant, and with $\Pi_3\gg 1$ we have
in section <ref> an equation of the form

$$ v_t = \Pi_3 v_{xx}, $$

since $w(u,v) =0$ has been imposed on (<ref>), while in section <ref> we have equations of one of the two forms

$$ v_t = \Pi_3 v_{xx} - \Pi_3\lambda , ~~~~~~~ v_t = \Pi_3 v_{xx}
-\Pi_3 \lambda v, $$
via Equations (<ref>), (<ref>) respectively.
These are accompanied by boundary conditions of the form

$$v(t,0)= \psi_1(t), ~~~~~~~v(t,1) = \psi_2(t), $$

and initial conditions $v(0,x)=v_0(x)$ are required to close the system. For simplicity, we assume the boundary conditions and initial conditions are consistent at $(t,x)=(0,0), (0,1)$.

To further proceed we firstly assume  $\Pi_3^{-1}\ll \lambda \ll \Pi_3$, so that $\lambda$ can be treated as unit order of magnitude in asymptotic methods based on $\varepsilon = \Pi_3^{-1}\ll 1$, and we also assume that $\psi_1(t),$ $\psi_2(t)$ have derivatives that are unit order of magnitude, or less. With these weak assumptions, our objective is show that the time derivative $v_t$ can be neglected at leading order, justifying the use of Eqs. (<ref>), (<ref>), (<ref>) in the main text and also justifying the neglect of the consideration of initial conditions in the main text on the grounds this only governs fast initial transients.

Below we work with the PDE

\begin{eqnarray} \label{vgen} v_t = \Pi_3 v_{xx} - \mu \Pi_3 \lambda v^\zeta, \end{eqnarray}

with  $\mu,~\zeta\in\{0,1\}$ so that $\mu=0$ gives one the above equations while $\mu=1,~\zeta=0$, gives another with the final possibility corresponding to $\mu=1,~\zeta=1$, allowing  the three cases to be considered  simultaneously.

We have an outer timescale of $t$ and an inner timescale of $\tau = \Pi_3 t=t/\varepsilon$. In the outer region, $t\gg\epsilon$, for the leading order outer solution $v^\mathrm{out}(t,x)$ one indeed has

$$ 0 =  v^\mathrm{out}_{xx} - \mu  \lambda (v^\mathrm{out})^\zeta, ~~~~~ v^\mathrm{out}(t,0)= \psi_1(t), ~~~~~~~v^\mathrm{out}(t,1) = \psi_2(t),$$

without consideration of the initial condition. Hence the solutions presented in the main text, for instance
Eqs. (<ref>), (<ref>),(<ref>),(<ref>), are the same as the  leading order outer solutions. With the leading order  inner solution $v^\mathrm{in}(\tau,x)$  we have

$$  v^\mathrm{in}_\tau =  v^\mathrm{in}_{xx} - \mu  \lambda ( v^\mathrm{in})^\zeta, ~~~~ v^\mathrm{in}(0,x)=v_0(x),$$

and the  leading order boundary conditions

$$ v^\mathrm{in}(\tau,0)= \psi_1(\varepsilon \tau) = \psi_1(0) , ~~~~~~
v^\mathrm{in}(\tau,1)= \psi_2(\varepsilon \tau) = \psi_2(0) ,$$

on noting $t=\varepsilon \tau$ and where $\mathcal{O}(\varepsilon)$ corrections are dropped in the final term for both boundary conditions.
Recalling   $\mu,~\zeta\in\{0,1\}$ and using the decomposition

$$ v^\mathrm{in}(\tau,x) = v^\mathrm{out}(0,x)  + q(\tau,x), $$

without loss of generality, we have

$$ q_\tau = q_{xx} -\mu\zeta \lambda q, ~~~ q(\tau,0)=q(\tau,1)=0, ~~~ q(0,x) = v_0(x)-v^\mathrm{out}(0,x).$$

With the Fourier decomposition of $q(\tau,x)$ and its initial condition

$$ q(x,\tau) = \sum_{n=1}^\infty q_n(\tau) \sin\left(n\pi x \right), ~~~ q(0,x) = \sum_{n=1}^\infty q_n^0 \sin\left(n\pi x \right),

we have
$$ \frac{\mathrm{d}q_n}{\mathrm{d}\tau} = -n^2\pi^2 q_n - \mu \zeta\lambda q_n, ~~~ q_n(0) = q_n^0, $$
so that

$$  q(x,\tau) = \sum_{n=1}^\infty q_n^0 \sin\left(n\pi x \right) \mathrm{e}^{-(n^2\pi^2+\mu\zeta\lambda)\tau} \approx 0 ~~ \mbox{for} ~~ \pi^2 \tau \gg 1.$$

Thus the leading order composite solution is given by

\begin{eqnarray} v(t,x) &=& v^\mathrm{out}(t,x) + v^\mathrm{in}(\tau(t),x) - \lim_{t \rightarrow 0}v^\mathrm{out}(t,x) =v^\mathrm{out}(t,x) + q(\tau(t),x)
\nonumber \\  &\approx & v^\mathrm{out}(t,x) ~~ \mbox{for} ~~ t  \gg \epsilon/\pi^2.
\end{eqnarray}

Thus, as implemented in the main text, working solely  with the outer solution and neglecting  the initial conditions is a rational asymptotic approximation at leading order with respect to $\varepsilon = \Pi_3^{-1} \ll 1$, once initial transients have decayed, that is for times satisfying $t\gg \epsilon/\pi^2 = \Pi_3^{-1}/\pi^2$. As the initial transients, and very short times, are not of interest for determining the behaviour of the invasive front of cells for the majority of its propagation we thus work only with the outer equations and solution in the main text.