(* ::Package:: *)

(* ::Subtitle:: *)
(*Custom Mathematica function for chemical reaction network analysis*)


(* ::Text:: *)
(*This Mathematica package contains the MakeCRN custom function used in several other notebooks. This package was originally developed by Robert L. Karp and modified as needed by Tathagata Dasgupta and Joseph P. Dexter. References: Karp et al., "Complex-linear invariants of biochemical reaction networks," J. Theor. Biol. 311 (2012) 130-138 and Dexter et al., "Invariants reveal multiple forms of robustness in bifunctional enzyme systems," Biophys. J., submitted. *)


(* ::Text:: *)
(*Preliminary functions: The following functions will be used in the main custom functions. *)


BeginPackage["CRN`"]

(* Creates a zero matrix *)
makeZero[n_] := Table[0, {i,n}, {j,n}]

(* Sets the off-diagonal components of A\[Kappa]. A\[Kappa] is the matrix of coefficients in the polynomial system.   *)
makeK[L_] := Module[{i=L[[1]], j=L[[2]], s},
 s = Subscript[\[Kappa],i->j];
 A\[Kappa][[i,j]] = s ]

(* Sets the diagonal components of A\[Kappa]. *)
diagSet[i_,NL_] := Module[{s=0},
 Do[s += A\[Kappa][[i,j]], {j,NL}];
 A\[Kappa][[i,i]] = -s]

(* Creates a vector (\[CapitalPsi]) of the species in the networks. *)
\[CapitalPsi]F[L_] := Module[{s=1},
 Do[ s *= species[[j]]^L[[j]], {j,SL}];
s ]

(* Converts a product into a sum. *)
ToSum[P_] := Module[{p,f=0},
 p = FactorList[P ];
 Do[f = f + p[[i]][[1]]
  ,{i,2,Length[p]}];
f ]

(* Converts a sum into a product. *)
ToProduct[P_] := Module[{p,f=1},
 p = MonomialList[P ];
 Do[ f = f * p[[i]]
  ,{i,Length[p]}];
 f ]

(* Converts a list into an arrow for plotting. *)
ToArrow[L_] := Module[{i=L[[1]], j=L[[2]]},
Return[ ToSum[ Comple[[i]] ]-i -> ToSum[ Comple[[j]] ]-j ];
(* i->j *) ]

(* Creates alternate notation for the rate constants. *)
RenumberRates[crn_]:= Module[{f={}},
Do[
f = Append[f,crn [[i]]-> Subscript[k, i]];
,{i,Length[crn]}]; f ]

(* Finds an element in a list and returns its position. *)
FindInList [x_,L_] := Module[{index=-1},
 Do[ If[x==L[[i]], index=i], {i,Length[L]} ];
 index]

(* Finds an expression in a list. *)
FindInList[x_, list_] := Module[{loc = -1},
 Do[ If[x===list[[i]], loc = i; Break[] ]
  ,{i,Length[list]}];
 Return[loc]]

(* Generates the Y matrix for later use. *)
makeYPart[Complex_,Vars_] := Module[{t=Vars},
 Do[t[[i]] = Exponent[Complex,Vars[[i]]]
 ,{i,Length[Vars]}]; 
t]
makeY[Complex_,Vars_] := Module[{L=Length[Complex],t},
 t = Table[0,{i,L}];
 Do[t[[i]] = makeYPart[Complex[[i]],Vars],{i,L}]; 
t]
(*Alternate numbering system. *)
RenumberRates[crn_]:= Module[{f={}},
Do[
f = Append[f,crn [[i]]-> Subscript[k, i]];
,{i,Length[crn]}]; f ]

numbered=RenumberRates[RateConsts];


(* ::Text:: *)
(*Basic network construction function: Simple function (makeCRN) for mass-action analysis and visualization of reaction networks. *)


