# Classification and stability analysis of polarising and depolarising
travelling wave solutions for a model of collective cell migration

Nizhum Rahman1,2,∗, Robert Marangell3 and Dietmar Oelz1 1 School of
Mathematics and Physics, The University of Queensland, QLD, 4072, Australia.

###### Abstract

We study travelling wave solutions of a 1D continuum model for collective cell
migration in which cells are characterised by position and polarity. Four
different types of travelling wave solutions are identified which represent
polarisation and depolarisation waves resulting from either colliding or
departing cell sheets as observed in model wound experiments. We study the
linear stability of the travelling wave solutions numerically and using
spectral theory. This involves the computation of the Evans function most of
which we are able to carry out explicitly, with one final step left to
numerical simulation.

Keywords: collective cell migration, travelling wave solution, stability
analysis, Evans function.

##  1 Introduction

Epithelial cells line the surfaces of our body. They are tightly packed and
found in various organ systems. They serve many roles such as secretion,
absorption, sensation, protection and transport. In the context of various
physiological and pathological processes including tissue repair, cancer and
wound healing [1, 2] they undergo collective cell migration which is an
important feature characterising the development and life-cycle of multi-
cellular organisms. The analysis of mathematical models for cell migration is
a powerful tool that allows to identify characteristic features of the
collective dynamics.

Cells can move as a group with shared responsibilities according to recent
studies on migratory epithelial tissues [3, 4, 5, 6, 7]. During collective
cell migration epithelial cells maintain stable cell-cell junctions [8, 9, 10,
11].

Directed cell migration is closely linked to the polarity of the cells [12],
which refers to the spatial differences between cells in shape, structure and
function.

How cell polarity is influenced by intercellular coupling is a fundamental
question. The best way of understanding its underlying biochemical signalling
mechanism is the planar cell polarity pathway coupling bistable intracellular
states among adjacent cells [13, 14, 15, 16, 17].

In cells, polarity propagates as a travelling wave from cell to cell [18, 19,
17]. The corresponding travelling wave solution of a 1D model for the
collective migration of epithelial cells has been computed explicitly in [17].
Preliminary numerical simulations in [17] have shown that the same model
exhibits other travelling wave solutions such as a depolarisation wave
entering the sheet of departing cells.

In the present article we identify various other types of travelling wave
solutions of this model. We analyse their stability with the help of spectral
theory and using tools from dynamical systems theory [20, 21, 22] such as the
Evans function which provides information about the point spectrum of the
linearised model. Note that we are able to compute most of the information to
setup the Evans function explicitly with only one last component being
evaluated numerically.

This article is structured as follows. Section 2 states the governing
equations of the 1d model for collective cell migration. Section 3 lists
various types of polarising and depolarising travelling wave solutions which
can be formulated building on the solution found in [17]. In section 4 we
provide numerical evidence for their stability and in section 5 we perform the
linear stability analysis. Finally we wrap up and discuss these results in
section 6.

##  2 Preliminary results

The mathematical model for collective migration of epithelial cells introduced
in [17] characterises each cell by its position along the real axis and by a
real-valued quantity called polarity $a\in\mathbb{R}$. It represents the
asymmetry of the cell and it is assumed to be linked to active migratory
velocity through a non-linear function $v_{i}=M(a_{i})$. This function may be
thought of as the following step function (or a smoothened version of it as in
Fig. 1A)

$M(a)=\begin{cases}0&a\leq\alpha\;,\\\ 1&a>\alpha\;,\end{cases}$ (1)

which models polarity-dependent motility with threshold polarity $\alpha>0$.

For a single cell the following differential equation for the polarity is
defined which models auto-depolarisation of the cell and adaption to actual
motion. It is given (omitting physical constants) by $\dot{a}=-a+v=-a+M(a)$
which - taken as a dynamical system - features two stable steady states, a
non-polarised one at $a=0$ and a polarised one at $a=1$ (see Fig. 1).

Figure 1: A: Visualisation of the steady states (arrows) of the single cell
model defined by $M(a)=a$. B: Sketch of 1d model for collective cell
migration.

The 1d model introduced in [17] treats the epithelial cell sheet as a 1d chain
of such cells connected by linearly elastic springs. As a consequence each
cell’s velocity is determined by the spring forces emanating from neighbouring
cells and by the active migratory force. The governing equations (after non-
dimensionalisation) are

$\left\\{\begin{aligned} &\dot{a}_{i}=-a_{i}+\dot{x}_{i}\;,\\\
&\dot{x}_{i}=M(a_{i})+\kappa(x_{i_{1}}-2x_{i}+x_{i-1})\;,\end{aligned}\right.$
(2)

where $\kappa$ is a phenomenological parameter representing cellular
contractility.

The continuum model associated to this particle model has been derived in
[17]. Its governing equations are the continuity equation for the cell density
$\rho=\rho(x,t)$ coupled to the momentum equation and a separate equation for
the cell polarity $a(x,t)$ which correspond to the system (2). After non-
dimensionalisation, these equations are given by

$\begin{cases}\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x}(\rho
v)=0,\\\ \frac{\partial a}{\partial t}+v\frac{\partial a}{\partial x}=-a+v,\\\
v=M(a)-\kappa\frac{1}{\rho^{3}}\frac{\partial\rho}{\partial x},\end{cases}$
(3)

which involves the velocity field as $v=v(x,t)$ and the cell polarity
$a=a(x,t)$.

In [17] the system (3) is reformulated using the travelling wave ansatz
$\rho(x,t)=R(z),a(x,t)=A(z)$, where $z=x-st$ parametrises the wave profiles
and $s$ is the wave speed. Then the reduced system of equations is coupled to
boundary conditions which correspond to one of the stable fixed points of the
single particle model (2) prescribed at $z=-\infty$, namely
$(R_{-\infty},A_{-\infty})=(1,0)$, where $R_{-\infty}$ is the cell density at
rest normalised to 1. This yields

$\begin{cases}R^{\prime}=\frac{R^{2}}{\kappa}((M(A)-s)\,R+s),\\\
A^{\prime}=1+(\frac{A-s}{s})\,R.\end{cases}$ (4)

Note that the second stationary point of (4) is
$(R_{\infty},A_{\infty})=(\frac{s}{s-1},1)$.

In [17] the authors identify an explicit solution to the travelling wave
problem for the specific active velocity function given in (1). The travelling
wave solution they single out models ongoing polarisation of cells initiated
by a departing cell sheet (Fig. 2). Its travelling wave speed and profiles $R$
and $A$ (all with subscript $1$) are

  1. S1:

Polarisation wave triggered by departing cell (Fig 2) sheet[17]:
$s_{1}=-\sqrt{\kappa\left(\frac{1}{\alpha}-1\right)}$ with travelling wave
profiles

$A_{1}(z)=\begin{cases}s_{1}\alpha\left(1-\frac{1}{g^{-1}(M-\frac{zs_{1}}{\kappa})}\right)\;,&z<\mbox{0}\\\
1+(\alpha-1)e^{\frac{z}{s_{1}-1}}\;,&z\geq\mbox{0}\end{cases}\quad\text{and}\quad
R_{1}(z)=\begin{cases}g^{-1}\left(M-\frac{zs_{1}}{\kappa}\right)\;,&z<\mbox{0}\\\
\frac{s_{1}}{s_{1}-1}\;,&z\geq\mbox{0}\end{cases}$

where $g(y)=\frac{1}{y}+\log\left(\frac{1}{y}-1\right)$ and
$M=g\left(\frac{s_{1}}{s_{1}-1}\right)$.

Figure 2: (A) shows a simulation of (2) corresponding to the polarisation
caused by a departing cell sheet. (B) shows the corresponding phase plane
diagram of (4) and (C) shows the travelling wave profiles of density R(z) and
polarity A(z).

##  3 Travelling wave solutions

In addition to the travelling wave solution S1 found in [17] we identify other
travelling wave solutions corresponding to polarization and depolarization
waves. The specific travelling wave speeds and profiles may be obtained
through computations which are analogous to those performed in [17].
Alternatively, one may also derive them directly through two different
transformations. These specific transformations of dependent and independent
variables allow us to recast the travelling wave solution S1 into other
travelling wave solutions.

We consider two different transformations:

  1. 1.

One which transforms a solution of system (4) into another solution of the
same system,

$T_{1}:(R(z),A(z),s)\longrightarrow(\bar{R}(\bar{z}),\bar{A}(\bar{z}),\bar{s}),$

where $\bar{A}=A$, $\frac{1}{R}+\frac{1}{\bar{R}}=2$,
$\bar{z}=\int_{0}^{z}\left(1-2\,R(\tilde{z})\right)d\tilde{z}$ and
$\bar{s}=-s$,

  2. 2.

and a transformation which transforms a solution to (3) to another solution of
(3),

$T_{2}:(\rho(x,t),a(x,t))\longrightarrow({\bar{\rho}}({\bar{x}},{\bar{t}}),{\bar{a}}({\bar{x}},{\bar{t}})),$

where ${\bar{t}}=t$, ${\bar{x}}=t-x$, ${\bar{v}}=1-v$, ${\bar{a}}=1-a$,
$\bar{\alpha}=1-\alpha$ and $\bar{\rho}=\rho$.

Specifically for travelling wave solutions the transformation $T_{2}$ is given
by

