#============================================================
#
#   Given a reaction scheme the following Maple procedures determine the 
#   system of differential equations, the associated conservation laws, and 
#   some of the species that have a zero steady state.
# 
#   This program works for chemical reaction schemes containing unary, 
#   binary or tertiary reactions.
# 
#   The ideas underlying the algorithms used here can be found in the Maple
#   Technical Newsletter article 'An application of Maple to chemical
#   kinetics,' by M. Holmes and J. Bell (number 7) and the article 'The 
#   application of symbolic computing to chemical kinetic reaction schemes,'
#   by the same authors in the Journal of Computational Chemistry (Dec, 1991).
#
#   Thanks to Jon Bell for helping develop the theory on which the algorithms
#   are based, to Aditya Nath for helping code and check parts of these 
#   procedures, and to Michael Monagan and Greg Fee for their suggestions.
#  
#		by 
#  
#		Mark H. Holmes
#		Dept of Mathematical Sciences, RPI, Troy, NY  12180
#		holmem@rpi.edu
#   
#		last changed 1/12/92
#
#============================================================

#============================================================
#                on-line help pages
#============================================================
 
`help/text/kinetics` := TEXT(
`FUNCTION: kinetics - procedures for analyzing chemical reaction schemes`,
``,
`SYNOPSIS: - The procedures included in this file are:`,
`  1. odes() - this procedure finds the ODE's from a list of reactions and`,
`              it must be called before the other two procedures`,
`  2. laws() - this procedure finds the independent conservation laws`,
`  3. steady() - procedure for determining some of the species which have a`,
`                zero steady state`,
``,
`COMMENTS:`,
`   1) There are no limitations on the number of reactions (other`,
`      than memory)`,
`   2) The letter  R  is reserved for the list of reactions. Others reserved`,
`      by the program are:`,
`      i) species[1], species[2], species[3], ...`,
`             [these are the species in the reactions]`,
`      ii) ode.1, ode.2, ode.3, ...`,
`             [these are the ODE's]`,
`      iii) law.1, law.2, law.3, ...`,
`             [these are the conservation laws]`,
`      iv) kappa0, kappa1, kappa2, kappa3, kappa4, kappa5`,
`      v) nreacts, nvars, nzeros, nlaws`,
`   3) As a reminder, Maple reserves the symbols E, I, O, and W. Others`,
`      are listed under initially-known functions and initially-known names.`,
`   4) There are certain restrictions on the individual reactions and the`,
`      Maple Newsletter article 'An application of Maple to chemical`,
`      kinetics' by M. Holmes and J. Bell (number 7) should be consulted`,
`      for specifics.`,
``,
`SEE ALSO:  odes, laws, steady`
):

`help/text/odes` := TEXT(
`FUNCTION: odes - procedure for finding the ODE's from a list of reactions`,
``,
`CALLING SEQUENCE:`,
`   odes();`,
``,
`SYNOPSIS: - this procedure finds the ODE's from a list of reactions.`,
`    It does this in four steps:`,
`    1) First determine the species in the reactions; these are denoted as:`,
`         species[i]   for  i = 1, ..., nvars`,
`    2) Determine the mapping matrix M[nreacts, nvars] where, nreacts is the`,
`       number of reactions and  nvars  is the number of variables`,
`    3) Determine the overall rate for each reaction.`,
`    4) Determine the ODE for each species. These are denoted as:`,
`         ode.i   for  i = 1, ..., nvars`,
` -  Individual differential equations can be printed using the command`,
`    print(ode.i)  where  i  is a specified integer (i = 1, 2, ..., nvars)`,
` -  The ith species can be displayed using the command  print(species[i])`,
` -  The input for this procedure is simply a list of reactions  R  in`,
`    the form`,
`	           R := [ a + b > c, b > d + w, ...];`,
``,
`EXAMPLES:`,
``,
`R:= [S + E > C, C > S + E, C > P + E];`,
`odes();`,
``,
`R:= [O2+hv>2*Ox, Ox+O2+M>O3+M, O3+hv>O2+Ox, O3+Ox>2*O2];`,
`odes();`,
``,
`SEE ALSO:  kinetics, laws, steady`
): 

