(* ::Package:: *)

(* Mathematica Package *)

(* Created by the Wolfram Workbench July 22, 2011 *)

BeginPackage["RedoxMech`"]

(* Main exported symbols and their usage *) 

RedoxComplex::usage = "RedoxComplex[{Dre, Dox}, {Are, Aox}, EnzymeName, paramList, Options] 
The main function of RedoxMech that produces the prototypical internal mechanism of an oxidoreductase complex in terms of 
xCellerator system of equations (Default Option: Output->xCellerator) or a reaction matrix map of all transitions between the microstates (Option: Output->Matrix).";

SSVelocityEq::usage = "SSVelocityEq[inputMatrix, Options] Determine the steady-state velocity equation from the reaction matrix";

ReactionMatrix::usage = "ReactionMatrix[inputMatrix, Options] Display the reaction table from the reaction matrix output of the EC1/RedoxComplex function";

(* Support Functions *)

ToMolecules::usage = "ToMolecules[uM] Convert concentration in uM to number of molecules based on the volume (global variable) of the system.";

ToConcentration::usage = "ToConcentration[molecules] Convert from number of molecules to concentration in uM based on the volume (global variable) of the system.";

TwoReactantsReaction::usage = "TwoReactantsReaction[species list1, species list2, reaction label, options] 
Automates the creation of all permutations of the forward and reverse binding reaction between two molecules; each molecule may contain multiple species encoded in a list. 
Each reaction is identified by the reaction label and an unique index number to that reaction label. User may specify an alternate indexing scheme with the ReverseIndex->True/False option,
which determines if the index is incremented by advancing the species list1 or species list2";

ReactionTable::usage = "ReactionTable[(TwoReactantsReaction) reaction list] Display the binding reaction list output of TwoReactantsReaction in a tabular format"; 

TransferTable::usage = "TransferTable[reaction list] Display the electron transfer reaction list in a tabular format";

PickReactions::usage = "PickReactions[bindingList, reactionList] Find the unique index for both forward and reverse elementary reactions from a list of binding reactions";

UnitaryConstant::usage = "UnitaryConstant[reactions list, first specie, second specie] 
Identifies the unitary constant of a particular reaction by the labels of the reactant species
";

TotalSubstrate::usage = "TotalSubstrate[simulation, substrate, time] Sum of all forms (bound or unbound) of a particular substrate in the system";

TotalSubstrateD::usage = "TotalSubstrateD[simulation, substrate, time]";

(* Options *)

ReverseIndex::usage = "Option for TwoReactantsReaction which determines if the index is incremented by advancing the species list1 or species list2";

ReactionType::usage = "Option for EC1/RedoxComplex. Choose between precoded mechanisms of EC1. Available selections are: Complete, PingPong, Ordered, and Random";

RemoveList::usage = "Option for EC1/RedoxComplex. Remove specific internal elementary reaction(s) from the prototypical EC1 representation";

Output::usage = "Option for EC1/RedoxComplex. Choose between the default xCellerator output or the Matrix output";

Report::usage = "Option for SSVelocityEq. A True value produces a brief summary of the enzyme distribution equations. Default: False";

ShowConstants::usage = "Option for SSVelocityEq and ReactionMatrix to switch the kinetic constants between unitary, steady-state, or numerical values. Default: Unitary";

(*xCellerator arrow-form function calls to RedoxComplex function*)
EC1[

\!\(\*OverscriptBox[\(S_List \[RightArrowLeftArrow] P_List\), \(EnzymeName_\)]\), paramList_, opts___Rule] :=
    RedoxComplex[S, P, EnzymeName, paramList, opts]

EC1[

\!\(\*OverscriptBox[\(S_List \[RightArrowLeftArrow] P_List\), \(EnzymeName_\)]\), paramList_, reactionType_String, opts___Rule] :=
    RedoxComplex[S, P, EnzymeName, paramList, ReactionType->reactionType, opts]
 
Options[RedoxComplex] = {ReactionType -> "Complete", RemoveList -> {}, Output -> "xCellerator"};
Options[ReactionMatrix] = {ShowConstants->"Unitary"};
Options[SSVelocityEq] = {Report->False, ShowConstants->"Unitary"};

Begin["`Private`"]
(* Implementation of the package *)

ToConcentration[molecules_] :=
    molecules 10^6/(Global`Volume Global`nA);

ToMolecules[uM_] :=
    Round[uM (Global`Volume Global`nA)/10^6];
    
