INTRODUCTION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Channeling is a hypothesis that states that pairs of enzymes that catalyze consecutive reactions in a biochemical pathway form protein-protein complexes in which a substrate intermediate is directly passed from the active site (a/s) of one enzyme to the a/s of the other without diffusing away into free solution; thus, the substrate will not equilibrate with the pool of substrates in the surrounding medium (Agius and Sherratt, 1997). The modern concept of channeling has emerged largely from experimental observations. Direct transfer of a metabolite might be possible because of close juxtapositioning of two enzymatic a/s in multi-enzyme complexes, which may be formed in vitro and in vivo. The main advantage of channeling would be that it furnishes a high quantity of metabolic flux, compared with what follows from a traditional view of intermediary metabolism (i.e., the free diffusion of intermediates between the enzymes involving an equilibrium with a metabolite pool in solution). The physiological role of metabolite channeling in different metabolic pathways has been widely discussed (Agius and Sherratt, 1997). There is clear evidence of metabolite channeling in the urea cycle (Cheung et al., 1989), the citrate cycle (Sumegi et al., 1993; Shatalin et al., 1999), glycolysis (Clegg and Jackson, 1990), and other metabolic pathways (Agius and Sherratt, 1997).

One of the pathways in which metabolite channeling is thought to occur is glycolysis. Ovadi and Orosz (1997) suggested a new concept for control of glycolysis based on the assumption that the glycolytic enzymes are not distributed homogeneously in the cytoplasm, but rather that they associate with elements of the cell cytoskeleton-forming multi-enzyme complexes. The catalytic activity of the enzymes in such complexes may vary significantly from what was observed in vitro, or in the cytoplasm. As a consequence of formation of multi-enzyme complexes, the direct transfer of metabolites becomes feasible. The changes in the enzymatic activity of glycolytic enzymes when they form complexes were shown experimentally. Vertessy and Ovadi (1987) observed a higher transient rate constant of the catalytic action for glycerol-3-phosphate dehydrogenase (GAPDH) on dihydroxyacetone phosphate in a coupled reaction involving fructose-1,6-bisphosphate aldolase (aldolase) than in the reaction with dehydrogenase alone. The analysis of the isotope dilution experiments led Orosz and Ovadi (1987) to conclude that there is a channeled transfer of glyceraldehyde-3-phosphate (GAP) from aldolase to GAPDH. The new analysis of transient-state kinetic experiments, however, showed that the kinetic behavior of the aldolase/GAPDH reactions is fully consistent with the free-diffusion mechanism of metabolite transfer (Petterson and Petterson, 1999). Such disagreements are common among metabolite channeling studies. Some agree that channeling exists and invoke this idea for explaining the experimental data, whereas others doubt the existence of channeling and explain the experimental results from a classical view of chemical reactions in solution. Theoretical studies of channeling based on structural information about the involved enzymes may be very useful in distinguishing such arguments; theoretical methods serve as a direct tool to show on microscopic level the possibility of channeling (if it exists) and its kinetic advantages.

Brownian dynamics (BD) is a powerful theoretical method to study intermolecular interactions in solution (Ermack and McCammon, 1978). BD models the relative translational and rotational diffusive motion of whole macromolecules under influence of the complicated electrostatic and excluded volume interactions present in biophysical systems (Northrup and Herbert, 1990). This method was applied to the study of electrostatic channeling of oxaloacetate in the fusion protein of malate dehydrogenase with citrate synthase (Elcock and McCammon, 1996) and the channeling of dihydrofolate in a bifunctional enzyme dihydrofolate reductase-thymidylate synthase (Elcock et al., 1996). High values for the efficiency of transfer (i.e., the relative number of BD trajectories for which the substrate reached the a/s of the target enzyme) were obtained. At zero ionic strength, the efficiency of transfer of dihydrofolate was ~95% (Elcock et al., 1996). The channeling efficiency of oxaloacetate was ~45% (Elcock and McCammon, 1996). The strong dependence of the channeling efficiency on the ionic strength was demonstrated in both papers (Elcock and McCammon, 1996; Elcock et al., 1996); as ionic strength increased channeling efficiency decreased (Elcock and McCammon, 1996; Elcock et al., 1996). In both cases, it was the electrostatic field of the enzymes that makes channeling possible (Elcock and McCammon, 1996; Elcock et al., 1996). BD simulations without electrostatic interactions resulted in very low efficiencies (Elcock and McCammon, 1996; Elcock et al., 1996). In both these pioneering BD studies, the substrate (oxaloacetate or dihydrofolate) was modeled as a charged sphere with a radius of 2 Å and charge of 2e. These theoretical BD results, when compared with experimental results, demonstrate that the transfer efficiencies in both cases (the bifunctional dihydrofolate reductase-thymidylate synthase enzyme and the fusion protein of malate dehydrogenase and citrate synthase) were consistent with experimental transient time measurements (Elcock et al., 1997).