$\tilde{T}_{2}:(R(z),A(z),s)\longrightarrow(\bar{R}(\bar{z}),\bar{A}(\bar{z}),\bar{s})\;,\quad\text{where}\quad\bar{A}=1-A\;,\quad\bar{R}=R\quad\text{and}\quad\bar{s}=1-s\;.$

Applying these transformations, we obtain three additional travelling wave
solutions in addition to S1 as illustrated in Fig. 3.

Figure 3: Scheme of travelling wave solutions (polarising/depolarising,
departing/colliding) obtain by applied the transfromations $T_{1}$ and
$T_{2}$.

Their specific speeds and profiles are:

  1. S2:

By applying the transformation $T_{1}$ to the travelling wave solution (3), we
obtain the polarisation wave caused by a colliding cell sheet. Its speed is
given by $s_{2}=\sqrt{\kappa\left(\frac{1}{\alpha}-1\right)}$ and the wave
profiles $A_{2}$ and $R_{2}$ (third column of Fig. 4 S2) are given by

$A_{2}(z)=\begin{cases}s_{2}\alpha\left(1-\frac{1}{h^{-1}(N-\frac{zs_{2}}{\kappa})}\right)\;,&z\geq\mbox{0}\\\
1+(\alpha-1)e^{\frac{z}{s_{2}-1}}\;,&z<\mbox{0}\end{cases}\quad\text{and}\quad
R_{2}(z)=\begin{cases}h^{-1}\left(N-\frac{zs_{2}}{\kappa}\right)\;,&z\geq\mbox{0}\\\
\frac{s_{2}}{s_{2}-1}\;,&z<\mbox{0}\end{cases}$

where $h(y)=\frac{1}{y}+\log(1-\frac{1}{y})$ and
$N=h\left(\frac{s_{2}}{s_{2}-1}\right)$. Note that the travelling wave
solution S2 can only be realised if the model parameters $\kappa$ and $\alpha$
are such that its wave speed $s_{2}>1$. If that is not the case the
mathematical solution is not physical and violates the impenetrability of
single cells (see supplementary section C).

  2. S3:

After applying $T_{2}$ to the solution S1, we find the depolarization wave due
to a departing cell sheet. The profiles $A_{3}$ and $R_{3}$ (third column of
Fig. 4 S3) are given by

$A_{3}(z)=\begin{cases}1+\left(1-s_{3}\right)\left(1-\alpha\right)\left(\frac{1}{g^{-1}\left(M+\frac{(1-s_{3})z}{\kappa}\right)}-1\right),&z\geq\mbox{0}\\\
\alpha
e^{\frac{z}{s_{3}}},&z<\mbox{0}\end{cases}\;\text{and}\;R_{3}(z)=\begin{cases}g^{-1}\left(M+\frac{(1-s_{3})z}{\kappa}\right),&z\geq\mbox{0}\\\
\frac{s_{3}-1}{s_{3}},&z<\mbox{0}\end{cases}$

where $M=g\left(\frac{s_{3}-1}{s_{3}}\right)$ and the wave speed is given by
$s_{3}=1+\sqrt{\kappa\left(\frac{1}{1-\alpha}-1\right)}$.

  3. S4:

After applying $T_{2}$ to the solution S2, we find the depolarization wave
caused by a colliding cell sheet. The profiles $A_{4}$ and $R_{4}$ (third
column of Fig. 4 S4) are given by

$A_{4}(z)=\begin{cases}1+\left(1-s_{4}\right)\left(1-\alpha\right)\left(\frac{1}{h^{-1}\left(N+\frac{(1-s_{4})z}{\kappa}\right)}-1\right),&z<\mbox{0}\\\
\alpha
e^{\frac{z}{s_{4}}},&z\geq\mbox{0}\end{cases}\;\text{and}\;R_{4}(z)=\begin{cases}h^{-1}\left(N+\frac{(1-s_{4})z}{\kappa}\right),&z<\mbox{0}\\\
\frac{s_{4}-1}{s_{4}},&z\geq\mbox{0}\end{cases}$

where $N=h\left(\frac{s_{4}-1}{s_{4}}\right)$ and the wave speed is given by
$s_{4}=1-\sqrt{\kappa\left(\frac{1}{1-\alpha}-1\right)}$. Note that similar to
S2 the travelling wave solution S4 violates the impenetrability of cells if
the model parameters $\kappa$ and $\alpha$ are such that the wave speed
satisfies $s_{4}\geq 0$. (see supplementary section C).

Figure 4: S2: Polarisation wave caused by a colliding cell sheet. S3:
Depolarisation wave due to a departing cell sheet. S4: Depolarisation wave
caused by a colliding cell sheet. The second and third columns of S1, S2 and
S3 show the corresponding phase plane diagrams of solutions to (4) and the
travelling wave profiles respectively.

In the rest of the paper we will be concerned with the stability of the
travelling wave solutions $S1$-$S4$. Note that the transformation $T_{2}$ acts
on solutions of the continuum model (3) and in a vicinity of the travelling
wave solution it is a smooth map. Therefore any perturbation of one the
travelling wave solutions S1 and S2, respectively, translates into a
perturbation of S3 and S4, respectively - and vice versa. Therefore we will
only investigate the spectral stability of the travelling wave solution S1,
which will immediately imply the spectral stability of S3.

The same argument will also apply to S2 and S4, respectively. Yet, since the
transformation T1 only applies to solutions of (4), but not of (3) we have to
investigate the stability of S2 separately along the same lines as for S1. For
this reason this will be added in the supplementary section B.

##  4 Numerical Simulation

To explore the stability of the travelling wave solution S1 we compute
numerical solutions of (3) numerically using a Lax-Friedrichs scheme [23].
Typically we use very fine spatial grids to minimise the approximation error.

We are interested in the question whether the choice of the threshold polarity
$0<\alpha<1$ which models the sensitivity of cells to polarisation affects the
stability of the travelling wave solutions. Running the simulations starting
with the initial condition given by the travelling wave solution S1 shows that
for both, small values of $\alpha$ and large values of $\alpha$, the
travelling wave solution is stable (Fig 5(A,B) and Fig 5(C,D)).

Figure 5: Travelling wave solution S1 for different parameter values $\alpha$
and $\kappa$.

We take this as an indication that that the travelling wave S1 (and therefore
S3) is stable for all parameter values. Further below we will investigate this
question using spectral analysis of the linearised operator.

Figure 6: (A): Travelling wave for small $\alpha$. (B): Depolarization
travelling wave due to large $\alpha$.

We also perform a numerical experiment in order to explore whether the size of
the attractive regions of the polarisation wave S1 and the respective
depolarisation wave S3 may depend on the sensitivity $\alpha$. To this end we
simulate (3) starting with a given initial condition that is distinct from the
two travelling wave profiles, namely a step function centred at $x=0$. The
simulations for small $\alpha$ (Fig 6(A)) illustrate the convergence of the
solution towards the polarisation wave. The same simulation for a large value
of $\alpha$, however, (Fig 6(B)) shows that the solution readjusts and
exhibits a depolarisation wave in the opposite direction. Taking into account
potential approximation errors due to the numerical discretisation, we find
that the threshold value $\bar{\alpha}$ (i.e. convergence to polarisation wave
if $\alpha<\bar{\alpha}$, otherwise convergence to depolarisation wave)
converges to about $\bar{\alpha}=0.79$ as the spatial grid is gradually tuned
finer (Fig 7).

Figure 7: The threshold value $\bar{\alpha}$ such that the solution converges
to the polarisation wave if $\alpha<\bar{\alpha}$, and to the depolarisation
wave otherwise.

##  5 Stability Analysis

We now investigate the spectral stability of these travelling wave solutions.
To this end we introduce the perturbations of density $\delta\rho$, polarity
$\delta a$ and velocity $\delta v$ and linearise the system (3) in the moving
coordinate frame $(t,z=x-st)$ at the travelling wave profile S1, i.e.
$s=s_{1}$ and in what follows we write $R=R_{1}$, $A=A_{1}$ and $V$ for the
associated velocity field. The linearised system of equations is given by

$\displaystyle\begin{pmatrix}\partial_{t}\delta\rho\\\ \partial_{t}\delta
a\end{pmatrix}$ $\displaystyle=\mathcal{L}\begin{pmatrix}\delta\rho\\\ \delta
a\end{pmatrix}=\begin{pmatrix}s\delta\rho^{\prime}-(V\delta\rho+R\delta
v)^{\prime}\\\ s\delta a^{\prime}-\delta v(A^{\prime}-1)-V\delta
a^{\prime}-\delta a\end{pmatrix}\;,$ (5) $\displaystyle\text{where}\quad\delta
v=M^{\prime}(A)\delta
a-\kappa\left(\frac{1}{R^{3}}\delta\rho\right)^{\prime}\;,$

and where we use the notation ${}^{\prime}=\partial_{z}$. To obtain the
associated eigenvalue problem $(\mathcal{L}-\lambda)(\delta\rho,\delta
a)^{T}=0$, consider the solution of the above system (5) to be
$\delta\rho(z,t)=e^{\lambda t}\delta\rho(z)$ and $\delta a(z,t)=e^{\lambda
t}\delta a(z).$ This converts the eigenvalue problem into a system of first
order ODEs,