makeCRN[species_,CRN_,YY_,Deficiency_, ODE_] := Module[
 {GL,NL,CG,CVects,Stoichiometry,CRNGl,\[CapitalPsi]ND,CGND},
 SL = Length[species];(* species *)
 GL = Length[CRN];    (* reactions *)
 NL = Length[YY];     (* complexes *)

 (* Checks for duplicates. *)
 If[
	Length[Comple] != Length[DeleteDuplicates[Comple] ],
	If[Length[Comple] != Length[DeleteDuplicates[Comple] ]+ Count[Comple,1]-1
		, Print["Wrong Complexes!"] ] ];

 (* Creates the A\[Kappa] matrix. *)
 A\[Kappa] = makeZero[NL];
 RateConsts = Map[makeK,CRN];
 Do[ diagSet[i,NL], {i,NL}];

(* Computes the deficiency of the network. *)
 If[Deficiency,
 Print["Deficiency: ", Length[NullSpace[Transpose[A\[Kappa].YY]]]-Length[NullSpace[Transpose[A\[Kappa]]]] ]; ];

 (* Constructs the various matricies. *)
 \[CapitalPsi]= Map[ \[CapitalPsi]F, YY];

 (* Generates the complex balanced ideal. *)
 CG = Simplify[ \[CapitalPsi].A\[Kappa] ];

(* Generates the steady-state ideal. *)
 SteadyState = FullSimplify[ CG.YY];

(* Display the polynomials in SteadyState as separate equations f_ 1,...,f_n. *)
	For[i=0,i<=Length[SteadyState], i++,Subscript[f, i]=SteadyState[[i]]];
	For[i=1,i<=Length[SteadyState],i++,Print[Subscript["f", i]," = ",Subscript[f, i]]];

 (* Generates and displays the reaction network diagram. *)
 CRNGl = Map[ToArrow, CRN];
 GraphPlot[CRNGl, VertexLabeling->True, DirectedEdges->True] ]



(* ::Text:: *)
(*Full network construction function: Full function (MakeCRN) for mass-action analysis and visualization of reaction networks that uses makeCRN. *)


(* Generates a list of species in the network. *)
MakeSpecies[L_] := Module[{s={}},
 Do[ s = Union[s, MonomialList[L[[i]]] ]
  ,{i,Length[L]}];
 Return[ Complement[ s, {1}] ] ]

(* Generates a list of all individual complexes. *)
MakeNetwork[L_] := Module[{l={},t},
 Do[
  Do[t=L[[i]][[j]];
	 (* If it is numeric complex (synthesis or degradation) then included it again *)
     If[ FindInList[t,l]==-1 || True==NumericQ[t], l=Append[l,t] ]
    ,{j,2}];
 ,{i,Length[L]}];
 Return[l] ]

(* Constructs the numbered chemical reaction network from the list of all the individual complexes. *)
MakeTheCRN[L_,C_] := Module[{crn={},t,loc={0,0},prev=0},
 Do[
  Do[ t = L[[i]][[j]];
   If[ True==NumericQ[t] ,
    loc[[j]] = prev + FindInList[ t,C[[prev+1;; ]]  ];
    prev = loc[[j]], loc[[j]] = FindInList[ t,C ] ]
   ,{j,2} ];
   crn = Append[crn, loc ];
   If[Length[ L[[i]] ]!=2, crn=Append[crn,{ loc[[2]], loc[[1]] }] ]
  ,{i,Length[L]}];
 Return[crn] ]

(* Displays reaction network and associated analysis. *)
MakeCRN[Reactions_,Deficiency_:True,ODE_:False] := Module[{},
 Comple = MakeNetwork[Reactions];
 CRN = MakeTheCRN[Reactions,Comple];
 species = MakeSpecies[Comple];
 YY = makeY[Comple,species];
 makeCRN[species,CRN,YY,Deficiency,ODE] ]



(* ::Text:: *)
(*Kernel function: Computes the kernel of the matrix of polynomial coefficients from the reaction network. *)


(* Function for generating the kernel. *)
TheKernel[] := Module[{tmp,tmz},
 tmp = Transpose[NullSpace[Transpose[A\[Kappa].YY]]];
 tmz = tmp;
 Do[ 
  Do[ 
   If[0===tmp[[i]][[j]], tmv=0, tmv="*"];
   tmz[[i]][[j]] = tmv
  ,{j, Length[ tmp[[1]] ]}]
 ,{i, Length[tmp]}];
 Print[Join[Transpose[{Comple}], tmz, 2]//MatrixForm];
 tmz=Join[Transpose[{Comple}], tmp, 2];
 tmz ]
EndPackage[];
