DEMO OF DETERMINANT OF    JACOBIAN OF CHEM REACTION NETWORK SOFTWARE

Here,  we consider the differential equation 
      (dx/dt)=f (x)=Sv(x) which acts on the non-\(negative\ orthant\ \(R\^d\)\_\( ≥ \_0\)\), where S is a d x d' Stoichiometric Matrix,  \(and\ v \((x)\)\ is\ a\ vector\ with\ d'\ components\ consisting\ of\ monomials\ multiplied\ by\ rate\ constants . \ 
              These\ are\ called\ fluxes . 
\ 
              This\ differential\ equation\ is\ associated\ to\) 
     \(a\ chemical\ network\ modeled\ \ with\ mass - 
          action\ kinetics . 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Runs\ in\ MMa\ 4.2  and\ \(\(higher\)\(.\)\)\)

Input the Differential Equation

Input the Stoichiometric matrix S

S={{-1, 0, 2}, {-1, -1, 
          0}, {0, -1, -1}, {1, 0, 0}, {0, 1, 
          0}}

S  // MatrixForm

mr3[m_] := 
          Length[Select[
              RowReduce[
                m], !  (Dot[#, #] === 0)  &]];

      mr3[S]

Make the flux vector and list of species concentrations automatically from the Stoiciometric matrix

(* Each column of S corresponds to a chemical reaction, modeled by a monomial multiplied by a rate constant (flux) .  
In particular, writing S = ((Subscript[S, ij])),x  =  ((((Subscript[x, 1]),  ... , (Subscript[x, d])))^T),  and v (x)  =  ((((Subscript[v, 1]) (x),  ... , (Subscript[v, d']) (x)))^T), we have 
(Subscript[v, j]) (x) = 
                (Subscript[k, j]) (∏_(i = 1)^d) ((Subscript[x, i])^(|min {0, (Subscript[S, ij])}|)) .  

The function ' 
                      makeMonomial' will construct v (x), and take S as input.*)

makeMonomial[Smatrix_] := 
      Module[{tmp, lista, pow, mat}, 
        refine[S_] := S  /. Map[# -> 1 &, Variables[S]]; 
        mat = refine[Smatrix]; 
pow = 
          Table[0, {Length[mat[[1]]]}]; 
        tmp = Table[
            If[mat[[i, j]] ≥ 0, 1, 
              If[pow[[j]] > 0, 
                a[i]^(-Smatrix[[i, j]]), pow[[j]]++;  
                k[j] a[i]^(-Smatrix[[i, j]])]], {i, 1, 
              Length[mat]}, {j, 1, 
              Length[mat[[1]]]}]; 
lista = 
          Map[Apply[Times, #] &, Transpose[tmp]];  lista]

makeMonomialNMA[
        Smatrix_] := 
Module[{reactantlist, tmp}, 
        tmp = Smatrix;  tmp = Map[If[# < 0, 1, 0] &, tmp, {2}]; 
        tmp = Transpose[tmp]; 
        reactantlist = 
          Table[tmp[[i, j]]*a[j], {i, 
              1, Dimensions[tmp][[1]]}, {j, 
              1, Dimensions[
                  tmp][[2]]}]; 
reactantlist = 
          DeleteCases[reactantlist, 0, {2}]; 
Table[
          Apply[k[i], reactantlist[[i]]], {i, 1, 
            Length[reactantlist]}]]

(* 
        The function ' 
          variables' takes a vector of monomials (fluxes) as input,  and returns a vector of all the variables (species) appearing in the vector of monomials .   WARNING :  
              if the stoichiometric matrix is not fully reversible,  then not all species may appear as variables in the monomial vector,  
        so this command may return too few variables. *)

variables[monom_] := 
      Module[{list},   list = 
            Select[Variables[monom], 
              Head[#] == 
                    Head[a[1]] ||  (Head[#] == Power  &&  
                      Head[Level[#, 1][[1]]] == 
                        Head[a[1]]) &];  
Table[ 
            If[ Length[list[[i]]] == 1,  
              list[[i]], ( Level[list[[i]], 1])[[1]] ], {i, 
              1, Length[list]}]] // Union

svars[Smatrix_] := Table[a[i], {i, 1, Length[Smatrix]}];

makeReversible[S_] := 
      Module[{tmp}, tmp = Transpose[S];  
        Flatten[ 
Table[ {tmp[[i]], -tmp[[i]]}, {i, 1, Length[tmp]}], 1] // Transpose]

rS = makeReversible[S]

monomialvec  =  makeMonomial[rS]; 
    monomialvec // MatrixForm

vars  =  variables[monomialvec]

Right - hand side of the ODE

fns = rS . monomialvec;
fns // MatrixForm

jac[Fns_, Vars_] := 
        Table[D[Fns[[i]], Vars[[j]]], {i, 1, Length[Fns]}, {j, 1, 
            Length[Vars]}]

Computations start here

The Jacobian of the RHS of the ODE and their Determinants

J = jac[fns, vars];

    J // MatrixForm

cfDet[Jac_] :=  
        Det[Jac - IdentityMatrix[Length[Jac]]] // Expand;

cf = cfDet[J]

coeffsList[expansion_] := 
      DeleteCases[ 
        Flatten[CoefficientList[expansion,  Variables[expansion]]], 
        0]

coeffs[expansion_] := 
      Module[{tmp}, 
        tmp = coeffsList[
            expansion]; 
Print["The number of terms in \
the det expansion is ", 
          Length[tmp], ","]; 
Print["and (a,b) says \
that the number of terms with coeff. a is b:"]; 
Map[{#, 
              Count[tmp, #]} &, Union[tmp]]]

coeffs[cf]

The Core Determinant

The core determinant   is the determinant of  the Jacobian, fns'(x), restricted to range S.
It is natural for analyzing chem reactions with no outflows, (see Helton-Klep-Gomez).

coreDet[Jac_,  Smatrix_]  := 
        Module[{d,  cdet}, 
          d =  Length[Jac]  -  
              mr3[Smatrix]; 
cdet  =  (Det[
                    Jac -  t IdentityMatrix[Length[Jac]]])/t^d // 
              Expand; 
cdet  /. t -> 0]

core  =  coreDet[J,  rS]

coeffs[core]