INTRODUCTION

Top Abstract Introduction Materials & Methods Results Discussion References

Lipases, a class of enzymes belonging to the serine hydrolases, have proven to be efficient for resolving racemic alcohols and esters. Many lipases are heat-stable and they do not need co-factors. They can often be used in organic solvents. It is of general interest to use molecular modeling for predicting the ability of these enzymes to discriminate between the enantiomers of a substrate. With the help of such a method, a suitable lipase can be chosen for a particular substrate. Molecular modeling can also be used for guidance on redesigning the enzyme by site-directed mutagenesis in order to improve or change the capability of the enzyme to accomplish a specific enantioselective transformation.

Several lipase structures have been determined during the last few years (Brady et al., 1990; Winkler et al., 1990; Schrag et al., 1991; Martinez et al., 1992, 1993; van Tilbeurgh et al., 1992, 1993; Grochulski et al., 1993; Schrag and Cygler, 1993; Uppenberg et al., 1994; Hermoso et al., 1996; Lang et al., 1996). Structure-function relationship studies by means of molecular graphics and molecular modeling are thus possible. Several of the lipases have been covalently inhibited by different phosphonates that mimic the tetrahedral intermediate along the reaction pathway. Chiral phosphonate inhibitors have been used in studies on the orientations of the enantiomers of a substrate in the active site. These structures can be used to align tetrahedral intermediates of various substrates for use in molecular modeling.

Attempts to calculate the enantioselectivity using force field methods have been made (DeTar, 1981; Bemis et al., 1992; Norin et al., 1993, 1994; Faber et al., 1994), with good prediction of the fast-reacting enantiomer. To be able to predict which of the two enantiomers is the fast-reacting one, we must understand the detailed reaction mechanism of the enzyme-catalyzed reaction. Furthermore, the binding of the two enantiomers in the active site during the reaction must be predicted. Lipase-catalyzed reactions proceed via the bi-bi ping-pong mechanism (Kraut, 1977). Various strategies can be used to determine the binding of the enantiomers. One strategy is to use chiral inhibitors that mimic the activated complexes of the enantiomers in the active site of the enzyme (Lee et al., 1996). Separate crystallizations of the enzyme with the two enantiomeric inhibitor complexes, followed by 3-D structure determination, will reveal how the enantiomers are oriented in the active site (Ahmed et al., 1994; Grochulski et al., 1994). Another way is to combine enzyme kinetics with molecular modeling to elucidate the binding orientation of the enantiomers (Norin et al., 1993; Kazlauskas, 1994; Holmquist et al., 1996; Orrenius et al., 1998).

Candida antarctica lipase B consists of 317 amino acid residues and has an / hydrolase fold (Uppenberg et al., 1994). It has been used successfully for resolving secondary alcohols and secondary amines (Öhrner et al., 1996). The combination of enzyme kinetics and molecular modeling was used to construct a model of how the enantiomers of secondary alcohols bind in the active site of this enzyme during transesterfication reactions (Orrenius et al., 1998). The two enantiomers of a secondary alcohol are shown to bind in two different orientations in the active site. The large group of the fast reacting R-enantiomer points out from the active site (mode I) (Fig. 1) while the large group of the slow reacting S-enantiomer points into the interior of the active site (mode II). When the large group is larger than ethyl, the selectivity increases dramatically due to the lack of space in the active site to accommodate it. The rate of conversion of the fast-reacting enantiomer is not significantly affected by an increase of the length of the large group. This is because it points out from the active site without unfavorable steric interactions with the enzyme. With these two modes of binding the enantiomers, all hydrogen bonds essential for the catalysis are present. These different modes of binding the enantiomers in C. antarctica lipase B result in a high enantioselectivity for secondary alcohols.

View larger version (47K): [in this window] [in a new window]   FIGURE 1   The four binding modes of the enantiomers. Productive bindings are shown at the top left and at the bottom right; nonproductive bindings at the top right and at the bottom left.

View larger version (20K): [in this window] [in a new window]   FIGURE 2   The fluctuation of the potential energy in a 1.5-ns-long MD simulation is shown with the transition state model comprising (R)-2-butanol in mode I orientation.