$\begin{pmatrix}\delta\rho^{\prime}\\\ \delta v^{\prime}\\\ \delta
a^{\prime}\end{pmatrix}=A(z,\lambda)\begin{pmatrix}\delta\rho\\\ \delta v\\\
\delta a\end{pmatrix}\;,\quad\text{where}\quad
A(z,\lambda)=\begin{pmatrix}\frac{3\,R^{\prime}}{R}&\frac{-R^{3}}{\kappa}&\frac{M^{\prime}\,R^{3}}{\kappa}\\\
\frac{2\,R^{\prime}s-R^{2}\,\lambda}{R^{3}}&\frac{-R^{\prime}}{R}-\frac{R\,s}{\kappa}&\frac{M^{\prime}\,R\,s}{\kappa}\\\
0&\frac{(A^{\prime}-1)\,R}{s}&\frac{R\,(\lambda+1)}{s}\end{pmatrix}\;,$ (6)

where we used that the velocity field (third equation in (3)) for the
travelling wave solution S1 satisfies

$V=s(1-1/R)\quad\text{and}\quad V^{\prime}=s\,R^{\prime}/R^{2}\;.$ (7)

We then denote the corresponding linear operator by $\mathcal{T}$, where
$\mathcal{T}(y)(z)=(\frac{\partial}{\partial z}-A(z,\lambda))y$ and
$y=(\delta\rho,\delta v,\delta a)^{T}$.

It is our goal to find the spectrum of the linearized operator $\mathcal{L}$
which is an operator on $H^{2}\times H^{1}$ mapping into $L^{2}\times L^{2}$
[20, 24]. If $(\mathcal{L}-\lambda)^{-1}$ is unbounded or does not exist for
$\lambda\in\mathbb{C}$, then $\lambda$ is in the spectrum
$\sigma(\mathcal{L})$. The Fredholm index of $\mathcal{L}$ is
$ind(\mathcal{L})=dim[ker(\mathcal{L})]-codim[R(\mathcal{L})]$, where
$R(\mathcal{L})$ and $ker(\mathcal{L})$ denote the range and kernel of
$\mathcal{L}$ respectively [24]. The spectrum of a Fredholm operator
$\mathcal{L}$ is composed of two disjoint sets, namely point spectrum and
essential spectrum. The point spectrum will consist of values
$\lambda\in\sigma(\mathcal{L})$ such that $(\mathcal{L}-\lambda)$ is a
Fredholm operator of index zero. The essential spectrum is the complement of
the point spectrum in $\sigma(\mathcal{L})$.

Figure 8: The essential spectrum of $\mathcal{L}$ (for travelling wave S1) is
given by $\lambda$ along the Fredholm borders, in this case the union of a
straight line ($\lambda_{1}$) and a parabola ($\lambda_{2}$) in the complex
plane.

###  5.1 Essential spectrum

The essential spectrum is the spectrum up to relatively compact perturbations.
To this end we consider the asymptotic matrices
$A^{\pm}(\lambda)=\lim_{z\to\pm\infty}A(z,\lambda)$ given by

$A^{+}(\lambda)=\begin{pmatrix}0&&\frac{-s^{3}}{\kappa\,(s-1)^{3}}&&0\\\
(\frac{1}{s}-1)\,\lambda&&\frac{s^{2}}{\kappa\,(1-s)}&&0\\\
0&&\frac{1}{1-s}&&\frac{\lambda+1}{s-1}\end{pmatrix}\quad\text{and}\quad
A^{-}(\lambda)=\begin{pmatrix}0&&\frac{-1}{\kappa}&&0\\\
-\lambda&&-\frac{s}{\kappa}&&0\\\
0&&-\frac{1}{s}&&\frac{\lambda+1}{s}\end{pmatrix}.$ (8)

The asymptotic operator of $\mathcal{T}$ is given by

$\mathcal{T}_{\infty}\begin{pmatrix}\delta\rho\\\ \delta v\\\ \delta
a\end{pmatrix}=\begin{pmatrix}\delta\rho^{\prime}\\\ \delta v^{\prime}\\\
\delta a^{\prime}\end{pmatrix}-A_{\infty}\begin{pmatrix}\delta\rho\\\ \delta
v\\\ \delta a\end{pmatrix}\;,\quad\text{where}\quad
A_{\infty}=\begin{cases}A_{-}(\lambda)&z<0\;,\\\ A_{+}(\lambda)&z\geq
0\;.\end{cases}$ (9)

Note that the opterator $\mathcal{L}$ is a relatively compact perturbation of
the asymptotic operator $\mathcal{L}_{\infty}$ which is defined as the limit
of $\mathcal{L}$ as $z\to\pm\infty$ and which is equivalent to the operator
$\mathcal{T_{\infty}}$ ([24], Theorem 3.1.11), see also [25]). According to
Weyl’s Essential Spectrum Theorem ([24], Theorem 2.2.6), the operators
$\mathcal{L}$ and $\mathcal{L}_{\infty}$ have the same spectra, i.e.,
$\sigma_{ess}(\mathcal{L})=\sigma_{ess}(\mathcal{L}_{\infty})$, or
equivalently $\sigma_{ess}(\mathcal{T})=\sigma_{ess}(\mathcal{T}_{\infty})$.

Exponential dichotomies can be used to characterise the spectrum of an
operator. According to this concept each solution to (5) decays exponentially
either for $z\rightarrow\infty$ or for $z\rightarrow-\infty$ [25]. For
spatially constant matrices the presence of an exponential dichotomy implies
that the matrix is hyperbolic. The Morse indices of the constant matrices
$A_{\pm}(\lambda)$ are defined as the dimension of their unstable subspaces
written as $i_{\pm}(\lambda)$, respectively. For $\lambda\in\mathbb{C}$ such
that $\mathcal{T_{\infty}}$ is Fredholm, we have
ind$(\mathcal{T_{\infty}}-\lambda)=i_{-}(\lambda)-i_{+}(\lambda)$ ([24], Lemma
3.1.10). As a consequence we can define the essential spectrum of
$\mathcal{L_{\infty}}$ as ([24])

$\sigma_{ess}(\mathcal{L_{\infty}})=\\{\lambda\in\mathbb{C}\mid
i_{-}(\lambda)\neq i_{+}(\lambda)\\}\cup\\{\lambda\in\mathbb{C}\mid
dim\;\mathbb{E}^{c}(A_{\pm}(\lambda))\neq 0\\}\;,$ (10)

where $\mathbb{E}^{c}$ denotes the centre subspace associated to the
asymptotic linearised system.

The dispersion relations of $A^{-}$ and $A^{+}$ which characterise
$dim\;\mathbb{E}^{c}(A_{\pm}(\lambda))\neq 0$ are defined by

$\det(A^{-}-i\,\mu\,\rm{Id})=0\quad\text{and}\quad\det(A^{+}-i\,\mu\,\rm{Id})=0\;,\quad\text{where}\quad\mu\in\mathbb{R}\;.$
(11)

These relations of both $A^{-}$ and $A^{+}$, after rescaling the spatial
eigenvalue $\mu$, coincide and are given by

$\lambda^{1}=-1+i\mu s\quad\text{and}\quad\lambda^{2}=-\mu^{2}\kappa+i\mu
s\;.$ (12)

Note that the matrix eigenvalues of $A^{-}(\lambda)$ defined in (8) are given
by

$\mu_{1}(\lambda)=\frac{\lambda+1}{s}\;,\quad\mu_{2}(\lambda)=\frac{-s}{2\kappa}+\frac{\sqrt{s^{2}+4\kappa\lambda}}{2\kappa}\;,\quad\mu_{3}(\lambda)=\frac{-s}{2\kappa}-\frac{\sqrt{s^{2}+4\kappa\lambda}}{2\kappa}\;.$
(13)

and those of $A^{+}(\lambda)$ are

$\displaystyle\mu_{1}(\lambda)$
$\displaystyle=\frac{\lambda+1}{s-1}\quad\text{and}\quad\mu_{2,3}(\lambda)=\frac{s^{2}}{2\kappa(1-s)}\mp
s\,\frac{\sqrt{s^{2}+4\kappa\lambda}}{2\kappa(1-s)}\;,$ (14)

which are multiples of the matrix eigenvalues of $A^{-}(\lambda)$ by
$\frac{s}{s-1}>0$. As a consequence the Morse indices coincide for all
$\lambda$, i.e. $i_{-}(\lambda)=i_{+}(\lambda)$. Therefore the essential
spectrum only consists of the the Fredholm borders (12) which are visualised
in Fig. 8.

###  5.2 Absolute spectrum

While the absolute spectrum is not spectrum [26], it gives information about
the stability of the operator $\mathcal{L}$ in exponentially weighted spaces.
Most notably it tells us how far the essential spectrum may be shifted to the
left by considering weighted function spaces with exponential weights at
$\pm\infty$. If the absolute spectrum lies completely in the open left half-
plane of the complex plane, then we say $\mathcal{L}$ is absolutely stable,
otherwise absolutely unstable.

We note that common value of the of the Morse index for $\lambda\gg 1$ by
$i_{\infty}$, i.e. $i_{\infty}=i_{\pm}(\lambda)$. We define the matrix
(spatial) eigenvalues $\mu_{\pm}^{j}\,,\,j=1,...,n\,$ of the asymptotic
matrices $A_{\pm}(\lambda)$ ordered according to the size of their real
parts[24],

