Abstract (original article): 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 reaction–advection–diffusion PDE models that are linear in the unobserved quantities. We show that the differential algebra approach can always, in theory, be applied to such models. Moreover, despite the perceived complexity introduced by the addition of advection and diffusion terms, consideration of spatial analogues of non-spatial models cannot exacerbate structural identifiability. We conclude by discussing future possibilities and the computational cost of performing structural identifiability analysis on more general PDE models. Citation (original article): Browning AP, Tașcǎ M, Falcó C, Baker RE. 2024. Structural identifiability analysis of linear reaction-advection-diffusion processes in mathematical biology. Proceedings of the Royal Society A 480:20230911. https://doi.org/10.1098/rspa.2023.0911

Identifiability of a reaction-advection-diffusion model

In this notebook, we reproduce the results in Section 2.2 of the original article to establish the structural identifiability of the reaction-advection-diffusion model
∂u
∂t
=
D
u
2
∂
u
∂
2
x
+
α
u
∂u
∂x
+
p
1
u+
p
2
v+
p
3
,​​
∂v
∂t
=
D
v
2
∂
v
∂
2
x
+
α
v
∂v
∂x
+
p
4
u+
p
5
v+
p
6
,
subject to non-informative (i.e., Neumann) boundary conditions, and observations of only the sum
n=u+v
.

Problem definition