UnitaryConstant[reactions_, first_, second_] :=
    Module[ {c, pos, length1, length2, exp},
        length1 = Length[first];
        length2 = Length[second];
        Which[
         length2 > length1,
         pos = Position[reactions, _ + first -> second];
         If[ pos == {},
             Return[0]
         ];
         c = DeleteCases[reactions[[pos[[1, 1]], 1, 1]] /. {Plus -> List}, 
            first][[1]];
         ,
         length1 > length2,
         pos = Position[reactions, first -> second + _];
         If[ pos == {},
             Return[0]
         ];
         c = 1;
         ,
         length1 == length2,
         pos = Position[reactions, first -> second];
         If[ pos == {},
             Return[0]
         ];
         c = 1;
         ];
        exp = c reactions[[pos[[1, 1]]]][[2]];
        Return[exp];
    ];

TwoReactantsReaction[list1_, list2_, type_, opts___Rule] :=
    Module[ {length1 = Length[list1], length2 = Length[list2], 
      reactionList = {}, newProductList = {}, ty = type, index, 
      reactant1, reactant2, newProduct, newForward, 
      newReverse},
        index = 1;
        If[ ReverseIndex /. {opts} /. {ReverseIndex -> False},
            For[j = 1, j <= length2, j++,
              For[i = 1, i <= length1, i++,
               reactant1 = list1[[i]];
               reactant2 = list2[[j]];
               newProduct = Diamond[reactant1, reactant2];
               If[ Length[reactant1] > 1,
                   newProduct = Append[reactant1, reactant2],
                   If[ Length[reactant2] > 1,
                       newProduct = Prepend[reactant2, reactant1]
                   ]
               ];
               newForward = {reactant1 + reactant2 -> newProduct, 
                 Subscript[Global`kf, ty][index]};
               newReverse = {newProduct -> reactant1 + reactant2, 
                 Subscript[Global`kr, ty][index]};
               reactionList = Append[reactionList, {newForward, newReverse}];
               newProductList = Append[newProductList, newProduct];
               index++;
               ]
              ];,
            For[i = 1, i <= length1, i++,
              For[j = 1, j <= length2, j++,
               reactant1 = list1[[i]];
               reactant2 = list2[[j]];
               newProduct = Diamond[reactant1, reactant2];
               If[ Length[reactant1] > 1,
                   newProduct = Append[reactant1, reactant2],
                   If[ Length[reactant2] > 1,
                       newProduct = Prepend[reactant2, reactant1]
                   ]
               ];
               newForward = {reactant1 + reactant2 -> newProduct, 
                 Subscript[Global`kf, ty][index]};
               newReverse = {newProduct -> reactant1 + reactant2, 
                 Subscript[Global`kr, ty][index]};
               reactionList = Append[reactionList, {newForward, newReverse}];
               newProductList = Append[newProductList, newProduct];
               index++;
               ]
              ];
        ];
        Return[{reactionList, newProductList}]
    ]
  
  ReactionTable[reactL_] :=
      TableForm[reactL[[1]],
      TableHeadings -> {reactL[[2]], {"Forward", "Reverse"}}
      ]
  
  TransferTable[reactL_] :=
      TableForm[reactL[[1 ;;, 1 ;;, 2]],
      TableHeadings -> {reactL[[1 ;;, 1, 1]], {"Forward", "Reverse"}}
      ]
  
  PickReactions[bindingList_, reactionList_] :=
      Module[ {returnList},
          returnList = 
           bindingList[[Flatten[Function[x, {2 x - 1, 2 x}] /@ reactionList]]];
          returnList
      ]
  
  TotalSubstrate[sim_, substrate_, time_] :=
      Module[ {subList, positions},
          subList = sim[[1, 1 ;;, 1]];
          positions = Position[subList, substrate][[1 ;;, 1]];
          Plus @@ ((subList[[#]][time]) & /@ positions) /. sim
      ]
      
  TotalSubstrateD[sim_, substrate_, time_] :=
      Module[ {subList, positions},
          subList = sim[[1, 1 ;;, 1]];
          positions = Position[subList, substrate][[1 ;;, 1]];
          Plus @@ ((subList[[#]]'[time]) & /@ positions) /. sim
      ]
      
 
      
RedoxComplex[{Dre_, Dox_}, {Are_, Aox_}, EnzymeName_, paramList_List: {}, opts___Rule] :=
    Module[ {Donor = {Dre, Dox}, Acceptor = {Aox, Are}, Epp, Enn, Epn, 
      Enp, Enz = {ToExpression[StringJoin[ToString[EnzymeName], "pp"]], 
        ToExpression[StringJoin[ToString[EnzymeName], "np"]], 
        ToExpression[StringJoin[ToString[EnzymeName], "pn"]], 
        ToExpression[StringJoin[ToString[EnzymeName], "nn"]]}, 
      BindingReactions,
      donorFirst, acceptorFirst, acceptorSecond, donorSecond, 
      TernaryTransfers, donorToEnzTri, EnzRhToRlTri, EnzToAcceptorTri, 
      BinaryTransfers, donorToEnzBi, EnzRhToRlBi, EnzToAcceptorBi, 
      UninaryTransfers, EnzRhToRlUni, BindingHeadings, donorFirstParam, 
      acceptorFirstParam, acceptorSecondParam, donorSecondParam, 
      donorToEnzTriParam, EnzRhToRlTriParam, EnzToAcceptorTriParam, 
      donorToEnzBiParam, EnzRhToRlBiParam, EnzToAcceptorBiParam, 
      EnzRhToRlUniParam, completeParam, nA = 6.022 10^23, F = 96.485, 
      R = 8.31, Emrd, Emra, Keq, removeTable, removeNull, matrixList, reactions, outputForm, removeList, reactionType},
        Epp = Enz[[1]];
        Enp = Enz[[2]];
        Epn = Enz[[3]];
        Enn = Enz[[4]];
        donorFirst = 
         TwoReactantsReaction[Donor, Enz, Global`D1, ReverseIndex -> False];
        acceptorFirst = 
         TwoReactantsReaction[Enz, Acceptor, Global`A1, ReverseIndex -> True];
        acceptorSecond = 
         TwoReactantsReaction[donorFirst[[2]], Acceptor, Global`A2, 
          ReverseIndex -> True];
        donorSecond = 
         TwoReactantsReaction[Donor, acceptorFirst[[2]], Global`D2, 
          ReverseIndex -> False];
        donorToEnzTri = {
          {{Dre\[Diamond]Epp\[Diamond]Aox -> Dox\[Diamond]Enp\[Diamond]Aox, 
            Subscript[Global`kf, Global`TD][1]}, {Dox\[Diamond]Enp\[Diamond]Aox -> 
             Dre\[Diamond]Epp\[Diamond]Aox, Subscript[Global`kr, Global`TD][1]}},
          {{Dre\[Diamond]Epn\[Diamond]Aox -> Dox\[Diamond]Enn\[Diamond]Aox, 
            Subscript[Global`kf, Global`TD][2]}, {Dox\[Diamond]Enn\[Diamond]Aox -> 
             Dre\[Diamond]Epn\[Diamond]Aox, Subscript[Global`kr, Global`TD][2]}},
          {{Dre\[Diamond]Epp\[Diamond]Are -> Dox\[Diamond]Enp\[Diamond]Are, 
            Subscript[Global`kf, Global`TD][3]}, {Dox\[Diamond]Enp\[Diamond]Are -> 
             Dre\[Diamond]Epp\[Diamond]Are, Subscript[Global`kr, Global`TD][3]}},
          {{Dre\[Diamond]Epn\[Diamond]Are -> Dox\[Diamond]Enn\[Diamond]Are, 
            Subscript[Global`kf, Global`TD][4]}, {Dox\[Diamond]Enn\[Diamond]Are -> 
             Dre\[Diamond]Epn\[Diamond]Are, Subscript[Global`kr, Global`TD][4]}}};
        EnzRhToRlTri =
         {{{Dre\[Diamond]Enp\[Diamond]Aox -> Dre\[Diamond]Epn\[Diamond]Aox, 
            Subscript[Global`kf, Global`TT][1]}, {Dre\[Diamond]Epn\[Diamond]Aox -> 
             Dre\[Diamond]Enp\[Diamond]Aox, Subscript[Global`kr, Global`TT][1]}},
          {{Dox\[Diamond]Enp\[Diamond]Aox -> Dox\[Diamond]Epn\[Diamond]Aox, 
            Subscript[Global`kf, Global`TT][2]}, {Dox\[Diamond]Epn\[Diamond]Aox -> 
             Dox\[Diamond]Enp\[Diamond]Aox, Subscript[Global`kr, Global`TT][2]}},
          {{Dre\[Diamond]Enp\[Diamond]Are -> Dre\[Diamond]Epn\[Diamond]Are, 
            Subscript[Global`kf, Global`TT][3]}, {Dre\[Diamond]Epn\[Diamond]Are -> 
             Dre\[Diamond]Enp\[Diamond]Are, Subscript[Global`kr, Global`TT][3]}},
          {{Dox\[Diamond]Enp\[Diamond]Are -> Dox\[Diamond]Epn\[Diamond]Are, 
            Subscript[Global`kf, Global`TT][4]}, {Dox\[Diamond]Epn\[Diamond]Are -> 
             Dox\[Diamond]Enp\[Diamond]Are, Subscript[Global`kr, Global`TT][4]}}};
        EnzToAcceptorTri =
         {{{Dre\[Diamond]Enn\[Diamond]Aox -> Dre\[Diamond]Enp\[Diamond]Are, 
            Subscript[Global`kf, Global`TA][1]}, {Dre\[Diamond]Enp\[Diamond]Are -> 
             Dre\[Diamond]Enn\[Diamond]Aox, Subscript[Global`kr, Global`TA][1]}},
          {{Dox\[Diamond]Enn\[Diamond]Aox -> Dox\[Diamond]Enp\[Diamond]Are, 
            Subscript[Global`kf, Global`TA][2]}, {Dox\[Diamond]Enp\[Diamond]Are -> 
             Dox\[Diamond]Enn\[Diamond]Aox, Subscript[Global`kr, Global`TA][2]}},
          {{Dre\[Diamond]Epn\[Diamond]Aox -> Dre\[Diamond]Epp\[Diamond]Are, 
            Subscript[Global`kf, Global`TA][3]}, {Dre\[Diamond]Epp\[Diamond]Are -> 
             Dre\[Diamond]Epn\[Diamond]Aox, Subscript[Global`kr, Global`TA][3]}},
          {{Dox\[Diamond]Epn\[Diamond]Aox -> Dox\[Diamond]Epp\[Diamond]Are, 
            Subscript[Global`kf, Global`TA][4]}, {Dox\[Diamond]Epp\[Diamond]Are -> 
             Dox\[Diamond]Epn\[Diamond]Aox, Subscript[Global`kr, Global`TA][4]}}};
        donorToEnzBi =
         {{{Dre\[Diamond]Epp -> Dox\[Diamond]Enp, 
            Subscript[Global`kf, Global`BD][1]}, {Dox\[Diamond]Enp -> Dre\[Diamond]Epp, 
            Subscript[Global`kr, Global`BD][1]}},
          {{Dre\[Diamond]Epn -> Dox\[Diamond]Enn, 
            Subscript[Global`kf, Global`BD][2]}, {Dox\[Diamond]Enn -> Dre\[Diamond]Epn, 
            Subscript[Global`kr, Global`BD][2]}}};
        EnzRhToRlBi =
         {{{Dre\[Diamond]Enp -> Dre\[Diamond]Epn, 
            Subscript[Global`kf, Global`BT][1]}, {Dre\[Diamond]Epn -> Dre\[Diamond]Enp, 
            Subscript[Global`kr, Global`BT][1]}},
          {{Dox\[Diamond]Enp -> Dox\[Diamond]Epn, 
            Subscript[Global`kf, Global`BT][2]}, {Dox\[Diamond]Epn -> Dox\[Diamond]Enp, 
            Subscript[Global`kr, Global`BT][2]}},
          {{Enp\[Diamond]Aox -> Epn\[Diamond]Aox, 
            Subscript[Global`kf, Global`BT][3]}, {Epn\[Diamond]Aox -> Enp\[Diamond]Aox, 
            Subscript[Global`kr, Global`BT][3]}},
          {{Enp\[Diamond]Are -> Epn\[Diamond]Are, 
            Subscript[Global`kf, Global`BT][4]}, {Epn\[Diamond]Are -> Enp\[Diamond]Are, 
            Subscript[Global`kr, Global`BT][4]}}
          };
        EnzToAcceptorBi =
         {{{Epn\[Diamond]Aox -> Epp\[Diamond]Are, 
            Subscript[Global`kf, Global`BA][1]}, {Epp\[Diamond]Are -> Epn\[Diamond]Aox, 
            Subscript[Global`kr, Global`BA][1]}},
          {{Enn\[Diamond]Aox -> Enp\[Diamond]Are, 
            Subscript[Global`kf, Global`BA][2]}, {Enp\[Diamond]Are -> Enn\[Diamond]Aox, 
            Subscript[Global`kr, Global`BA][2]}}};
        EnzRhToRlUni =
         {{{Enp -> Epn, Subscript[Global`kf, Global`UT][1]}, {Epn -> Enp, 
            Subscript[Global`kr, Global`UT][1]}}
          };
        BindingReactions = 
         Flatten[{donorFirst[[1]], acceptorFirst[[1]], acceptorSecond[[1]], 
           donorSecond[[1]]}, 2];
        BindingHeadings = Join[
          Flatten[
           Table[Flatten[{"D First (molecule/s)", 
              "  Dissociation (1/s)"}], {i, 1, 8}]],
          Flatten[
           Table[Flatten[{"A First (molecule/s)", 
              "  Dissociation (1/s)"}], {i, 1, 8}]],
          Flatten[
           Table[Flatten[{"A Second(molecule/s)", 
              "  Dissociation (1/s)"}], {i, 1, 16}]],
          Flatten[
           Table[Flatten[{"D Second(molecule/s)", 
              "  Dissociation (1/s)"}], {i, 1, 16}]]
          ];
        TernaryTransfers = 
         Flatten[{donorToEnzTri, EnzRhToRlTri, EnzToAcceptorTri}, 2];
        BinaryTransfers = 
         Flatten[{donorToEnzBi, EnzRhToRlBi, EnzToAcceptorBi}, 2];
        UninaryTransfers = Flatten[{EnzRhToRlUni}, 2];
        donorFirstParam = Flatten[Table[{
            Subscript[Global`kf, Global`D1][i] -> N[Global`BindingOrder/(nA Global`Volume)],
            Subscript[Global`kr, Global`D1][i] -> N[Global`BindingOrder Global`Ksd]}, {i, 1, 8}]];
        acceptorFirstParam = Flatten[Table[{
            Subscript[Global`kf, Global`A1][i] -> N[Global`BindingOrder/(nA Global`Volume)],
            Subscript[Global`kr, Global`A1][i] -> N[Global`BindingOrder Global`Ksa]}, {i, 1, 8}]];
        acceptorSecondParam = Flatten[Table[{
            Subscript[Global`kf, Global`A2][i] -> N[Global`BindingOrder/(nA Global`Volume)],
            Subscript[Global`kr, Global`A2][i] -> N[Global`BindingOrder Global`Ksda]}, {i, 1, 16}]];
        donorSecondParam = Flatten[Table[{
            Subscript[Global`kf, Global`D2][i] -> N[Global`BindingOrder/(nA Global`Volume)],
            Subscript[Global`kr, Global`D2][i] -> N[Global`BindingOrder Global`Ksad]}, {i, 1, 16}]];
        donorToEnzTriParam = Flatten[Table[{
            Subscript[Global`kf, Global`TD][i] -> Global`boundaryTransfer,
            Subscript[Global`kr, Global`TD][i] -> Global`boundaryTransfer Global`Krd}, {i, 1, 4}]];
        EnzRhToRlTriParam = Flatten[Table[{
            Subscript[Global`kf, Global`TT][i] -> Global`kCat,
            Subscript[Global`kr, Global`TT][i] -> N[Global`kCat/Keq]}, {i, 1, 4}]];
        EnzToAcceptorTriParam = Flatten[Table[{
            Subscript[Global`kf, Global`TA][i] -> Global`boundaryTransfer,
            Subscript[Global`kr, Global`TA][i] -> Global`boundaryTransfer Global`Kra}, {i, 1, 4}]];
        donorToEnzBiParam = Flatten[Table[{
            Subscript[Global`kf, Global`BD][i] -> Global`boundaryTransfer,
            Subscript[Global`kr, Global`BD][i] -> Global`boundaryTransfer Global`Krdb}, {i, 1, 2}]];
        EnzRhToRlBiParam = Flatten[Table[{
            Subscript[Global`kf, Global`BT][i] -> Global`kCatb,
            Subscript[Global`kr, Global`BT][i] -> N[Global`kCatb/Keq]}, {i, 1, 4}]];
        EnzToAcceptorBiParam = Flatten[Table[{
            Subscript[Global`kf, Global`BA][i] -> Global`boundaryTransfer,
            Subscript[Global`kr, Global`BA][i] -> Global`boundaryTransfer Global`Krab}, {i, 1, 2}]];
        EnzRhToRlUniParam = Flatten[Table[{
            Subscript[Global`kf, Global`UT][i] -> Global`kCatu,
            Subscript[Global`kr, Global`UT][i] -> N[Global`kCatu/Keq]}, {i, 1, 1}]];
        completeParam = 
         Flatten[{donorFirstParam, acceptorFirstParam, acceptorSecondParam, 
           donorSecondParam, donorToEnzTriParam, EnzRhToRlTriParam, 
           EnzToAcceptorTriParam, donorToEnzBiParam, EnzRhToRlBiParam, 
           EnzToAcceptorBiParam, EnzRhToRlUniParam}];
        Emrd = -Log[Global`Krd] (R Global`T)/(Global`n F)(*mV*);
        Emra = -Log[Global`Kra] (R Global`T)/(Global`n F)(*mV*);
        Keq = Exp[(Global`n F)/(R Global`T) (Global`\[CapitalDelta]Em - Emrd - Emra)];
        reactionType = ReactionType/.{opts}/.Options[RedoxComplex];
        outputForm = Output/.{opts}/.Options[RedoxComplex];
        removeList = RemoveList/.{opts}/.Options[RedoxComplex];
        If[ reactionType=="Ordered",
            reactions = Join[PickReactions[Flatten[donorFirst[[1]], 1], {1}], 
            PickReactions[Flatten[acceptorFirst[[1]], 1], {5}], 
            PickReactions[Flatten[acceptorSecond[[1]], 1], {1}], 
            PickReactions[Flatten[donorSecond[[1]], 1], {13}], 
            PickReactions[
             TernaryTransfers, {1, 6, 12}], {{Enn -> Global`\[EmptySet], 
              0}, {Enp -> Global`\[EmptySet], 0}, {Epn -> Global`\[EmptySet], 0}}
            ];,
            If[ reactionType=="Random",
                reactions = Join[PickReactions[Flatten[donorFirst[[1]], 1], {1, 5}], 
                PickReactions[Flatten[acceptorFirst[[1]], 1], {1, 5}], 
                PickReactions[Flatten[acceptorSecond[[1]], 1], {1, 13}], 
                PickReactions[Flatten[donorSecond[[1]], 1], {1, 13}], 
                PickReactions[
                 TernaryTransfers, {1, 6, 12}], {{Enn -> Global`\[EmptySet], 
                  0}, {Enp -> Global`\[EmptySet], 0}, {Epn -> Global`\[EmptySet], 0}}
                ];,
                If[ reactionType=="PingPong",
                    reactions = Join[PickReactions[Flatten[donorFirst[[1]], 1], {1, 6}],
                    PickReactions[Flatten[acceptorFirst[[1]], 1], {2, 5}],
                    PickReactions[
                     BinaryTransfers, {1, 5, 7}], {{Enn -> Global`\[EmptySet], 
                      0}, {Epn -> Global`\[EmptySet], 0}}
                    ];,
                    reactions = Join[BindingReactions, TernaryTransfers, BinaryTransfers, 
                    UninaryTransfers];
                ];
            ];
        ];
        Which[outputForm == "Matrix",
        removeNull = Select[reactions[[1 ;;, 1]], MatchQ[#, b_ -> Global`\[EmptySet]]&];
        If[ Length[removeNull]>0,
            removeTable = 
                 Table[{Position[reactions, removeNull[[i]]]}, {i, 1, Length[removeNull]}];
            reactions = Delete[reactions, Transpose[{Flatten[removeTable, 2][[1 ;;, 1]]}]];
            matrixList = Union[Complement[Enz, removeNull[[1 ;;, 1]]], DeleteCases[reactions[[1 ;;, 1, 1]], Plus[b_, c_]], 
            DeleteCases[reactions[[1 ;;, 1, 2]], Plus[b_, c_]]];,
            matrixList = Union[Enz, DeleteCases[reactions[[1 ;;, 1, 1]], Plus[b_, c_]], 
            DeleteCases[reactions[[1 ;;, 1, 2]], Plus[b_, c_]]];
        ];
        Return[{matrixList, Table[UnitaryConstant[reactions, matrixList[[i]], matrixList[[j]]], 
            {i, Length[matrixList]}, {j, Length[matrixList]}],completeParam,paramList}];
        , 
        outputForm == "xCellerator"
        ,
        removeTable = 
        Table[{Position[reactions, removeList[[i]]]}, {i, 1, 
        Length[removeList]}];
        If[ Length[removeTable]>0,
            reactions = Delete[reactions, Transpose[{Flatten[removeTable, 2][[1 ;;, 1]]}]
                ];
        ];
        Return[reactions/. completeParam /. paramList];
        ];
    ]

(*SSVelocityEq*)  
SSVelocityEq[inputVector_List, opts___Rule] :=
    Module[ {rcm, iel, el, sl, rm, xx, yy, rate, bm, cm,nm2, nm3, p5, 
    cpa, i, j, l, n, o, X, cpm, nmf, npm, en, denom1, denom2, num, report, show, values, substrates, products, vmaxf, vmaxr, pos, donorBoundaryTransferPos, acceptorBoundaryTransferPos, ssflux, nb},
        rcm = inputVector[[2]];
        iel = inputVector[[1]];
        If[ iel == {},
            el = Table[e[i], {i, Length[rcm]}],
            el = iel
        ];
        Off[Solve::ratnz];
        donorBoundaryTransferPos = Position[inputVector[[2]], Subscript[Global`kf, Global`BD][b_]];
        If[ Length[donorBoundaryTransferPos]>0,
            acceptorBoundaryTransferPos = Position[inputVector[[2]], Subscript[Global`kf, Global`BA][b_]];,
            donorBoundaryTransferPos = Position[inputVector[[2]], Subscript[Global`kf, Global`TD][b_]];
            acceptorBoundaryTransferPos = Position[inputVector[[2]], Subscript[Global`kf, Global`TA][b_]];
        ];
        substrates = {inputVector[[1]][[donorBoundaryTransferPos[[1, 1]]]][[1]],inputVector[[1]][[acceptorBoundaryTransferPos[[1, 1]]]][[-1]]};
        products = {inputVector[[1]][[donorBoundaryTransferPos[[1, 2]]]][[1]],inputVector[[1]][[acceptorBoundaryTransferPos[[1, 2]]]][[-1]]};
        sl = Select[Flatten[rcm /. Times -> List], ! MatchQ[#, 0] && ! MatchQ[#, x_Subscript] && !NumericQ[#] && ! MatchQ[Head[#], x_Subscript] && !MatchQ[#, x_k] &];
        sl = DeleteDuplicates[sl];
        rm = Table[rcm[[i]]*el[[i]], {i, Length[rcm]}];
        xx[1] = Table[Position[rm, sl[[i]]*b_], {i, Length[sl]}];
        yy[1] = Table[Reverse[xx[1][[i, j]]],
          {i, Length[xx[1]]}, {j, Length[xx[1][[i]]]}];
        xx[2] = 
         Table[Extract[rm, xx[1][[i, j]]], {i, Length[xx[1]]}, {j, 
           Length[xx[1][[i]]]}];
        yy[2] = 
         Table[Extract[rm, yy[1][[i, j]]], {i, Length[yy[1]]}, {j, 
           Length[xx[1][[i]]]}];
        xx[3] = Table[xx[2][[i]] /. List -> Plus, {i, Length[xx[2]]}];
        yy[3] = Table[yy[2][[i]] /. List -> Plus, {i, Length[yy[2]]}];
        Do[
         rate[sl[[i]]] = yy[3][[i]] - xx[3][[i]], {i, Length[sl]}];
        n = Length[rcm];
        bm = Table[Delete[rcm, i], {i, n}];
        cm = Table[DeleteCases[bm[[i, j]], 0], {i, Dimensions[bm][[1]]}, {j,Dimensions[bm][[2]]}];
        nm2 = Table[Select[Tuples[cm[[i, ;;]]], Position[#, 0] == {} &], {i, n}];
        p1 = Table[Position[nm2[[i]], rcm[[j, i]]], {i, n}, {j, n}];
        p2 = Table[Flatten[p1[[i]], 1], {i, n}];
        p3 = Table[p2[[j, i, 1]], {j, n}, {i, Length[p2[[j]]]}];
        p4 = Table[Union[p3[[i]]], {i, n}];
        nm3 = Table[nm2[[i, p4[[i, l]]]], {i, n}, {l, Length[p4[[i]]]}];
        For[i = 1, i < n + 1, i++,
         For[j = 1, j < Length[nm3[[i]]] + 1, j++,
           For[m = 1, m < n, m++,
             p5[i, j, m, 1] = Position[rcm, nm3[[i, j, m]]][[1, 2]];
             If[ p5[i, j, m, 1] == i,
                 Continue[]
             ];
             (o = 1;
              Label[xxx];
              X = o + 1;
              p6 = Flatten[Table[Position[
                  rcm[[p5[i, j, m, o]]], nm3[[i, j, t]]], {t, n - 1}]];
              If[ p6 === {},
                  cpa[i, j] = nm3[[i, j]];
                  m = n;
                  Continue[]
              ];
              p5[i, j, m, X] = p6[[1]];
              If[ p5[i, j, m, X] == i,
                  Continue[]
              ];
              o++;
              If[ o < n,
                  Goto[xxx]
              ]);
             cpa[i, j] = nm3[[i, j]];
             ];
      ];
    ];
        cpm = DeleteCases[Table[cpa[i, j], 
            {i, n}, {j, Length[nm3[[i]]]}], 
          cpa[x_, y_], Infinity];
        nmf = Table[Complement[nm3[[i]], cpm[[i]]], {i, n}];
        npm = Table[Apply[Times, nmf[[i, j]]],
            {i, n}, {j, Length[nmf[[i]]]}];
        Do[en[i] = Apply[Plus, npm[[i]]], {i, n}];
        denom1 = Apply[Plus, Flatten[npm]];
        denom2 = Collect[denom1, sl];
        num = Simplify[Table[rate[sl[[i]]], {i, Length[sl]}] /. 
           Table[el[[i]]->en[i], {i, Length[el]}]];
        report = Report/.{opts}/.Options[RateEquation];
        show = ShowConstants/.{opts}/.Options[ReactionMatrix];
        values = {};
        Which[show=="SteadyState", values = inputVector[[3]];,
          show=="Values", values = inputVector[[3]]/.inputVector[[4]];
      ];
        pos = Flatten[Position[sl, products[[1]]]];
        vmaxf = Limit[Limit[num[[pos[[1]]]]/denom2/.values, substrates[[1]] -> Infinity], 
        substrates[[2]] -> Infinity]/.{products[[1]]->0, products[[2]]->0};
        vmaxr = Limit[Limit[num[[pos[[1]]]]/denom2/.values, products[[1]] -> Infinity], 
        products[[2]] -> Infinity]/.{substrates[[1]]->0, substrates[[2]]->0};
        ssflux = Join[
                {Global`flux -> ((Global`e0 num[[pos[[1]]]])/denom2)/.values},
                {Global`vmaxf->vmaxf},
                {Global`vmaxr->vmaxr},
                {Global`Km[ToString[substrates[[1]]]]->Limit[substrates[[1]] /. 
            Flatten[Solve[(num[[pos[[1]]]]/denom2/.values) == vmaxf/2, 
            substrates[[1]]]/.{products[[1]] -> 0, products[[2]] -> 0}], substrates[[2]] -> Infinity]},
                {Global`Km[ToString[substrates[[2]]]]->Limit[substrates[[2]] /. 
            Flatten[Solve[(num[[pos[[1]]]]/denom2/.values) == vmaxf/2, 
            substrates[[2]]]/.{products[[1]] -> 0, products[[2]] -> 0}], substrates[[1]] -> Infinity]},
                {Global`Km[ToString[products[[1]]]]->Limit[products[[1]] /. 
            Flatten[Solve[(num[[pos[[1]]]]/denom2/.values) == vmaxr/2, 
            products[[1]]]/.{substrates[[1]] -> 0, substrates[[2]] -> 0}], products[[2]] -> Infinity]},
                {Global`Km[ToString[products[[2]]]]->Limit[products[[2]] /. 
            Flatten[Solve[(num[[pos[[1]]]]/denom2/.values) == vmaxr/2, 
            products[[2]]]/.{substrates[[1]] -> 0, substrates[[2]] -> 0}], products[[1]] -> Infinity]}    
            ];
        If[ report,
            nb = DocumentNotebook[
                Join[
                    {TextCell["Enzyme Distribution Functions", "Section"]},
                    Table[TextCell[Row[{el[[i]]/Global`eo," = ", en[i] /. values}], "Text"], {i, n}],
                    {TextCell["Steady-State Differential Equations", "Section"]},
                    Table[TextCell[Row[{"d", sl[[i]],"/dt = ", Global`e0 num[[i]]/Global`Denom /. values}], "Text"],{i,Length[sl]}],
                    {TextCell["Denominator"]},
                    {TextCell[denom2 /. values, "Text"],
                    TextCell["Steady-State Flux Equation and Kinetic Parameters", "Section"],
            		ExpressionCell[ssflux // TraditionalForm, "Output"]}                      
                ]
                ];
            Return[Join[ssflux,{Global`SSReport->nb}]],
            Return[ssflux]
        ];
    ];
   
 ReactionMatrix[inputVector_List, opts___Rule] :=
     Module[ {matrix,show,values},
         show = ShowConstants/.{opts}/.Options[ReactionMatrix];
         values = {};
         Which[show=="SteadyState", values = inputVector[[3]];,
          show=="Values", values = inputVector[[3]]/.inputVector[[4]];
      ];
         matrix = Join[Transpose[{inputVector[[1]]}], inputVector[[2]], 2];
         Join[{Join[{"Reaction Matrix"}, inputVector[[1]]]}, 
         matrix/.values] // MatrixForm
     ]  
   
    
End[]

EndPackage[]

