# Structural identifiability analysis of linear reaction-advection-diffusion
processes in mathematical biology

Alexander P Browning  Mathematical Institute, University of Oxford, Oxford,
United Kingdom  Maria Tască  Mathematical Institute, University of Oxford,
Oxford, United Kingdom  Somerville College, University of Oxford, Oxford,
United Kingdom  Carles Falcó  Mathematical Institute, University of Oxford,
Oxford, United Kingdom  Ruth E Baker  Mathematical Institute, University of
Oxford, Oxford, United Kingdom

###### Abstract

Effective application of mathematical models to interpret biological data and
make accurate predictions often requires that model parameters are
identifiable. Approaches to assess the so-called structural identifiability of
models are well-established for ordinary differential equation models, yet
there are no commonly adopted approaches that can be applied to assess the
structural identifiability of the partial differential equation (PDE) models
that are requisite to capture spatial features inherent to many phenomena. The
differential algebra approach to structural identifiability has recently been
demonstrated to be applicable to several specific PDE models. In this brief
article, we present general methodology for performing structural
identifiability analysis on partially observed linear reaction-advection-
diffusion (RAD) PDE models. We show that the differential algebra approach can
always, in theory, be applied to linear RAD models. Moreover, despite the
perceived complexity introduced by the addition of advection and diffusion
terms, identifiability of spatial analogues of non-spatial models cannot
decrease structural identifiability. Finally, we show that our approach can
also be applied to a class of non-linear PDE models that are linear in the
unobserved variables, and conclude by discussing future possibilities and
computational cost of performing structural identifiability analysis on more
general PDE models in mathematical biology.

###### Keywords

structural identifiability, reaction-diffusion, advection, cell migration,
chemotaxis, partial differential equations

##  1 Introduction

Mathematical models play an irreplacable role in the interpretation of
biological data. Model parameters are now routinely used to objectively
quantify observed behaviour and characterise behaviours that cannot be
directly measured [1, 2]. Whether it is possible, in theory, for model
parameters to be uniquely identified given a particular mathematical model,
and a particular set of observed quantities, is referred to as structural
identifiability [3, 4, 5, 6, 7]. Assessing the structural identifiability of a
model and observation process can provide vital insights that guide data
collection before experiments have been conducted [8]; establish identifiable
parameter combinations that resolve issues of non-identifiability either
through reparameterisation or additional experimentation [9]; and provide
confidence in predictions drawn from calibrated mathematical models [10].

Tools for assessing structural identifiability are long-established for
deterministic, ordinary differential equation (ODE) models [11, 12, 13, 14].
The advent of open-source and even web-based software [15] to automate
otherwise the tedious analysis has ingrained questions related to structural
identifiability into the inference process for ODE models [14]. Methods
originate with differential-algebra-based approaches [16], which for linear or
polynomial systems assess identifiability through a so-called input-output
relation: a set of monic polynomials in the derivatives of observed variables,
the coefficients of which form the set of identifiable parameter combinations.
Differential algebra approaches are trivial if all variables in a system are
observed, although questions of structural identifiability very often relate
to partially observed systems, where only measurements of a subset of state
variables are obtainable [3, 8, 17]. Here, the algorithmic complexity of
differential algebra approaches comes from the reduction of a high-dimensional
(in the number of state variables), lower-derivative-order, system to a set of
higher-derivative-order polynomials that include only the observed quantities
[11]. Several other approaches for structural identifiability have been since
been established; notably those based on Taylor series, generating series and
Lie derivatives [7, 18], and similarity transforms [19, 7]. Many of these more
modern approaches are more broadly applicable to analytic systems, and
significantly more computationally efficient than those based on differential-
algebra.

However, many forms of biological data are inherently spatial, and therefore
not well-described by ODE models [20]; data relating to cell migration [21,
22], or diffusive processes [23, 24], for example. Yet, tools for assessing
the structural identifiability of the partial differential equation (PDE)
models that capture spatial heterogeneity remain almost entirely undeveloped.
Recent work has demonstrated the application of the differential algebra
method to assess structural identifiability of a class of first-order age
structured PDE models [25], and a system comprising a single diffusive species
[26]. However, structural identifiability of more generic systems that contain
both first and second order spatial derivatives (hereafter referred to as
reaction-advection-diffusion (RAD) systems) have not, to the best of our
knowledge, been explored. In particular, it has remained unclear whether the
differential algebra approach can always be applied to PDE models, or whether
such analysis of PDE models is more restrictive. Even more basic questions,
including whether the apparent additional complexity of including a diffusion
term exacerbates or alleviates issues relating to structural identifiability,
remain unanswered. Given the prevalence of spatial data in biology, brought in
no small part due to advancements in microscopy and imaging technologies, it
is imperative that tools for structural identifiability are developed for the
PDE models that account for inherently spatial behaviours.

In this brief article, we present a methodology and a series of general
results for the structural identifiability linear and semi-linear RAD PDE
models subject to observation processes that do not include all state
variables (i.e., models that are partially). For a more thorough and formal
introduction to structural identifiability and the differential algebra
approach, we direct the reader to [27, 25]. We present our results through
example, first presenting a didactic guide to the procedure for a two-state
model of cell proliferation and migration inspired by common scratch assay
experiments fig. 1 [21]. Next, in Section 2.2, we present more general
methodology and general results for a generic two-state RAD model; these
results are then compared to the corresponding ODE model before the role of
the initial condition—which differs significantly to the ODE case—is discussed
in Section 2.3. In Section 2.4 we formally present results for generic linear
systems comprising an arbitrary number of states. Finally, in Section 2.5 we
present an example set of results for two semi-linear systems: a model of
bacterial chemotaxis, and a three-state nonlinear analogue of the scratch
assay model. Avenues for future work and key features that distinguish the
structural identifiability problem for ODE models from those of PDE models are
discussed in Section 3.

