Notebook for constructing and manipulating polynomials generated by the Matrix-Tree Theorem, as applied to the Laplacian matrix of a labelled, directed graph.

Background

This notebook provides capabilities for the linear framework, which is developed in a series of papers

[1] J. Gunawardena, “A linear framework for time-scale separation in nonlinear biochemical systems”, PLoS ONE 7:e36321, 2012. 
[2] I. Mirzaev, J. Gunawardena, “Laplacian dynamics on general graphs”, Bull Math Biol 75:2118-49 2013. 
[3] J. Gunawardena, “Time-scale separation: Michaelis and Menten’s old idea, still bearing fruit”, FEBS J, invited review, in press, 2013
[4] T. Ahsendorf, F. Wong, R. Eils, J. Gunawardena, “A framework for modelling gene regulation that accommodates non-equilibrium mechanisms”, submitted, 2013. 

The linear framework is based on a labelled, directed graph, G, which gives rise to a dynamics, dx/dt = ℒ(G).x, where ℒ(G) is the Laplacian matrix of G [1]. This linear dynamics can be thought of as as first-order chemistry or as the corresponding stochastic master equation [2] and has wide applications in different biological contexts [3,4]. The steady-states of this dynamics or, equivalently, the elements of the kernel of the Laplacian, are the subject of interest. For a strongly-connected graph, the kernel of the Laplacian has dimension 1 and a basis element can be constructed by the Matrix-Tree Theorem [1].

Usage

Before using this notebook, a text file with the extension .txt should be created (for example, graph.txt) in which a graph is described as a list of triples, (i,a,j), signifying a directed edge from vertex i to vertex j with label a. The graph is assumed to be strongly connected. The python script MTT.py (see Further Information below) should then be executed. This reads graph.txt and produces, for each vertex i in the graph, a file graph-i.txt, which enumerates the spanning trees rooted at i, as prescribed by the Matrix-Tree Theorem. The script also produces a cleaned description of the graph in the filegraph-parsed.txt. This notebook reads these files and expects to be executed in the directory which contains them. 

Set the variable name GraphName below to the name of the graph in the graph description file (ie: GraphName = “graph” if the file is graph.txt) and then execute the functions below. The notebook builds the basis vector ρ (called here “rho”) in the kernel of the Laplacian matrix of the graph. For each vertex i, the i-th component of the basis vector is created as a polynomial in the edge labels

rhoi[(Subscript[a, 1]), ... , (Subscript[a, n])]

with the labels ordered as they appear in graph.txt.These polynomials can then be interrogated using Mathematica functions, as required, depending on the context.

Further information

The script MTT.py and further details, including contact information, are available at vcp.med.harvard.edu/software.html#linear. This notebook was developed by Felix Wong for the paper [4].

(* kill kernel *)
Quit[]

(* specify graph *)
GraphName = "";
graph=Import[NotebookDirectory[]<>GraphName<>"-parsed.txt","Table"];
MathematicaGraph=ToExpression[StringReplace[ToString[graph],{"{{"->"{","}}"->"}","("->"{",")"->"}","'"->"\""}]]/.{f_Integer,m_,l_Integer}:>{f->l,ToString[m]};
(* get labels *)
listoflabels = DeleteDuplicates[MathematicaGraph[[All,2]]];
lstring = "{";
For[m=1,m≤Length[listoflabels],m++,
lstring = lstring <> "{"<>ToString[m]<>",\""<>listoflabels[[m]]<>"\"},";]
flstring = StringDrop[lstring ,-1]<> "}";
ToExpression["labelrule="<>flstring];
(* function that comes in handy when plotting *)HashString[f_String]:=Mean[ToCharacterCode[f]]+ToCharacterCode[f].ToCharacterCode[f]

Check to see what the graph looks like