`help/text/laws` := TEXT(
`FUNCTION: laws - procedure for finding the conservation laws`,
``,
`CALLING SEQUENCE:`,
`   laws();`,
``,
`SYNOPSIS: - This procedure finds the independent conservation laws of`,
`     a reaction scheme using the results obtained from the procedure odes().`,
`     The number of independent conservation laws is denoted as  nlaws`,
`     and the individual laws are denoted as:`,
`		law.i   for  i = 1, ..., nlaws`,
`     This procedure finds the laws that are independent of the rate`,
`     constants.`,
`  -  Individual conservation laws can be printed using the command`,
`     print(law.i)  where  i = 1, 2, ..., nlaws.`,
``,
`EXAMPLES:`,
``,
`R:= [S + E > C, C > S + E, C > P + E];`,
`odes();`,
`laws();`,
``,
`R:= [O2+hv>2*Ox, Ox+O2+M>O3+M, O3+hv>O2+Ox, O3+Ox>2*O2];`,
`odes();`,
`laws();`,
`print(law.2);`,
``,
`SEE ALSO:  kinetics, odes, steady`
): 

`help/text/steady` := TEXT(
`FUNCTION: steady - procedure for finding species with a zero steady state`,
``,
`CALLING SEQUENCE:`,
`   steady();`,
``,
`SYNOPSIS:`,
`   -  This determines some of the species which have a zero steady state.`,
`      The number of such species is denoted as  nzeros  and the species are`,
`  	     species[zeros[i]]   for  i = 1, ..., nzeros`,
`   -  Individual species can be printed using the command`,
`      print(species[zeros[i]])  where  i = 1, 2, ..., nzeros.`,
`   -  This procedure uses the results of the procedure odes()`,
``,
`EXAMPLES:`,
``,
`R:= [S + E > C, C > S + E, C > P + E];`,
`odes();`,
`steady();`,
`print(species[zeros[1]]);`,
``,
`R:= [O2+hv>2*Ox, Ox+O2+M>O3+M, O3+hv>O2+Ox, O3+Ox>2*O2];`,
`odes();`,
`steady();`,
``,
`SEE ALSO:  kinetics, odes, laws`
): 

#============================================================
#                procedure odes
#============================================================

macro( var=`kinetics/var`, M=`kinetics/M`, reduce=`kinetics/reduce`,
       Mentry=`kinetics/Mentry`  );

odes:=proc()
local i,j,k,l,ll,part,pp,r,rr,sparse,t:
options `Copyright 1992 by Mark H. Holmes`;

#
# STEP 1:  Determine the variables (i.e., the species in the reaction scheme)
#
nreacts:=nops(R):

# Determine species[1] from the first reaction.

nvars:=1:
l:=op(2,R[1]):
if nops(l)=1 then
    species[1]:=l 
elif nops(l)=2 then
    ll:=op(1,l):
    if  (ll=2) or (ll=3)  then
        species[1]:=op(2,l)  
    elif nops(ll)=2 then
        species[1]:=op(2,ll)
    else
        species[1]:=ll
    fi:
elif nops(l)=3 then
    species[1]:=op(1,l)
fi:
#
# Check the LHS & RHS for new variables.
# For every reaction, each of the reactants and products are
# singled out and are sent to procedure  var()  to update the
# list of variables.
#
# loop over each of the reactions and loop over each side of a reaction
#
for i from 1 to nreacts do
    for pp from 1 to 2 do
        l:=op(pp,R[i]):
        if nops(l)=1 then
            var(l)  
        elif nops(l)=2 then
            ll:=op(1,l): 
            if (ll=2) or (ll=3) then
                var(op(2,l))
            else
                for part from 1 to 2 do
                    ll:=op(part,l):
                    if nops(ll)=2 then
                        var(op(2,ll))
                    else
                        var(op(1,ll))
                    fi:
                od:
            fi:
        elif nops(l)=3 then
            for k to 3 do
                var(op(k,l))
            od:
        fi:
    od:
od:

# STEP 2:  Construct the Mapping Matrix  M[nreacts, nvars]
# Determine which species are on the LHS of the reactions
# (i.e., the reactants) and temporarily fill the  M  matrix with the
# coefficient corresponding to that species.

M:=array(sparse,1..nreacts,1..nvars):
for i from 1 to nreacts do
    l:=op(2,R[i]):
    if nops(l)=1 then
        for j to nvars do
            if l=species[j] then break fi 
        od:
        M[i,j]:=-1:
    elif nops(l)=2 then
        ll:=op(1,l):
        if (ll=2) or (ll=3) then
            for j to nvars do
                if op(2,l)=species[j] then break fi
            od:
            M[i,j]:=-1*ll:
        else
            for part from 1 to 2 do
                ll:=op(part,l):
                if nops(ll)=2 then
                    for j to nvars do
                        if op(2,ll)=species[j] then break fi
                    od:
                    M[i,j]:=-2:
                else
                    for j to nvars do
                        if ll=species[j] then break fi
                    od:
                    M[i,j]:=-1:
                fi:
            od:
        fi:
    else
        for k to 3 do
            for j to nvars do
                if op(k,l)=species[j] then break fi
            od:
            M[i,j]:=-1:
        od:
    fi:
od:

# Check the variables on the RHS of reactions (the products) and
# fill up the  M  matrix keeping the entries due to the reactants
# in consideration. The procedure  'Mentry'  is used to place the
# appropriate value in the  M  matrix.

for i from 1 to nreacts do
    r:=op(1,R[i]):
    if nops(r)=1 then
        for j to nvars do
            if r=species[j] then break fi:
        od:
        Mentry(1,i,j):
    elif nops(r)=2 then
        rr:=op(1,r):
        if (rr=2) or (rr=3) then
            for j to nvars do
                if op(2,r)=species[j] then break fi
            od:
            if rr=2 then
                    Mentry(2,i,j)
                else
                    Mentry(3,i,j)
            fi:
        else 
            for part from 1 to 2 do
                rr:=op(part,r):
                if nops(rr)=2 then 
                    for j to nvars do
                        if op(2,rr)=species[j] then break fi
                    od:
                    Mentry(2,i,j):
                else
                    for j to nvars do
                        if rr=species[j] then break fi
                    od:
                    Mentry(1,i,j):
                fi:
            od:
        fi:
    elif nops(r)=3 then
        for k to 3 do
            for j to nvars do
                if op(k,r)=species[j] then break fi
            od:
            Mentry(1,i,j):
        od:
    fi:
od:
#
# STEP 3:  Determine the rate of each reaction;
# the rate constant for the ith reaction is:  k.i 
#
for i from 1 to nreacts do
    r.i:=k.i:
    for j from 1 to nvars do
        if M[i,j]=-1 or M[i,j]=kappa0 or
            M[i,j]=kappa1 or M[i,j]=kappa2 then
                r.i:=r.i * species[j]
        elif M[i,j]=-2 or M[i,j]=-kappa1 or
            M[i,j]=kappa4 or M[i,j]=kappa3 then
                r.i:=r.i * (species[j]**2)
        elif M[i,j]=-3 or M[i,j]=-kappa2 or
            M[i,j]=-kappa3 or M[i,j]=kappa5 then
                r.i:=r.i * (species[j]**3)
        fi:
    od:
od:
#
# STEP 4: Determine the ODE for Each Species
# The ODE for the species is computed based upon the entries of the
# mapping matrix and the rate of the reactions.
#
# loop over all the species and over all the reactions
#
for j from 1 to nvars do
    rhs.j:=0:
    for i from 1 to nreacts do
        if M[i,j]=-1 or M[i,j]=-kappa1
            or M[i,j]=-kappa3 then
                rhs.j:= rhs.j - r.i
        elif M[i,j]=-2 or M[i,j]=-kappa2  then
                rhs.j:= rhs.j - r.i*2
        elif M[i,j]=1 or M[i,j]=kappa1
            or M[i,j]=kappa3 then
                rhs.j:= rhs.j + r.i
        elif M[i,j]=2 or M[i,j]=kappa2 then
            rhs.j:= rhs.j + r.i*2
        elif M[i,j]=3 then
            rhs.j:= rhs.j + r.i*3
        elif M[i,j]=-3 then
            rhs.j:= rhs.j - r.i*3
        else 
            rhs.j:= rhs.j + 0
        fi:
    od:
od:
#
for i from 1 to nvars do
    ode.i:=diff(species[i](t),t)=rhs.i:
od:
#
# print out the differential equations and the CPU time taken
# to compute the equations
#
if printlevel > 0 then
print();
lprint(`The differential equations are:`);
print();
for i from 1 to nvars do
    print(ode.i)
od:
print();
fi;
#
# end of procedure odes()
#
end:
 
# The input to this procedure is a species name.
# The existing list of species names is checked and if the species
# does not appear in the list then it is added on.

var:=proc(side)
local it,j:
options `Copyright 1992 by Mark H. Holmes`;
it:=0:
for j to nvars do
    if side=species[j] then it:=1 fi
od:
if it=0 then
    nvars:=nvars + 1:
    species[nvars]:=side:
fi:
end:

# Procedure ' Mentry'  determines the entries of the mapping matrix.
# Before this procedure is called, the mapping matrix has entries
# corresponding to the LHS of the reactions only.  Based upon these
# entries and the species of the RHS, new values are entered into the
# mapping matrix.
# The input parameters to this procedure are the position of the
# mapping matrix entry (i and j) and the coefficient of the 
# species on RHS.

Mentry:=proc(x,i,j)
options `Copyright 1992 by Mark H. Holmes`;
if x=1 then
    if M[i,j]=-1 then M[i,j]:=kappa0
    elif M[i,j]=-2 then M[i,j]:=-kappa1
    elif M[i,j]=-3 then M[i,j]:=-kappa2
    else M[i,j]:=M[i,j] + x
    fi:
elif x=2 then
    if M[i,j]=-1 then M[i,j]:=kappa1
    elif M[i,j]=-2 then M[i,j]:=kappa4
    elif M[i,j]=-3 then M[i,j]:=-kappa3
    else M[i,j]:=M[i,j] + x
    fi:
elif x=3 then
    if M[i,j]=-1 then M[i,j]:=kappa2
    elif M[i,j]=-2 then M[i,j]:=kappa3
    elif M[i,j]=-3 then M[i,j]:=kappa5
    else M[i,j]:=M[i,j] + x
    fi:
fi:
end:
 
#============================================================
#             procedure laws
#============================================================

# This procedure computes the conservation laws based upon the 
# properties of the mapping matrix.
# Use is made of the linear algebra functions of Maple to compute
# the kernel of the mapping matrix.

laws:=proc()
local i,ii,kk,jj,vec:
options `Copyright 1992 by Mark H. Holmes`;

# The symbols are replaced by their corresponding stoichiometric number.

kappa0:=0:kappa1:=1:kappa4:=0:
kappa2:=2:kappa3:=1:kappa5:=0:
kk:=linalg[kernel](M,'nlaws'):

for i to nlaws do
    vec:=kk[i]:
    law.i:=0:
    for jj to nvars do
        law.i:=law.i+vec[jj]*species[jj]
    od:
od:

# Having finished with the linear algebra operations on M
# those entries which were replaced by numbers earlier are
# replaced with symbols again

kappa0:='kappa0':kappa1:='kappa1':kappa2:='kappa2':
kappa4:='kappa4':kappa3:='kappa3':kappa5:='kappa5':

# printout the conservation laws and the CPU time elapsed

if printlevel > 0 then
print();
if nlaws=0 then
    lprint(` There are no conservation laws.`)
else
    lprint(` The conservation laws are:`)
fi:
print();
for i to nlaws do
    ii:=i*1.:
    lprint(`   `,ii,`  constant = `,law.i);
od:
print();
fi;

# end of procedure laws()

end:
 

#============================================================
#             procedure steady
#============================================================

steady:=proc()
local i,j,ii,jj,i0,i1,i2,i3,j0,j1,j2,j3,jr,
in0,in1,MM,nzeros0,nzeros1:
options `Copyright 1992 by Mark H. Holmes`;
 
