Structural Identifiability

The StructuralIdentifiability.jl package is used to provide local and global identifiability of ODE models. It has been branched to the Filtration.jl package to provide identifiability of models of membrane filtration systems. To do so, let's first define a problem.

using Filtration

Jv = 2.71e-6; JBW = 5.66e-6;
C = 2184e-3; β = 1e-6; η = 4.9e-3
k1f = 2.0769; k2f = 0.0185; k1BW = 9.0551; k2BW = 0.0870
p = [Jv, JBW, C, β, η, k1f, k2f, k1BW, k2BW]

f₁(m,p) = p[3] * p[1] - p[3] * p[1] * p[4] * m
f₂(m,p) = -p[5] * m
g₁(m,p) = p[1]
g₂(m,p) = -p[2]
e₁(m,p) = p[6] * m + p[7]
e₂(m,p) = p[8] * m + p[9]

model = Model(e₁, e₂, f₁, f₂, g₁, g₂, p, :min; initial_guess = 20.);
nothing; #hide

To check the identifiability of parameters and states, you can use the assess_identifiability function on the constructed model. Here's an example:

assess_identifiability(model)

OrderedCollections.OrderedDict{SymbolicUtils.BasicSymbolic{Real}, Symbol} with 12 entries:
  x1(t) => :globally
  x2(t) => :nonidentifiable
  x3(t) => :globally
  p[1]  => :globally
  p[2]  => :globally
  p[3]  => :nonidentifiable
  p[4]  => :nonidentifiable
  p[5]  => :globally
  p[6]  => :nonidentifiable
  p[7]  => :globally
  p[8]  => :nonidentifiable
  p[9]  => :globally

To attempt to reparametrize the model in the case that states and parameters are not identifiable, you can use the reparametrize_global function.

reparam = reparametrize_global(model)

(new_ode = X3'(t) = (-1//2*X2(t)*a6*u1(t) + 1//2*X2(t)*a6 + 1//2*X2(t)*a7*u1(t) + 1//2*X2(t)*a7 - 1//2*a1*a6*u1(t) + 1//2*a1*a6 + 1//2*a2*a6*u1(t) + 1//2*a2*a6)//a6
X1'(t) = 1//2*a4*u1(t) - 1//2*a4 + 1//2*a5*u1(t) + 1//2*a5
X2'(t) = 1//2*X2(t)*a3*u1(t) - 1//2*X2(t)*a3 - 1//2*X2(t)*a5*a8*u1(t) - 1//2*X2(t)*a5*a8 + 1//2*a5*a6*u1(t) + 1//2*a5*a6
y2(t) = X1(t)
y1(t) = X3(t)
, new_vars = Dict{Nemo.QQMPolyRingElem, AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}(a6 => p[3]*p[8], X3(t) => x1, a2 => p[7], X1(t) => x3, a4 => p[2], y1(t) => y1, a1 => p[9], a7 => p[3]*p[6], u1(t) => u, y2(t) => y2…), implicit_relations = Nemo.QQMPolyRingElem[])

To print the new proposed model, you can use the new_ode field of the reparam dictionary.

reparam[:new_ode]

X3'(t) = (-1//2*X2(t)*a6*u1(t) + 1//2*X2(t)*a6 + 1//2*X2(t)*a7*u1(t) + 1//2*X2(t)*a7 - 1//2*a1*a6*u1(t) + 1//2*a1*a6 + 1//2*a2*a6*u1(t) + 1//2*a2*a6)//a6
X1'(t) = 1//2*a4*u1(t) - 1//2*a4 + 1//2*a5*u1(t) + 1//2*a5
X2'(t) = 1//2*X2(t)*a3*u1(t) - 1//2*X2(t)*a3 - 1//2*X2(t)*a5*a8*u1(t) - 1//2*X2(t)*a5*a8 + 1//2*a5*a6*u1(t) + 1//2*a5*a6
y2(t) = X1(t)
y1(t) = X3(t)

In the same way, you can use the new_vars field of the reparam dictionary to prind the relation between odd and new parameters and states.

reparam[:new_vars]

Dict{Nemo.QQMPolyRingElem, AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}} with 14 entries:
  a6    => p[3]*p[8]
  X3(t) => x1
  a2    => p[7]
  X1(t) => x3
  a4    => p[2]
  y1(t) => y1
  a1    => p[9]
  a7    => p[3]*p[6]
  u1(t) => u
  y2(t) => y2
  a5    => p[1]
  a8    => p[3]*p[4]
  X2(t) => x2*p[8]
  a3    => p[5]

You can easily check if the new model is identifiable

assess_identifiability(reparam[:new_ode])

OrderedCollections.OrderedDict{Any, Symbol} with 11 entries:
  X3(t) => :globally
  X1(t) => :globally
  X2(t) => :globally
  a1    => :globally
  a2    => :globally
  a3    => :globally
  a4    => :globally
  a5    => :globally
  a6    => :globally
  a7    => :globally
  a8    => :globally