BD methods should also be useful for predicting the interactions between enzymes involved in the glycolytic pathway. BD has predicted the interaction between aldolase with an actin filament (Ouporov et al., 1999) and GAPDH with an actin filament (Ouporov et al., 2001). The possible binding modes of the enzymes to an actin filament is based on the electrostatic attraction between the molecules. BD should also be able to predict the interaction between aldolase and GAPDH. With well-predicted aldolase/GAPDH complexes, BD can then follow the transfer of GAP from the a/s of one enzyme to any of the a/s's of any other enzyme. Furthermore, simulations of GAP diffusing out of solution to any of the a/s of one enzyme can be compared with the channeling simulations to see if there is an observable difference between the channeling and diffusion mechanisms.

	COMPUTATIONAL METHODS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Protein models

The x-ray structures of aldolase and GAPDH monomers from human muscle were retrieved from the RCSB Protein Data Bank (entries 1ALD and 3GPD) (Berman et al., 2000). Both enzymes are homotetramers. The complete tetramers were built using the crystallographic symmetry information in the RCSB Protein Data Bank file and the molecular modeling package Insight II (MSI, San Diego, CA). The GAPDH subunit contained nicotinamide adenine dinucleotide (NAD⁺) as a coenzyme, and two molecules of SO₄ (which were replaced by two molecules of inorganic phosphate). Only heavy atoms were included in the simulations (11,056 aldolase atoms and 10,312 GAPDH atoms).

Charge assignments

The program package MacroDox (version 3.0)(Northrup et al., 1997) was used to assign the titratable charges on the proteins, solve the linearized Poisson-Boltzmann equation, and run the various BD simulations; the BD algorithm for this package is detailed in Northrup et al. (1987), and is overviewed along with the charge assignment and Poisson-Boltzmann algorithms in Northrup et al. (1993). Using the atomic coordinates of each model, the charges of titratable amino acids were assigned by applying the Tanford-Kirkwood method with static accessibility modification (Tanford and Kirkwood, 1957; Tanford and Roxby, 1972; Shire et al., 1974; Matthew, 1985) with the MacroDox charge set (Northrup et al., 1997) at pH 7.0, a temperature of 298.15 K, and an ionic strength of 0.1 M. The partial charges for NAD⁺, the inorganic phosphates, and GAP were calculated using the CHARMM force field (Brooks et al., 1983) within the QUANTA program (MSI). The total charges of proteins were +17.9e for aldolase and 17.3e for GAPDH; the total charge for GAP was 2.2e.

Electrostatic potential

After charge assignments, the electrostatic potential about the proteins was determined by numerically solving the linearized Poisson-Boltzmann equation as implemented in the MacroDox program (Northrup et al., 1997). The electrostatic potential around each protein was determined on two lattices of 81 × 81 × 81 nodes each. The resolution of the coarse lattice was 5.25 Å, and the resolution the fine lattice was 1.75 Å.

Brownian dynamics of aldolase-binding GAPDH

The predominant binding mode of aldolase and GAPDH was estimated by 7,300 BD trajectories, each beginning with the center of mass (COM) of aldolase placed randomly on a sphere with a radius radius 110 Å around GAPDH COM. The orientation and angular position of aldolase was chosen randomly. Both molecules were allowed to rotate and translate. When COM of aldolase reached the surface of a sphere 200 Å away from the GAPDH COM, the trajectory was terminated. During each trajectory, the structure of the complex between aldolase and GAPDH with the lowest value of electrostatic interaction energy was saved if this value was <6 kT, where k is the Boltzmann constant and T = 298 K. The simulation took almost 288 CPU hours on a SGI Indy R4400 workstation.