(* plot graph *)
GraphPlot[MathematicaGraph, VertexLabeling->True, DirectedEdges->True,VertexCoordinateRules->Auto,VertexLabeling->True,
VertexRenderingFunction->({White,EdgeForm[Black],Disk[#,.1],Black,Style[Text[#2,#1],Large,Italic,FontFamily->"Helvetica"]}&),
EdgeRenderingFunction->(
{Hue[HashString[#3],1,0.7]
,Arrow[#1,0.1],If[#3=!=None,Text[Style[#3,Medium,Hue[HashString[#3],1,0.7],Italic,FontFamily->"Helvetica"],
If[Length[#1]>2,Mean[#1]+(#1[[1]]-#1[[2]])*0.5,Mean[#1]]
,Background->White],{}]}&),
ImageSize->700
]

(* plot graph in 3D *)
GraphPlot3D[MathematicaGraph, VertexLabeling->True, DirectedEdges->True,VertexCoordinateRules->Auto,VertexLabeling->True,
VertexRenderingFunction->({White,EdgeForm[Black],Sphere[#,.1],Black,Style[Text[#2,#1],Large,Italic,FontFamily->"Helvetica"]}&),
EdgeRenderingFunction->(
{Hue[HashString[#3],1,0.7],Arrow[#1,0.1],If[#3=!=None,Text[Style[#3,Medium,Hue[HashString[#3],1,0.7],Italic,FontFamily->"Helvetica"],If[Length[#1]>2,Mean[#1]+(#1[[1]]-#1[[2]])*0.7,Mean[#1]],Background->White],{}]}&),
ImageSize->700
]

Make the components of the basis vector rho, in the kernel of the Laplacian

(* initializing variables and functions *)
(* find number of nodes in graph *)
numberofnodes = Max[MathematicaGraph/.{f_Integer->l_Integer,m_String}:>{f,l}];

rhofnarg= "";
identityelement = "";
For[m=1,m≤Length[listoflabels],m++,
rhofnarg = rhofnarg <> labelrule[[m]][[2]]<>"_,";
identityelement = identityelement <> "1,";
]
rhoargument = StringDrop[rhofnarg,-1];
identity = StringDrop[identityelement,-1];

(* build the rho's from the imported data *)
For[i=1,i≤numberofnodes,i++,
ToExpression[ToString["rhofn"<>ToString[i]<>"=Import[NotebookDirectory[]<>"<>ToString["\""<>GraphName<>"-"<>ToString[i]<>".txt"]<>"\",\"Table\"]"]];
totalarg = "";

For[m=1,m≤Length[listoflabels],m++,
ToExpression[ToString["N"<>ToString[m]<>ToString[i]<>"=rhofn"<>ToString[i]<>"[[All,"<>ToString[m]<>"]]"]];
totalarg = StringReplace[totalarg <> "*"<>labelrule[[m]][[2]]<>"^N"<>ToString[m]<>ToString[i],{"*"<>labelrule[[1]][[2]]->labelrule[[1]][[2]]}];
]

ToExpression[ToString["rhoexpression"<>ToString[i]<>"=Total["<>totalarg<>"]"]]; 
ToExpression[ToString["rho"<>ToString[i]<>"["<>rhoargument<>"]=rhoexpression"<>ToString[i]]];
]

sumstring = "";
For[i=1,i≤numberofnodes,i++,
sumstring = sumstring <>"+"<>"rhoexpression"<>ToString[i];
]
sumstring=StringDrop[sumstring,1];
ToExpression["rhotot["<>rhoargument<>"]="<>sumstring];

Work out the total number of spanning trees in the graph

(* number of spanning trees is sum of rho_i's at identity *) 
rhosum = 0;
For[i=1,i≤numberofnodes,i++,
rhosum = rhosum + ToExpression["rho"<>ToString[i]<>"["<>identity<>"]"]
]

rhosum

Print out the rhoi. Note that rhotot gives rho1 + ... + rhon (ie: the partition function)

(* print out rhos *) 
For[i=1,i≤numberofnodes,i++,
Print["rho"<>ToString[i]<>"["<>StringReplace[rhoargument,{"_"->""}]<>"]="ToExpression["rho"<>ToString[i]<>"["<>StringReplace[rhoargument,{"_"->""}]<>"]"]]
]