$Re\,\mu_{\pm}^{1}(\lambda)\geq.....\geq\,Re\,\mu_{\pm}^{l}(\lambda)\geq.....\geq\,Re\,\mu_{\pm}^{n}(\lambda)\;,$

and introduce the stable and unstable extrema

$\mu_{\pm}^{u}(\lambda)=Re\,\mu_{\pm}^{i_{\infty}}(\lambda)\;,\quad\mu_{\pm}^{s}(\lambda)=Re\,\mu_{\pm}^{i_{\infty}+1}(\lambda)\;.$

Hence $\mu_{\pm}^{u}(\lambda)$ denotes the smallest (positive) real part of
any of the matrix eigenvalues, and $\mu_{\pm}^{s}(\lambda)$ denotes the
largest(negative) real part of any of the matrix eigenvalues, of the
asymptotic matrices $A_{\pm}(\lambda)$. Here we have
$\mu_{\pm}^{u}(\lambda)>0>\mu_{\pm}^{s}(\lambda)$ for $Re\,\lambda\gg 1$. As
we move the eigenvalue $\lambda$ towards the Fredholm border coming from the
very far right of the complex plane, then at least one of the matrix
eigenvalues will become close to the imaginary axis. Moreover, for the
weighted spaces, when the distance between the spatial eigenvalues,
$\mu_{\pm}^{u}\,(\lambda)-\mu_{\pm}^{s}\,(\lambda)$, becomes zero, we can not
choose a weight that renders the operator $\mathcal{L}-\lambda$ Fredholm with
index zero. This motivates the following definition.

The subset $\sum_{abs}^{+}$ of $\mathbb{C}$ consists exactly of those
$\lambda$ for which
$Re\,\mu_{+}^{i_{\infty}}(\lambda)\,=\,Re\,\mu_{+}^{i_{\infty}+1}(\lambda)\,$.
Analogously, $\lambda$ is in $\sum_{abs}^{-}$ if, and only if,
$Re\,\mu_{-}^{i_{\infty}}(\lambda)\,=\,Re\,\mu_{-}^{i_{\infty}+1}(\lambda)\,$.
Finally, we say that $\lambda$ is in the absolute spectrum $\sum_{abs}$ of an
operator if $\lambda$ is in $\sum_{abs}^{+}$ or in $\sum_{abs}^{-}$ (or in
both).

Figure 9: Absolute spectrum shown in red (for $s=-2$, $\kappa=1$).

For $Re(\lambda)\gg 1$ the matrix eigenvalues of $A^{-}(\lambda)$ given in
(13) satisfy $Re(\mu_{2})>0$ as well as $Re(\mu_{3})<0$ and $Re(\mu_{1})<0$
(since $s<0$). In order to find which of the two asymptotically negative
spatial eigenvalues is larger (real part) for a given $\lambda$, we solve the
inequality $Re(\mu_{1})\leq Re(\mu_{3})$, i.e.,

$\displaystyle Re\left(\frac{1+\lambda_{1}+i\lambda_{2}}{s}\right)$
$\displaystyle\leq
Re\left(\frac{-s}{2\kappa}-\frac{\sqrt{s^{2}+4\kappa(\lambda_{1}+i\lambda_{2})}}{2\kappa}\right).$

In this case we get, writing
$S(\lambda_{1})=\frac{(s^{2}+2\kappa\,(1+\lambda_{1}))\,\sqrt{(s^{2}+\kappa\,(1+\lambda_{1})^{2})}}{s^{2}\sqrt{\kappa}}$,
that

$\Omega_{1}=\left\\{\lambda=\lambda_{1}+i\,\lambda_{2}\,\,\text{such
that}\,\,\lambda_{1}\in\left[-\frac{\,s^{2}+2\kappa}{2\kappa},\infty\right)\,\,\text{and}\,\,\lambda_{2}\in\left[-S(\lambda_{1}),S(\lambda_{1})\right]\right\\}$

(hatched area in Fig. 9) and $Re(\mu_{1})>Re(\mu_{3})$ is satisfied by all
$\lambda$ in the complement $\Omega_{2}=(\mathbb{C}\setminus\Omega_{1})$ of
this set. As a consequence the absolute spectrum is the set of all $\lambda$
such that $Re(\mu_{2})=Re(\mu_{3})$ in $\Omega_{1}$ and
$Re(\mu_{2})=Re(\mu_{1})$ in $\Omega_{2}$.

Now for the absolute spectrum in the set $\Omega_{1}$, we solve the equation
$Re(\mu_{2})=Re(\mu_{3})$, i.e.

$\displaystyle
Re\left(\frac{-s}{2\kappa}+\frac{\sqrt{s^{2}+4\kappa\lambda}}{2\kappa}\right)$
$\displaystyle=Re\left(\frac{-s}{2\kappa}-\frac{\sqrt{s^{2}+4\kappa\lambda}}{2\kappa}\right),$

which implies that
$Re(\lambda)\in\left[-\frac{s^{2}+2\kappa}{2\kappa},-\frac{s^{2}}{4\kappa}\right]$
and $Im(\lambda)=0$.

To obtain the absolute spectrum in the set $\Omega_{2}$, we solve the equation
$Re(\mu_{1})=Re(\mu_{3})$, i.e.

$\displaystyle Re\left(\frac{1+\lambda_{1}+i\lambda_{2}}{s}\right)$
$\displaystyle=Re\left(\frac{-s}{2\kappa}+\frac{\sqrt{s^{2}+4\kappa(\lambda_{1}+i\lambda_{2})}}{2\kappa}\right),$

which is the case for all $\lambda=\lambda_{1}+i\,\lambda_{2}$ such that
$\lambda_{1}\,<-\frac{\,s^{2}+2\kappa}{2\kappa}\,\,\text{and}\,\,\lambda_{2}=\pm\frac{(s^{2}+2\kappa\,(1+\lambda_{1}))\,\sqrt{(s^{2}+\kappa\,(1+\lambda_{1})^{2})}}{s^{2}\sqrt{\kappa}}$.
Thus the absolute spectrum of $A^{-}(\lambda)$ is given by

$\sigma_{abs}^{-}=\left[-\frac{s^{2}+2\kappa}{2\kappa},-\frac{s^{2}}{4\kappa}\right]\,\cup\Bigg{\\{}\lambda=\lambda_{1}+i\,\lambda_{2}\,\mid\lambda_{1}\,<-\frac{\,s^{2}+2\kappa}{2\kappa}\,\,\text{and}\\\
\lambda_{2}=\pm\frac{(s^{2}+2\kappa\,(1+\lambda_{1}))\,\sqrt{(s^{2}+\kappa\,(1+\lambda_{1})^{2})}}{s^{2}\sqrt{\kappa}}\Bigg{\\}}.$
(15)

Finally the matrix eigenvalues of $A^{+}(\lambda)$ given in (14) are a
multiple of the matrix eigenvalues of $A^{-}(\lambda)$, namely by
$\frac{s}{s-1}>0$. So the absolute spectrum of $A^{+}(\lambda)$ is equal to
the absolute spectrum of $A^{-}(\lambda)$ and the absolute spectrum is given
by (15) (see Fig. 9).

This implies in a conveniently chosen weighted space, and without the presence
of any point spectrum with non-negative real part, save for at $\lambda=0$
(see Section 5.3), the wave is spectrally stable. To determine the weights for
which we potentially have spectral stability, we follow [27, 26] looking for a
so-called ideal weight. The ideal weight will be the (two-sided) weight which
maximises the resolvent set. We find the ideal weight $\eta^{*}_{+}$ to
satisfy $\textrm{Re}(\eta^{*}_{+})=\frac{s^{2}}{2\kappa(1-s)}$, for the
operator corresponding to $A^{+}(\lambda)$, and
$\textrm{Re}(\eta^{*}_{-})=\frac{-s}{2\kappa}$ for the operator corresponding
to $A^{-}(\lambda)$. We note that these also differ by a factor of
$\frac{s}{s-1}$. By considering perturbations that decay at least like
$\exp{(-\eta z)}$ as $z\to\infty$ with $\eta\in(0,\frac{s^{2}}{\kappa(1-s)}]$,
and like $\exp{\eta z}$ as $z\to-\infty$, with $\eta\in(0,-\frac{s}{\kappa}]$,
we have that the essential spectrum will be contained in the (open) left half
plane. This is worth noting because the derivative of the S1 wave decays as
$z\to+\infty$ like $\exp{\left(\frac{z}{s-1}\right)}$, and as $z\to-\infty$,
like $\exp{\left(\frac{-3sz}{\kappa}\right)}$ and so will remain in the
weighted space for the weights we want to consider. This means that
$\lambda=0$ will still be an eigenvalue both in $L^{2}$ and in the weighted
space. As we shall see in the next section, $\lambda=0$ is the only element of
the point spectrum that we can numerically find.

###  5.3 Point Spectrum

In this section we investigate the point spectrum of $\mathcal{L}$. We compute
the Evans function whose zeros in the complex plan characterise the point
spectrum. In 1972, J. W. Evans used this technique to investigate the
stability of the solution to equations modelling the nerve axon [28, 29]. He
showed that the Evans function $D(\lambda)$ is always analytic to the right
side of the essential spectrum. In addition it has an analytic extension up to
the absolute spectrum (Fig. 9).

