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

,


Introduction
Mathematical models play an irreplaceable 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].The question of whether it is possible, in theory, for model parameters to be uniquely identified given a specified mathematical model and a specified 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 well 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 differentialalgebra-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 available [3,8,17].For such partially observed systems, 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 alternative 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 [7,19].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].Significantly, Renardy et al. [25] demonstrate that the differential-algebra approach may be applicable more generally to systems of linear first-order PDE models.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, 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 behaviour.
In this brief article, we present a methodology and a series of general results for the structural identifiability of linear RAD PDE models subject to partial observation.We also demonstrate that our methodology can be applied to any non-linear RAD PDE model that is linear in the unobserved quantities: we refer to such models as semi-linear.For a more thorough and formal introduction to structural identifiability and the differential algebra approach, we direct the reader to [25,27].
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.

Linear cell-cycle model
To demonstrate the differential algebra approach for linear PDE models, we first consider a twostate cell cycle model of cell migration subject to exponential growth [21,31].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 ) for t > 0 and x ∈ [0, L], and 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; and k 2 the mitotic rate of green cells.For compactness, Snapshots collected at 0 h and 48 h from a typical scratch assay experiment conducted using PC3 prostate cancer cells [28].Averaging in the direction of the scratch provides onedimensional 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.
from this point on we denote derivatives using superscripts: r (i,j) = ∂ i+j r/∂x i ∂t j and denote the undifferentiated variables interchangeably without r (0,0) ≡ r.We prescribe homogeneous Neumann (no-flux) boundary conditions on a domain x ∈ [0, L], such that no additional information is available from the equations that govern the behaviour at the boundaries.The initial condition is assumed to be general: r(x, 0) = r 0 (x) and g(x, 0) = g 0 (x).
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 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); i.e., by measuring only the total cell density?
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 to eliminate r.At this point, it is clear that the analysis is far simpler in the case that 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.(1b) and (3) simultaneously for g (0,1) and g (2,0) in terms of n and its derivatives to obtain 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 both the previous and all following steps may be carried out using only multiplication and elimination to avoid division by zero in determining 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, 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.
A pathological example that illustrates why we must do this is to consider that we could otherwise multiply eq. ( 5) through by an arbitrarily chosen parameter.This would give the impression that the chosen parameter is identifiable through the coefficient of n (0,2) .We take a more practical approach to ensure uniqueness and divide the resultant polynomial through by an arbitrarily chosen coefficient.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 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, . ..), 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.

Generic two-state reaction-advection-diffusion system
Having established the general procedure for applying the differential algebra framework to reactiondiffusion models, we now consider the general two-state model for t > 0 and subject to observations of the form n(x, t) = u(x, t)+v(x, t).The initial conditions are fully prescribed by n(x, 0) = n 0 (x) and v(x, 0) = v 0 (x).We do not consider a domain or boundary conditions, only assuming that no information is available through observation of the system at the boundaries.
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 αu(x, t) + βv(x, t) where α and β 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 We proceed by solving eqs.(7b) and ( 8) for v (2,0) , v (0,1) 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 ith order 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 where v (i) ∈ R (i+1)(i+2)/2 denotes a vector of all possible ith order derivatives.We then proceed by performing Gaussian elimination to reduce the augmented matrix (A|b) into row echelon form.
This provides a set of expressions, based upon linear combinations of the elements of 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: the number of polynomial equations we require for linear and semi-linear problems will be given by the difference in rank of A and the rank of the augmented matrix (A|b) from the overdetermined system.For the two-state RAD system, we thus require a single polynomial equation.
It happens in our case that we are able to obtain this single polynomial equation from a subset of third order equations.Algebraic manipulations are performed in Mathematica [34] and available as supplementary code with algorithmic details provided in Appendix A.
For the system defined by eq. ( 7), we arrive at the set of polynomial coefficients eq10Equation (10) provides an exhaustive set of identifiable parameter combinations.We then reduce this set further to establish a fully-reduced set of identifiable parameters and parameter combimnations.For example, denoting two combinations by c 6 = D u + D v and c 8 = D u D v we can solve for the parameters D u and D v in terms of the identifiable parameter combinations c 6 and c 8 , thus establishing that D u and D v are identifiable.For eq. ( 7), the reduced set of identifiable parameter combinations is given by By setting α u = α v = 0 in eq. ( 10), we see that these results are identical to the no-advection model.However, results differ slightly for the non-spatial model.Setting to remove spatial derivatives, we see that the reduced set of identifiable parameter combinations is now so we can no longer identify any rate parameters individually.This result demonstrates that we are potentially able to learn more about the process by introducing spatial heterogeneity into the system.

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 0 (x).A key point of distinction between such partially observed PDE models and their ODE counterparts is that the initial condition is an unknown function 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 can consider that the right-hand-side of eq. ( 8) is "observed" such that the functional 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.
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 → 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.