nzeros:=0:
zeros:=array(1..nvars):
 
#  MM is a working matrix identifying which species are
#    consumed in the scheme
kappa0:=11:kappa1:=11:kappa2:=11:kappa3:=11:kappa4:=11:kappa5:=11:
MM:=array(1..nreacts,1..nvars):
for i to nreacts do
    for j to nvars do
        MM[i,j]:=M[i,j]:
        if MM[i,j]=11 or MM[i,j]=-11 then  MM[i,j]:=-1 fi:
    od:
od:
kappa0:='kappa0':kappa1:='kappa1':kappa2:='kappa2':
kappa4:='kappa4':kappa3:='kappa3':kappa5:='kappa5':
 
#  apply conditions (1) and (2) repeatedly until  nzero 
#  doesn't change
 
for jr to nvars do
    nzeros0:=nzeros:
# apply condition (1)
    for j0 to nvars do
        i1:=0:
        i2:=0:
        for i to nreacts do
            if MM[i,j0]>0 then i1:=1 fi: 
            if MM[i,j0]<0 then i2:=1 fi:
        od:
        if i2=0 then next elif i1=1 then next fi:
        for i0 to nreacts do
            if MM[i0,j0]>=0 then next fi:
            i3:=0:
            for jj to nvars do
                if jj=j0 then next
                elif MM[i0,jj]<0 then i3:=1
                fi:
            od:
            if i3=1 then next else 
                in0:=0:
            for ii to nzeros do
                if j0=zeros[ii] then in0:=1 fi
            od:
            if in0=0 then nzeros:=nzeros+1: zeros[nzeros]:=j0 fi:
            fi:
        od:
    od:
 
# remove any species with a zero steady state using procedure reduce

    if nzeros>nzeros0 then reduce(MM) fi:
    nzeros1:=nzeros:
 
# apply condition (2)
    for j1 to nvars do
        j2:=0:
        i1:=0:
        for i to nreacts do
            if MM[i,j1]<0 then j2:=1 fi:
            if MM[i,j1]>0 then i1:=1 fi:
        od:
        if i1=0 or j2=1 then next fi:
        for i0 to nreacts do
            if MM[i0,j1]<=0 then next fi:        
            for j0 to nvars do
                if MM[i0,j0]>=0 then next else j3:=0 fi:
                for jj to nvars do
                    if jj=j0 then next
                    elif MM[i0,jj]<0 then j3:=1
                    fi:
                od:
                if j3=1 then next else 
                    in1:=0:
                    for ii to nzeros do
                    if j0=zeros[ii] then in1:=1 fi
                    od:
                if in1=0 then nzeros:=nzeros+1:
                zeros[nzeros]:=j0 fi:
                fi:
            od:
        od:
    od:
 
# remove any species with a zero steady state
    if nzeros>nzeros1 then reduce(MM) fi:
 
# if the number hasn't changed after applying both conditions then exit
    if nzeros0=nzeros then break fi:
 
od:

if printlevel > 0 then
print();
if nzeros=0 then
    lprint(`No species were found that have a zero steady state`)
else
    lprint(`There are`,nzeros,
    `species that were found to have a zero steady state and they are:`);
    for j to nzeros do
        jj:=zeros[j]:
        lprint(`   `,species[jj]);
    od:
fi:
print();
fi;
 
end:

 
#  procedure to set to zero those elements in the 
#  mapping matrix associated with the zero steady states
 
reduce:=proc(MM)
local i,j,jj,jn:
options `Copyright 1992 by Mark H. Holmes`;
 
for jn to nzeros do
    jj:=zeros[jn]:
    for i to nreacts do
        if MM[i,jj]<0 then
            for j to nvars do
                MM[i,j]:=0
            od:
        fi:
    od:
od:
for jn to nzeros do
    jj:=zeros[jn]:
    for i to nreacts do
        MM[i,jj]:=0
    od:
od:
end:

macro( var=var, M=M, reduce=reduce, Mentry=Mentry );
#============================================================
#				The End
#============================================================

# save `kinetics.m`;
# quit