To compute the Evans function, we start by identifying the stable and unstable
eigenvectors of the asymptotic matrices $A^{\pm}(\lambda)$ respectively in an
attempt to construct an integrable function which plays the role of an
eigenvector for the linearised system.

The unstable eigenvalue of $A^{-}$ and the stable eigenvalues of $A^{+}$ are
given by

$\mu_{2}^{-}(\lambda)=\frac{-s+\sqrt{s^{2}+4\kappa\lambda}}{2\kappa}\;,\quad\mu_{1}^{+}(\lambda)=\frac{\lambda+1}{s-1}\quad\text{and}\quad\mu_{3}^{+}(\lambda)=\frac{s^{2}+s\sqrt{s^{2}+4\kappa\lambda}}{2\kappa(1-s)}\;,$
(16)

which satisfy $\mu_{2}^{-}(\lambda)>0$ and
$\mu_{1}^{+}(\lambda)\,,\mu_{3}^{+}(\lambda)<0$ for ${\rm Re}\,\lambda\gg 1$
since that travelling wave speed is negative, $s=s_{1}<0$.

We introduce the eigenvectors $X_{0}^{-}(\lambda)$ associated to the
eigenvalue $\mu_{2}^{-}(\lambda)$ of $A^{-}$ as well as $X_{0}^{+}(\lambda)$
and $Y_{0}^{+}(\lambda)$ being eigenvectors of $A^{+}$ associated to
$\mu_{1}^{+}(\lambda)$ and $\mu_{3}^{+}(\lambda)$.

To construct an eigenfunction of the linearised system we solve (6) using
these vectors as initial, respectively terminal conditions until $z=0$. Then
the Evans function is defined as the Wronskian

$D(\lambda)={\rm
det}\;[X^{-}(z=0,\lambda),X^{+}(z=0,\lambda),Y^{+}(z=0,\lambda)]\;,$ (17)

which vanishes if the stable and unstable eigenvectors propagated to $z=0$ are
linearly dependent and can be combined into a smooth eigenfunction.

We start by computing $X^{+}(z=0,\lambda)$ and $Y^{+}(z=0,\lambda)$ which can
be done in terms of a closed-form expression. Only for the computation of
$X^{-}(z=0,\lambda)$ we will resort to numerical results.

For $z>0$ it holds that $R\equiv s/(s-1)$ and the matrix $A(z,\lambda)$ from
(6) is given by

$A(z,\lambda)=\begin{pmatrix}0&&\frac{-s^{3}}{\kappa\,(s-1)^{3}}&&0\\\
(\frac{1}{s}-1)\,\lambda&&\frac{s^{2}}{\kappa\,(1-s)}&&0\\\
0&&\frac{A^{\prime}-1}{s-1}&&\frac{\lambda+1}{s-1}\end{pmatrix}\,.$ (18)

Note that due to the zeros in the third column the equations for $\delta\rho$
and $\delta v$ are not coupled to $\delta a$. They satisfy

$\begin{pmatrix}\delta\rho^{\prime}\\\ \delta
v^{\prime}\end{pmatrix}=\begin{pmatrix}0&&\frac{-s^{3}}{\kappa\,(s-1)^{3}}\\\
(\frac{1}{s}-1)\,\lambda&&\frac{s^{2}}{\kappa\,(1-s)}\end{pmatrix}\begin{pmatrix}\delta\rho\\\
\delta v\end{pmatrix}\;,$ (19)

which is a system of two linear, constant-coefficient equations. It admits one
fundamental solution with negative eigenvalue (we omit the unstable
fundamental solution) given by

$\begin{pmatrix}\delta\rho\\\ \delta
v\end{pmatrix}=C_{3}\,\exp(\mu_{3}^{+}(\lambda)\,z)\begin{pmatrix}\frac{s^{2}\,(-s+\sqrt{s^{2}+4\,\kappa\,\lambda})}{(s-1)^{2}\,\kappa}\\\
2\,\lambda\end{pmatrix}$ (20)

With this information we can rewrite the equation for $\delta a$ which is
contained in the third row of (18),

$\displaystyle\delta a^{\prime}$
$\displaystyle=\frac{A^{\prime}-1}{s-1}\,\delta
v+\frac{\lambda+1}{s-1}\,\delta a$
$\displaystyle=\frac{A^{\prime}-1}{s-1}\,2\,\lambda\,C_{3}\,e^{\mu_{3}\,z}+\mu_{1}^{+}\,\delta
a\;.$

The general solution is given by

$\delta
a=C_{1}e^{\mu_{1}^{+}z}+\frac{2\lambda}{s-1}C_{3}\int_{0}^{z}e^{(\mu_{3}^{+}-\mu_{1}^{+})z}\left(\frac{\alpha-1}{s-1}e^{\frac{z}{s-1}}-1\right)\,d\tilde{z}$

This implies that

$\lim_{z\to 0^{+}}\begin{pmatrix}\delta\rho\\\ \delta v\\\ \delta
a\end{pmatrix}=\begin{pmatrix}C_{3}\,\frac{s^{2}\,(-s+\sqrt{s^{2}+4\,\kappa\,\lambda})}{(s-1)^{2}\,\kappa}\\\
C_{3}\,2\,\lambda\\\ C_{1}\end{pmatrix}\;.$

While most components of $A(z,\lambda)$ in (6) are functions, the derivative
of the active speed of migration (1) is a $\delta$-distribution,
$M^{\prime}(a)=\delta(a-\alpha)$. The change of variables between $z$ and the
monotone function $a=A(z)=A_{1}(z)$ shows that when integrating with respect
to $z$ the following expression is a $\delta$-distribution centred at $z=0$,
$M^{\prime}(A(z))A^{\prime}(z)=\delta(z)$. Therefore $\delta\rho$ and $\delta
v$ in the solution of (6) undergo a jump at $z=0$ which involves the factor
$1/A^{\prime}(0)$ (for details see appendix A). The left limit as
$z\rightarrow 0$ of the solution vector is given by

$\displaystyle\lim_{z\to 0^{-}}\begin{pmatrix}\delta\rho\\\ \delta v\\\ \delta
a\end{pmatrix}=$
$\displaystyle\begin{pmatrix}C_{3}\,\frac{s^{2}\,(-s+\sqrt{s^{2}+4\,\kappa\,\lambda})}{(s-1)^{2}\,\kappa}\\\
C_{3}\,2\,\lambda\\\
C_{1}\end{pmatrix}-\frac{C_{1}}{A^{\prime}(0)}\frac{1}{\kappa}\begin{pmatrix}R^{3}(0)\\\
R(0)\,s\\\ 0\end{pmatrix}$ $\displaystyle=$ $\displaystyle
C_{3}\,\begin{pmatrix}\frac{s^{2}\,(-s+\sqrt{s^{2}+4\,\kappa\,\lambda})}{(s-1)^{2}\,\kappa}\\\
2\,\lambda\\\
0\end{pmatrix}+\,C_{1}\,\begin{pmatrix}-\frac{s^{3}}{(s-1)^{2}\,(\alpha-1)\,\kappa}\\\
-\frac{s^{2}}{(\alpha-1)\,\kappa}\\\ 1\end{pmatrix}.$

where we take the values of $R(0)$ and $A^{\prime}(0)$ from the travelling
wave profile S1.

We obtain the closed-form solutions $X^{+}(0,\lambda)$ and $Y^{+}(0,\lambda)$
corresponding to the stable eigenvalues $\mu_{1}^{+}(\lambda)$ and
$\mu_{3}^{+}(\lambda)$,

$\displaystyle
X^{+}(0,\lambda)=\begin{pmatrix}-\frac{s^{3}}{(s-1)^{2}\,(\alpha-1)\,k}\\\
-\frac{s^{2}}{(\alpha-1)\,k}\\\ 1\end{pmatrix}\quad\text{and}\quad
Y^{+}(0,\lambda)=\begin{pmatrix}\frac{s^{2}\,(-s+\sqrt{s^{2}+4\,k\,\lambda})}{(s-1)^{2}\,k}\\\
2\,\lambda\\\ 0\end{pmatrix}.$ (21)

Finally, in order to evaluate the Evans function (17) for a given
$\lambda\in\mathbb{C}$, we compute $X^{-}(0,\lambda)$ propagating the
eigenvector associated to the spatial eigenvalue $\mu_{2}$ given by

$v_{2}=\begin{pmatrix}\frac{-s^{2}+\sqrt{s^{4}+4s^{2}\kappa\lambda}}{2\lambda
s\kappa}\\\ 1\\\
\frac{2\kappa}{s^{2}+2\kappa(1+\lambda)+\sqrt{s^{4}+4s^{2}\kappa\lambda}}\end{pmatrix}\;\text{for}\;\lambda\neq
0\quad\text{and}\quad v_{2}=\begin{pmatrix}\frac{1}{s}\\\ 1\\\
\frac{\kappa}{s^{2}+\kappa}\end{pmatrix}\;\text{for}\;\lambda=0$

from an arbitrary small value of $z$ (we choose $z=-20$) until $z=0$ solving
the system (6) numerically.

We apply the Argument principle to identify zeros of the Evans function [25].
The argument principle states that

