Supplementary notebook for article : "Tree Graphs and Parameter Estimation in Metabolic Networks"
This notebook contains an example time series data, example adjacency matrix to model the network, a function to fit natural splines to data and a function to estimate the dynamic kinetic rates.
l.j.astola@tue.nl

Example data.

7 example  time series of measured metabolites (these can be replaced, e.g., by the averages of the daily measurements under varying conditions in the supplementary Excel file: Metabolite_time_series_data.xlsx)

timepoints={5.,5.5,6.,6.5,7.,7.5,8.,8.5,9.};

measurements={{12.3222,15.8784,16.0821,14.2798,11.8181,10.9275,10.3022,9.52015,8.15967},{31.9721,46.3335,54.1569,57.2115,57.2661,56.7307,55.4519,53.9173,52.6142},{28.515,26.8844,24.0293,20.9113,18.4919,18.9783,19.5958,18.815,15.1067},{1.49297,1.70505,1.75752,1.71758,1.65242,1.66866,1.71517,1.78028,1.85229},{0.682998,0.974363,1.03428,0.971729,0.895717,1.01166,1.13927,1.19471,1.09413},{3.98829,4.80188,4.84432,4.50029,4.15447,4.56483,4.9962,5.08667,4.47436},{2.25139,2.15267,1.99159,1.8065,1.63578,1.56848,1.49088,1.33995,1.05267}};