Figure 1: Partially observed spatial data in a scratch assay experiment. (a)
Scratch assays involve growing a monolayer of cells (purple) in a well before
making an artificial scratch or wound (white) and imaging a central region.
(b) Snapshots collected at $0\text{\,}\mathrm{h}$ and $48\text{\,}\mathrm{h}$
from a typical scratch assay experiment conducted using PC3 prostate cancer
cells [28]. Averaging in the direction of the scratch provides one-dimensional
spatio-temporal information relating to the total cell density. (c) Scratch
assay experiment with cell cycle information conducted using WM983C melanoma
cells [29]. Encoding cells with FUCCI technology allows additional
visualisation of the cell cycle [30], shown schematically in (d). FUCCI
scratch assay experiments provide spatio-temporal information relating to the
density of cells in each stage of the cell cycle. Panels in (b) and (c) are
reprinted from [28] and [29], respectively, under a CC-BY and CC-attribution
licence.

##  2 Methods and Results

###  2.1 Linear cell-cycle model

To demonstrate the differential algebra approach for linear PDE models, we
first consider a two-state cell cycle model of cell migration subject to
exponential growth [31, 21]. We divide the cell cycle (fig. 1d) into two
subspecies that approximately correspond to the fluorescent markers in a FUCCI
[30] scratch assay (fig. 1c): cells in G1 phase fluoresce red, and cells in
S/G2/M phase fluoresce green [21]. During M-phase (mitosis), a green cell
proliferates into two daughter cells that eventually fluoresce red. Finally,
we model cell motility using linear diffusion. The model is given by

$\displaystyle\partialderivative{r}{t}$
$\displaystyle=D_{1}\partialderivative[2]{r}{x}-k_{1}r+2k_{2}g$ (1a)
$\displaystyle\partialderivative{g}{t}$
$\displaystyle=D_{2}\partialderivative[2]{g}{x}+k_{1}r-k_{2}g,$ (1b)

where $r(x,t)$ and $g(x,t)$ denote the density of cells in G1 (red) and G2
(green), respectively; $D_{1}$ and $D_{2}$ are the diffusion coefficients
associated with each subspecies; $k_{1}$ is the G1 (red) to G2 (green)
transition rate; $k_{2}$ the mitotic rate of green cells. For compactness,
from this point on we denote derivatives using superscripts:
$r^{(i,j)}=\partial^{i+j}r/\partial x^{i}\partial t^{j}$. We prescribe
homogeneous Neumann (no-flux) boundary conditions, such that no additional
information is available from the equations that govern the behaviour at the
boundaries.

While the model is motivated by FUCCI experiments that provide information on
the cell-cycle status of individual cells (fig. 1b), more typical experimental
data is of the form

$n(x,t)=r(x,t)+g(x,t),$ (2)

as shown in fig. 1c [32]. A key question of interest is: can the model
parameters, including the transition rates between cell-cycle states that are
not individually observed, be inferred from an observation process of the form
given by eq. 2?

Application of the differential algebra approach to assess this question of
structural identifiability requires first that we write the system eq. 1 as a
differential-algebraic equation in terms only of the variable that we observe,
$n$, and its derivatives. To do so, we first write

$n^{(0,1)}=D_{1}n^{(2,0)}+(D_{2}-D_{1})g^{(2,0)}+k_{2}g,$ (3)

to eliminate $r$. At this point, it is clear that the analysis is far simpler
in the case that $D_{2}=D_{1}=D$: all that remains in this case is to solve
eq. 3 for $g$ in terms of $n^{(0,1)}$ and $n^{(2,0)}$, and substitute into eq.
1b. For distinct diffusivities, progress is made by solving eqs. 3 and 1b
simultaneously for $g^{(0,1)}$ and $g^{(2,0)}$ to obtain

$\displaystyle g^{(0,1)}$
$\displaystyle=k_{1}n+\dfrac{\big{[}D_{2}(k_{1}+2k_{2})-D_{1}(k_{1}+k_{2})\big{]}g-D_{2}n^{(0,1)}+D_{1}D_{2}n^{(2,0)}}{D_{1}-D_{2}},$
(4a) $\displaystyle g^{(2,0)}$
$\displaystyle=\dfrac{k_{2}g-n^{(0,1)}+D_{1}n^{(2,0)}}{D_{1}-D_{2}}.$ (4b)

This step is advantageous as yields a first order equation, giving greater
flexibility in subsequent equations that can be produced up to a specified
order. Importantly, we note that this and following steps can be completed
without division but multiplication and elimination to derive the requisite
first-order equation.

Next, we differentiate eqs. 4a and 4b with respect to $x$ twice and $t$,
respectively, to obtain two expressions for $g^{(2,1)}$ that can be equated to
give, after some simplification,

$0=-k_{1}k_{2}n+(k_{1}+k_{2})n^{(0,1)}+n^{(0,2)}-(D_{2}k_{1}+D_{1}k_{2})n^{(2,0)}-(D_{1}+D_{2})n^{(2,1)}+D_{1}D_{2}n^{(4,0)}.$
(5)

At this point, many authors divide through by the coefficient of the “highest-
order” term, following some ordering, to ensure that the resultant polynomial
system is monic, and therefore unique. We take a simpler approach and divide
the resultant polynomial through by any coefficient: this still ensures
uniqueness, however simplies calculations. For eq. 5, we trivially divide
through by the coefficient of $n^{(0,2)}$. Thus, the following set of
polynomial coefficients (multiplied through by $-1$ where appropriate), are
structurally identifiable

$\Big{\\{}k_{1}k_{2},k_{1}+k_{2},D_{2}k_{1}+D_{1}k_{2},D_{1}+D_{2},D_{1}D_{2}\Big{\\}}.$
(6)