Analysis of the BD identified complexes included: (1) a visual examination of the complexes to determine the general orientations and (2) a statistical analysis to determine how often certain salt bridges formed between the two enzymes. The statistical analysis was performed using MacroDox (Northrup et al., 1997). For purposes of this analysis, the salt bridge was defined to be intermolecular distances between charged residues that were <6 Å.

Brownian dynamics of substrate channeling

One of the complexes representing the dominant binding mode of aldolase to GAPDH was chosen for BD simulations of GAP channeling between the a/s's of the enzyme. According to the reaction mechanism of GAPDH, an inorganic phosphate binds to GAPDH a/s after GAP. For this reason, all inorganic phosphates were removed from GAPDH a/s's. The electrostatic field around the complex was calculated same way as it was for aldolase and GAPDH, but two 101 × 101 × 101 grids of size 3.36 and 1.4 Å were used for the complex. In the channeling BD simulations, GAP started its diffusive motion from one of the a/s of aldolase (subunit A or D) with a constant timestep of 1 ps. The Protein Data Bank file for aldolase (1ALD) contains the model coordinates of GAP heavy atoms in an aldolase monomer a/s. These coordinates were chosen as the initial position. The conformation of GAP was held rigid during the BD simulations. GAP was allowed to diffuse around the aldolase/GAPDH complex. During a BD trajectory, the distances between the carbonyl carbon (C₁) of GAP and the sulfur atoms of Cys 151 from any GAPDH subunit were monitored. If one of these distances became less than 10 Å, binding was assumed to have occurred, and the GAP trajectory and the time passed since the beginning of trajectory were saved for analysis. If GAP never approached any of the GAPDH a/s's, the trajectory was terminated after the substrate moved 300 Å away from the center of mass of the complex. Four hundred thousand BD trajectories for GAP, starting from the a/s of aldolase subunit A, took 178 single CPU hours on a SGI R12000 OCTANE workstation. The same number of trajectories for GAP starting from the a/s of aldolase D subunit took 134 single CPU hours on the same computer. [Note: to achieve 400,000 trajectories it was necessary to combine four 100,000 trajectory simulations each with a different random number seed; this was necessary to insure the randomness of the random force term at the beginning of each trajectory. A simulation with 400,000 trajectories would exceed the random numbers possible from the random number generator available on the workstation. This practice of breaking down the total number of trajectories into a group of simulations was done throughout.]

Similar BD simulations of GAP transference from GAPDH a/s's to aldolase a/s's were performed. The initial position of GAP in each GAPDH a/s was chosen by visual analysis of the final positions of GAP in GAPDH a/s determined in BD simulations of GAP transfer from aldolase to GAPDH. The GAP structure was identical to one used in previous simulations. For the reverse reaction, the distances between the GAP carbonyl carbon and the epsilon amino nitrogen of Lys 229 from each aldolase subunit were monitored. Three hundred thousand BD trajectories of GAP starting from the R and G active sites of GAPDH were performed. Each simulation took approximately 80 single CPU hours on a SGI R12000 OCTANE workstation. The binding criteria, the condition of trajectory termination, and the time step were similar to that described above.

Brownian dynamics of substrate binding from solution

To compare the efficiency of GAP channeling with GAP binding an aldolase/GAPDH complex from solution, a new series BD simulations were performed to mimic GAP binding from solution. These series of solution simulations began with GAP randomly positioned and oriented on the surface of a 3D sphere with a 150 Å radius; the sphere was centered about the center of mass of the aldolase/GAPDH complex. The 150 Å radius was chosen because the electrostatic potential of the aldolase/GAPDH complex has essentially returned to zero this far away from the center of mass of the complex. For comparison, a second similar simulation with a radius of 200 Å was performed; this radius placed the GAP far out into bulk solution. One million, six hundred thousand BD trajectories followed GAP binding to the a/s of GAPDH; 1,200,000 BD trajectories followed GAP binding to the a/s of aldolase. The number of trajectories was chosen to provide reasonable statistical results within a minimal usage of CPU time. Each trajectory began in a different random position. Because each trajectory represents a single GAP molecule diffusing in three dimensions around an aldolase/GAPDH complex; it is not a concentration of GAP in solution. The number of trajectories chosen allowed for enough two-body trajectories to observe a pattern statistically in the results.