Generic linear systems
A question naturally arising from 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 ℓ 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., nth 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 : is straightforward to see that n * = 4(m − 1), as required.
Theorem 1 highlights the potential computational cost of performing structural identifiability for PDE models with large numbers of states: 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.While analysis of the two-state system required a walltime of approximately 1 s, attempting to analyse a corresponding generic three-state model subject to a single observation did not yield results within a 12 h runtime (Mathematica 13.1 [34], Apple M1 Pro with 16GB RAM).
Importantly, however, theorem 1 demonstrates that it is always theoretically possible to perform structural identifiability analysis on linear RAD models.Furthermore, the Gaussian-eliminationbased reduction procedure requires only that the system be linear in the unobserved quantities, thus Theorem 1 will 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 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, Q, from the set of polynomial coefficients.Next, we can recover the set of identifiable parameter combinations in the spatially homogeneous model, Q * , by considering that all terms involving spatial derivatives vanish.Thus, Q * ⊆ 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 consequentially Corollary 1 will also hold for semilinear systems.Furthermore, we propose that one cannot increase the set of identifiable parameter combinations by imposing additional constraints 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.

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.

Two-state model of bacteria chemotaxis
First, we consider a model of bacteria chemotaxis.Bacteria, ρ(x, t), secrete a factor, c(x, t), at rate k > 0. 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 χ (χ > 0 and χ < 0 correspond to negative and regular chemotaxis, respectively).The chemotactic factor degrades at constant rate α > 0. The dynamics are described by the non-linear system [20] ρ (0,1) = D ρ ρ (2,0) + Chemotaxis χ(ρc (0,1) ) (0,1) , c (0,1) = D c c (2,0) + kρ − αc, for t > 0 and x ∈ [0, L].We consider a realistic experimental setup where only information relating to the density of bacteria is available (i.e., the concentration of the 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 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 Thus, D ρ , D c , and k are structurally identifiable as is the product αχ, while the individual constituents α and χ are not.This analysis agrees with results that one would arrive through the scaling c = αc to eliminate the disparate appearance of α and χ in eq. ( 14).
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 α and χ 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 α and χ can compensate for each other without varying the initial condition (fig.3).

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,35].The model is given by 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 = −λ 3 (g (0,0) ) 2 + λ 3 g (0,0) (K − r (0,0) ) + K g (0,1) − D 3 g (2,0) Kλ 2 + λ 3 g (0,0) , and then substitute the resultant expression directly into eqs.(16a) and (16b).This approach is advantageous in that it can be applied to any non-linear partial differential equation model in which we can solve for the unobserved quantities in terms of the observed quantities and respective derivatives; 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,36].

Discussion and Conclusion
Many 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 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-and three-state systems in one spatial dimension.Results in Theorem 1 show that the derivative order required grows rapidly as the number of state variables increases; we expect this to be exacerbated for PDE models with more than one spatial dimension.This presents not only computational but practical issues with interpretation: the resultant set of polynomial coefficients will be verbous and potentially complicated, making it difficult to establish a reduced set of 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 [37].Common systems of PDE models that we can analyse include models of chemotaxis [20], Turing patterning and morphogenesis [38], 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 [39]; 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,40].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 [41].
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 [42]) 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 [43].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 [44].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.

A Algorithm of implementation
In algorithm 1 we provide an algorithm that outlines our implementation of the differential algebra approach to assess the structural identifiability of m-state RAD PDE subject to ℓ linearly independent observations.An implementation of the analysis in Mathematica [34] is available as supplementary material on GitHub (https://github.com/ap-browning/pde_structural_identifiability).
Algorithm 1 Differential algebra for linear and semi-linear RAD PDE models 1: Given a system of m equations, denoted Y, and a set of ℓ < m linearly-independent observed quantities.2: Rewrite the system Y in terms of only the ℓ observed quantities and a remaining set of k = m−ℓ unobserved quantities that are to be eliminated.The rewritten system of m equations is denoted by Y 1 .3: Reduce Y 1 such that the first equation is first-order in the unobserved quantities.The reduced system is denoted by Y 2 .This can be done by, for example, solving the last m − 1 equations simultaneously for the second-order spatial derivatives of the unobserved quantities and substituting the resultant expressions into the first equation.4: Apply Theorem 1 to determine that the required order of the expanded system is at most n * = 4(m − 1).5: Expand Y 2 up to order n * by taking all possible order n * − 1 partial derivatives of the first equation, and all possible order n * − 2 partial derivatives of the remaining equations.The expanded system, denoted by Y 3 will comprise m equations, linear in a total of ñ < m partial derivatives of the unobserved variables.6: Write the expanded system as the linear system Av (n * ) = b, where A ∈ R ( m,ñ) and b ∈ R ( m) .7: Perform Gaussian elimination to reduce the symbolic augmented matrix (A|b) into row-echelon form, denoted (A re |b re ).8: The set of polynomial equations, denoted by R, is given by the non-zero elements of b re corresponding to rows of A re that are identically zero.Elements of R are polynomial in the set of partial derivatives of the observed ℓ quantities and are also considered to be observed.9: Normalise each polynomial to ensure uniqueness by dividing through by a chosen non-zero coefficient.The set of monic polynomials is denoted by R 1 .10: The set of identifiable parameter combinations, Q, is given by union of the the sets of polynomial coefficients of each element of R 1 .Coefficients that do not depend on the unknown quantities should be removed before further reduction.

Figure 1 .
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 h and 48 h from a typical scratch assay experiment conducted using PC3 prostate cancer cells[28].Averaging in the direction of the scratch provides onedimensional 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.

Figure 2 .
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 → 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, α 1 = −5, α 2 = 5, p i = 0.1 (i ̸ = 2) and p 2 = 0.2.
(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 → 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