Clearly, since both the product $k_{1}k_{2}$ and sum $k_{1}+k_{2}$ are
structurally identifiable, so too are the individual parameters $k_{1}$ and
$k_{2}$. A similar observation can be made for the diffusivities, and so all
model parameters are structurally identifiable. As we expect, removing spatial
information (for example, by initialising the system with a spatially
homogeneous initial condition such that $n^{(i,0)}=0$ for all $i=1,2,\ldots$),
the polynomial system reduces to terms that do not include $D_{1}$ nor
$D_{2}$, however the rate constants $k_{1}$ and $k_{2}$ remain identifiable.

###  2.2 Generic two-state reaction-advection-diffusion system

Having established the general procedure for applying the differential algebra
framework to reaction-diffusion models, we now consider the general two-state
model

$\displaystyle u^{(0,1)}$
$\displaystyle=D_{u}u^{(2,0)}+\alpha_{u}u^{(1,0)}+p_{1}u+p_{2}v+p_{3},$ (7a)
$\displaystyle v^{(0,1)}$
$\displaystyle=D_{v}v^{(2,0)}+\alpha_{v}v^{(1,0)}+p_{4}u+p_{5}v+p_{6},$ (7b)

subject to observations of the form $n(x,t)=u(x,t)+v(x,t)$.

Analysis of such a system captures both extensions to the linear cell-cycle
model that incorporates advection, and canonical models analysed in the
structural identifiability literature such as the so-called “two-pool model”
[33]. The following results are trivially applicable to other linear
combinations of the states (for example, observations of the form $\alpha
u(x,t)+\beta v(x,t)$ where $\alpha$ and $\beta$ are unknown parameters)
through a rescaling.

The procedure for determining the required polynomial equation is now
complicated by the presence of first order spatial derivatives $v^{(1,0)}$ in
the system. To proceed, we first eliminate $u$ to obtain a two-state system in
$n$ and $v$, with

$\displaystyle n^{(0,1)}$
$\displaystyle=D_{u}n^{(2,0)}+(D_{v}-D_{u})v^{(2,0)}+\alpha_{u}n^{(1,0)}+(\alpha_{v}-\alpha_{u})v^{(1,0)}$
(8)
$\displaystyle\qquad+(p_{1}+p_{4})n+(p_{2}+p_{5}-p_{1}-p_{4})v+p_{3}+p_{6}.$

Solving eqs. 7b and 8 for $\left(v^{(2,0)},v^{(0,1)}\right)$ so as to obtain a
system of one first order and one second order equation. We can then expand
the resultant system by taking appropriate derivatives up to $i$th order so as
to end up with a linear system with more equations than unknown variables. In
this case, this is achieved by considering derivatives up to order $i=4$,
which yields 16 equations in all 15 possible fourth order derivatives of the
unobserved variable $v$. The resultant system will always be linear for all
linear and what we term “semi-linear” systems, the latter defined as systems
that are linear in the variables we wish to eliminate (a logistic term of the
form $v(1-n)$ would be semi-linear, since $n$ is observed).

Thus, we arrive at the over-determined system

$\mathbf{A}\mathbf{v}^{(4)}=\mathbf{b},$ (9)

where $\mathbf{v}^{(i)}\in\mathbb{R}^{(i+1)(i+2)/2}$ denotes a vector of all
possible $i$th order derivatives. We proceed by performing Gaussian
elimination to reduce the augmented matrix $(\mathbf{A}|\mathbf{b})$ into row
echelon form. This provides a set of expressions, based upon linear
combinations of the elements of $\mathbf{b}$, that necessarily vanish; thus
providing the required set of polynomial equations.

It is not necessarily the case that we must work with the full set of fourth
order equations, just that we are required to do so to ensure that the system
is appropriately determined. This is important as it is possible for the
system to be characterised by more than one polynomial equation: indeed, the
number of polynomial equations we expect to require for linear and semi-linear
problems will be given by the difference in rank of $\mathbf{A}$ and the
augmented matrix $(\mathbf{A}|\mathbf{b})$ from the over-determined system. In
our case, we expect a single polynomial equation. It also happens in our case
that we can obtain this polynomial equation by working with a subset of third
order equations. Algebraic manipulations are performed in Mathematica and
available as supplementary code.

For the system defined by eq. 7, we arrive at the set of polynomial
coefficients

$\displaystyle\Big{\\{}p_{1}p_{5}-$ $\displaystyle
p_{2}p_{4},p_{1}+p_{5},p_{5}\alpha_{u}+p_{1}\alpha_{v},\alpha_{u}+\alpha_{v},D_{v}p_{1}+D_{u}p_{5}+\alpha_{u}\alpha_{v},$
(10) $\displaystyle-
D_{u}-D_{v},D_{v}\alpha_{u}+D_{u}\alpha_{v},D_{u}D_{v},(p_{5}-p_{4})p_{3}+(p_{1}-p_{2})p_{6}\Big{\\}},$

thus identifiable parameter combinations are given by

$\displaystyle\Big{\\{}D_{u},D_{v},\alpha_{u},\alpha_{v},p_{1},p_{5},p_{2}p_{4},(p_{5}-p_{4})p_{3}+(p_{1}-p_{2})p_{6}\Big{\\}}.$
(11)

By setting $\alpha_{u}=\alpha_{v}=0$, we see that these results are identical
to the no-advection model. However, results differ slightly for the non-
spatial model. Setting $D_{u}=D_{v}=\alpha_{u}=\alpha_{v}=0$ to remove spatial
derivatives, we see that the set of identifiable parameter combinations is now

$\displaystyle\Big{\\{}p_{1}p_{5}-p_{2}p_{4},p_{1}+p_{5},(p_{5}-p_{4})p_{3}+(p_{1}-p_{2})p_{6}\Big{\\}},$
(12)