The enantioselectivity is caused by the difference in the free energy gaps between the ground states of the enantiomers and the corresponding transition states, which the enantiomers form with the enzyme. The thermodynamic ground states of the enantiomers are the same. Therefore, the enantioselectivity is determined by the difference in free energy between the two transition states. This difference in free energy (G^#) is related to the enantioselectivity, expressed as the enantiomeric ratio (E), as follows:

(1)

The tetrahedral intermediate is, at a first approximation, a good model of the transition state (Warshel et al., 1989). According to Burkert, the translational and rotational terms of the Gibbs free energy cancel to a very good approximation when the energy difference between two isomers is computed. The difference in the vibrational energy part is small compared to the difference in the potential energy part of two isomers (Burkert et al., 1982). The potential energy difference is calculated by the force field. Therefore, the force field energy (steric energy) well approximates the enthalpic part of the Gibbs free energy. The difference in steric energy (V^#) will be close to the free energy difference (G^#) in systems where the entropy contribution is of minor importance to the enantioselectivity.

Earlier attempts to calculate enantioselectivity fixed most of the atoms in the protein to their crystallographic coordinates, leaving only the substrate and various numbers of amino acid residues close to it free to move during the simulations (DeTar, 1981; Bemis et al., 1992; Norin et al., 1993, 1994). With systems where the lock and key concept (Fischer, 1894) works, this should be a successful way to calculate the difference in free energy. However, this approach is not appropriate when the enzyme must rearrange its geometry in order to accommodate the substrate in accordance with the induced fit concept (Haldane, 1930; Black et al., 1996). Releasing all geometrical constraints means more costly simulations in addition to larger energy fluctuations.

Different ways of forming subsets have been investigated in this study for the purpose of distinguishing the energy causing the enantioselectivity from the large energy fluctuations in the enzyme. Two fundamentally different approaches toward defining subsets for subsequent use in energy evaluations were tested. The first approach used structure-based subsets including different parts of the enzyme and the tetrahedral intermediate. In this approach the subsets could not change with changes of substrates. The second approach defined subsets on a basis of energy differences so that they could vary with the substrates.

	MATERIALS AND METHODS

Top Abstract Introduction Materials & Methods Results Discussion References

Preparation of the enzyme starting structure

The enzyme starting structure to be used in all computations was prepared as follows. The coordinates of C. antarctica lipase B were collected from the structure with entry number 1TCA in the Brookhaven Protein Data Bank (Bernstein et al., 1977). The two sugar units (NAG) were removed, since these are far from the active site and should have no influence on the enantioselectivity. The coordinates were read into the CHARMM program (Brooks et al., 1983) (version c24b2), which was used for all of the subsequent force field calculations. Hydrogen atoms were added using the ADDH command, except Asp-134, which was explicitly protonated. The cysteine residues 293 and 311, 22 and 64, and 216 and 258 were connected to cystine residues. The parameter set PAR  ALL22  PROT with all hydrogen atoms explicitly treated was used in all calculations. The water molecules were represented by the TIP3 model. In order to find all loosely bound water molecules we performed a 50-ps molecular dynamics (MD) simulation. The time step used was 1 fs and the temperature was 300 K. The system was heated to this temperature at the rate of 50 K/ps. The algorithm used was the Verlet algorithm. A distance-dependent dielectric function was used and a nonbonded cutoff of 8 Å was applied. The nonbonded interactions were updated every 25 steps. The SHAKE routine (Ryckaert et al., 1977) was applied to all bonds connected to hydrogens. This way of stripping the protein of water molecules resulted in the loss of three water molecules (HOH61, HOH126, and HOH242). The active site of the x-ray structure contained three water molecules. However, during the 50-ps MD simulation several water molecules entered the active site. In order to give space for the substrates, four water molecules in the active site (HOH130, HOH149, HOH238, and HOH285), in addition to the three that escaped the surface during MD, were removed and not present in the subsequent calculations. The remaining 279 water molecules were included in the following simulations.

The tetrahedral intermediate

As a model for the transition state we used the tetrahedral intermediate formed in the gas phase reaction between methyl acetate and a methoxide ion. The geometry was optimized by means of the restricted Hartree-Fock method using Gaussian 94. The basis set used was 6-31+G*. To account for the charged anion, diffuse functions were used in the basis set. Two local minima were found and optimized. The energy difference of these conformations was 2.9 kcal/mol (Fig. 3). Gaussian 94 (Frisch et al., 1995) was used in all quantum mechanical calculations.

View larger version (19K): [in this window] [in a new window]   FIGURE 3   The dimethyl orthoacetate anion geometry optimized using RHF/6-31+G*. Two conformations were found; conformation A was lower in energy than conformation B by 2.9 kcal/mol.

Force field parameters for bonding terms

Force field parameter values for bonding terms, needed in the calculations but missing in the CHARMM parameter file, were selected as follows. Equilibrium lengths and angles were selected from the ab initio optimized geometry of the dimethyl orthoacetate anion in conformation A (see Fig. 3). Force constants for stretching terms and bending terms were set at high values (see Table 2).

Missing torsional parameter values were selected by analogy with CHARMM-standard parameter values (X-C'-C"-Y) V_N = 0.20 kcal/mol and (X-C-O-Y) V_N = 0.14 kcal/mol, where X and Y refer to arbitrary atoms. The atoms labeled C, C', and C" refer to any carbon atom and the symbol O refers to any oxygen atom. The van der Waals parameters R* = 1.9 Å and * = 0.01 kcal/mol for the atom type CIN were selected by analogy with CHARMM-standard parameter values.

Force field parameters for nonbonded terms

In order to allow the assignment of accurate partial charges to the tetrahedral intermediates, these were geometry optimized using RHF/STO-3G. The electrostatic potential was then computed and partial charges derived (see Table 2), from these geometries with the basis set 6-31+G* using the option CHelpG (Breneman and Wiberg, 1990). To adjust the overall charge of the covalently linked Ser-105 and the tetrahedral intermediate to 1, the charge of C and its attached hydrogens were adjusted. However, their values were chosen to be as close as possible to the CHARMM partial charges of a Ser residue. All other charges needed in the calculations were chosen from the CHARMM PAR_ALL22_PROT parameter set.

Building the enzyme-substrate complexes

The x-ray structure of C. antarctica lipase B prepared as described above was docked with the tetrahedral intermediates of the octanoates of 2-butanol, 3-methyl-2-butanol, and 3,3-dimethyl-2-butanol. The acyl part of the intermediate (octanoyl) was modeled into the acyl pocket (Uppenberg et al., 1995) and the alcohol part was modeled in such a way that all essential hydrogen bonds were formed to His-224 and Thr-40. The side chain oxygen of Ser-105 was covalently linked to the central carbon in the tetrahedral intermediate. The R-and S-enantiomers were docked in both mode I and mode II in accordance with previous results (Orrenius et al., 1998) (Fig. 1).

Molecular dynamics simulations

At first, all of the water molecules were allowed to relax through 4000 steps of steepest-descent minimization, while all of the amino acid residues and the tetrahedral intermediate were kept fixed. All atoms were then allowed to move during 6000 steps of minimization using the steepest-descent method. After this, the structure was subjected to 2000 steps of minimization using the Powell algorithm. This was followed by a 200-ps molecular dynamics simulation using the Verlet algorithm. The system was heated by 10 K/ps to the equilibrium temperature 300 K. The nonbonded list was updated every 25 steps. A nonbonded cutoff of 8 Å was applied and a distance-dependent dielectric function was used. A time step of 1 fs was applied. The SHAKE routine (Ryckaert et al., 1977) was applied to all bonds connected to hydrogens. In all, 12 MD simulations were performed with a simulation time of 200 ps each. Samples of the trajectories were stored every 32 fs.

Low energy conformations

The 120 conformations collected from the MD simulations (chosen from all conformations stored, 75,000) in the time interval 70-200 ps with the lowest potential energy and separated by at least 0.5 ps (Karplus and Petsko, 1990) were subjected to subsequent refinement. The structures were first minimized with 8000 steps using the steepest-descent method and then with 6000 steps using the Powell method. The largest energy gradient in any of the structures was <0.08 kcal/mol Å.

Structural comparisons of the enzyme-substrate complexes with the x-ray structure

All minimized structures with the tetrahedral intermediates of 2-butanol and 3,3-dimethyl-2-butanol including R-and S-enantiomers and different binding modes were compared with the x-ray structure as follows. Three different sets of atoms were selected from the last collected and minimized structure from the MD simulation of (R)-3,3-dimethyl-2-butanol in mode I. All atoms of amino acids with at least one atom within a radius of 5, 8, or 10 Å of any atom of the tetrahedral intermediate were included in the sets. The x-ray structure was compared with the 40 structures of 2-butanol and the 40 structures of 3,3-dimethyl-2-butanol based on these three sets of atoms.

Structure-based subsets

Various subsets of atoms were used in the evaluation of the energies responsible for the enantioselectivity. The largest subset (subset 1) included all atoms. The second largest one (subset 2) contained all atoms but excluded the water molecules. The energy of the interaction between the water molecules and subset 2 was excluded.

The two subsets based on the radii 3 Å and 6 Å (subsets 3 and 4) were defined using the last collected and minimized structure from the MD simulation of (R)-2-butanol in mode I. All amino acid residues in this structure with at least one atom within a distance of 3 Å or 6 Å from the central carbon atom of the tetrahedral intermediate (the atom with the atom type CIN, see Fig. 4) were included in the subsets. The bonded and nonbonded interaction energies were calculated on these selected atoms and included the energies of nonbonded interactions with the surrounding atoms.

View larger version (15K): [in this window] [in a new window]   FIGURE 4   Schematic picture of the tetrahedral intermediate showing atom types used in the calculations.

The smallest subset (subset 5) included only the modified Ser-105 (tetrahedral intermediate) with its intermolecular interactions with the neighboring atoms.

Energy-based subsets

The selection of these subsets was based on the energy difference between the averaged interactions of enantiomers with a specific amino acid residue or water molecule in the enzyme (Eq. 2).

(2)

Mean values of the energies of interactions between any amino acid residue (_res,sub^xR) or water molecule and one enantiomer of the substrate were calculated over the 10 minimized structures. The difference in the mean energy value between the two enantiomer-residue pairs (_res^x,sub) was used to include or exclude that particular residue in the subset. The procedure was repeated with every substrate. Thus, the subset could be different for every substrate, including only those residues, which differed in the energies of interaction with the tetrahedral intermediates formed from the enantiomers. Energies were calculated on these energy-based subsets for different absolute values of _res^x,sub.

Nonbonded energies of interactions between the enantiomers, in different modes, of the three substrates and the 317 amino acid residues as well as the 279 water molecules were calculated. Four subsets were defined using different absolute values of _res^x,sub (>0.5, >0.3, >0.1, and >0.01 kcal/mol). The selection of amino acid residues and water molecules was based only on the R-enantiomer in mode I and the S-enantiomer in mode II. Intra and intermolecular energies were calculated on the selected atoms. Intermolecular energies with surrounding atoms were included.

	RESULTS

Top Abstract Introduction Materials & Methods Results Discussion References

Molecular dynamics

An initial MD simulation with a simulation time of 50 ps was made in order to find loosely bound water molecules on the protein surface. This resulted in three water molecules leaving the protein structure. These were not included in the subsequent calculations. Four water molecules occupying the active site were also excluded. The remaining 279 water molecules were included in the simulations. The tetrahedral intermediates formed by the enzyme and the octanoates of the secondary alcohols (2-butanol, 3-methyl-2-butanol, and 3,3-dimethyl-2-butanol) were manually docked into the active site as described elsewhere (Uppenberg et al., 1995).

A 1.5-ns-long MD simulation was performed with the transition state model comprising (R)-2-butanol in the mode I orientation. The potential energy fluctuation over time is presented in Fig. 2. No large conformational changes in the enzyme structure were observed in the MD simulation. The largest deviations from the x-ray structure were located near the C-terminus (Phe-304 to Pro-317), and in two other loop regions (Leu-140 to Val-149) and (Pro-260 to Ala-279). The overall root mean square (RMS) deviation (the backbone atoms were compared) between the structure sampled at 1.5 ns from the MD trajectory and the x-ray structure was 0.96 Å. Based on the results from this simulation, 12 MD simulations of 200 ps each were performed for the enzyme-tetrahedral intermediates of the three secondary alcohols with their enantiomers and two docking modes. The energy gradient remaining in any of the 120 energy minimized structures, described in Materials and Methods, was <0.08 kcal/mol Å.

The influence of the substrate on the lipase structure was studied by means of calculation of the RMS deviations between the x-ray structure and the refined structures of the enzyme-tetrahedral intermediates including the small substrate 2-butanol and the bulky 3,3-dimethyl-2-butanol. The values are presented in Table 1.

The tetrahedral intermediate

Two conformers of the dimethyl orthoacetate anion [CH₃CO(OCH₃)₂] (Fig. 3) were found (Ewig and Van Wazer, 1986) and the geometries were optimized using RHF/6-31+G*. Only conformation A with the lowest energy (lower by 2.9 kcal/mol than that of conformation B) could be fitted into the active site without unfavorable van der Waals interactions and with all of those hydrogen bonds formed, which are essential for the catalysis. The calculated equilibrium distances and angles of the central part of A were used as force field parameters in the subsequent molecular mechanics and molecular dynamics calculations (Table 2). To restrict the distances and angles to these values, large force constant values were applied (Table 2). After MD simulation and subsequent molecular mechanics refinement, the central part geometry of the tetrahedral intermediate in the active site of the enzyme remained very close to the set ab initio equilibrium geometry of the dimethyl orthoacetate ion.

Partial charges

All atoms of the tetrahedral intermediates of the three substrates studied were assigned partial charges from electrostatic potential calculations using ab initio methods (see Materials and Methods for details). The partial charges derived and subsequently used in the MM and MD calculations are presented in Table 2.

Catalytically important hydrogen bonds

The catalytic activity of serine hydrolases depends on a stabilization of the transition state caused by several hydrogen bonds. With all three substrates studied, three hydrogen bonds were formed in the oxyanion hole. The hydrogen bond between Ser-105 (O) and His-224 (N) was formed in all of the minimized structures. However, the hydrogen bond between the oxygen of the secondary alcohol and His-224 (N) was not formed in all of the minimized structures containing 3-methyl-2-butanol or 3,3-dimethyl-2-butanol (Table 3).

Structure-based and energy-based subsets

When all atoms of the refined structures were included in the energy evaluations, large negative energy values were obtained. When all the water molecules (279), which remained after stripping the protein of seven water molecules, were excluded from the calculations, the negative energy values increased by ~80% (7000 to 1000 kcal/mol). Table 4 presents the calculated force field energies of these two structure-based subsets and three smaller ones, which were predefined parts of the enzyme-substrate complexes that did not vary with the substrates. A number of subsets were constructed on the basis of pairwise interaction energies between a specific residue (or water molecule) and the tetrahedral intermediate. We will call these subsets energy-based subsets, as they varied with the substrate and the interaction energies. The energy differences (_res^x,sub) (Eq. 2) of these pairwise interactions between the R-configuration in mode I and the S-configuration in mode II were used to construct four subsets of amino acid residues and water molecules. The cutoff values |_res^x,sub| were chosen to be >0.01, >0.1, >0.3, and >0.5 kcal/mol (see Table 5). The calculated force field energies of these subsets are presented in Table 6.

The residues giving the largest contributions to the enantioselectivity (see Table 5) were Thr-40, Gly-41, Thr-42, Trp-104, Gln-106, Gly-108, Leu-109, Phe-131, Leu-140, His-224, Leu-278, and Ile-285. One water molecule (W1), hydrogen-bonded to the backbone of the catalytic Ser-105, the backbone of Phe-131, and the side chain of Thr-103, also made a significant contribution. These residues contributed differently to the enantioselectivities of the three substrates.

	DISCUSSION

Top Abstract Introduction Materials & Methods Results Discussion References

Calculating the enantioselectivity

The calculation of the enantioselectivity as the energy difference between the diastereomeric enzyme-substrate complexes at the transition state is complicated by the high energy values of the enzyme-substrate complexes and their fluctuations during the simulations. The difference in free energy between the transition states of the enzyme-substrate complexes is often not much larger than 1 kcal/mol. The CHARMM force field energy of the C. antarctica lipase B is ~1000 kcal/mol and when all crystal water molecules are included the energy lowers to ~7000 kcal/mol. For this reason, a selective way of determining the small energy differences responsible for the enantioselectivity must be found. Several attempts to do this with various enzymes and substrates have been reported (Norin et al., 1993, 1994). In most cases the fast-reacting enantiomer has been correctly predicted, though the calculated value of the potential energy difference at the transition state has not agreed quantitatively with the experimental value. The reasons for this can be several. The main assumption, using only the difference in force field energy to predict the free energy difference at the transition state, is that the entropy contribution can be negligible. The entropy contribution can be approximated at <1 kcal/mol (DeTar, 1981), provided that the enzyme can be considered rigid during the catalysis when bound to the two enantiomers, and that solvent molecules do not interact differently with the two enzyme-enantiomer complexes. However, it has been demonstrated that entropy can contribute significantly to the enantioselectivity (Jenny Ottosson and Karl Hult, personal communication). This calls for calculating the entropy contribution to the free energy as well. Åqvist and Warshel have developed a method to directly calculate the change in free energy when a ligand becomes bound to a protein. However, that method needs to be tested and calibrated against long free energy integrations (Åqvist and Warshel, 1993).

If carefully parametrized, the force field energy can well approximate the enthalpic part of the free energy. The structures of enzymes are distributed over a huge number of conformations close in energy and separated by low energy barriers (Karplus and Petsko, 1990). One conformation is obtained via the crystallographic structure. If the substrate is small and fits well in the active site according to the lock and key concept, this conformation should be maintained during the simulations. One way of keeping the enzyme from deviating from its crystallographic coordinates during simulations is to fix large parts of it, preferably the backbone atoms. However, if the substrate cannot easily be fitted into the active site without a rearrangement of the enzyme structure, the results will be unreliable. This was the case with the enantioselectivity of a bulky substrate, which was predicted wrong in an earlier study (Orrenius et al., 1998). In order for the enzyme to accommodate the substrate properly by adopting a different conformation, all atoms of the enzyme should be free to move during molecular dynamics simulations.

The tetrahedral intermediate approximation and the hydrogen bonding pattern

It is not possible to calculate the geometry and energy of a transition state using force field methods. This can be done only by using quantum mechanical methods. However, the size of enzymes still keeps the use of these methods out of reach. It has been proposed that the tetrahedral intermediate closely resembles the transition state and therefore could be used as a substitute for it (Warshel et al., 1989). We approximated the central part of the enzymatic tetrahedral intermediate with a small model system. The gas phase geometry of the dimethyl orthoacetate anion was used to construct the central part of the enzymatic tetrahedral intermediate. Using ab initio calculations, two conformations were found and optimized. Only one of these could be fitted into the active site without unfavorable van der Waals interactions and with all hydrogen bonds essential for catalysis formed (conformation A, Fig. 3). This conformation was the lowest in energy.

In order to study how the enzyme can accomodate and stabilize the gas phase structure of this simple model system instead of the tetrahedral intermediate, we applied large force constant values for stretching and bending energy terms to the central part of it. This was done in order to keep the geometry of that part fixed during the simulations. The geometry closely resembled the ab initio geometry throughout the molecular dynamics and molecular mechanics calculations.

To stabilize the transition state and thereby increase the rate of catalysis, a number of hydrogen bonds between the enzyme and the substrate must be formed. In C. antarctica lipase B the oxyanion hole is assumed to be built from three hydrogen bond donors (one from the backbone of Gln-106 and two from the backbone and side chain of Thr-40), the construction being similar to that of the oxyanion hole in cutinase (Nicolas et al., 1996). The proton transfer between His-224, Ser-105, and the substrate is supported by two hydrogen bonds.

Hydrogen bonds were formed between the oxyanion hole and the oxyanion in all refined structures collected from the MD trajectory. This indicates that these hydrogen bonds contribute considerably to the stabilization of the transition state. The hydrogen bond between Ser-105 (O) and His-224 (N) was formed in all structures. This demonstrated that the gas phase structure of the model system was representative of the enzymatic transition state. The hydrogen bond between the oxygen atom of the alcohol and His-224 (N) was formed by all structures of 2-butanol with R-and S-enantiomers in mode I and mode II. With the larger and more bulky substrates 3-methyl-2-butanol and 3,3-dimethyl- 2-butanol, the hydrogen bond of all structures except one (3,3-dimethyl-2-butanol) was formed with the S-enantiomer in mode II (Table 3). However, with the fast-reacting enantiomer in the productive mode I orientation, this hydrogen bond was formed only in seven of the 10 refined structures of 3-methyl-2-butanol. With 3,3-dimethyl-2-butanol the hydrogen bond was formed in six of the 10 structures. The difficulty of forming this essential hydrogen bond in the catalysis of the fast reacting enantiomer might be a reason for the lower rate of conversion of bulky substrates (Orrenius et al., 1998).

Sampling the conformational space

It is not possible to search the total conformational space of an enzyme/substrate system of this size, but one is restricted to making a local search. The feasible methods to use are either of the Monte Carlo (MC) or of the MD type. We chose MD for searching local minima. After a heating period of 30 ps, the temperature was kept constant at 300 K during the simulations. Cooling and heating the system between a low and a high temperature is a strategy called simulated annealing, in which low energy conformations are supposed to be trapped during the cooling period. The enzyme studied had no geometrical constraints. Hence, heating the system above 300 K did not seem feasible. In order fully to use the relatively short MD simulation to search local minima via barrier crossing, the temperature was kept constant at 300 K. Low energy conformations were collected from the MD trajectory with a time interval of at least 0.5 ps to avoid structures belonging to the same local minimum. The structures were collected within the time interval 70-200 ps, when the system had equilibrated, and exhibited its lowest potential energy.

To evaluate the influence of the size of the substrate on the lipase structure, we calculated the RMS deviations between the x-ray structure and the enzyme-tetrahedral intermediates. The steric influences exerted by the tetrahedral intermediates on the surrounding enzyme were investigated by varying the volume of the enzyme included in the comparisons. The fast enantiomers of 2-butanol and 3,3-dimethyl- 2-butanol, in mode I orientation, caused smaller deviations of the enzyme structure than the slow reacting enantiomers in mode II orientation, when compared with the x-ray structure. In the low energy reactive mode, the bulky 3,3-dimethyl-2-butanol caused a larger RMS deviation than 2-butanol when compared with the x-ray structure. The R-enantiomer and the S-enantiomer of 2-butanol showed larger RMS deviations in the sterically more hindered mode II than in mode I. The S-enantiomer of 3,3-dimethyl-2-butanol showed the same behavior as 2-butanol with a higher RMS deviation in mode II than mode I. The R-enantiomer of 3,3-dimethyl-2-butanol, on the other hand, showed smaller RMS deviations in mode II than in mode I. However, none of the 10 refined structures of mode II had a hydrogen bond between the alcohol oxygen atom and His-224 (N). Although the R-enantiomer of 3,3-dimethyl-2-butanol is more sterically hindered in mode I than in mode II, mode I must be the productive and mode II the nonproductive one, because the essential hydrogen bond is lacking in mode II (Table 3). This demonstrates the importance of avoiding geometrical constraints during MD and MM. Moderate movements in the tertiary structure occurred during the simulations. The loop near the C-terminus moved in relation to that of the x-ray structure.

In an earlier study, the S-enantiomer of 3,3-dimethyl-2-butanol was predicted to be the fast-reacting one (Orrenius et al., 1998). This erroneous prediction was probably due to the geometrical constraints applied in those simulations. The necessity of binding the enantiomers in different orientations in order to go through catalysis, and the unnatural constraints applied on some parts of the enzyme, could have resulted in strong repulsive interactions between the R-enantiomer and the constrained surrounding enzyme.

Structure-based and energy-based subsets

The most straightforward way of evaluating the energies involved in the enantioselectivity would be to include the entire enzyme, the tetrahedral intermediate, and all of the water molecules (subset 1) in the energy calculations. It was demonstrated that this led to predictions not compatible with the experimental results (Table 4). The greatest uncertainty in this procedure was probably the water molecules with very high mobilities causing considerable energy fluctuations. Another problem with these large systems, when no geometrical constraints were applied, is that very extensive refinement must be carried out to reduce the energy gradient and find the correct energy of a specific local minimum. Thus, the inclusion of the energy from the whole system led to a greater error than the inclusion of only a limited part of the system. In a second approach all of the water molecules were excluded in the energy evaluations (subset 2). This method could not reproduce the experimental results. The reason was probably the strong hydrogen bonding between the water molecules and the enzyme. By just excluding all energy contributions caused by these hydrogen bonds, considerable errors were introduced in the calculations. Two smaller subsets were then tested (subsets 3 and 4), which contained the same residues at the studies of all of the substrates. The larger one of these subsets predicted the wrong enantiomer of 2-butanol. When only the tetrahedral intermediate with its bonded and nonbonded interactions with the surrounding enzyme was included in the energy evaluations (subset 5), the magnitude of the enantioselectivity went down, but the wrong enantiomer was predicted for 2-butanol and 3,3-dimethyl-2-butanol. The calculation of the enantioselectivity using structure-based subsets (subsets 1-5) was clearly demonstrated not to be appropriate. Only subset 4 could qualitatively predict the enantioselectivities toward all three substrates correctly. However, this structure-based subset predicted C. antarctica lipase B to be more enantioselective toward 2-butanol than toward 3-methyl-2-butanol.

In order to increase the accuracy of the energy evaluations but also to find the residues, which were representative of the enantioselectivity, we examined the influence of every amino acid residue and water molecule on the enantioselectivity (forming energy-based subsets) using Eq. 2. If the difference in mean value |_res^x,sub|, between a specific residue or water molecule, and the R- and the S-enantiomers was larger than 0.01, 0.1, 0.3, or 0.5 kcal/mol, that residue or water molecule was included in the subset. Using this concept of forming energy-based subsets it was possible to find all major contributors to the enantioselectivity. Several residues caused the enantioselectivity, and they did not surround the tetrahedral intermediate uniformly. Therefore, the subsets should not be selected on a structural basis by using, for example, a sphere of residues to be included around the tetrahedral intermediate (Fig. 5).

View larger version (76K): [in this window] [in a new window]   FIGURE 5   Stereoscopic view of Candida antarctica lipase B (gray ribbon representation) with the transition state comprising 3,3-dimethyl-2-butanol in mode I orientation (black). The energy-based subset with a cutoff = 0.3 kcal/mol is shown at the top (dark gray), and subset 3 is shown at the bottom (dark gray).

The fast enantiomer of all the substrates was correctly predicted with the strategy of forming energy-based subsets. The R-enantiomer of all the substrates was the lowest in energy when modeled in mode I orientation. The S-enantiomer of 2-butanol was modeled in mode I orientation, but during molecular dynamics it rotated into an orientation close to mode II. For 3-methyl-2-butanol the S-enantiomer was lowest in energy in mode I orientation. However, only two of the 10 refined structures showed an intact hydrogen bond pattern. The S-enantiomer of 3,3-dimethyl-2-butanol was lowest in energy in mode II orientation.

The above results show that by using too large subsets, energy-based as well as structure-based ones, large fluctuations result that overrule the small energy difference causing the enantioselectivity. It is also clear that if too little of the surrounding enzyme is included, by constructing small subsets, important interactions with the tetrahedral intermediates are not considered. This too leads to a poor prediction of the enantioselectivity. Comparisons with experimental data (Table 7) show that the calculated values of the enantioselectivity of the various substrates are systematically larger than the experimental ones. However, this approach does predict the fast-reacting enantiomer correctly for all three substrates studied.

View this table: [in this window] [in a new window]   TABLE 7   Differences in potential energies (V#) between the transition state complexes of the R- and S-enantiomers of the substrates and the corresponding experimentally determined values of (H#)

We suggest that the calculation of enantioselectivity should be performed in three steps: generation of structures and energies, forming energy-based subsets, and energy evaluations based on these energy-based subsets. The energy values obtained from the calculations can then be used qualitatively to predict the fast-reacting enantiomer. It has been clearly demonstrated in this work that by using this approach and applying no constraints on the enzyme structure, the enantioselectivity of small as well as bulky substrates can be predicted correctly.