SIMULATION OF KINETICS, GIVEN: MECHANISM RATE CONSTANTS AND STARTING REACTANT CONCENTRATIONS aims: 1. tutorial on how to simulate reaction kinetics given mechanism; numerical integration of rate equations. 2. illustration of kinetic properties of several typical mechanisms: Michaelis-Menten Briggs-Haldane one-step equilibrium binding two-step equilibrium binding LDH (steady state) lysozyme chymotrypsin (amide and ester substrates) login at gilmour c tmp # go to a scratch directory cp ~jrupley/kinetics/run.template/xrun . # copy program # so you can make changes # note: if you are logged into gilmour through a # non-X terminal, then copy "run" instead of "xrun" vi xrun # edit at places marked "MAKE CHANGES": # mechanism and rate constants # starting concentrations # time range and time incrementation xrun # run with changes # and display results c ~jrupley/kinetics/run.examples # go to directory # with various results files xgobi basename # display a file of results; # there are several, for various mechanisms # with basename lys_1, ldh_50_1, etc; # to see how they were generated, see the files # run* where * is lys, ldh_50, etc # or see below in this file, for a # minimal description MICHAELIS-MENTEN ANALYSIS OF INITIAL RATE DATA DERIVED FROM SIMULATION aims: 1. illustration of use of simulation of reaction kinetics to obtain values for use in other analyses. 2. illustration of fit of pairs of values of rate and concentration to a model mechanism, here the Michaelis-Menten 3. illustration of the appropriate weighting of data in a reciprocal fit login at gilmour c tmp # go to a scratch directory cp ~jrupley/kinetics/run.rate/rate . # copy programs so you cp ~jrupley/kinetics/run.rate/mk_conc_rate . # can edit your copies kS # start S # the program "rate" controls the reaction simulation and # evaluates the reaction rate, v=v_o(0-10%)/E_o; # if you may want to change reaction rate constants # or otherwise revise the mechanism, # then you need to edit rate: !vi rate # (ignore this step for demo) # the program "mk_conc_rate" runs the program "rate" # for various E,S concentrations and massages the output # to extract values for the matrix of [S],v pairs; # if you want to change the list of E,S concentrations, # then you need to edit mk_conc_rate: !vi mk_conc_rate # (ignore this step for demo) # run the simulation (this takes a moment) # and read the [S],v pairs into an S matrix xmat<-matrix(as.numeric(unix("mk_conc_rate")),ncol=2,byrow=T, dimnames=list(NULL,c("S","v"))) # display results and calculate MM parameters; # NOTE: unweighted least squares fitting of data to # the Lineweaver-Burke equation (reciprocal plot), # as done below to evaluate the MM parameters, # is quick and easy, but a no-no x11() # startup the graphics driver # x11() or iris4d() plot(xmat) # plot all simulation data plot(xmat[xmat[,"S"]<1.e-4,]) # plot subset of data == omit high [S] plot(1/xmat[xmat[,"S"]<1.e-4,],xlab="1/S",ylab="1/v") # reciprocal plot # and label axes # for subset of data plot(1/xmat,xlab="1/S",ylab="1/v") # same for all data xdf<-as.data.frame(xmat) # make data frame from data matrix aa<-lm(1/v ~ I(1/S), data=xdf) # linear least squares fit #aa<-lsfit(1/xmat[,"S"], 1/xmat[,"v"]) # old style lsq fit command cbind(1/xdf,aa$residuals) # display results of fit # note substantial deviations # in residuals # and of kcat, Km from theory cat("k_cat =",1/aa$coef[[1]],"\nK_m =",aa$coef[[2]]/aa$coef[[1]],"\n") bb<-lm(1/v ~ I(1/S), data=xdf, weights=xdf$v^4) # weighted lsq fit cbind(1/xdf,bb$residuals) # display results of fit # note deviations smaller with wghting cat("k_cat =",1/bb$coef[[1]],"\nK_m =",bb$coef[[2]]/bb$coef[[1]],"\n") xxdf<-xdf[xdf[,"S"]>=1.e-6,] # subset of data == remove low [S] cc<-lm(1/v ~ I(1/S), data=xxdf) # linear least squares fit cbind(1/xxdf,cc$residuals) # display fit of subset of data # note smaller deviations for # unweighted fit cat("k_cat =",1/cc$coef[[1]],"\nK_m =",cc$coef[[2]]/cc$coef[[1]],"\n") dd<-lm(1/v ~ I(1/S), data=xxdf, weights=xxdf$v^4) # weighted subset fit cbind(1/xxdf,dd$residuals) # display wghted fit of subset of data # note deviations about same as for # wghted fit of full set of data cat("k_cat =",1/dd$coef[[1]],"\nK_m =",dd$coef[[2]]/dd$coef[[1]],"\n") graphics.off() LOGIC OF COMPUTATIONS AND SOFTWARE xrun or run: shell script; front end for react; sets up input files for react; runs react; massages output of react, with prep or mk_xfiles; displays massaged output, with ttygraph or xgobi; prep (run or xrun): react output filtered with prep is displayed with ttygraph (for dumb terminal); mk_xfiles (xrun): react output filtered with mk_xfiles is displayed with xgobi (requires X terminal) react: executable from C source; front end for lsoda; parses input files that contain, in a moderately user-friendly format: the reaction mechanism values for rate constants values for starting concentrations of reactants and products time range and time incrementation and miscellaneous control constants creates output file with time profile of concentrations of reactants and products and intermediates; there is no formal documentation for react; the comments in the front end, xrun, should be sufficient; comments by the author are in gilmour:~jrupley/kinetics/Docs, in the form of postings; if you want more, please read the sources lsoda: executable from C source; applies numerical integration to solve ordinary differential equations; documentation for lsoda is in gilmour:~jrupley/kinetics/Docs rate: shell script; a fancy version of run; does what run does (runs react and displays results with ttygraph); also uses the output of react to calculate the halftime and initial rate of reaction (appearance of product); command line options control what is printed and allow setting of E and S initial conentrations; for usage, do "rate -?" mk_concn_rate: shell script; runs rate for a set of E,S concentrations and massages output xgobi: X-based package - display multivariate data & look for correlations; quick viewing of graphs, Ramachandran plot, etc. S: statistical analysis; graphics, for display and for plots suitable for publication; computation, especially matrix algebra; programming language; UNIX interface. DESCRIPTION OF FILES IN kinetics/run.examples cht_amide.r #values for chymotrypsin from hammes, enzyme catalysis and regulation #for hydrolysis of N-furylacryloyl-L-Trp amide (page 129) e = s = 1.e-4 M E + S -> X1 <--> X2 --> X3 + P1 -> E + P2 rate constant reaction 1 6.000000e+07 E + S -> X1 2 1.000000e+04 X1 -> E + S 3 1.500000e+00 X1 -> X2 4 3.000000e+01 X2 -> X1 5 4.300000e-02 X2 -> X3 + P1 6 5.000000e+01 X3 -> E + P2 cht_ester.r cht_ester1.r #values for chymotrypsin from hammes, enzyme catalysis and regulation #for hydrolysis of N-acetyl-L-Trp ethyl ester (Table 7-3, page 129) #kon for X1 is estimated cht_ester: e = s = 1.e-4 M cht_ester1: e = 1.e-3 M; s = 1.e-2 M E + S -> X1 --> X2 + P1 -> E + P2 rate constant reaction 1 1.000000e+08 E + S -> X1 2 2.100000e+05 X1 -> E + S 3 3.500000e+01 X1 -> X2 + P1 4 8.400000e-01 X2 -> E + P2 ldh_0_1.r ldh_0_2.r #values for ldh from borgmann, moon, and laidler, B 13, 5152 (1974) #for beef heart ldh at 0 C #A=NAD+, B=lactate, X=pyruvate, Y=NADH _1: e = a = 1.e-3 M; b = 1.e-2 M _2: e = a = b = 1.e-3 M rate constant reaction 1 4.000000e+05 E + A -> EA 2 3.100000e+01 EA -> E + A 3 1.510000e+03 EA + B -> EAB 4 8.930000e-01 EAB -> EA + B 5 3.600000e+01 EAB -> EY + X 6 9.260000e+07 EY + X -> EAB 7 2.530000e+00 EY -> E + Y 8 2.980000e+06 E + Y -> EY ldh_50_1.r ldh_50_2.r #values for ldh from borgmann, moon, and laidler, B 13, 5152 (1974) #for beef heart ldh at 50 C #A=NAD+, B=lactate, X=pyruvate, Y=NADH _1: e = a = b = 1.e-3 M _2: e = 1.e-4 M; a = b = 1.e-3 M rate constant reaction 1 1.450000e+06 E + A -> EA 2 1.270000e+03 EA -> E + A 3 2.590000e+06 EA + B -> EAB 4 3.220000e+04 EAB -> EA + B 5 2.200000e+03 EAB -> EY + X 6 1.750000e+06 EY + X -> EAB 7 4.330000e+02 EY -> E + Y 8 1.280000e+08 E + Y -> EY ldh_50_rev_1.r ldh_50_rev_2.r #values for ldh from borgmann, moon, and laidler, B 13, 5152 (1974) #for beef heart ldh at 50 C #A=NAD+, B=lactate, X=pyruvate, Y=NADH _1: e = x = y = 1.e-3 M _2: e = 1.e-4 M; x = y = 1.e-3 M rate constant reaction 1 1.450000e+06 E + A -> EA 2 1.270000e+03 EA -> E + A 3 2.590000e+06 EA + B -> EAB 4 3.220000e+04 EAB -> EA + B 5 2.200000e+03 EAB -> EY + X 6 1.750000e+06 EY + X -> EAB 7 4.330000e+02 EY -> E + Y 8 1.280000e+08 E + Y -> EY lys_1.r lys_2.r #values for HEW lysozyme, from banerjee et al, jbc 250, 4355 (1975) #for hydrolysis of beta(1->4)-linked hexamer of N-acetylglucsamine #at pH = 6.3, 25 deg #kon for EU1 and ESB1 are estimated _1: e = s = 1.e-5 M _2: e = 1.e-5 M; s = 1.e-4 M rate constant reaction 1 2.000000e+01 EU2 -> EU1 2 2.000000e+02 EU1 -> EU2 3 6.250000e+03 EU1 -> E + S 4 1.000000e+08 E + S -> EU1 5 1.000000e+08 E + S -> ESB1 6 5.000000e+03 ESB1 -> E + S 7 1.400000e+01 ESB1 -> ESB2 8 5.000000e+00 ESB2 -> ESB1 9 4.600000e-01 ESB2 -> ESG 10 1.300000e-01 ESG -> ESB2 11 3.800000e-02 ESG -> E + P1 + P2 one_step.r #single step equilibrium between E + S and EU1 complex e = s = 1.e-3 M rate constant reaction 1 1.000000e+05 EU1 -> E + S 2 1.000000e+08 E + S -> EU1 two_step.r #two step equilibrium between E + S and EU1 intermediate and EU2 complex e = s = 1.e-3 M rate constant reaction 1 1.000000e+00 EU2 -> EU1 2 1.500000e+02 EU1 -> EU2 3 1.000000e+05 EU1 -> E + S 4 1.000000e+08 E + S -> EU1 MM.r #MM mechanism e = 1.e-5 M; s= 1.e-4 M rate constant reaction 1 1.000000e+08 E + S -> ES 2 1.000000e+02 ES -> E + S 3 1.000000e+00 ES -> E + P Haldane.r #Briggs-Haldane mechanism e = 1.e-5 M; s = 1.e-4 M rate constant reaction 1 1.000000e+08 E + S -> ES 2 1.000000e+00 ES -> E + S 3 1.000000e+02 ES -> E + P