$W=\frac{1}{2\pi
i}\,\oint_{C}\,\frac{D^{\prime}(\lambda)}{D(\lambda)}\,d\lambda$

corresponds to the winding number around the origin of the image of $C$ under
the map $D$. Here $C$ is a closed curve, oriented in a counterclockwise
direction. Since we choose the contour $C$ to the right of the absolution
spectrum the Evans function $D(\lambda)$ is analytic on and inside $C$ and the
the winding number $W$ corresponds to the number of zeros of $D(\lambda)$
inside $C$. Note that $\lambda=0$ is always an eigenvalue corresponding to the
propagation of the travelling wave profile [20].

Figure 10: (A) and (B) show the contours $C_{1}$ and $C_{2}$ including and
excluding the origin, respectively.

We take the contours $C_{1}$ and $C_{2}$ such that these are always to the
right of the absolute spectrum. For the contour $C_{1}$, we take the line $l$,
parallel to the imaginary axis and right to the branch point with distance
$d_{l}$ from the axis and the semicircle right to the origin has radius $r$.
In the right hand side of the complex plane, we draw the contour $C_{2}$ with
the radii $r_{i}$ and $r_{o}$ of the inner and outer semicircles respectively
excluding the origin.

Figure 11: $D(C_{1})$ for $\alpha=0.2\,,0.4\,,0.5\,\text{and}\,0.7$.

In figure 11, we draw the image of $C_{1}$ (where $d_{l}=-0.05$ and $r=0.1$)
under $D(\lambda)$. According to the argument principle, the winding number
for different $\alpha$’s is one, so there is only one zero which is at
$\lambda=0$. Hence we see the $0$ is the only eigenvalue close to the origin.

Now we investigate whether there are eigenvalues on the right side of the
complex plane. To do this we define a semicircle excluding the origin. We draw
the contour $C_{2}$ (see Fig. 10 (B)).

Figure 12: $D(C_{2})$ for $\alpha=0.2\,,0.4\,,0.5$ and $0.7$.

In figure 12, we see that the winding number about the origin is zero for the
contour $C_{2}$ with $r_{i}=0.1$ and $r_{o}=5$.

Figure 13: $D(C_{1})$ for $\kappa=5$. Figure 14: $D(C_{2})$ for $\kappa=5$.

This illustrates that there is no point spectrum on the right side of the
complex plane other than the eigenvalue at $\lambda=0$ which corresponds to
the translocation of the travelling wave. This indicates that the travelling
wave solution S1 is spectrally stable.

##  6 Conclusion

We study the travelling wave solutions of a 1D model for collective cell
migration in an epithelial layer. We identify four different travelling waves,
two of them corresponding to gradual polarisation of the cells and the other
two corresponding to depolarisation. These travelling wave solutions are
related by the two different transformations $T1$ and $T2$.

We apply the transformation $T1$ obtaining the polarisation, respectively
depolarisation waves associated to departing cell sheets from the
(de)polarisation waves of colliding cell sheets and vice versa. The
transformation $T2$ can be used to turn polarisation waves into depolarisation
waves and vice versa.

A preliminary test of the stability of these travelling waves computing
numerical solutions of the underlaying PDE suggest that they are stable for
all parameter values. Yet, we also illustrate that the threshold polarity
$\alpha$ which represents the sensitivity to polarisation has a significant
impact on the respective sizes of the domains of convergence of the
polarisation, respectively depolarisation waves.

Using spectral theory and the Evans function we investigate the stability of
travelling waves. This analysis confirms the spectral stability of the
travelling wave solutions. While the essential spectrum touches the origin, we
show that the absolution spectrum is on the left hand side of the complex
plane. As a consequence in appropriately chosen weighted spaces the essential
spectrum is contained in the left hand side of the complex plane. The
evaluation of the Evans function is done in a way which combines explicit
computations and numerics, and illustrates that the only eigenvalue is the
simple eigenvalue $\lambda=0$ associated with the propagation of the
travelling wave.

We perform this analysis in the context of the polarisation wave associated to
departing cell sheets. Due to the smoothness of the map $T_{2}$ the stability
of the polarisation wave also translates into stability of the associated
depolarisation wave S3. We add the analogous analysis for the polarisation
wave associated to colliding cell sheets in the supplementary material, which
through the map $T_{2}$ also applies to S4.

In conclusion, we find that the travelling waves of cell polarisation (or
depolarisation) arising in cell migration due to departing or colliding cell
sheets are always stable. For biological tissues this implies that for both,
sensitivity to polarisation ($\alpha$) and strength of intercellular
mechanical interaction ($\kappa$), there is no absolute threshold such that
gradual recruitment of cells into migration, respectively the gradual
transition into a non-migratory rest state comes to a halt. Yet, the
preliminary numerical experiments (Fig. 6) indicate that the sensitivity to
polarisation $\alpha$ characterises behaviour which is reminiscent of a domain
of attraction, i.e. lower sensitivity to polarisation (high $\alpha$) renders
polarisation waves more prone to disruption.

## Acknowledgments

DO was supported by ARC Discovery Project DP180102956. RM was supported by ARC
Discovery Project DP200102130. NR was supported by an RTP scholarship funded
by the University of Queensland (UQ).

The authors are very grateful to Zoltan Neufeld (UQ) and Hamid Khataee (UQ)
for useful conversations and suggestions.

##  Appendix A Solution of (6) around $z=0$.

The first component of Eq. (6) is given by

$\displaystyle\delta\rho^{{}^{\prime}}$
$\displaystyle=\frac{3R^{{}^{\prime}}}{R}\delta\rho-\frac{R^{3}}{k}\delta
v+\frac{M^{\prime}(A)R^{3}}{k}\delta a\;,$

which we write as the 1st order ODE

$\displaystyle\delta\rho^{\prime}-C(z)\delta\rho$ $\displaystyle=D(z)\;,$

where $C(z)=\frac{3R^{{}^{\prime}}}{R}$ and $D(z)=-\frac{R^{3}}{k}\delta
v+\frac{M^{\prime}(A)R^{3}}{k}\delta a$. Note that the definition of $M=M(A)$
in (1) implies that $M^{\prime}(A)=\delta(A-\alpha)$. Multiplying both sides
by the integrating factor $e^{-\int_{0}^{z}C(\tilde{z})d\tilde{z}}$ we find

$\displaystyle(e^{-\int_{0}^{z}C(\tilde{z})d\tilde{z}}\,\delta\rho)^{\prime}$
$\displaystyle=D(z)\,e^{-\int_{0}^{z}C(\tilde{z})d\tilde{z}}\;.$

Integrating on $[-\varepsilon$, $\varepsilon]$ for a small value
$\varepsilon>0$ we obtain

$\displaystyle
e^{-\int_{0}^{\varepsilon}C(\tilde{z})d\tilde{z}}\,\delta\rho(\varepsilon)-e^{-\int_{0}^{-\varepsilon}C(\tilde{z})d\tilde{z}}\,\delta\rho(-\varepsilon)$
$\displaystyle=\int_{-\varepsilon}^{\varepsilon}D(\hat{z})\,e^{-\int_{0}^{\hat{z}}C(\tilde{z})d\tilde{z}}d\hat{z}$
$\displaystyle=\int_{-\varepsilon}^{\varepsilon}\left(-\frac{R^{3}}{k}\delta
v+\frac{R^{3}}{k}\,\delta(A-\alpha)\,\delta
a\right)\,e^{-\int_{0}^{\hat{z}}C(\tilde{z})d\tilde{z}}d\hat{z}\;.$

The change of variables $\hat{a}=A(\hat{z})$ implies

$e^{-\int_{0}^{\varepsilon}C(\tilde{z})d\tilde{z}}\,\delta\rho(\varepsilon)-e^{-\int_{0}^{-\varepsilon}C(\tilde{z})d\tilde{z}}\,\delta\rho(-\varepsilon)=\\\
=\int_{A(-\varepsilon)}^{A(\varepsilon)}\left(-\frac{R^{3}}{k}\delta
v+\frac{R^{3}}{k}\,\delta(\hat{a}-\alpha)\,\delta
a\right)\,\frac{1}{A^{\prime}(A^{-1}(\hat{a}))}e^{-\int_{0}^{A^{-1}(\hat{a})}C(\tilde{z})d\tilde{z}}d\hat{a}\;.$

Note that $A$ is monotonically increasing with $A(0)=\alpha$. Therefore
$\int_{A(-\varepsilon)}^{A(\varepsilon)}\delta(\hat{a}-\alpha)d\hat{a}=1$ and
we obtain in the limit as $\varepsilon\to 0$ that

$\lim_{\varepsilon\to 0^{+}}\delta\rho(\varepsilon)-\lim_{\varepsilon\to
0^{-}}\delta\rho(\varepsilon)=\frac{R^{3}}{k}\,\delta
a\frac{1}{A^{\prime}\,(A^{-1}(\alpha))}=\frac{R^{3}}{k}\,\delta
a\frac{1}{A^{\prime}\,(0))}\;.$

##  Appendix B Stability Analysis of polarisation wave S2