To compare GAP binding the aldolase/GAPDH complex with GAP binding the isolated proteins, two more BD solution simulations were performed. One million BD trajectories followed GAP binding to a single GAPDH enzyme, and another 1,000,000 BD trajectories followed GAP binding to a single aldolase enzyme from solution. For these calculations, the electrostatic potentials around GAPDH (with the inorganic phosphates removed) and aldolase were calculated by numerical solving the linearized Poisson-Boltzmann equation using two 101 × 101 × 101 grids, 3.36 and 1.4 Å each. For each of these simulations, GAP started from a sphere with a radius of 150 Å around COM of the isolated protein and at random orientations. The binding criteria and termination condition were identical to BD simulations of GAP transfer from aldolase to GAPDH.

BD simulations provided information about the efficiency of binding (the relative number of successful trajectories, those that reached any a/s of GAPDH or aldolase) and the distribution of transfer times. The saved successful trajectories were analyzed to determine possible pathways of GAP moving between aldolase and GAPDH.

	RESULTS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Aldolase/GAPDH complexes

Seven thousand, three hundred BD trajectories identified 847 complexes (with interaction electrostatic energy less than 6 kT) indicating that the proteins do interact. Statistical analysis revealed that the amino acids most frequently involved in intermolecular contacts are aldolase Lys 341, 288, 317, 321, 316, GAPDH Glu 105, 82, and, 78, and Asp 80 and 63. These residues form salt bridges with oppositely charged amino acids contributing to the stability of complexes. Visual analysis of the 100 most stable complexes revealed two classes of complexes between aldolase and GAPDH. In class I, two adjacent subunits of GAPDH (subunits G/R or H/S) form salt bridges with two subunits of aldolase (A/D or B/C) (Fig. 1). As in the case of rabbit muscle aldolase (Ouporov et al., 1999), positively charged Lys and Arg of human aldolase form positively charged grooves between subunits A/D or B/C. The clusters of negatively charged Asp and Glu on the GAPDH surface are attracted to these grooves. The size of both proteins and their quaternary structures allow for contacts of two subunits of aldolase A/D or B/C with two subunits of GAPDH G/R or H/S. Complexes of class I are typically stabilized by 3 to 4 salt bridges, and the value of the electrostatic interaction energy for these complexes is around 12 kcal/mol. In class II complexes, the proteins contact each other by the charged residues located at the outer edges of subunits exposed to the solvent (Fig. 1). For these complexes, charged amino acids from any GAPDH subunit contact residues from any subunit of aldolase; i.e., one subunit of one protein contacts one subunit of the other protein. Class I complexes are more compact than class II. The a/s's of GAPDH and aldolase are closer for class I complexes; the typical distance between a/s's of aldolase and GAPDH that face each other is ~45 to 60 Å. For these reasons complexes of class I are the potential candidates for the channeling of GAP from an aldolase a/s to GAPDH a/s and back.

View larger version (56K): [in this window] [in a new window]   FIGURE 1   Two classes of complexes between aldolase and GAPDH. Aldolase and GAPDH are presented as ribbons of the C carbon trace. Aldolase appears on the right and GAPDH appears on the left in both classes.

GAP channeling

One of the complexes from class I (Fig. 1) was chosen to model GAP channeling. In this complex, the R and G GAPDH subunits face A and D aldolase subunits. The distances between a/s of the proteins for this complex are: aldolase A to GAPDH G, 47.0 Å; from A to R, 48.0 Å; from D to G, 54.2 Å; from D to R, 56.5 Å. The electrostatic field around the complex shows a mix of positive and negative patches on the surface (Fig. 2); there is no obvious electrostatic channel visible. In spite of no obvious visible electrostatic channel, a significant number of trajectories found their way into any GAPDH a/s. The number of successful trajectories (i.e., those that reached any GAPDH a/s), the efficiency of GAP transfer from aldolase to GAPDH, and the number of successful trajectories are presented in Table 1. Similar information for the reverse reaction (GAP transfer from GAPDH to aldolase) is presented in Table 2. Fig. 3 shows a more detailed structure of the chosen complex and the pathway of the most successful trajectories among the GAPDH a/s's.