First, we need to define the set of equations. These are defined as a set of expressions that satisfy the PDE when vanished.
In[]:=
sys1={​​D[u[x,t],t]-(Du*D[u[x,t],{x,2}]+αu*D[u[x,t],x]+p1*u[x,t]+p2*v[x,t]+p3),​​D[v[x,t],t]-(Dv*D[v[x,t],{x,2}]+αv*D[v[x,t],x]+p4*u[x,t]+p5*v[x,t]+p6)​​};
We’ll also extract the parameters (i.e., everything not a state variable or a derivative of a state variable) in the above expression.
In[]:=
params=Quiet[Complement[Variables[sys1],Join[{u[x,t],v[x,t]},Select[Variables[sys1],#[[0]][[1]]===u||#[[0]][[1]]===v&]]]]
Out[]=
{Du,Dv,p1,p2,p3,p4,p5,p6,αu,αv}
Finally, we need to define our observation function. We do this in terms of the substitution that we need to make to eliminate u from the system of PDEs; thus arriving at an equivalent set involving only the observed state variable n and the variable that we wish to eliminate, v.
In[]:=
obs=u[x,t]->n[x,t]-v[x,t];

Reduction to a system with at least one first order equation

Next, we make the substitution. Note that we must also substitute all required derivatives of the observation substitution.
In[]:=
sys2=sys1/.Table[D[obs,{t,i[[1]]},{x,i[[2]]}],{i,{{0,0},{1,0},{0,1},{0,2}}}];
Now, we manipulate the system such that at least one equation is first order in the variable that we wish to eliminate. This reduces the number of derivatives we must take when expanding the system later. To do this, we solve the system for
v
xx
and
v
t
. The corresponding expressions for each become the new system (the expression for
v
t
will not contain
v
xx
and will, therefore, be of first order.
In[]:=
sol1=Solve[Table[expr==0,{expr,sys2}],{D[v[x,t],{x,2}],D[v[x,t],t]}][[1]];​​sys3={​​D[v[x,t],{x,2}]-(D[v[x,t],{x,2}]/.sol1),​​D[v[x,t],t]-(D[v[x,t],t]/.sol1)​​}
Out[]=

1
Du-Dv
-p3-p6-p1n[x,t]-p4n[x,t]+p1v[x,t]-p2v[x,t]+p4v[x,t]-p5v[x,t]+
(0,1)
n
[x,t]-αu
(1,0)
n
[x,t]+αu
(1,0)
v
[x,t]-αv
(1,0)
v
[x,t]-Du
(2,0)
n
[x,t]+
(2,0)
v
[x,t],
(0,1)
v
[x,t]+
1
Du-Dv
-Dvp3-Dup6-Dvp1n[x,t]-Dup4n[x,t]+Dvp1v[x,t]-Dvp2v[x,t]+Dup4v[x,t]-Dup5v[x,t]+Dv
(0,1)
n
[x,t]-Dvαu
(1,0)
n
[x,t]+Dvαu
(1,0)
v
[x,t]-Duαv
(1,0)
v
[x,t]-DuDv
(2,0)
n
[x,t]

Expansion to a fully determined system

We must now expand the system by taking repeated derivatives to ensure that the system is fully determined (i.e., there are more equations than there are derivatives of v that we wish to eliminate. Generally, one must expand fully up to a given derivative (see Theorem 1 of the original article): here, we are guaranteed that the system will be fully determined if we expand to derivatives of order 4. However, here we find that we can expand to a much smaller system, provided we check that the system is still fully determined.
In[]:=
sys4=Join[sys3,{​​D[sys3[[1]],x],​​D[sys3[[1]],t],​​D[sys3[[2]],x],​​D[sys3[[2]],{x,2}]​​}];
Let us now convert the system of expressions into the matrix system Av + b = 0, where v is a vector of the state v and it’s derivatives (which we wish to eliminate), and b contains the residual terms (i.e., terms involving only n, which is observed).
In[]:=
vvars=Join[{v[x,t]},Quiet[Select[Variables[sys4],#[[0]][[1]]===v&]]];​​A=Table[Coefficient[Expand[expr],var],{expr,sys4},{var,vvars}];​​b=Table[ConstantArray[0,Length[vvars]]/.CoefficientRules[expr,vvars],{expr,sys4}];
The system will be fully determined if A is not of full rank, which is, indeed, the case here:
In[]:=
Length[A]-MatrixRank[A]
Out[]=
1
The nullity of A gives us the number of polynomial equations that we will obtain. Note that we should still validate that we expect only one polynomial equation for this model by examining the nullity of the fully expanded system (for full details, see Appendix A of the original article). To do this, we fully expand the system to order four, and construct the appropriate augmented matrix system:
sysfull=Join[​​Catenate[Table[Table[D[sys3[[1]],{x,i},{t,k-i}],{i,0,k}],{k,0,2}]],​​Catenate[Table[Table[D[sys3[[2]],{x,i},{t,k-i}],{i,0,k}],{k,0,3}]]​​];​​vvarsfull=Join[{v[x,t]},Quiet[Select[Variables[sysfull],#[[0]][[1]]===v&]]];​​A2=Table[Coefficient[Expand[expr],var],{expr,sysfull},{var,vvarsfull}];​​b2=Table[ConstantArray[0,Length[vvarsfull]]/.CoefficientRules[expr,vvarsfull],{expr,sysfull}];​​M2=Transpose[Insert[Transpose[A2],b2,Length[vvarsfull]+1]];​​MatrixRank[M2]-MatrixRank[A2]
Out[]=
1
Since the difference rank is one, we do, indeed, expect only a single polynomial equation.

Reduction to a polynomial equation

We’ll now perform row reduction on the augmented matrix A. We expect the rank of A to be at least one less than it’s height: this is indicative of the single polynomial expression we expect. However, the actual rank will possibly be less than this. We’ll augment A with an identity matrix so that we can apply the same operations to b (since Mathematica outputs only the reduced row echelon form).
In[]:=
op=RowReduce[Join[A,IdentityMatrix[Length[sys4]],2]][[;;,-Length[sys4];;]];
The nullity of A will correspond to the number of elements in the row-reduced b that must vanish for the system to be satisfied. These expressions give the set of polynomial equations. Note that this may be more than the number of polynomial equations we expect as these elements may not be linearly independent.
In[]:=
expr=Expand[FullSimplify[(op.b)[[-1]]]]
Out[]=
-
p3p4
Dv(p1-p2)+Du(p4-p5)
+
p3p5
Dv(p1-p2)+Du(p4-p5)
+
p1p6
Dv(p1-p2)+Du(p4-p5)
-
p2p6
Dv(p1-p2)+Du(p4-p5)
-
p2p4n[x,t]
Dv(p1-p2)+Du(p4-p5)
+
p1p5n[x,t]
Dv(p1-p2)+Du(p4-p5)
-
p1
(0,1)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
-
p5
(0,1)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
(0,2)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
p5αu
(1,0)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
p1αv
(1,0)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
-
αu
(1,1)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
-
αv
(1,1)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
Dvp1
(2,0)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
Dup5
(2,0)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
αuαv
(2,0)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
-
Du
(2,1)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
-
Dv
(2,1)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
Dvαu
(3,0)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
Duαv
(3,0)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)
+
DuDv
(4,0)
n
[x,t]
Dv(p1-p2)+Du(p4-p5)

Extracting identifiable coefficients

The full set of identifiable coefficients will be proportional to the coefficients of monomials in n and it’s derivatives.
In[]:=
nvars=Join[{n[x,t]},Quiet[Select[Variables[sys4],#[[0]][[1]]===n&]]];​​coef1=Values@CoefficientRules[expr,nvars]
Out[]=
-
p2p4
Dvp1-Dvp2+Dup4-Dup5
+
p1p5
Dvp1-Dvp2+Dup4-Dup5
,-
p1
Dvp1-Dvp2+Dup4-Dup5
-
p5
Dvp1-Dvp2+Dup4-Dup5
,
1
Dvp1-Dvp2+Dup4-Dup5
,
p5αu
Dvp1-Dvp2+Dup4-Dup5
+
p1αv
Dvp1-Dvp2+Dup4-Dup5
,-
αu
Dvp1-Dvp2+Dup4-Dup5
-
αv
Dvp1-Dvp2+Dup4-Dup5
,
Dvp1
Dvp1-Dvp2+Dup4-Dup5
+
Dup5
Dvp1-Dvp2+Dup4-Dup5
+
αuαv
Dvp1-Dvp2+Dup4-Dup5
,-
Du
Dvp1-Dvp2+Dup4-Dup5
-
Dv
Dvp1-Dvp2+Dup4-Dup5
,
Dvαu
Dvp1-Dvp2+Dup4-Dup5
+
Duαv
Dvp1-Dvp2+Dup4-Dup5
,
DuDv
Dvp1-Dvp2+Dup4-Dup5
,-
p3p4
Dvp1-Dvp2+Dup4-Dup5
+
p3p5
Dvp1-Dvp2+Dup4-Dup5
+
p1p6
Dvp1-Dvp2+Dup4-Dup5
-
p2p6
Dvp1-Dvp2+Dup4-Dup5

We’ll normalise so that one coefficient is unity.
In[]:=
coef2=FullSimplify[coef1/coef1[[3]]]
Out[]=
{-p2p4+p1p5,-p1-p5,1,p5αu+p1αv,-αu-αv,Dvp1+Dup5+αuαv,-Du-Dv,Dvαu+Duαv,DuDv,p3(-p4+p5)+(p1-p2)p6}
To establish the identifiable parameter combinations, we can (manually) solve for the individual parameters in terms of the above coefficients. In this case, since
D
u
+
D
v
and
D
u
D
v
are identifiable, so to are the constituents
D
u
and
D
v
. A similar observation can be made for
α
u
and
α
v
.
In conclusion, the set of identifiable parameter combinations is given by
In[]:=
coef3={p1,p5,αu,αv,Du,Dv,p2p4,p3(-p4+p5)+(p1-p2)p6}
Out[]=
{p1,p5,αu,αv,Du,Dv,p2p4,p3(-p4+p5)+(p1-p2)p6}

CITE THIS NOTEBOOK

Structural identifiability analysis of linear reaction–advection–diffusion processes in mathematical biology​
by Alexander Browning, Maria Taşcă, Carles Falcó and Ruth Baker
Wolfram Community, STAFF PICKS, April 11, 2024
​https://community.wolfram.com/groups/-/m/t/3156653