[R-gregmisc-users] SF.net SVN: r-gregmisc: [950] trunk/PathwayModeling/thesispaper
Brought to you by:
warnes
Menu
▾
▴
[R-gregmisc-users] SF.net SVN: r-gregmisc: [950] trunk/PathwayModeling/thesispaper
|
From: <r_b...@us...> - 2006-03-29 00:37:37
|
Revision: 950 Author: r_burrows Date: 2006-03-28 16:37:26 -0800 (Tue, 28 Mar 2006) ViewCVS: http://svn.sourceforge.net/r-gregmisc/?rev=950&view=rev Log Message: ----------- incorporated Greg's edits Modified Paths: -------------- trunk/PathwayModeling/thesispaper/R/plotMSD16.R trunk/PathwayModeling/thesispaper/discussion.Snw trunk/PathwayModeling/thesispaper/introduction.Snw trunk/PathwayModeling/thesispaper/methods.Snw trunk/PathwayModeling/thesispaper/paper.Snw trunk/PathwayModeling/thesispaper/results.Snw Modified: trunk/PathwayModeling/thesispaper/R/plotMSD16.R =================================================================== --- trunk/PathwayModeling/thesispaper/R/plotMSD16.R 2006-03-25 14:42:14 UTC (rev 949) +++ trunk/PathwayModeling/thesispaper/R/plotMSD16.R 2006-03-29 00:37:26 UTC (rev 950) @@ -9,11 +9,14 @@ } temp3 <- temp3/10 ymax <- max(temp1[,2],temp2[,2],temp3[,2]) - plot(log10(temp1[,1]),temp1[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),xlab="log10(#evals)",ylab="mean squared residual") +# plot(log10(temp1[,1]),temp1[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),xlab="log10(#evals)",ylab="mean squared residual") + plot(temp1[,1],temp1[,2],type='l',log="x",ylim=c(5000,ymax),xlim=c(1,10^7),xlab="#evals",ylab="mean squared residual") par(new=T) - plot(log10(temp2[,1]),temp2[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=2,lty=2,xlab="",ylab="") +# plot(log10(temp2[,1]),temp2[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=2,lty=2,xlab="",ylab="") + plot(temp2[,1],temp2[,2],type='l',log="x",ylim=c(5000,ymax),xlim=c(1,10^7),col=2,lty=2,xlab="",ylab="") par(new=T) - plot(log10(temp3[,1]),temp3[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=3,lty=3,xlab="",ylab="") - legend(4.5,ymax,c("1-comp","all-comp","NKC"),lty=1:3,col=1:3) +# plot(log10(temp3[,1]),temp3[,2],type='l',ylim=c(5000,ymax),xlim=c(0,7),col=3,lty=3,xlab="",ylab="") + plot(temp3[,1],temp3[,2],type='l',log="x",ylim=c(5000,ymax),xlim=c(1,10^7),col=3,lty=3,xlab="",ylab="") + legend(10^4.5,ymax,c("1-comp","all-comp","NKC"),lty=1:3,col=1:3) } Modified: trunk/PathwayModeling/thesispaper/discussion.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/discussion.Snw 2006-03-25 14:42:14 UTC (rev 949) +++ trunk/PathwayModeling/thesispaper/discussion.Snw 2006-03-29 00:37:26 UTC (rev 950) @@ -6,9 +6,9 @@ We are interested in assessing the usefulness of Markov chain Monte Carlo (MCMC) methods for the fitting of models of metabolic pathways. Fitting -these models can be difficult because the pathways are described by sets +metabolic pathway models can be difficult because pathways are described by sets of reaction rate equations with overlapping sets of parameters. We have -found that the MCMC algorithms handled this situation well, and produced +found that the MCMC algorithms handle this situation well, and produce reasonable joint probability densities for the model parameters. This output can be used for the estimation of confidence intervals for the parameters and the @@ -37,32 +37,22 @@ All three algorithms produced very similar results but the execution times varied considerably depending on the size of the model and the algorithm. -The all-component Metropolis and NKC algorithms converged fairly easily -but the 1-component algorithm was quite slow. This is undoubtably due to +The all-component Metropolis and NKC algorithms converged quickly +but the 1-component algorithm was quite slow to converge. This is undoubtably due to the high correlations between some parameters which is evident in Figure~\ref{scatter}. The 1-component algorithm has trouble sampling from this density because it is restricted to movements parallel to the axes. This is less of a problem for the all-components and NKC algorithms since they can move diagonally. -The similarity of the results obtained with the different algorithms is -apparent from Table~\ref{MSq}. The mean squared residuals do not vary +Each of the algorithms produced essentially identical results. The mean squared residuals do not vary much, nor do the calculated values for $R_{adj}^2$. These last values are indicative of a useful fit of the models to the data. These results also demonstrate the usefulness of MCMC simulation. The derived probability densities can be used to find good estimates for the parameters as well as reveal the correlations between the parameters. -The result is a well-fitted model with realistic confidence intervals (Figure~\ref{fitted}). +The result is a well-fitted model with realistic confidence intervals (Figure~\ref{fits}). -It was of interest to see if the calculated parameter values could be used -to infer values for the underlying rate constants by fitting a simple -linear regression model. For this, expressions for the parameters as -polynomial functions of the rate constants were found by systematically -varying the rate constants according to a face-centered cubic experimental -design. It was found that the parameters are rather complex functions of -the rate constants, too complex to allow useful inferences concerning the -rate constants to be made given the parameter values using this method. - In this paper we have developed models of linear sequences of irreversible biochemical reactions and shown that several Markov chain methods can be used to generate useful fits of the models to simulated data. We also have Modified: trunk/PathwayModeling/thesispaper/introduction.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/introduction.Snw 2006-03-25 14:42:14 UTC (rev 949) +++ trunk/PathwayModeling/thesispaper/introduction.Snw 2006-03-29 00:37:26 UTC (rev 950) @@ -33,32 +33,31 @@ to be the result of reactions with glycolytic intermediates such as methylglyoxal. Nishikawa et al. \cite{Nishikawa00} have hypothesized that this is caused by the increased catabolism of glucose in diabetes -which leads to an increase of superoxide. Superoxide is known to inhibit the enzyme +which leads to an increase of superoxide anion ($O_2^-$). Superoxide is known to inhibit the enzyme glyceraldehyde-3-phosphate dehydrogenase (G3PDH) by about 50\% under these conditions. The inhibition of this enzyme then increases the concentration of glyceraldehyde-3-phosphate (G3P) and dihydroxyacetone phosphate (DHAP) which breaks down to form methylglyoxal, an active aldehyde which readily reacts with and inactivates proteins. To test this hypothesis we need to be able to answer -questions such as: Is the observed inhibition of this enzyme by 50\% sufficient +questions such as: Is the observed inhibition of glyceraldehyde-3-phosphate dehydrogenase by 50\% sufficient to explain the observed increase in methylglyoxal concentration? If not, then what other mechanisms might be operative? What is the most effective way to modulate these reactions? We address these questions by developing analytical methods that will allow us to characterize these types of biochemical systems at the molecular level. -We will use a Bayesian modeling approach with fitting via Markov chain -Monte Carlo (MCMC) simulation +We have used a Bayesian modeling approach with fitting via Markov chain +Monte Carlo (MCMC) methods \cite{Metropolis53,Gilks95} to estimate the kinetic parameters of the enzymes in -a metabolic pathway. With this method we can combine known information -about the pathway with newly-generated experimental data to produce +a metabolic pathway. With this approach we can combine known information +about the pathway with newly-generated experimental data to produce the joint -probability densities for all of the parameters. This joint +probability density of all of the parameters. This joint probability density can then be used to answer various questions about the behavior of the system. -We will explore the utility of the the Bayesian approach using a -simulation technique: First we will generate data for a model with -known parameters, and then we will fit this data using Bayesian +We have explored the utility of the the Bayesian approach using a +simulation technique: First we generated data for hypothetical sequences of enzymatic reactions, and then we fit models to this data using Bayesian methods. The biochemical data is generated using Gillespie's stochastic simulation algorithm \cite{Gillespie77}. This algorithm is @@ -75,13 +74,12 @@ systems. We simulate perturbation experiments. In these -experiments, the concentration of a reactant in a pathway is +experiments, the concentration of the first reactant in a pathway is transiently increased and the concentrations of all the reactants is monitored as they return to a steady state. This type of experiment is quite feasible with cultured animal cells or suspensions of microbial cells such as yeast, and it yields information on the -behavior of the pathway as a single system. The MCMC methods described -here are very useful for analyzing this type of data. +behavior of the pathway as a single system. MCMC simulation was found to be a very useful method for fitting systems of equations to data. It is capable of fitting several Modified: trunk/PathwayModeling/thesispaper/methods.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/methods.Snw 2006-03-25 14:42:14 UTC (rev 949) +++ trunk/PathwayModeling/thesispaper/methods.Snw 2006-03-29 00:37:26 UTC (rev 950) @@ -4,11 +4,9 @@ $Id$ \end{verbatim} - \subsection{The Data} + \subsection{The Pathway} - The data is generated with Gillespie's - algorithm for stochastic simulation of chemical reactions - \cite{Gillespie77}. In this paper we present the results for a + While we examined several pathway structures, in this paper we present the results for a 5-reaction sequence that is sampled at 12, 16, or 25 time points with 3 replicates at each time point. <<data, echo=F, eval=F>>=1 @@ -22,7 +20,7 @@ times16 <- c(16:19,20.5,21,seq(22,40,2)) # times25 <- c(16:19,20.5,21:40) # generating the MCMC input files, e.g., input16.dat -data <- getPaperData(rawdata,20,times16,npoints=3, varVector=c(50,70,60,40,40)) +data <- getPaperData(rawdata,20,times16,npoints=3, varVector=c(50,50,50,50,50)) rates <- getPaperRates(data,1:12,13:48) temp <- matrix(0,nrow=16,ncol=11) for (i in 1:16) temp[i,] <- (rates[3*i,]+rates[3*i-1,]+rates[3*i-2,])/3 @@ -43,18 +41,22 @@ R5 &\yields^{k_{14}}& sink \end{chemarray*} + \subsection{The Experiment} + To estimate the kinetic parameters of the enzymes in a sequence of reactions, we need to find expressions of the form \[ \text{reaction velocity} = f(\mathit{[substrate],[enzyme]}) \] - These expressions can be found by perturbing a reaction network + These expressions are found by perturbing a reaction network and observing the metabolite concentrations as the system relaxes - back to a steady state. The metabolite concentrations are - measured, and the velocities are calculated as the - slopes of the curves of metabolite concentration as functions of + back to a steady state. Specifically, the Gillespie simulation was run with the initial values for the reactant concentrations until a steady state was achieved. The steady state appeared to have been reached by \textit{time} = 15. At \textit{time} = 20 the concentration of $R1$ was increased to 10000. At the indicated time points the metabolite concentrations are + measured, and each reaction velocity is calculated as the + slope of the curve of metabolite concentration as a function of time. + \subsection{The data} + For the 5-reaction sequence of reactions the perturbation of $R1$ at $time = 20$ results in the time courses plotted in Figure~\ref{pulse}. @@ -79,12 +81,12 @@ concentrations and 5 values for the estimated reaction rates. The data at each time point was generated by sampling 3 values for each reactant from a normal distribution centered on the Gillespie - output and with a standard deviation of 50. Three sets of time + output and with a standard deviation of 50, i.e., the measurement error is 50. Three sets of time points containing 12,16, and 25 points were used. The reaction velocities were estimated with the \textit{smooth.spline()} function in \textit{R}~\cite{R}. - \subsection{Biochemical Models} + \subsection{The Biochemical Model} Individual reactions were fit using the Michaelis-Menten equation \cite{Michaelis13} for individual enzymes. This is a @@ -183,10 +185,9 @@ \] where the $\mu_i$ are the estimates of the parameter values before any data are collected. With this distribution, the $d_i$s will - have a mode of $\mu_i$ and large values will be unlikely. In this situation, $\mu_i \equiv 50$ for + have a mode of $\mu_i$ and large values will be unlikely. For this situation, we selected $\mu_i \equiv 50$ for the rate constants and $\sigma \equiv 2$ for the data standard - deviation seem appropriate, and are used in the - priors. + deviation. The likelihood is the probability of the estimated reaction velocities $v_j$ given the model. @@ -217,7 +218,7 @@ \subsubsection{MCMC sampling algorithms} - The performance of a MCMC simulation + The efficiency of a MCMC simulation \cite{Metropolis53,Hastings70} depends heavily on the method used to find candidate points to add to the Markov chain [$q(\cdot |\cdot)$]. @@ -230,7 +231,7 @@ The component-wise Metropolis and the all-components Metropolis algorithms both operate on single Markov chains. The component-wise algorithm (Figure~\ref{1comp}) generates candidate points by - changing the value of only one component (dimension) of the current point at + changing the value of only one component (dimension) of the current state at each iteration by sampling from a univariate normal distribution centered at the current point. The all-components algorithm (Figure~\ref{allComp}) changes all the components @@ -257,7 +258,7 @@ The component-wise Metropolis algorithm has the advantage of simplicity but may move very slowly if the components are highly correlated. The all-components Metropolis avoids the problems with - correlation by using a covariance matrix and updating two or more + correlation, if an appropriate covariance matrix is supplied, by updating all components simultaneously but may perform poorly in high dimensions. A problem with both of these algorithms is that if the component distributions have two or more modes @@ -268,7 +269,7 @@ The Normal Kernel Coupler (NKC) algorithm operates on multiple chains, so that there are several ``current'' states. The NKC uses a normal kernel density estimate as the proposal distribution for choosing candidate - states. The kernel density estimate is generated using the entire set of + states. The estimate is generated using the entire set of current points. Since the algorithm uses the entire set of current points it can move over areas of low probability in the parameter space, especially if the user takes @@ -277,10 +278,10 @@ \section{Convergence} The MCMC simulations are run until the Markov chains have reached - stable distributions as assessed by \textit{mcgibbsit()} + stable distributions as assessed by the \textit{mcgibbsit()} algorithm \cite{Warnes00}. \textit{mcgibbsit()} calculates the number of iterations necessary to estimate a user-specfied quantile $q$ to within $\pm r$ with probability $s$, i.e., \textit{mcgibbsit()} indicates when the MCMC sampler has run long enough to provide good confidence interval estimates. The defaults, which are used - in this paper, are $q=0.025$, $r=0.0125$, and $s=0.95$. + in this paper, are $q=0.025$, $r=0.0125$, and $s=0.95$. These values generate estimates of the 2.5\% and 97.5\% quantiles to $\pm$ 0.0125 quantiles with probability 0.95. Modified: trunk/PathwayModeling/thesispaper/paper.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/paper.Snw 2006-03-25 14:42:14 UTC (rev 949) +++ trunk/PathwayModeling/thesispaper/paper.Snw 2006-03-29 00:37:26 UTC (rev 950) @@ -34,16 +34,16 @@ \address{University of Rhode Island, Kingston, RI} \begin{abstract} -\textit{We have examined} the usefulness of \textit{Bayesian statistical methods} for the modeling +We have examined the usefulness of Bayesian statistical methods for the modeling of biochemical reactions. With simulated data, it is -shown that mechanistic models can \textit{effectively} be fit to sequences of enzymatic -reactions. These methods have the advantages of being relatively easy +shown that these methods can effectively fit mechanistic models of sequences of enzymatic +reactions to experimental data. These methods have the advantages of being relatively easy to use and producing probability distributions for the model -parameters rather than point estimates, \textit{allowing more informative statistical inferences to be drawn}. +parameters rather than point estimates, allowing more informative inferences to be drawn. Three Markov chain Monte Carlo algorithms are used to fit models to data from a -sequence of 4 enzymatic reactions. These algorithms +sequence of 4 enzymatic reactions. The algorithms are evaluated with respect to the goodness-of-fit of the fitted models and the time to completion. It is shown that the algorithms produce essentially the same Modified: trunk/PathwayModeling/thesispaper/results.Snw =================================================================== --- trunk/PathwayModeling/thesispaper/results.Snw 2006-03-25 14:42:14 UTC (rev 949) +++ trunk/PathwayModeling/thesispaper/results.Snw 2006-03-29 00:37:26 UTC (rev 950) @@ -4,21 +4,10 @@ $Id$ \end{verbatim} -The joint probability density for the parameters of a model of a -metabolic pathway has been estimated using MCMC simulation. Three -algorithms were used for the simulations and they converged to similar -distributions. - +All three algorithms converged to similar distributions and produced essentially the same metabolic pathway model. The marginal distributions for the parameters were mound-shaped -(Figure~\ref{densities}). They -generally were skewed to the right and had thicker -tails than normal distributions. The modes of the distributions -turned out to be poor estimates of the maximum likelihood values for -the parameters. The maximum likelihood values were all larger than -the modes, and the difference increased as the position of the reaction -in the sequence increased. Correlated pairs of parameters -(e.g., $d_1$ --$d_2$ and $d_3$ --$d_4$, Figure~\ref{scatter}) were similar relative distances from their -modes. The significance of all this is not clear at this point. +(Figure~\ref{densities}), were skewed to the right and had thicker +tails than normal distributions. <<fig4,echo=F,eval=F>>=3 source("R/paperSSQ.R") @@ -84,7 +73,7 @@ points} \label{converged} \end{figure} -There is some narrowing of the distributions as the number of points +We see some improved precision as the number of points increases from 16 to 25 but it is not pronounced. Overall, the reduction in width varied from 18-fold for $d_9$ and 9-fold for $d_1$ to 1.5-fold for $d_4$. @@ -114,11 +103,11 @@ \caption[Bivariate scatter plots of the parameter distributions]{Bivariate scatter plots of the parameter distributions for the 5-reaction model found with the all-components Metropolis algorithm (upper triangle); correlation coefficients (lower - triangle).} + triangle). The red lines indicate the maximum likelihood estimates of the parameters.} \label{scatter} \end{figure} -The usefulness of the probability density for inference was assessed +The value of the probability density for inference was assessed graphically. The probability density was used to find the maximum likelihood estimate for the model parameters and the 95\% confidence intervals. The fits of the resulting models to the 16-point data set @@ -203,12 +192,7 @@ @ \begin{figure} \centering -% \subfloat[12 points]{ -% \includegraphics[scale=0.5]{figures/MSR12}} -% \subfloat[16 points]{ \includegraphics[scale=0.5]{figures/tempDir/MSR16} -% \subfloat[25 points]{ -% \includegraphics[scale=0.5]{figures/MSR25}} \caption[Mean squared residuals vs. number of likelihood evaluations]{Mean squared residuals vs. number of likelihood evaluations for the 5-reaction model} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
