Supplementary notebook for article : "Tree Graphs and Parameter Estimation in Metabolic Networks"

l.j.astola@tue.nl

Sloppiness/identifiability of the parameters

Network parameters can be estimated from time series data using the other supplementary notebook “DynamicKineticRateSolver.nb”.  This notebook is for testing to which extent the parameters can be identified.

After the parameters are estimated and the concentration functions reconstructed, we use the parameter matrix (automatically) obtained in the process to study the identifiability. The results are listed below.

1. Estimate the parameters

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

2. Choose a parameter

pos=Position[networkmatrix,1](*this network model is from notebook DynamicKineticRateSolver.nb, you may substitute with an alternative tree structured network matrix*)

Use finite differences to compute Hessian eigenvalues

parameter wise variation (this function is used in computing the Hessian)

vary[parameterposition_,matrixlist_,model_,sampleset_]:=Module[{matlist,perturbed,timeser,odes},funs=Join[{1},Table[ToExpression["x"<>ToString[i]<>"[t]"],{i,Dimensions[sampleset][[1]]}]];
funns=Table[ToExpression["x"<>ToString[i]],{i,Dimensions[sampleset][[1]]}];
perturbed=Table[ReplacePart[matrixlist[[i]],{parameterposition->(1+0.2*RandomVariate[NormalDistribution[]])*matrixlist[[i,parameterposition[[1]],parameterposition[[2]]]]}],{j,1,10},{i,Length[matrixlist]}];timeser=Table[{sampleset[[1,All,1]][[j]],perturbed[[i,j,k,l]]},{i,Dimensions[perturbed][[1]]},{k,Dimensions[perturbed][[3]]},{l,Dimensions[perturbed][[4]]},{j,Dimensions[perturbed][[2]]}];
ippolsysmats=Map[(model+IdentityMatrix[Dimensions[model][[1]]])*#&,Table[Interpolation[timeser[[i,j,k]],InterpolationOrder->1][t],{i,Dimensions[timeser][[1]]},{j,Dimensions[timeser][[2]]},{k,Dimensions[timeser][[3]]}]];
odes=Table[Join[Table[D[funs[[i]],t]==Sum[ippolsysmats[[k,i,j]]*funs[[j]],{j,dim+1}],
{i,2,dim+1}],
Table[funns[[i]][init]==sampleset[[i,1,2]],{i,dim}]],{k,Length[ippolsysmats]}];
sols=Table[NDSolve[odes[[k]],funns,{t,init,end}],{k,Length[odes]}];
functions=Table[funns/.sols[[k,1]],{k,Length[sols]}];
residuals=Table[Norm[Table[functions[[k,i]][t],{i,dim},{t,sampleset[[1,All,1]]}]-sampleset[[All,All,2]],"Frobenius"],{k,Length[functions]}];
{residuals,functions,ippolsysmats}]

sloppy parameter wise variation (this function is used in computing the Hessian)

take each matrix given by timevariant separately to obtain derivatives

δ1=0.000001

δ2=0.000001

parder[parameterposition_,matrix_,δ_]:=ReplacePart[matrix,{parameterposition->(1+δ)*matrix[[parameterposition[[1]],parameterposition[[2]]]]}]

parcur[parameterposition1_,parameterposition2_,matrix_,δ1_,δ2_]:=ReplacePart[ReplacePart[matrix,{parameterposition1->(1+δ1)*matrix[[parameterposition1[[1]],parameterposition1[[2]]]]}],{parameterposition2->(1+δ2)*matrix[[parameterposition2[[1]],parameterposition2[[2]]]]}]

the Hessian

hessian[m_,l_]:=CompoundExpression[parpos1=pos[[m]];parpos2=pos[[l]];
(*perturb the optimal solution matrix to =-directions*)
mlist1=Table[parder[parpos1,matrixlist[[i]],δ1],{i,5}];
mlist2=Table[parder[parpos1,matrixlist[[i]],-δ1],{i,5}];
mlist11=Table[parcur[parpos1,parpos2,matrixlist[[i]],δ1,δ2],{i,5}];
mlist12=Table[parcur[parpos1,parpos2,matrixlist[[i]],-δ1,δ2],{i,5}];mlist21=Table[parcur[parpos1,parpos2,matrixlist[[i]],δ1,-δ2],{i,5}];
mlist22=Table[parcur[parpos1,parpos2,matrixlist[[i]],-δ1,-δ2],{i,5}];
(*get new time series*)

tseries1=Table[{sampleset[[1,i,1]],mlist1[[i,j,k]]},{j,Length[matrixlist[[1]]]},{k,Length[Transpose[matrixlist[[1]]]]},{i,Length[sampleset[[1,All,1]]]}];
tseries2=Table[{sampleset[[1,i,1]],mlist2[[i,j,k]]},{j,Length[matrixlist[[1]]]},{k,Length[Transpose[matrixlist[[1]]]]},{i,Length[sampleset[[1,All,1]]]}];tseries11=Table[{sampleset[[1,i,1]],mlist11[[i,j,k]]},{j,Length[matrixlist[[1]]]},{k,Length[Transpose[matrixlist[[1]]]]},{i,Length[sampleset[[1,All,1]]]}];
tseries12=Table[{sampleset[[1,i,1]],mlist12[[i,j,k]]},{j,Length[matrixlist[[1]]]},{k,Length[Transpose[matrixlist[[1]]]]},{i,Length[sampleset[[1,All,1]]]}];tseries21=Table[{sampleset[[1,i,1]],mlist21[[i,j,k]]},{j,Length[matrixlist[[1]]]},{k,Length[Transpose[matrixlist[[1]]]]},{i,Length[sampleset[[1,All,1]]]}];
tseries22=Table[{sampleset[[1,i,1]],mlist22[[i,j,k]]},{j,Length[matrixlist[[1]]]},{k,Length[Transpose[matrixlist[[1]]]]},{i,Length[sampleset[[1,All,1]]]}];
(*interpolate both time series*)

tmat1=model*Map[Interpolation[#,InterpolationOrder->1][t]&,tseries1,{2}];tmat2=model*Map[Interpolation[#,InterpolationOrder->1][t]&,tseries2,{2}];tmat11=model*Map[Interpolation[#,InterpolationOrder->1][t]&,tseries11,{2}];tmat12=model*Map[Interpolation[#,InterpolationOrder->1][t]&,tseries12,{2}];tmat21=model*Map[Interpolation[#,InterpolationOrder->1][t]&,tseries21,{2}];tmat22=model*Map[Interpolation[#,InterpolationOrder->1][t]&,tseries22,{2}];(*solve ode's to get solutions corresponding to perturbation of a parameter*)

odes1=Join[Table[D[funs[[i]],t]==Sum[tmat1[[i,j]]*funs[[j]],{j,dim+1}],
{i,2,dim+1}],
Table[funns[[i]][init]==sampleset[[i,1,2]],{i,dim}]];odes2=Join[Table[D[funs[[i]],t]==Sum[tmat2[[i,j]]*funs[[j]],{j,dim+1}],
{i,2,dim+1}],
Table[funns[[i]][init]==sampleset[[i,1,2]],{i,dim}]];odes11=Join[Table[D[funs[[i]],t]==Sum[tmat11[[i,j]]*funs[[j]],{j,dim+1}],
{i,2,dim+1}],
Table[funns[[i]][init]==sampleset[[i,1,2]],{i,dim}]];odes12=Join[Table[D[funs[[i]],t]==Sum[tmat12[[i,j]]*funs[[j]],{j,dim+1}],
{i,2,dim+1}],
Table[funns[[i]][init]==sampleset[[i,1,2]],{i,dim}]];
odes21=Join[Table[D[funs[[i]],t]==Sum[tmat21[[i,j]]*funs[[j]],{j,dim+1}],
{i,2,dim+1}],
Table[funns[[i]][init]==sampleset[[i,1,2]],{i,dim}]];odes22=Join[Table[D[funs[[i]],t]==Sum[tmat22[[i,j]]*funs[[j]],{j,dim+1}],
{i,2,dim+1}],
Table[funns[[i]][init]==sampleset[[i,1,2]],{i,dim}]];(*compute the residuals w.r.t. real data points*)

sols1=Table[NDSolve[odes1,funns,{t,init,end}],{k,Length[odes1]}];
functions1=funns/.sols1[[1]];
residual1=Norm[Table[functions1[[1,i]][t],{i,dim},{t,sampleset[[1,All,1]]}]-sampleset[[All,All,2]],"Frobenius"];sols2=Table[NDSolve[odes2,funns,{t,init,end}],{k,Length[odes2]}];
functions2=funns/.sols2[[1]];
residual2=Norm[Table[functions2[[1,i]][t],{i,dim},{t,sampleset[[1,All,1]]}]-sampleset[[All,All,2]],"Frobenius"];sols11=Table[NDSolve[odes11,funns,{t,init,end}],{k,Length[odes11]}];
functions11=funns/.sols11[[1]];
residual11=Norm[Table[functions11[[1,i]][t],{i,dim},{t,sampleset[[1,All,1]]}]-sampleset[[All,All,2]],"Frobenius"];sols12=Table[NDSolve[odes12,funns,{t,init,end}],{k,Length[odes12]}];
functions12=funns/.sols12[[1]];
residual12=Norm[Table[functions12[[1,i]][t],{i,dim},{t,sampleset[[1,All,1]]}]-sampleset[[All,All,2]],"Frobenius"];
sols21=Table[NDSolve[odes21,funns,{t,init,end}],{k,Length[odes21]}];
functions21=funns/.sols21[[1]];
residual21=Norm[Table[functions21[[1,i]][t],{i,dim},{t,sampleset[[1,All,1]]}]-sampleset[[All,All,2]],"Frobenius"];sols22=Table[NDSolve[odes22,funns,{t,init,end}],{k,Length[odes22]}];
functions22=funns/.sols22[[1]];
residual22=Norm[Table[functions22[[1,i]][t],{i,dim},{t,sampleset[[1,All,1]]}]-sampleset[[All,All,2]],"Frobenius"];(*take the average as a derivative w.r.t. a parameter*){1/7*1/5*(residual1-residual2)/(2*δ1),1/7*1/5*(residual11-residual12-residual21+residual22)/(4*δ1*δ2)}]

Eigenvalues for the current data

hess=Table[hessian[m,l][[2]],{m,7},{l,7}](**)

Eigenvectors[hess]//MatrixForm(*cotyledon control*)

ListPlot[Eigenvalues[hess],PlotRange->All,AspectRatio->1/2]

perturbations

news=Table[vary[pos[[i]],matrixlist,model,sampleset],{i,Length[pos]}];

pic1=Plot[news[[1,2,1,1]][t],{t,5,9},PlotRange->All,PlotStyle->{ColorData["Rainbow"][1/Length[sampleset]]},PlotLabel->text["parameter varied",Medium],Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},FrameTicksStyle->Directive[10],FrameLabel->{text["days",Medium],text["concentration",Medium]}];

pic2=Show[Table[newplotmean[Range[5,9],bigdata[[i]],ColorData["Rainbow"][i/Dimensions[data][[1]]],0.015,Opacity[0.],ColorData["Rainbow"][i/Dimensions[data][[1]]],.06,""],{i,dim}]];

pic3=Show[Table[Plot[result[[2,k]][t],{t,5,9},PlotRange->All,PlotStyle->{ColorData["Rainbow"][k/Length[sampleset]]}],{k,Length[sampleset]}],Axes->None];

pic41=Show[Table[Plot[news[[1,2,k,1]][t],{t,5,9},PlotRange->All,Axes->None,PlotStyle->{ColorData["Rainbow"][1/Length[sampleset]]}],{k,10}]];

pic42=Show[Table[Plot[news[[1,2,k,2]][t],{t,5,9},PlotRange->All,Axes->None,PlotStyle->{ColorData["Rainbow"][2/Length[sampleset]]}],{k,10}]];

pic43=Show[Table[Plot[news[[1,2,k,3]][t],{t,5,9},PlotRange->All,Axes->None,PlotStyle->{ColorData["Rainbow"][3/Length[sampleset]]}],{k,10}]];

pic44=Show[Table[Plot[news[[1,2,k,4]][t],{t,5,9},PlotRange->All,Axes->None,PlotStyle->{ColorData["Rainbow"][4/Length[sampleset]]}],{k,10}]];

pic45=Show[Table[Plot[news[[1,2,k,5]][t],{t,5,9},PlotRange->All,Axes->None,PlotStyle->{ColorData["Rainbow"][5/Length[sampleset]]}],{k,10}]];

pic46=Show[Table[Plot[news[[1,2,k,6]][t],{t,5,9},PlotRange->All,Axes->None,PlotStyle->{ColorData["Rainbow"][6/Length[sampleset]]}],{k,10}]];

pic47=Show[Table[Plot[news[[1,2,k,7]][t],{t,5,9},PlotRange->All,Axes->None,PlotStyle->{ColorData["Rainbow"][7/Length[sampleset]]}],{k,10}]];

fig22=Show[pic1,pic2,pic3,pic41,pic42,pic43,pic44,pic45,pic46,pic47]

Export["D:\\My_submissions_2012\\bmc_article\\images\\Revised_Images\\Figure9.eps",Show[GraphicsRow[{fig22,fig11}],ImageSize->800]]

Additional study: what if the data was perturbed?

noisysamplesets=Table[Map[{#[[1]],(1+0.1*RandomVariate[NormalDistribution[]])*#[[2]]}&,sampleset,{2}],{i,10}];

perturbations=Table[timevariant[noisysamplesets[[i]],model1],{i,10}];

pic1=Plot[perturbations[[1,2,1]][t],{t,5,9},PlotRange->All,PlotStyle->{ColorData["Rainbow"][1/Length[sampleset]]},PlotLabel->text["data varied",Medium],Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},FrameTicksStyle->Directive[10],FrameLabel->{text["days",Medium],text["concentration",Medium]}];

pic2=Show[Table[newplotmean[Range[5,9],bigdata[[i]],ColorData["Rainbow"][i/Dimensions[data][[1]]],0.015,Opacity[0.],ColorData["Rainbow"][i/Dimensions[data][[1]]],.06,""],{i,dim}]];

pic3=Show[Table[Plot[result[[2,k]][t],{t,5,9},PlotRange->All,PlotStyle->{ColorData["Rainbow"][k/Length[sampleset]]}],{k,Length[sampleset]}],Axes->None];

pic4=Table[Plot[perturbations[[i,2,k]][t],{t,5,9},PlotRange->All,PlotStyle->{ColorData["Rainbow"][k/Length[sampleset]]}],{i,10},{k,7}];

fig11=Show[pic1,pic2,pic3,pic4]

Dimensions[sampleset]