Here we investigate the linear stability of the travelling wave solution S2.
Linearisation of the system (3) in the moving coordinate frame $(t,z=x-st)$ at
the travelling wave profile S2 lead to (5) where $s=s_{2}$ as well as
$R=R_{2}$ and $A=A_{2}$. The associated eigenvalue problem
$(\mathcal{L}-\lambda)(\delta\rho,\delta a)^{T}=0$ is given by (6), (7).

At the far left and right ends the travelling wave profiles $R_{2}$ and
$A_{2}$ converge to $(R_{-\infty},A_{-\infty})=(\frac{s_{2}}{s_{2}-1},1)$ and
$(R_{\infty},A_{\infty})=(1,0)$ which corresponds to the limit values of R1
and A1 ”swapped around”. This implies that the asymptotic matrices
$B^{\pm}(\lambda)=\lim_{z\to\pm\infty}A(z,\lambda)$ for S2 are the same as for
S1 given in (22), however with $+$ and $-$ inverted.

$B^{-}(\lambda)=\begin{pmatrix}0&&\frac{-s^{3}}{\kappa\,(s-1)^{3}}&&0\\\
(\frac{1}{s}-1)\,\lambda&&\frac{s^{2}}{\kappa\,(1-s)}&&0\\\
0&&\frac{1}{1-s}&&\frac{\lambda+1}{s-1}\end{pmatrix}\quad\text{and}\quad
B^{+}(\lambda)=\begin{pmatrix}0&&\frac{-1}{\kappa}&&0\\\
-\lambda&&-\frac{s}{\kappa}&&0\\\
0&&-\frac{1}{s}&&\frac{\lambda+1}{s}\end{pmatrix},$ (22)

where $s=s_{2}$. As a consequence the essential spectra as well as the
absolute spectra of S2 and S1 coincide.

Note that the spatial eigenvalues of $B^{-}$ are given by those of $A^{+}$
given in (14) and those of $B^{+}$ are given by those of $A^{-}$ listed in
(13), again where $s=s_{2}(=-s_{1})$.

To compute the Evans function, we start by identifying the stable and unstable
eigenvectors of the asymptotic matrices $B^{\pm}(\lambda)$ respectively in an
attempt to construct an integrable function which plays the role of an
eigenvector for the linearised system. The spatial eigenvalues (13) and (14)
are also spatial eigenvalues of $B^{\pm}$. Since here $s=s_{2}>1$ to avoid
violations of the impenetrability constraint (see section C), the unstable
eigenvalue of $B^{-}$ and the stable eigenvalues of $B^{+}$ are given by

$\mu_{1}^{-}(\lambda)=\frac{\lambda+1}{s-1}\;,\quad\mu_{2}^{-}(\lambda)=\frac{s^{2}-s\sqrt{s^{2}+4\kappa\lambda}}{2\kappa(1-s)}\quad\text{and}\quad\mu_{3}^{+}(\lambda)=\frac{-s-\sqrt{s^{2}+4\kappa\lambda}}{2\kappa}\;,$
(23)

which satisfy $\mu_{3}^{+}(\lambda)<0$ and
$\mu_{1}^{-}(\lambda)\,,\;\mu_{2}^{-}(\lambda)>0$ for ${\rm Re}\,\lambda\gg
1$.

We introduce the eigenvectors $X_{0}^{+}(\lambda)$ associated to the
eigenvalue $\mu_{2}^{+}(\lambda)$ of $A^{+}$ as well as $X_{0}^{-}(\lambda)$
and $Y_{0}^{-}(\lambda)$ being eigenvectors of $A^{-}$ associated to
$\mu_{1}^{-}(\lambda)$ and $\mu_{3}^{-}(\lambda)$.

To construct an eigenfunction of the linearised system we solve (6) using
these vectors as initial, respectively terminal condition until $z=0$. Then
the Evans function is defined as the Wronskian

$D(\lambda)={\rm
det}\;[X^{+}(z=0,\lambda),X^{-}(z=0,\lambda),Y^{-}(z=0,\lambda)]\;,$ (24)

which vanishes if the stable and unstable eigenvectors propagated to $z=0$ are
linearly dependent and can be combined into a smooth eigenfunction. We start
by computing $X^{-}(z=0,\lambda)$ and $Y^{-}(z=0,\lambda)$ which can be in
terms of a closed-form expression. Only for the computation of
$X^{+}(z=0,\lambda)$ we will resort to numerical results.

For $z<0$ it holds that $R\equiv s_{2}/(s_{2}-1)$ and the matrix (6) is given
by (18). Again, due to the zeros in the third column the equations for
$\delta\rho$ and $\delta v$ are not coupled to $\delta a$. They satisfy (19)
which is a system of two linear, constant-coefficient equations. It admits one
fundamental solution with positive eigenvalue (we omit the stable fundamental
solution) given by

$\begin{pmatrix}\delta\rho\\\ \delta
v\end{pmatrix}=C_{2}\,\exp(\mu_{2}^{-}(\lambda)\,z)\begin{pmatrix}-\,\frac{s^{2}\,(s+\sqrt{s^{2}+4\,\kappa\,\lambda})}{(s-1)^{2}\,\kappa}\\\
2\,\lambda\end{pmatrix}$ (25)

With this information we can rewrite the equation for $\delta a$ which is
contained in the third row of (18),

$\displaystyle\delta a^{\prime}$
$\displaystyle=\frac{A^{\prime}-1}{s-1}\,\delta
v+\frac{\lambda+1}{s-1}\,\delta a$
$\displaystyle=\frac{A^{\prime}-1}{s-1}\,2\,\lambda\,C_{2}\,e^{\mu_{2}^{-}\,z}+\mu_{1}^{-}\,\delta
a\;.$

The general solution is given by

$\delta
a=C_{1}e^{\mu_{1}^{-}z}+\frac{2\lambda}{s-1}C_{2}\int_{0}^{z}e^{(\mu_{2}^{-}-\mu_{1}^{-})z}\left(\frac{\alpha-1}{s-1}e^{\frac{z}{s-1}}-1\right)\,d\tilde{z}\;.$

This implies that

$\lim_{z\to 0^{-}}\begin{pmatrix}\delta\rho\\\ \delta v\\\ \delta
a\end{pmatrix}=\begin{pmatrix}C_{2}\,\left(\,-\,\frac{s^{2}\,(s+\sqrt{s^{2}+4\,\kappa\,\lambda})}{(s-1)^{2}\,\kappa}\right)\\\
C_{2}\,2\,\lambda\\\ C_{1}\end{pmatrix}\;.$

While most components of $A(z,\lambda)$ in (6) are functions, the derivative
of the active speed of migration (1) is a $\delta$-distribution,
$M^{\prime}(a)=\delta(a-\alpha)$. The change of variables between $z$ and the
monotonically decreasing function $a=A(z)=A_{2}(z)$ shows that when
integrating with respect to $z$ the following expression is a
$\delta$-distribution centred at $z=0$,
$-M^{\prime}(A(z))A^{\prime}(z)=\delta(z)$. Therefore $\delta\rho$ and $\delta
v$ in the solution of (6) undergo a jump at $z=0$ which involves the factor
$-1/A^{\prime}(0)$. The right limit as $z\rightarrow 0$ of the solution vector
is given by

$\displaystyle\lim_{z\to 0^{+}}\begin{pmatrix}\delta\rho\\\ \delta v\\\ \delta
a\end{pmatrix}=$
$\displaystyle\begin{pmatrix}C_{2}\,\left(\,-\frac{s^{2}\,(s+\sqrt{s^{2}+4\,\kappa\,\lambda})}{(s-1)^{2}\,\kappa}\right)\\\
C_{2}\,2\,\lambda\\\
C_{1}\end{pmatrix}+\frac{(-1)C_{1}}{A^{\prime}(0)}\frac{1}{\kappa}\begin{pmatrix}R^{3}(0)\\\
R(0)\,s\\\ 0\end{pmatrix}$ $\displaystyle=$ $\displaystyle
C_{2}\,\begin{pmatrix}-\,\frac{s^{2}\,(s+\sqrt{s^{2}+4\,\kappa\,\lambda})}{(s-1)^{2}\,\kappa}\\\
2\,\lambda\\\
0\end{pmatrix}+\,C_{1}\,\begin{pmatrix}-\frac{s^{3}}{(s-1)^{2}\,(\alpha-1)\,\kappa}\\\
-\frac{s^{2}}{(\alpha-1)\,\kappa}\\\ 1\end{pmatrix}.$

where we take the values of $R(0)$ and $A^{\prime}(0)$ from the travelling
wave profile S2.

We obtain the closed-form solutions $X^{-}(0,\lambda)$ and $Y^{-}(0,\lambda)$
corresponding to the stable eigenvalues $\mu_{1}^{-}(\lambda)$ and
$\mu_{2}^{-}(\lambda)$,

$\displaystyle
X^{-}(0,\lambda)=\begin{pmatrix}\frac{s^{3}}{(s-1)^{2}\,(\alpha-1)\,k}\\\
\frac{s^{2}}{(\alpha-1)\,k}\\\ 1\end{pmatrix}\quad\text{and}\quad
Y^{-}(0,\lambda)=\begin{pmatrix}-\,\frac{s^{2}\,(s+\sqrt{s^{2}+4\,k\,\lambda})}{(s-1)^{2}\,k}\\\
2\,\lambda\\\ 0\end{pmatrix}.$ (26)

