we have an ODE system as following

   S'[t] = \[Mu]*N - \[Beta]*S[t]*I[t] - \[Mu]*S[t];
   I'[t] = \[Beta]*S[t]*I[t] - (\[Mu] - \[Gamma])*I[t];
   R'[t] = \[Gamma]*I[t] - \[Mu]*R[t];
   y[t] = k*I[t];

I would like to find the characteristics of this system, I use Groebner basis to do that, but sometimes my answer is different with Ritt-s pseudo-division method.

Is there any function in Mathematica that I could find characteristics of a system?

  • 1
    $\begingroup$ (1) What are "characteristics" for a system of ODEs? I am not familiar with the term here. (2) Is the last equation meant to involve y'[t]? (3) For purposes at hand I am pretty sure you want to use Equal infix, which is ==, rather than the single = which is infix for Set. $\endgroup$ – Daniel Lichtblau Apr 22 '18 at 14:14
  • $\begingroup$ Also, you cannot use I as a function name since I == Sqrt[-1]. Also N is a built-in function and cannot be used as a variable. $\endgroup$ – Bob Hanlon Apr 22 '18 at 17:14
  • $\begingroup$ @DanielLichtblau A lowest rank autoreduced set that can be formed with polynomials from a given set S of differential polynomials, is called a characteristic set of S. $\endgroup$ – ali Apr 22 '18 at 22:43
  • $\begingroup$ @DanielLichtblau Actually, I am looking for something like this library.wolfram.com/infocenter/MathSource/5716 $\endgroup$ – ali Apr 23 '18 at 6:16
  • $\begingroup$ My limited understanding is that one has observable input and output variables, nonobservable state variables, and parameters. It is impossible to tell from the ODE system in this post what is what, so it is not obvious which variables are to be retained and which to eliminate. $\endgroup$ – Daniel Lichtblau Apr 23 '18 at 17:34

This is a bit of a wild guess, based on two articles sent by the poster (references below).

First rewrite the system using lower case names, making the N into a time-dependent "input" variable, and calling Y, renamed to yy, the "output" variable. All other time-dependent variables are regarded as non-observable state variables.

odes = {-ss'[t] + mu*nn[t] - beta*ss[t]*ii[t] - mu*ss[t],
   -ii'[t] + beta*ss[t]*ii[t] - (mu - gamma)*ii[t],
   -rr'[t] + gamma*ii[t] - mu*rr[t],
   -yy[t] + k*ii[t]};

We prolong the system by taking two derivatives. I did not try in any way to optimize on what gets differentiated, I simply take two derivatives of all inputs. To get some idea of how one might determine the actual number needed (other than by trial and error), I gave a quick-and-dirty counting approach here. But I am not sure it applies perfectly in the situation we have now.

dodes = D[odes, t];
ddodes = D[dodes, t];

Now we gather our full system of differential polynomials, and separate input/output from internal state variables from parameters.

diffpolys = Join[odes, dodes, ddodes];
allsyms = Variables[diffpolys];
params = Select[allsyms, FreeQ[#, t] &];
tvars = Select[allsyms, ! FreeQ[#, t] &];
inoutvars = Select[tvars, ! FreeQ[#, yy | nn] &];
statevars = Complement[tvars, inoutvars];

We use a Groebner basis rather than Ritt-Wu pseudodivision in order to eliminate the state variables and their derivatives.

 gb1 = GroebnerBasis[diffpolys, inoutvars, statevars, 
   CoefficientDomain -> RationalFunctions, 
   MonomialOrder -> EliminationOrder]]

(* Out[232]= {0.683641, {(gamma k mu - k mu^2) yy[t]^2 + 
   beta k mu nn[t] yy[t]^2 + (beta gamma - beta mu) yy[t]^3 - 
   k mu yy[t] Derivative[1][yy][t] - 
   beta yy[t]^2 Derivative[1][yy][t] + k Derivative[1][yy][t]^2 - 
   k yy[t] (yy^\[Prime]\[Prime])[t]}} *)

For assessing identifiability the idea is as follows. Is it possible to get the same differential polynomial(s) in the input/output variables, using more than one set of parameters? In other words, do the coefficients of the i/o variables uniquely determine the parameters? This can be answered in the affirmative in this case, simply by gathering those coefficients and determining how many solutions there are to a "random" substitution of values. If the solution is unique 9that is, only solution is exactly the random substitution), then the system is globally identifiable.

We gather all coefficients using an internal function (could be done with System context functions but this is a shortcut).

coeffs = 
 Flatten[GroebnerBasis`DistributedTermsList[gb, inoutvars, 
    CoefficientDomain -> RationalFunctions][[1, All, All, 2]]]

(* Out[206]= {beta k mu, beta gamma - beta mu, -beta, 
 gamma k mu - k mu^2, -k mu, -k, k} *)

Now see if there are multiple solutions if given a (sorta) "random" substitution of values.

 coeffs == (coeffs /. Thread[params -> {3, -7, -2, 11}]), params]

(* Out[241]= {{beta -> 3, gamma -> -7, k -> -2, mu -> 11}} *)

Unique solution, so globally identifiable. This of course is subject to the choices made above as to how to separate into input, output, state and parameter variables.

References: (1) Nicolette Meshkat, Chris Anderson, and Joseph DiStefano Jr. Alternative to Ritt's pseudodivision for finding input-output equations in algebraic structural identifiability analysis. Math Biosci. 2012 Sep;239(1):117-23. doi: 10.1016/j.mbs.2012.04.008. Epub 2012 May 22.



(2) Giuseppina Bellu, Maria Pia Saccomani, Stefania Audoly, and Leontina D’Angiò. DAISY: a new software tool to test global identifiability of biological and physiological systems. Comput Methods Programs Biomed. 2007 Oct;88(1):52-61. Epub 2007 Aug 20.



| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.