pointpairs=Map[Transpose[{timepoints,#}]&,measurements];

Show[Map[ListPlot[#,Joined->True,PlotRange->{{First[timepoints],Last[timepoints]},{Min[measurements],Max[measurements]}}]&,pointpairs],Frame->True,FrameTicksStyle->Directive[20]]

Network topology

The networkmatrix has number one in element (i, j), when there is "matter flow" from node j to node i. Since in our case there is influx from outside of the network (node 0), the matrix rows are numbered 0, 1, 2, ... and the columns are numbered 1, 2, 3, ...

networkmatrix={{0,0,0,0,0,0,0,0},{1,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0},{0,1,0,0,0,0,0,0},{0,0,1,0,0,0,0,0},{0,0,1,0,0,0,0,0},{0,0,0,1,0,0,0,0},{0,0,0,1,0,0,0,0}};

MatrixForm[networkmatrix]

Functions.

Fit natural cubic splines to the given data. This function is used in the subsequent kinetic rate solving function.

natspline::usage="produce cubic spline polynomial coefficients, for a given set of datapoints. Datapoints should be in the following form: {{x0,y0},{x1,y1},...}"

natspline[timepoints_,values_]:=
Module[{numberofpoints,mat,bvecx,bvecy,dx,dy,data},
data=Transpose[{timepoints,values}];
   numberofpoints=Length[data];
   
mat=1.*Inverse[DiagonalMatrix[Flatten[{2,Table[4,{i,numberofpoints-2}],2}]]+
DiagonalMatrix[Table[1,{i,numberofpoints-1}],1]+DiagonalMatrix[Table[1,{i,numberofpoints-1}],-1]];

bvecx=Flatten[{3*(data[[2,1]]-data[[1,1]]),Table[3*(data[[i,1]]-data[[i-2,1]]),
{i,3,numberofpoints}],3*(data[[numberofpoints,1]]-data[[numberofpoints-1,1]])}];

bvecy=Flatten[{3*(data[[2,2]]-data[[1,2]]),Table[3*(data[[i,2]]-data[[i-2,2]]),
{i,3,numberofpoints}],3*(data[[numberofpoints,2]]-data[[numberofpoints-1,2]])}];

dx=mat.bvecx;
dy=mat.bvecy;

Transpose[{Table[{data[[i,1]],dx[[i]],3*(data[[i+1,1]]-data[[i,1]])
-2*dx[[i]]-dx[[i+1]],2*(data[[i,1]]-data[[i+1,1]])+dx[[i]]+dx[[i+1]]},{i,1,numberofpoints-1}],
Table[{data[[i,2]],dy[[i]],3*(data[[i+1,2]]-data[[i,2]])-2*dy[[i]]-dy[[i+1]],2*(data[[i,2]]
-data[[i+1,2]])+dy[[i]]+dy[[i+1]]},{i,1,numberofpoints-1}]}]]

coeftopoly::usage="convert polynomial coefficients to polynomials of x"

coeftopoly[coefmat_,x_]:=
Partition[Table[Apply[#1+#2*x+#3*x^2+#4*x^3&,Flatten[coefmat,1][[i]]],
{i,1,Dimensions[Flatten[coefmat,1]][[1]]}],2]

Solve the time varying kinetic rates.

KineticRates::usage="given a sampleset of points in 2D, and a network topology(adjacency)
matrix, solves timevariant kinetic constants, computes the residual and the
solution curves to timevariant ODEs as well as the timevariant discrete (pointwise) kinetic parameter matrix.
Output is given as a vector containing the residual, the reconstructed concentration curves,
2. order approximation of estimated time varying kinetic functions and the system matrices evaluated at data points."

KineticRates[timepoints_,measurements_,networkmatrix_]:=
Module[{edges,dimension,adjacencymatrix,systemmatrix,parameters,polynomials,derivatives,solutions,funs,odes,continuous,discrete,curves,parameterpolynomials},
pointpairs=Map[Transpose[{timepoints,#}]&,measurements];
init=First[timepoints];
end=Last[timepoints];
δt=(end-init)/(Length[timepoints]-1);

edges=Position[networkmatrix,1];
dimension=Length[networkmatrix]-1;
adjacencymatrix=ReplacePart[
Table[0,{i,dimension+1},{j,dimension+1}],
Table[edges[[i]]->ToExpression["k"<>ToString[edges[[i,2]]]<>ToString[edges[[i,1]]]],
{i,Length[edges]}]];

systemmatrix=ReplacePart[adjacencymatrix,Table[{i,i}->-Total[Transpose[adjacencymatrix][[i]]],{i,dimension+1}]];
parameters=Cases[Flatten[adjacencymatrix],Except[0]];

polynomials=Table[coeftopoly[natspline[timepoints,measurements[[k]]],x],{k,Length[measurements]}];
derivatives=Table[D[polynomials[[l,k,2]],x]/D[polynomials[[l,k,1]],x],{l,Dimensions[polynomials][[1]]},
{k,Dimensions[polynomials][[2]]}];

solutions=Join[Table[NMinimize[{((derivatives[[All,m]]/.{x->0})-
((Rest[systemmatrix].Prepend[polynomials[[All,All,2]],
Table[1,{l,Dimensions[polynomials][[2]]}]])[[All,m]]/.{x->0})).((derivatives[[All,m]]/.{x->0})
-((Rest[systemmatrix].Prepend[polynomials[[All,All,2]],
Table[1,{l,Dimensions[polynomials][[2]]}]])[[All,m]]/.{x->0})),
ToExpression[StringReplace[ToString[Map[#>=0.0&,Flatten[parameters]]],","->"&&"]]},parameters],
{m,Dimensions[derivatives][[2]]}],{NMinimize[{((derivatives[[All,Dimensions[derivatives][[2]]]]/.{x->1})
-((Rest[systemmatrix].Prepend[polynomials[[All,All,2]],
Table[1,{l,Dimensions[polynomials][[2]]}]])[[All,Dimensions[derivatives][[2]]]]/.{x->1})).
((derivatives[[All,Dimensions[derivatives][[2]]]]/.{x->1})
-((Rest[systemmatrix].Prepend[polynomials[[All,All,2]],
Table[1,{l,Dimensions[polynomials][[2]]}]])[[All,Dimensions[derivatives][[2]]]]/.{x->1})),
ToExpression[StringReplace[ToString[Map[#>0.&,Flatten[parameters]]],","->"&&"]]},parameters]}];

discrete=Map[systemmatrix/.#[[2]]&,solutions];

curves=Table[ToExpression["x"<>ToString[i]],{i,dimension}];

funs=Join[{1},Table[ToExpression["x"<>ToString[i[t]]],{i,dimension}]];

continuous=Table[Fit[Table[{5+(i-1)*δt,discrete[[i,j,k]]},
{i,Dimensions[discrete][[1]]}],{1,t,t^2},t],
{j,Dimensions[discrete][[2]]},
{k,Dimensions[discrete][[3]]}];

odes=Join[Table[D[funs[[i]],t]==Sum[continuous[[i,j]]*funs[[j]],{j,dimension+1}],
{i,2,dimension+1}],
Table[curves[[i]][init]==pointpairs[[i,1,2]],{i,dimension}]];

sols=NDSolve[odes,curves,{t,init,end}];

curves=curves/.sols[[1]];
parameterpolynomials=Cases[Flatten[ReplacePart[continuous,Table[{i,i}->0.,{i,Length[continuous]}]]],Except[0.]];

residual=Norm[Table[curves[[i]][t],{i,dimension},{t,init,end,δt}]-pointpairs[[All,All,2]],"Frobenius"];

Print["residual =", residual ];
Do[Print["estimation of parameter "<>ToString[parameters[[i]]]<>" = ",parameterpolynomials[[i]]],{i,1,Length[parameters]}];{solutions,residual,curves,discrete}]

Example of usage and visualization of results.

Plot the network topology as a graph.

GraphPlot[Transpose[networkmatrix],DirectedEdges->True]

Plot the reconstructed concentrations vs.original data points.

Clear[solutions,residual,curves,discrete]

{solutions,residual,curves,discrete}=KineticRates[timepoints,measurements,networkmatrix];

Show[Table[Plot[curves[[i]][t],{t,5,9},PlotStyle->{Thick,Dashed,ColorData["Rainbow"][i/Length[measurements]]}],{i,Length[measurements]}],ListPlot[pointpairs,PlotStyle->Table[{ColorData["Rainbow"][i/Length[measurements]],PointSize[0.02]},{i,Length[measurements]}],PlotRange->All],AxesLabel->{Text[Style["time",Medium]],Text[Style["concentration",Medium]]},PlotRange->All]