so we can no longer identify any rate parameters individually. This results
demonstrates that we are potentially able to learn more about the process
through spatial heterogeneity in the system.

###  2.3 Role of the initial condition in identifiability

Clearly, if we observe only $n(x,t)=u(x,t)+v(x,t)$ we cannot always fully
ascertain the initial condition $v(x,0)$. A key point of distriction between
such partially observed PDE models and their ODE counterparts is that the
initial condition is an unknown function, denoted here by $v(x,0):=v_{0}(x)$,
whereas for an ODE model the initial condition is merely an additional unknown
scalar parameter. In the case of structural identifiability, we assume that we
perfectly observe not only $n(x,0):=n_{0}(x)$ but also its derivatives. Thus,
we expect to observe right-hand-side of eq. 8 such that the functional

$\displaystyle f(x):=n_{0}^{(0,1)}(x)$
$\displaystyle=D_{u}n_{0}^{(2,0)}(x)+(D_{v}-D_{u})v_{0}^{(2,0)}(x)+\alpha_{u}n_{0}^{(1,0)}(x)+(\alpha_{v}-\alpha_{u})v_{0}^{(1,0)}(x)$
(13)
$\displaystyle\qquad+(p_{1}+p_{4})n(x)+(p_{2}+p_{5}-p_{1}-p_{4})v_{0}(x)+p_{3}+p_{6},$

is also “identifiable” in a similar sense to the other identifiable parameter
combinations in eq. 11. This result shows a clear interdependence between
model parameters and the initial condition: specifiying a functional form for
the initial condition has a similar role to fixing a non-identifiable
parameter that appears in combination with other model parameters. Thus, it is
possible for a practitioner to unknowingly render structurally non-
identifiable parameters identifiable through an assumption related to the
initial condition, highlighting the importance of considering structural
identifiability when working with partially observed PDE models.

Figure 2: Non-identifiability of the linear RAD model. (a) Individual
(unobserved) states, $u(x,t)$ and $v(x,t)$ at $t=0$, 2, 4 and 6. The arrow
indicates the direction of increasing $t$. Shown are curves for the original
set of parameter values (green and blue solid curves, top row) and curves for
the modified set of parameter values where $p_{2}\mapsto 2p_{2}$ (red dashed,
bottom row). Note the difference in initial condition between parameter sets.
(b) Observed variable, $n(x,t)=u(x,t)+v(x,t)$ at the original set of parameter
values (grey) and modified set of parameter values (red dashed). Note that the
parameter sets are indistinguishable. (c) Equivalent dynamics for the ODE
model, showing $n(t)$ and $v(t)$ at the original parameter set (solid) and
modified parameter set (dashed). Original parameters are $D_{u}=20$,
$D_{v}=10$, $\alpha_{1}=-5$, $\alpha_{2}=5$, $p_{i}=0.1$ $(i\neq 2)$ and
$p_{2}=0.2$.

We illustrate these results in fig. 2 by noting that eq. 11 gives combinations
of parameters that can produce identical model outputs. For example, we can
set $p_{2}\mapsto 2p_{2}$ and see no change to $n(x,t)$ provided other
parameters and the initial condition are adjusted accordingly. Solving the
model with an initial parameter set and initial condition for $n(x,0)$
provides $f(x)$ to determine the new initial condition for $v(x,t)$ in terms
of eq. 13, a steady-state reaction-advection equation (fig. 2a). Solving the
full system again (fig. 2b) demonstrates structural non-identifiability: model
outputs are indistinguishable. We show the equivalent set of results for the
corresponding ODE model in fig. 2c.

###  2.4 Generic linear systems

A question that naturally arises out of our analysis of the generic two-state
system is whether the procedure is applicable to linear systems with an
arbitrary number of states and outputs; a second is whether the reduction in
the number of identifiable parameter combinations when reducing to a
corresponding ODE model is generally true of RAD models.

To determine whether it is, in theory, possible to find a polynomial system
for all $m$-state linear RAD models subject to observations of $\ell$
variables, we need only consider whether it is possible to close such a system
through repeated differentiation. That is, is there an order of derivative at
which we are guaranteed that we will have more determining equations than
variables? We arrive at the following result.

######  Theorem 1.

All linear RAD models of $m$ states can be reduced to a set of polynomial
relations involving derivatives of order no more than $4(m-1)$.

###### Proof.

We can rewrite the system of $m$ equations to include at least one first-order
equation. Expanding the system to include derivatives up to order $n$ yields
at least $q(n)=(mn^{2}-mn+2n)/2$ equations, with no more than
$v(n)=(n+1)(n+2)(m-1)/2$ unknown variables (i.e., $n$th order and lower
derivatives of the unobserved quantities). The number of excess variables is
no more than $d(n)=q(n)-v(n)$, which we require to be positive. Defining
$n^{*}=\min n:d(n)\geq 0,n\in\mathbb{Z}$, it is straightforward to see that
$n^{*}=4(m-1)$, as required. ∎

The result of Theorem 1 highlights the computational cost of performing
structural identifiability for PDE models: application of our method for a
linear system of only three states would potentially require us to perform
Gaussian elimination symbolically on a 92 equation system. Importantly,
however, the result does demonstrate that it is always possible to perform
structural identifiability analysis on linear RAD models. Furthermore, the
Gaussian-elimination-based reduction procedure requires only that the system
be linear in the unobserved quantities, thus we also expect Theorem 1 to hold
for all semi-linear systems, including, for example, analogues of eq. 1 with
logistic growth terms and observation processes that include the total
population.

Secondly, we observe for the linear two-state system that the PDE model
contains a greater number of identifiable parameter combinations than the
corresponding ODE (i.e., spatially homogeneous) model. We present the
following result.