View larger version (114K): [in this window] [in a new window]   FIGURE 2   Electrostatic potential for the aldolase/GAPDH complex used for GAP channeling simulations. The red contours represent an isopotential surface where a 1e charge possesses an electrostatic energy equal to 0.5 kcal/mol. The blue isopotential surfaces are for an electrostatic energy of +0.5 kcal/mol. The black lines are the C carbon traces of the two proteins. The complex is in the exact same orientation as Fig. 3 with GAPDH on the left and aldolase on the right.

View larger version (68K): [in this window] [in a new window]   FIGURE 3   Aldolase/GAPDH complex used for GAP channeling simulations. Aldolase (blue) and GAPDH (red) are presented as ribbons of their C carbon traces. The subunits of both proteins are labeled. The green spheres present the trace for shortest trajectory (3.5 ns) from aldolase A a/s to G a/s of GAPDH. The brown spheres are the trace for shortest trajectory (9 ns) from aldolase D a/s to GAPDH R a/s.

First, only in 1.3% of the cases where GAP started from solution did GAP reach any a/s of isolated GAPDH. Second, GAP binding to isolated GAPDH from solution does not show a preference for a particular GAPDH a/s; the number of trajectories that reached any a/s of GAPDH are approximately the same (Table 1, row 1). The analysis of the distribution of the trajectory duration which reached a particular active site of GAPDH (data not shown), unveiled that the distributions are almost identical, indicating that similar kinetics of GAP binding from solution to any a/s of GAPDH. The distribution of successful trajectories among the aldolase a/s's (Table 2, row 1) shows that a/s's of A and B subunits are twice as attractive for GAP than a/s's of C and D subunits. The kinetics of GAP binding to aldolase a/s's was indistinguishable among them (data not shown). The preference for GAP binding of A and B aldolase a/s's compared to C and D a/s's is probably caused by the quaternary structure of this protein. The binding of GAP to the isolated proteins confirmed the absence of bias to particular a/s of GAPDH or aldolase caused by the BD simulation set up (or in other words, all aldolase and GAPDH a/s's were nearly equally attractive for GAP). Thus, simulations of GAP transfer in a complex of aldolase/GAPDH should not show preference for any a/s's unless the quinary structure provides the preference.

For a complex between aldolase and GAPDH, the values for the efficiency of transfer depend on the initial location of GAP in aldolase. If GAP was initially located in aldolase subunit A a/s, 2.6% of the BD trajectories reached GAPDH a/s's. When GAP started its motion from aldolase subunit D a/s, only 1.5% of the trajectories were successful. The efficiency of GAP binding GAPDH beginning from 150 Å away in solution was the same (1.3%), as in a case of binding to single GAPDH molecule (Table 1); when GAP began farther out in the solution (e.g., 200 Å away), the efficiency dropped to 0.87% (data not shown). The drop in efficiency as GAP moves farther out into the solution is expected because GAP is undergoing random diffusion in the bulk of the solution, and only simple chance will steer it in the direction where it would encounter the electrostatic field of the aldolase/GAPDH complex. Thus, there is a slightly higher preference for GAP channeling from aldolase to GAPDH compared to GAP binding from solution. The numbers for successful trajectories among different GAPDH a/s's vary significantly. Almost half of the successful trajectories ended in GAPDH G a/s (Table 1; Fig. 3); whereas, the R a/s of GAPDH was less attractive. The number of trajectories that ended in the GAPDH R a/s increased when GAP started from the aldolase D a/s. Such division of successful trajectories among GAPDH a/s's is definitely caused by the overall geometry of the complex (i.e., quinary structure). For this particular complex, the GAPDH a/s that belongs to the R subunit is hindered by the nearby subunits of aldolase. The pathways that lead to the R subunit a/s are not completely blocked (some trajectories found this a/s), but the electrostatic field of the nearby aldolase subunits and the excluded volume interaction between them make these paths less likely to finish in an a/s. The effect of the geometry of a complex was even stronger than the distance factor between aldolase and GAPDH a/s's; the H and S GAPDH a/s's are farther from either A or D aldolase a/s's, but they attracted more BD trajectories than the R a/s of GAPDH.

For each BD simulation, a distribution function of transfer time (the duration of successful trajectory) was calculated (Fig. 4). These distributions represent the relative number of trajectories that have a duration within particular limits. Fig. 4 shows that the shape and position of maximum for all three distribution functions (for each BD simulation) are approximately the same. On average, a little more than 100 ns are needed for GAP to reach an a/s of GAPDH, whether it starts from aldolase a/s or from solution. Fig. 4 shows that G and R a/s's were reached more quickly than H or S a/s's when GAP started from aldolase A a/s. The fast GAP transfer from aldolase D a/s to GAPDH R a/s is almost an order of magnitude faster than the average transfer time from D to G. 

View larger version (39K): [in this window] [in a new window]   FIGURE 4   Distribution of transfer times for ending in particular active sites of GAPDH. The top represents GAP beginning in aldolase a/s A. The middle represents GAP beginning from aldolase a/s D, and the bottom represents GAP beginning in solution on the surface of a sphere 150 Å away from the COM of the complex.

The results of BD simulations of reverse GAP transfer (from GAPDH to aldolase) are presented in Table 2. The values of the transfer efficiency are in the range from 1 to 1.5%, and are close to the values of efficiency for aldolase to GAPDH transfer. Both GAPDH a/s's showed a high percentage of trajectories ending in aldolase A a/s (72% when GAP started from G a/s of GAPDH and 61% when starting from R a/s). The number of trajectories that ended in other aldolase a/s's was approximately the same. The distribution transfer time (data not shown) demonstrate a bell curve with a maximum around decimal logarithm 5.25 to 5.5. This corresponds to the average time needed for GAP to reach an aldolase a/s, ~200 to 300 ns. Thus, the reverse GAP transfer can occur with a similar efficiency and on a similar time scale as the forward GAP transfer from aldolase to GAPDH.

	DISCUSSION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

One substantial methodological difference between these BD simulations and earlier simulations (Elcock and McCammon, 1996; Elcock et al., 1996) is that the current simulations use an all heavy atoms presentation of the substrate, GAP (in the conformation published in Protein Data Bank entity 1ALD). The earlier channeling simulations treated the substrate as a structureless charged sphere with a radius of 2 Å (Elcock and McCammon, 1996; Elcock et al., 1996). For this reason, the channeling efficiencies calculated in our work are relatively low. The all-atom approach, however, does have several advantages. For example, BD simulations in which the GAP reactive atom set consisted of carbonyl carbon (C₁) and the two other backbone carbons (C₂ and C₃) provided a 10% increase in the number of successful trajectories, indicating that the carbonyl atom of GAP and the sulfur atom of Cys 151 are the closest atoms for ~90% of the successful trajectories. This means that the orientation of GAP in GAPDH a/s's corresponds to the enzymatic mechanism of GAPDH, in which the carbonyl carbon of GAP and the sulfur atom of Cys 151 bind covalently to form thiohemiacetal (Voet and Voet, 1995). In other words, the current BD simulations provided a proper orientation of GAP in GAPDH a/s's.

To compare the effect of the substrate representation (either charged sphere or all heavy atoms), BD simulations of GAP transfer from aldolase to GAPDH without any electrostatic interaction present a pure random diffusion around the enzyme-enzyme complex. Without electrostatic forces the efficiency of GAP transfer was 0.05% (for GAP starting from either aldolase A or D a/s's). The efficiency of GAP binding to GAPDH a/s's from solution without electrostatic forces was 0.03%. These numbers are significantly smaller than the efficiencies without electrostatic interactions published for channeling in the malate dehydrogenase/citrate synthase fusion protein or the dihydrofolate reductase-thymidylate synthase system (Elcock and McCammon, 1996; Elcock et al., 1996) (1% for oxaloacetate and 6% for dihydrofolate). Comparing our results for GAP channeling with the results for oxaloacetate channeling (Elcock and McCammon, 1996) (both molecules are approximately the same size and have similar structures), one can conclude that the all-atom substrate presentation for BD simulations brings at least a 20-fold decrease in the efficiency of transfer. This factor might explain the difference between the GAP transfer efficiency obtained in our work (2.6%) and the efficiency of oxaloacetate transfer (9%; approximated from the work of Elcock and McCammon, 1996) at an ionic strength of 0.1 M.

The high values of transfer efficiencies calculated by Elcock et al. (1996) and Elcock and McCammon (1996) may be explained by the structure of their protein systems. The electrostatic potential around the fusion protein of malate dehydrogenase and citrate synthase calculated at zero solution ionic strength (Elcock and McCammon, 1996, Fig. 4) shows a wide uninterrupted region of positive potential that covers a/s's of both proteins and guides oxaloacetate (which is negatively charged) from one a/s to another. At higher values of ionic strength this positive region of electrostatic potential disappears, but there are many local positive islands that keep oxaloacetate from wandering into solution and facilitate the transfer between a/s's of the studied proteins. For the aldolase/GAPDH complex there is no such steering region of positive potential (GAP is negatively charged) from aldolase to GAPDH a/s's. The local positive patches are separated on the complex surface and do not provide an uninterrupted directed diffusion path from aldolase to GAPDH a/s's or back (Fig. 2).

Just as for the systems studied by Elcock et al. (1996) and Elcock and McCammon (1996) our simulations showed that GAP transfer from GAPDH (R or G GAPDH a/s) to any aldolase a/s (reverse transfer) has similar efficiencies and dynamical characteristics of channeling as observed for GAP transfer from aldolase to GAPDH. More than half of the successful trajectories that began in GAPDH R or G a/s's ended in aldolase A a/s indicating a preference of aldolase A a/s for channeling in this complex. The difference of GAPDH or aldolase a/s's as considered for GAP binding depends on the general geometry of the complex (quinary structure). Another series of BD simulations for a different class I complex produced efficiencies of GAP transfer that were similar to those presented herein; however, the affinity of enzyme a/s's was different.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

BD results show a small preference of GAP channeling from aldolase to GAPDH when compared to binding from solution. Substrate transfer between neighboring a/s's is more favorable when compared to solution binding from the perspective of transfer time and transfer efficiency. Twice as many trajectories reacted when they started from aldolase A a/s; the average transfer time for GAP to reach GAPDH R or G a/s is approximately two times less than the time necessary to reach GAPDH H or S a/s. Unfortunately, this difference is not large enough to be considered as solid theoretical evidence of channeling between these proteins.

The applied method is quite detailed, and represents the physical reality well (an all-atom presentation of substrate, for example). Future improvements that will make results more credible are (1) including hydrophobic interactions, (2) including conformational changes of the proteins, and (3) including the internal mobility of the substrate. Although hydrophobic forces have been shown to be more important for protein folding than for complex formation (Sheinerman et al. 2000, Larsen et al. 1998, Tsai et al. 1997), the spread around average observations of protein-protein interfaces is large (Lo Conte et al., 1999), indicating that there are exceptions where hydrophobic forces may contribute significantly to complex formation. Thus, the inclusion of hydrophobic interactions may allow for the prediction the structure of protein-protein complexes with a higher accuracy. Future studies will explore the possibility of including hydrophobic interactions based in part on the surface area lost as a complex forms. The early stages of interactions between macromolecules in solution is governed by electrostatic interactions which are represented with high degree of accuracy in these BD simulations. Moreover, BD considers all molecules as solid objects, and their internal dynamics are not included. This factor does not allow conformational changes during complex formation. Conformational rearrangements in participating proteins can provide a great deal of complex stability and changes of complex geometry, which may affect the substrate transfer. It may be possible to develop a traditional molecular dynamics protocol for the BD-predicted protein-protein complexes that would allow for minor conformational flexibility without loosing the openness of the active sites. Finally, the internal mobility of the substrate may facilitate its movement across channeling pathways, and speed up the leaving and the entering of active sites. The inclusion of substrate internal mobility would definitely affect the results of the simulation. The only way for the current BD method to include such substrate mobility would be to perform separate BD simulations for different conformations of the substrate and compare the results. In spite of these shortcomings, this initial study does prove that BD is a significant method for following the path of a substrate from one active site to another and the time factors associated with the path, and from these there is a real albeit small chance of channeling.