Finally, in order to evaluate the Evans function (17) for a given
$\lambda\in\mathbb{C}$, we compute $X^{+}(0,\lambda)$ propagating the
eigenvector associated to the spatial eigenvalue $\mu_{3}^{+}$ given by

$v_{3}=\begin{pmatrix}\frac{-s^{2}+\sqrt{s^{4}+4s^{2}\kappa\lambda}}{2\lambda
s\kappa}\\\ 1\\\
\frac{2\kappa}{s^{2}+2\kappa(1+\lambda)+\sqrt{s^{4}+4s^{2}\kappa\lambda}}\end{pmatrix}\;\text{for}\;\lambda\neq
0\quad\text{and}\quad v_{3}=\begin{pmatrix}\frac{1}{s}\\\ 1\\\
\frac{\kappa}{s^{2}+\kappa}\end{pmatrix}\;\text{for}\;\lambda=0$

from an arbitrary small value of $z$ (we choose $z=20$) until $z=0$ solving
the system (6) numerically.

Figure 15: $D(C_{1})$ for $\alpha=0.2,\,0.4,\,0.5\quad\text{and}\quad 0.6$.
Figure 16: $D(C_{2})$ for $\alpha=0.2,\,0.4,\,0.5\quad\text{and}\quad 0.6$.

Finally we compute the image of both contours $C_{1}$ and $C_{2}$ (Fig.10)
under the Evans function defined by (24). For various parameter values
$\alpha$ and $\kappa$ we find that for $C_{1}$ the winding number around the
origin is $1$ (Fig. 15) and for the contour $C_{2}$ which does not enclose the
origin, the winding number is 0 (Fig. 16). As for S1 this indicates that the
point spectrum on the right of the complex plane only consists of $\lambda=0$
which is expected for a travelling wave solution. This suggests that the
travelling wave solutions S2 and – through to transformation $T_{2}$ – S4 are
linearly stable.

##  Appendix C Unphysical travelling wave solutions

Note that the travelling wave solution S2 can only be realised if the model
parameters $\kappa$ and $\alpha$ are such that $s_{2}>1$. If that is not the
case the mathematical solution is not physical and violates the
impenetrability of single cells as illustrated in Fig. 17.

Figure 17: Un-physical travelling wave solutions in violation of the cells’
impenetrablity. (A) shows a simulation of the polarisation wave S2 for a set
of parameters for which $s_{2}<1$. (B) shows a simulation of the polarisation
wave S4 for a set of parameters for which $s_{4}>0$.

## References

  * [1] D. Montell, “Morphogenetic cell movements: Diversity from modular mechanical properties,” _Science_ , vol. 322, no. 5907, pp. 1502–1505, 2008.
  * [2] P. Friedl and D. Gilmour, “Collective cell migration in morphogenesis, regeneration and cancer,” _Nat. Rev. Mol. Cell Biol._ , vol. 10, pp. 445–457, 2009.
  * [3] R. McMinn and F. Johnson, “Mitosis in migrating epithelial cells,” _Nature_ , vol. 178, p. 212, 1956.
  * [4] N. Koshikawa, G. Giannelli, V. Cirulli, K. Miyazaki, and V. Quaranta, “Role of cell surface metalloprotease mt1-mmp in epithelial cell migration over laminin-5.” _The Journal of cell biology_ , vol. 148, no. 3, pp. 615–624, 2000.
  * [5] M. Vishwakarma, J. D. Russo, D. Probst, U. S. Schwarz, T. Das, and J. P. Spatz, “Mechanical interactions among followers determine the emergence of leaders in migrating epithelial cell collectives,” _Nature Communications_ , vol. 9, no. 1, pp. 1–12, 2018.
  * [6] M. L. Zorn, A.-K. Marel, F. J. Segerer, and J. O. Rädler, “Phenomenological approaches to collective behavior in epithelial cell migration,” _BBA - Molecular Cell Research_ , vol. 1853, no. 11, pp. 3143–3152, 2015.
  * [7] B. Dalton and J. Steele, “Migration mechanisms: Corneal epithelial tissue and dissociated cells,” _Experimental Eye Research_ , vol. 73(6), pp. 797–814, 2001.
  * [8] K. Ebnet, “Organization of multiprotein complexes at cell–cell junctions,” _Histochemistry and Cell Biology_ , vol. 130(1), pp. 1–20, 2008.
  * [9] M. Cavey and T. Lecuit, “Molecular bases of cell–cell junctions stability and dynamics,” _Cold Spring Harbor perspectives in biology_ , vol. 1(5), pp. 1–20, 2009.
  * [10] A. S. Yap, M. S. Crampton, and J. Hardin, “Making and breaking contacts: the cellular biology of cadherin regulation,” _Current Opinion in Cell Biology_ , vol. 19, no. 5, pp. 508–514, 2007.
  * [11] T. Lecuit and A. S. Yap, “E-cadherin junctions as active mechanical integrators in tissue dynamics,” _Nature Cell Biology_ , vol. 17(5), pp. 533–539, 2015.
  * [12] Y. Mori, A. Jilkine, and L. Edelstein-Keshet, “Wave-pinning and cell polarity from a bistable reaction-diffusion system,” _Biophysical Journal_ , vol. 94, no. 9, pp. 3684–3697, 2008.
  * [13] R. s. Gray, I. Roszko, and L. Solnica-Krezel, “Planar cell polarity: Coordinating morphogenetic cell behaviors with embryonic polarity,” _Developmental Cell_ , vol. 21, no. 1, pp. 120–133, 2011.
  * [14] F. Bosveld, I. Bonnet, B. Guirao, S. Tlili, Z. Wang, A. Petitalot, R. Marchand, P.-L. Bardet, P. Marcq, F. Graner, and Y. Bellaiche, “Mechanical control of morphogenesis by fat/dachsous/four-jointed planar cell polarity,” _Science_ , vol. 336, no. 6082, p. 724, 2012\.
  * [15] Y. Burak and B. I. Shraiman, “Order and stochastic dynamics in drosophila planar cell polarity,” _PLoS Computational Biology_ , vol. 5(12), pp. 1–10, 2009.
  * [16] A. Asnacios and O. Hamant, “The mechanics behind cell polarity,” _Trends in Cell Biology_ , vol. 22(11), pp. 584–591, 2012.
  * [17] D. Oelz, H. Khataee, A. Czirok, and Z. Neufeld, “Polarization wave at the onset of collective cell migration,” _Phys. Rev. E_ , vol. 100, p. 032403, Sep 2019.
  * [18] X. Serra-Picamal, V. Conte, R. Vincent, E. Anon, D. T. Tambe, E. Bazellieres, J. P. Butler, J. J. Fredberg, and X. Trepat, “Mechanical waves during tissue expansion,” _Nature Physics_ , vol. 8, no. 8, p. 628, 2012.
  * [19] Y. Zhang, G. Xu, R. Lee, Z. Zhu, J. Wu, S. Liao, G. Zhang, Y. Sun, A. Mogilner, W. Losert, T. Pan, F. Lin, Z. Xu, and M. Zhao, “Collective cell migration has distinct directionality and speed dynamics,” _Cellular and Molecular Life Sciences_ , vol. 74, no. 20, pp. 3841–3850, 2017.
  * [20] B. Sandstede, “Chapter 18 - stability of travelling waves,” in _Handbook of Dynamical Systems_ , ser. Handbook of Dynamical Systems, B. Fiedler, Ed. Elsevier Science, 2002, vol. 2, pp. 983 – 1055.
  * [21] K. Zumbrun and P. Howard, “Pointwise semigroup methods and stability of viscous shock waves,” _Indiana University Mathematics Journal_ , vol. 47(3), pp. 741–871, 1998.
  * [22] C. Mascia and K. Zumbrun, “Pointwise green’s function bounds and stability of relaxation shocks,” _Indiana University Mathematics Journal_ , vol. 51(4), pp. 773–904, 2002.
  * [23] J. Thomas, _Numerical Partial Differential Equations: Finite Difference Methods by J.W. Thomas._ , 1st ed., ser. Texts in Applied Mathematics, 22. New York, NY: Springer New York : Imprint: Springer, 1995.
  * [24] T. Kapitula and K. Promislow, _Spectral and Dynamical Stability of Nonlinear Waves_. New York: Springer, 2013\.
  * [25] M. Chan, P. Kim, and R. Marangell, “Stability of travelling waves in a wolbachia invasion,” _Discrete and Continuous Dynamical Systems - B_ , vol. 23, p. 609, 2018.
  * [26] B. Sandstede and A. Scheel, “Absolute and convective instabilities of waves on unbounded and large bounded domains,” _Physica D: Nonlinear Phenomena_ , vol. 145, no. 3, pp. 233–277, 2000.
  * [27] P. N. Davis, P. van Heijster, and R. Marangell, “Absolute instabilities of travelling wave solutions in a Keller–Segel model,” _Nonlinearity_ , vol. 30, no. 11, p. 4029, 2017.
  * [28] J. Evans, “Nerve axon equations: II stability at rest,” _Indiana University Mathematics Journal_ , vol. 22, no. 1, pp. 75–90, 1972.
  * [29] ——, “Nerve axon equations: III stability of the nerve impulse,” _Indiana University Mathematics Journal_ , vol. 22(6), pp. 577–593, 1972\.