######  Corollary 1.

Identifiable parameter combinations in a linear spatially homogeneous system
are a subset of the identifiable parameter combinations in all corresponding
spatially heterogeneous models subject to linear advection and/or diffusion
with general unknown initial conditions.

###### Proof.

By Theorem 1 we are guaranteed that the spatially heterogeneous system may be
written as set of polynomial relations, and, therefore, that we can express
the system as a set of polynomial relations involving only observed variables
and their derivatives. Thus, we can obtain a set of identifiable parameter
combinations, $\mathcal{Q}$, from the set of polynomial coefficients. Next, we
can recover the set of identifiable parameter combinations in the spatially
homogeneous model, $\mathcal{Q}^{*}$, by considering that all terms involving
spatial derivatives vanish. Thus, $\mathcal{Q}^{*}\subseteq\mathcal{Q}$ as
required. ∎

As stated previously, we will also be able to write semi-linear systems as
over-determined linear system of differential-algebraic equations, and thus we
expect Corollary 1 to also hold true for semi-linear systems. Furthermore, we
propose that one cannot increase the set of identifiable parameter
combinations by imposing additional constrains on the problem (for example, by
prescribing a particular form for the initial condition). Thus, we also expect
Corollary 1 to hold for systems with non-linear source terms by considering
identifiability of systems with no-flux boundary conditions and homogeneous
initial conditions, which reduce to the corresponding set of ODEs.

###  2.5 Semi-linear systems

We now briefly present and discuss analysis of two specific semi-linear
systems: a two-state model of bacteria migration due to chemotaxis, and a
three-state cell-cycle model subject to logistic growth.

####  2.5.1 Two-state model of bacteria chemotaxis

First, we consider a model of bacteria chemotaxis. Bacteria, $\rho$, secrete a
chemotactic factor, $c$, at rate $k$. Both bacteria and the chemotactic factor
diffuse, however bacteria additionally undergo directed motion due to the
spatial gradients in the chemotactic factor, with strength and direction given
by $\chi$ ($\chi>0$ and $\chi<0$ correspond to negative and regular
chemotaxis, respectively). The chemotactic factor degrades at constant rate
$\alpha$. The dynamics are described by the non-linear system [20]

$\displaystyle\rho^{(0,1)}$
$\displaystyle=D_{\rho}\rho^{(2,0)}+\overbrace{\chi(\rho
c^{(0,1)})^{(0,1)}\vphantom{\Big{(}}}^{\text{Chemotaxis}},$ $\displaystyle
c^{(0,1)}$ $\displaystyle=D_{c}c^{(2,0)}+k\rho-\alpha c.$

We consider a realistic experimental setup where only information relating to
the density of bacteria is available such that chemotactic factor is not
observed. The experiment is initiated without the chemotactic factor,
$c(x,0)=0$, and the field of view is such that no-flux boundary conditions are
appropriate. We consider the system to be semi-linear as it is linear the
unobserved quantity, $c$.

Results in Theorem 1 indicate that, after eliminating $c^{(2,0)}$ from one
equation, we can expand the system to include derivatives up to fourth order
to ensure that the system is fully determined. Performing row reduction on the
expanded system yields a polynomial expression containing 212 coefficients,
which determine the identifiable parameter combinations (full calculations
available as supplementary code) as

$\left\\{D_{\rho},D_{c},k,\alpha\chi\right\\}.$ (15)

Thus, $D_{\rho}$, $D_{c}$, and $k$ are structurally identifiable, while only
the product $\alpha\chi$ is identifiable while the constituents are not. This
analysis agrees with results that one would arrive through the scaling
$\tilde{c}=\alpha c$ to eliminate the disparate appearance of $\alpha$ and
$\chi$ in the resultant equations.

A key difference between analysis of the chemotaxis model and the general
linear model considered previously is the role of the initial condition. While
it may be true that $\alpha$ and $\chi$ would become individually identifiable
if information relating to a non-homogeneous non-zero initial condition were
known, the initial condition $c(x,0)$ conveys no information about the
parameters. Indeed, we can demonstrate that $\alpha$ and $\chi$ can compensate
for each other without varying the initial condition (fig. 3).

Figure 3: Non-identifiability of the chemotaxis model. Solutions show the
observed variable $\rho(x,t)$ (top row), and unobserved variable $c(x,t)$, at
the original parameter set (solid curves) and modified parameter set where
$\chi\mapsto 2\chi$ (red dashed). Original parameters are $D_{\rho}=20$,
$D_{c}=100$, $\chi=k=\alpha=1$. Note that the choice $\chi>0$ prescribes
negative chemotaxis: bacteria move down the chemotactic gradient.

####  2.5.2 Three-state cell-cycle model of cell migration subject to logistic
growth

Finally, we consider a three-state cell cycle model of cell migration that
extends on that presented in section 2.1 by capturing the intermediate stage
where cells appear to fluoresce yellow in FUCCI assays, and where contact
inhibition prevents cells in G2 from undergoing mitosis [31, 34]. The model is
given by

$\displaystyle r^{(0,1)}$
$\displaystyle=D_{1}r^{(2,0)}-\lambda_{1}r+2\lambda_{3}g(1-n/K),$ (16a)
$\displaystyle y^{(0,1)}$
$\displaystyle=D_{2}y^{(2,0)}-\lambda_{2}y+\lambda_{1}r,$ (16b) $\displaystyle
g^{(0,1)}$ $\displaystyle=D_{3}g^{(2,0)}-\lambda_{3}g(1-n/K)+\lambda_{2}y.$
(16c)

Here, the variable $y$ corresponds to cells that fluoresce both red and green
and hence appear to fluoresce yellow. We consider that FUCCI scratch assay
experiments are conducted (fig. 1b), but that we do not attempt to distinguish
the period of fluorescent overlap so that only measurements relating to $r$
and $g$ are available.

While this problem seems more complex—both due to the non-linearity introduced
by the logistic growth term and the additional state variable—it is relatively
straightforward to eliminate $y$ and its derivatives without expanding the
system through differentiation. To do this, we solve eq. 16b explicitly for
$y$ to obtain

$y=\dfrac{-\lambda_{3}g^{2}+\lambda_{3}g(K-r)+K\big{(}g^{(0,1)}-D_{3}g^{(2,0)}\big{)}}{K\lambda_{2}+\lambda_{3}g}$

and then substitute the resultant expression directly into eqs. 16a and 16b.
Where possible, this approach is advantageous in that it can be applied to any
non-linear partial differential equation model; we could, for example, apply
such a technique to a model where growth inhibition is modelled through
alternative growth models such as Gompertz or Richards (for models with non-
polynomial terms the orthogonality of terms will also need to be established
to construct the set of coefficients). Expanding the resultant system reveals
that all parameters are structurally identifiable, in agreement with the
result in Corollary 1, where we find that the corresponding ODE model is also
structurally identifiable according to the online tool SIAN [15, 35].

##  3 Discussion and Conclusion

Many biological processes are inherently spatial, and not well described by
ODE models. Yet, tools for assessing the identifiability of the PDE models
that are increasingly used to interpret such spatial data are almost entirely
undeveloped. In this brief article, we illustrate how the well-established
differential algebra approach to structural identifiability of ODE models can
be transferred to analyse a large class of spatial models.

While we demonstrate that the differential-algebra approach can be applied to
any system of linear RAD equations, we have restricted our demonstration to
two-dimensional systems. Results in Theorem 1 show that the derivative order
required grows rapidly as the number of state variables increases. This
presents not only computational difficulties, but practical issues with
interpretation: the resultant set of polynomial coefficients will be verbous
and potentially complicated, making it difficult to establish identifiable
parameter combinations. We see this even for the relatively simple semi-linear
chemotaxis model (supplementary code). Such computational issues are a well-
established shortcoming of the differential algebra approach, even for ODE
models [13]. While this bottleneck presents a clear shortcoming, many PDE
models used in practice only possess a small number of state variables:
perhaps more so than for ODE models, where large systems are common [36].
Common systems of PDE models that we can analyse include models of chemotaxis
[20], Turing patterning and morphogenesis [37], age-structured epidemic models
[25], systems of heat and wave equations, and many more. For some systems, one
approach to deal with larger numbers of state variables is through a
transformation that decouples states [38]; identifiability can then be
established for species sequentially.

Another approach to alleviate the computational cost is to develop more
efficient algorithmic implementations of state-variable elimination that do
not involve expanding to a fully determined system. In the ODE literature,
Ritt’s algorithm is commonly employed alongside a Gröbner basis factorisation
to eliminate state variables automatically [12, 39]. An advantage of eventual
implementation of these algorithms for analysis of PDE models is their
application to non-linear problems: computing the Gröbner basis is equivalent
to Gaussian elimination for linear systems of equations, but can be applied
generally to any polynomial system. While it remains unclear whether this can,
given infinite computational time, always be applied to PDE models (as we have
shown is the case for linear systems), such software will enable analysis of a
much larger class of PDE models that can be written as polynomial in both
state variables and their corresponding derivatives.

A unique feature in the application of PDE models to interpret data is the
mode of measurements that may be collected. It is not always the case that
experimental observations correspond to measurements of a spatiotemporal
function: rather, observations could comprise scaler measurements at a single
point in space (for example, measurements of temperature by a probe in a heat
conduction experiment [40]) or scalar measurements of a spatial average (i.e.
of total cell count). Summary-statistic type measurements can be even more
complicated for models with two or more spatial dimensions [41]. The scratch
assay experiment is an example of this, as data collected from the full, two-
dimensional, process often comprises spatial averages in the direction
parallel to the scratch [32]. This example is trivial as, provided appropriate
constraints are placed on the initial condition, the extraneous spatial
dimension can be integrated out of the full two-dimensional model to yield the
one-dimensional model that we analyse. It is entirely unclear whether it is
possible, in the general case, to reduce a system of partial differential
equations to a polynomial system involving an observation function that does
not include space.

This question of lower-dimensional observation functions is particularly
relevant to structural identifiability analysis of stochastic differential
equation models [42]. It is unclear whether such equations can be reduced to
eliminate unobserved states, since the underlying stochastic process is
typically not differentiable. One approach is to analyse the corresponding
system of Fokker-Plank equations, a PDE in as many spatial dimensions as SDE
state variables. Partially observing the system, therefore, corresponds to
observing marginals of the Fokker-Plank equation solution; in effect,
integrating over unobserved spatial variables.

PDE models are widely used to characterise spatial processes and interpret
spatial biological data. It is, therefore, paramount that effective and
efficient tools to assess the structural identifiability of these models be
developed. We build on existing work [25, 26] to show that the differential
algebra approach to structural identifiability analysis can be applied to all
linear RAD PDE models, and some classes of non-linear PDE models. We
demonstrate the interdependence between structural identifiability and the
initial conditions in partially observed models, highlighting the importance
of assessing the structural identifiability before attempting to infer
parameters in PDE models from experimental data.

## Data availability

Mathematica and Julia code used to perform the symbolic and numerical
computations, respectively, available on GitHub at https://github.com/ap-
browning/pde_structural_identifiability.

## Acknowledgements

APB thanks the Mathematical Institute for a Hooke Fellowship. MT and APB thank
the EPSRC and Mathematical Institute for funding a vacation internship that
funded research conducted by MT. CF acknowledges support of a fellowship from
“la Caixa” Foundation (ID 100010434) with code LCF/BQ/EU21/11890128.

## References

  * [1] Liepe J, Kirk P, Filippi S, Toni T, Barnes CP, Stumpf MPH. 2014.  A framework for parameter estimation and model selection from experimental data in systems biology using approximate Bayesian computation.  _Nature Protocols_ 9:439 456. (doi:10.1038/nprot.2014.025).
  * [2] Gábor A, Banga JR. 2015.  Robust and efficient parameter estimation in dynamic models of biological systems.  _BMC Systems Biology_ 9:74. (doi:10.1186/s12918-015-0219-2).
  * [3] Bellman R, Åström KJ. 1970.  On structural identifiability.  _Mathematical Biosciences_ 7:329–339. (doi:10.1016/0025-5564(70)90132-x).
  * [4] Walter E, Lecourtier Y. 1981.  Unidentifiable compartmental models: what to do?  _Mathematical Biosciences_ 56:1–25. (doi:10.1016/0025-5564(81)90025-0).
  * [5] Walter E. 1987.  _Identifiability of Parametric Models_.  Elsevier Science & Technology, London, United Kingdom.  (doi:10.1016/C2013-0-03836-4).
  * [6] Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, KlingmÃ¼ller U, Timmer J. 2009.  Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood.  _Bioinformatics_ 25:1923–1929. (doi:10.1093/bioinformatics/btp358).
  * [7] Chis OT, Banga JR, Balsa-Canto E. 2011.  Structural identifiability of systems biology models: a critical comparison of methods.  _PLOS One_ 6:e27755. (doi:10.1371/journal.pone.0027755).
  * [8] Villaverde AF, Barreiro A, Papachristodoulou A. 2016.  Structural identifiability of dynamic systems biology models.  _PLOS Computational Biology_ 12:e1005153. (doi:10.1371/journal.pcbi.1005153).
  * [9] Meshkat N, Eisenberg M, DiStefano JJ. 2009.  An algorithm for finding globally identifiable parameter combinations of nonlinear ODE models using Gröbner Bases.  _Mathematical Biosciences_ 222:61–72. (doi:10.1016/j.mbs.2009.08.010).
  * [10] Eisenberg MC, Jain HV. 2017.  A confidence building exercise in data and identifiability: Modeling cancer chemotherapy as a case study.  _Journal of Theoretical Biology_ 431:63–78. (doi:10.1016/j.jtbi.2017.07.018).
  * [11] Ljung L, Glad T. 1994.  On global identifiability for arbitrary model parametrizations.  _Automatica_ 30:265–276. (doi:10.1016/0005-1098(94)90029-9).
  * [12] Bellu G, Saccomani MP, Audoly S, D’Angiò L. 2007.  DAISY: A new software tool to test global identifiability of biological and physiological systems.  _Computer Methods and Programs in Biomedicine_ 88:52–61. (doi:10.1016/j.cmpb.2007.07.002).
  * [13] Raue A, Karlsson J, Saccomani MP, Jirstrand M, Timmer J. 2014.  Comparison of approaches for parameter identifiability analysis of biological systems.  _Bioinformatics_ 30:1440–1448. (doi:10.1093/bioinformatics/btu006).
  * [14] Barreiro XR, Villaverde AF. 2023.  Benchmarking tools for a priori identifiability analysis.  _Bioinformatics_ 39:btad065. (doi:10.1093/bioinformatics/btad065).
  * [15] Hong H, Ovchinnikov A, Pogudin G, Yap C. 2019.  SIAN: software for structural identifiability analysis of ODE models.  _Bioinformatics_ 35:2873–2874. (doi:10.1093/bioinformatics/bty1069).
  * [16] Meshkat N, Kuo CEz, DiStefano J. 2014.  On finding and using identifiable parameter combinations in nonlinear dynamic systems biology models and COMBOS: a novel web implementation.  _PLOS One_ 9:e110261. (doi:10.1371/journal.pone.0110261).
  * [17] Wieland FG, Hauber AL, Rosenblatt M, TÃ¶nsing C, Timmer J. 2021.  On structural and practical identifiability.  _Current Opinion in Systems Biology_ 25:60–69. (doi:10.1016/j.coisb.2021.03.005).
  * [18] Ligon TS, Frohlich F, Chis OT, Banga JR, Balsa-Canto E, Hasenauer J. 2017.  GenSSI 2.0: multi-experiment structural identifiability analysis of SBML models.  _Bioinformatics_ 34:1421–1423. (doi:10.1093/bioinformatics/btx735).
  * [19] Vajda S, Godfrey KR, Rabitz H. 1989.  Similarity transformation approach to identifiability analysis of nonlinear compartmental models.  _Mathematical Biosciences_ 93:217–248. (doi:10.1016/0025-5564(89)90024-2).
  * [20] Murray JD. 2002.  _Mathematical Biology_.  Springer-Verlag, Berlin. 3 edition.  (doi:10.1007/b98868).
  * [21] Simpson MJ, Baker RE, Vittadello ST, Maclaren OJ. 2020.  Practical parameter identifiability for spatio-temporal models of cell invasion.  _Journal of The Royal Society Interface_ 17:20200055. (doi:10.1098/rsif.2020.0055).
  * [22] Liu Y, Suh K, Maini PK, Cohen DJ, Baker RE. 2023.  Parameter identifiability and model selection for partial differential equation models of cell invasion.  _arXiv_ .
  * [23] Kicheva A, Pantazis P, Bollenbach T, Kalaidzidis Y, Bittig T, Jülicher F, González-Gaitán M. 2007.  Kinetics of Morphogen Gradient Formation.  _Science_ 315:521–525. (doi:10.1126/science.1135774).
  * [24] Romanova-Michaelides M, Hadjivasiliou Z, Aguilar-Hidalgo D, Basagiannis D, Seum C, Dubois M, Jülicher F, Gonzalez-Gaitan M. 2022.  Morphogen gradient scaling by recycling of intracellular Dpp.  _Nature_ 602:287–293. (doi:10.1038/s41586-021-04346-w).
  * [25] Renardy M, Kirschner D, Eisenberg M. 2022.  Structural identifiability analysis of age-structured PDE epidemic models.  _Journal of Mathematical Biology_ 84:9. (doi:10.1007/s00285-021-01711-1).
  * [26] Ciocanel MV, Ding L, Mastromatteo L, Reichheld S, Cabral S, Mowry K, Sandstede B. 2023.  Parameter identifiability in PDE models of fluorescence recovery after photobleaching.  _arXiv_ .
  * [27] Margaria G, Riccomagno E, Chappell MJ, Wynn HP. 2001.  Differential algebra methods for the study of the structural identifiability of rational function state-space models in the biosciences.  _Mathematical Biosciences_ 174:1–26. (doi:10.1016/s0025-5564(01)00079-7).
  * [28] Browning AP, Jin W, Plank MJ, Simpson MJ. 2020.  Identifying density-dependent interactions in collective cell behaviour.  _Journal of The Royal Society Interface_ 17:20200143. (doi:10.1098/rsif.2020.0143).
  * [29] Carr MJ, Simpson MJ, Drovandi C. 2021.  Estimating parameters of a stochastic cell invasion model with fluorescent cell cycle labelling using approximate Bayesian computation.  _Journal of the Royal Society Interface_ 18:20210362. (doi:10.1098/rsif.2021.0362).
  * [30] Sakaue-Sawano A, Kurokawa H, Morimura T, Hanyu A, Hama H, Osawa H, Kashiwagi S, Fukami K, Miyata T, Miyoshi H, Imamura T, Ogawa M, Masai H, Miyawaki A. 2008.  Visualizing spatiotemporal dynamics of multicellular cell-cycle progression.  _Cell_ 132:487–498. (doi:10.1016/j.cell.2007.12.033).
  * [31] Vittadello ST, McCue SW, Gunasingh G, Haass NK, Simpson MJ. 2018.  Mathematical models for cell migration with real-time cell cycle dynamics.  _Biophysical Journal_ 114:1241–1253. (doi:10.1016/j.bpj.2017.12.041).
  * [32] Johnston ST, Shah ET, Chopin LK, McElwain DLS, Simpson MJ. 2015.  Estimating cell diffusivity and cell proliferation rate by interpreting IncuCyte ZOOMTM assay data using the Fisher-Kolmogorov model.  _BMC Systems Biology_ 9:38. (doi:10.1186/s12918-015-0182-y).
  * [33] Cobelli C, DiStefano JJ. 1980.  Parameter and structural identifiability concepts and ambiguities: a critical review and analysis.  _American Journal of Physiology-Regulatory, Integrative and Comparative Physiology_ 239:7—24. (doi:10.1152/ajpregu.1980.239.1.r7).
  * [34] Donker L, Houtekamer R, Vliem M, Sipieter F, Canever H, Gï¿½mez-Gonzï¿½lez M, Bosch-Padrï¿½s M, Pannekoek WJ, Trepat X, Borghi N, Gloerich M. 2022.  A mechanical G2 checkpoint controls epithelial cell division through E-cadherin-mediated regulation of Wee1-Cdk1.  _Cell Reports_ 41:111475. (doi:10.1016/j.celrep.2022.111475).
  * [35] Hong H, Ovchinnikov A, Pogudin G, Yap C. 2020.  Global identifiability of differential models.  _Communications on Pure and Applied Mathematics_ 73:1831–1879. (doi:10.1002/cpa.21921).
  * [36] Karlsson J, Anguelova M, Jirstrand M. 2012.  An efficient method for structural identifiability analysis of large dynamic systems.  _IFAC Proceedings Volumes_ 45:941–946. (doi:10.3182/20120711-3-be-2027.00381).
  * [37] Turing AM. 1952.  The chemical basis of morphogenesis.  _Philosophical Transactions of the Royal Society of London Series B, Biological Sciences_ 237:37–72. (doi:10.1098/rstb.1952.0012).  URL https://doi.org/10.1098/rstb.1952.0012.
  * [38] Sun Y, Petersen J, Clement T. 1999.  Analytical solutions for multiple species reactive transport in multiple dimensions.  _Journal of Contaminant Hydrology_ 35:429–440. (doi:10.1016/s0169-7722(98)00105-3).
  * [39] Meshkat N, Anderson C, III JJD. 2012.  Alternative to Ritt’s pseudodivision for finding the input-output equations of multi-output models.  _Mathematical Biosciences_ 239:117–123. (doi:10.1016/j.mbs.2012.04.008).
  * [40] McInerney S, Carr EJ, Simpson MJ. 2019.  Parameterising continuum models of heat transfer in heterogeneous living skin using experimental data.  _International Journal of Heat and Mass Transfer_ 128:964–975. (doi:10.1016/j.ijheatmasstransfer.2018.09.054).
  * [41] Browning AP, Maclaren OJ, Buenzli PR, Lanaro M, Allenby MC, Woodruff MA, Simpson MJ. 2021.  Model-based data analysis of tissue growth in thin 3D printed scaffolds.  _Journal of Theoretical Biology_ 528:110852. (doi:10.1016/j.jtbi.2021.110852).
  * [42] Browning AP, Warne DJ, Burrage K, Baker RE, Simpson MJ. 2020.  Identifiability analysis for stochastic differential equation models in systems biology.  _Journal of The Royal Society Interface_ 17:20200652. (doi:10.1098/rsif.2020.0652).