Simulating the catalytic effect of a designed mononuclear zinc metalloenzyme that catalyzes the hydrolysis of phosphate triesters.

One of the greatest challenges in biotechnology and in biochemistry is the ability to design efficient enzymes. In fact, such an ability would be one of the most convincing manifestations of a full understanding of the origin of enzyme catalysis. Despite some progress on this front, most of the advances have been made by placing the reacting fragments in the proper places rather than by optimizing the preorganization of the environment, which is the key factor in enzyme catalysis. A rational improvement of the preorganization and a consistent assessment of the effectiveness of different design options require approaches capable of evaluating reliably the actual catalytic effect. In this work we examine the ability of the empirical valence bond (EVB) to reproduce the results of directed evolution improvements of the catalysis of diethyl 7-hydroxycoumarinyl by a designed mononuclear zinc metalloenzyme. Encouragingly, our study reproduced the catalytic effect obtained by directed evolution and offers a good start for further studies of this system.

Introduction

Rational enzyme design has a wide scope ranging from general industrial applications to medicine.[1] In fact, designing of an enzyme with a novel function can be considered as the best manifestation of the understanding of enzyme catalysis and enzyme evolution. However, the present generation of designers enzymes are much less efficient than naturally evolved enzymes.[1,2] The problems with current efforts of rational design is most likely due to an incomplete modeling of the transition state (TS) in the enzyme active site, and in particular to the limited awareness to the key role of the reorganization energy.[3] Thus, an effective enzyme design methodology should be judged by its ability to determine the activation free energy, along with firm understanding of the factors governing the change in the TS energy in directed evolution experiments.

The challenges of modeling enzymatic transition states is far from trivial as it requires both, extensive sampling and reliable potential surfaces. Here perhaps the most effective option is the use of the empirical valence bond (EVB). The EVB is a semiempirical quantum mechanics/molecular mechanics (QM/MM) approach,[3b,4] where the QM part is represented by empirical approximations of the relevant valence bond integrals.[4] The EVB has been successfully used in reproducing and predicting mutational effects,[5] as well as in quantitative screening of design proposals and in reproducing observed effect of directed evolution refinement of Kemp eliminases.[6] In addition to the EVB, one can use molecular orbital–QM/MM (QM(MO)/MM)[7] methods. This type of approach is in principal effective, but at present it involves major difficulties in obtaining reliable free energies by sampling the surfaces obtained with high level ab initio methods. Some effective options like paradynamics method[8] can help in this respect.

In considering the EVB as an effective tool for computer-aided enzyme design, it is useful to note that this approach has reproduced reliably the observed activation barriers for different mutants of trypsin,[5a] dihydrofolate reductase[5b] and kemp eliminase.[6] Nevertheless, it is important to further validate the EVB approach with newer sets of designed enzyme and different types of active sites. In this work we will focus on a designed mononuclear zinc metalloenzyme, which catalyzes the hydrolysis of a model organophosphate.[9] The design of this metalloenzyme started from adenosine deaminase with was manipulated by a denovo methodology[10] with the aim of generating an enzyme that can catalyze the hydrolysis of an organophosphate.[9] As in other previous cases, the most effective steps in the refinement were achieved by directed evolution experiments that mimic natural evolution by selecting mutations which are beneficial to the overall catalytic activity of an enzyme.[11] Thus, studies of this designed enzyme give us both an opportunity to validate our approach on metalloenzymes, and provide (at least in principle) the opportunity to study an evolutionary trajectory where enzyme evolves to perform a completely new function.

Systems and Methods

Systems

As stated above, the enzyme chosen for this study is a designed mononuclear zinc metalloenzyme, which catalyzes hydrolysis of diethyl 7-hydroxycoumarinyl phosphate (DECP) (Figure 1a) (mimicking organophosphate nerve agents).[9] This enzyme was designed from adenosine deaminase which is a mononuclear zinc metalloenzyme, where metal ion is thought to be primarily acting as an activating agent for a hydroxyl ion nucleophile.[12] Directed evolution process leads to different mutants with different catalytic power. The first variant that was found to show detectable activity (k/Km) contains eight mutations (designated as PT3). Three other variants, PT3.1, PT3.2, and PT3.3, in the evolutionary trajectory were found to have activities of (k/Km, M–1 s–1) of 4, 154, 959, and 9750, respectively, and kcat (×10–3 s–1) of 5 × 10–5, 0.2, 4, 47, and 351, respectively.

Figure 1

(a). Chemical structure of diethyl 7-hydroxycoumarinyl phosphate (DECP). (b). Evolutionary trajectory of the DECP hydrolysis activity.

In order to verify our ability to reproduce the results of the directed evolution experiments, we have simulated the activation barriers for the hydrolysis of DECP by adenosine deaminase and its four variants (PT3, PT3.1, PT3.2 and PT3.3) (Figure 1b). The calculations used as starting points of the crystal structures of adenosine deaminase (PDB id: 1A4L) and PT3.1 (PDB id: 3T1G). The structures of PT3, PT3.2 and PT3.3 were created by generating the corresponding mutations in the crystal structure of PT3.1 (PDB id: 3T1G).

Dealing with a designed enzyme, which performs a new function, it is imperative to examine the binding mode and reaction mechanism of the given substrate. The structure of the substrate–enzyme complex was created with the help of the model of the enzyme–TS complex used in a previous study[9] and with the help of docking using Autodock 4.[13] Of course, the configurations generated by Autodock were subsequently subjected to extensive relaxation in the EVB calculations. Figure 2 depicts both the enzyme (with the bound substrate) and a representative schematic of the binding modes of the ligands around the Zn metal. As seen from the figure the activated hydroxyl and phosphoryl oxygen both coordinate with the Zn metal and occupy two out of six coordinating positions in octahedral structure around the metal center. All five systems were subjected to relaxation through 100 ps of molecular dynamics simulation prior to the EVB calculations.

Figure 2

(a) PT3 enzyme and the bound substrate. (b) Substrate binding mode around the Zn metal ion.

Ab Initio Potential Energy Surface

The first step in establishing the mechanism of the phosphotriester hydrolysis requires information about the free energy surface of the reference solution reaction. Such information is best obtained by ab initio calculations,[14] but performing such calculations in a fully consistent way can be very challenging (see for example[15]). Thus, considering the fact that our main interest is in mutational effect rather than the whole catalytic effect, we only performed qualitative ab initio mapping of the relevant surface. This was done using the Gaussian03 software package,[16] and modeling for simplicity a phosphotriester, which has only methoxy groups attached to the phosphorus center. All structural optimizations and energy evaluations of the ab initio potential energy surface were performed using the 6-31++G** basis set with the B3LYP hybrid density functional. The solvent was treated with the COSMO implicit solvent model[17] and the resulting free energy surface is given in Figure 3. Obviously, more systematic calculations are needed, including careful QM/MM studies such as those preformed in our recent studies[15] and a systematic examination of the effect of the Zn (that would require to include the ligands in the ab initio calculations) . However, our main point here is that the present surface and other studies,[18] as well as studies of related systems,[19] indicate that we have a high-energy reaction intermediate or a plateau that leads to a mechanism of the type described in Figure 4. This mechanism justifies the use of three states EVB surface, which will be discussed below. At any rate, since our effort is dedicated to exploring mutational effects, we left further exploration of the solution surface and its modification by the Zn ion to subsequent studies.

Figure 3

Ab initio energy surface for the hydrolysis of dimethyl coumarinphosphate in solution (Onuc and Olg designate, respectively, the nucleophilic oxygen and the leaving group oxygen).

Figure 4

Schematic description of the energetics of stepwise hydrolysis of DECP.

EVB Calculations

As in previous studies,[6] we performed our simulations using the EVB method, which is described in great details elsewhere.[3b,4] In setting out the EVB potential, we note that the stepwise hydrolysis of DECP can be studied by using a three state EVB description of the system[20] (see Figure 5). This behavior can be described by an analytical expressionwhere Δg‡ is the actual activation barrier, Max (X;Y) is the maximum of the two variables X and Y and Δg2‡ = Δg23‡ + ΔG12. Although we evaluate the activation barriers, Δg‡, by the full EVB calculations we note that they can also be estimated by our linear free energy relationship (LFER) expression[4]where λ, ΔG, and H are reorganization energy, reaction energy, and the off diagonal mixing term, respectively. The effect of the specific environment is incorporated by taking into account the changes in the corresponding by the reorganization energies and/or by changing the value of ΔG.

Figure 5

Schematic description of the three state EVB model used to describe the hydrolysis of DECP.

In calibrating the EVB potential, we did not try to use the observed energetics of the reaction of hydroxide attack in water[21] since it does not include the effect of the Zn ion. Furthermore, here we are interested in the mutational effects, rather than the catalytic effect relative to water. Thus, we calibrated the surface taking the reaction in 1A4L as the reference reaction, where we fixed the rate determining barrier around 27.5 kcal/mol, while assuming that in 1A4L this barrier is Δg‡2 (this assumption is based on the group contribution reported below). We also selected a value of 20 kcal/mol for ΔG12 in 1A4L (see below) and took for ΔG23 in 1A4L a value of −6 kcal/mol. We note in this respect that the results do not depend strongly on the values of ΔG12 and ΔG23. That is, as can be seen by using eq 2 the trend in the rate-determining barrier for different mutants does not depend the corresponding ΔG23 (with a proper adjustment of H23) and this is also partially true with regards to ΔG12, since we can adjust H12 to obtain the same trend. The real uncertainty is in determining whether the second barrier is rate determining and at what point the first barrier starts to be rate limiting (the change in the LFER). Resolving this issue requires LFER experiments or very careful PD calculations. Thus, the decision on the point of change in the LFER is somewhat arbitrary at the present case. At any rate, our EVB parameters are given in the Supporting Information.

The EVB calculations were performed with the MOLARIS program[22] in conjunction with ENZYMIX force field.[23] The EVB activation barriers were estimated at configurations selected by the same free energy perturbation umbrella sampling (FEP/US) approach described extensively elsewhere.[3b,4] The simulation systems were solvated by the surface constrained all atom solvent (SCAAS) model,[23] with a water sphere of 18 Å radius around the substrate and surrounded by 2 Å grid of Langevin dipoles followed by a bulk solvent. The long-range electrostatic effects were treated by the local reaction field (LRF) method.[23] The EVB region consisted of the substrate molecule and the hydroxide group. The FEP mapping was evaluated by 21 frames of 20 ps each for moving along the reaction coordinate using SCAAS model. All the simulations were performed at 300 K with a time step of 1 fs for integration. In order to obtain converged results, the calculations were repeated 5 times with different initial conditions.

Estimating Group Contributions

The contributions from each residue to the activation barrier (the group contributions) were estimated by calculating the effect of change of substrate charges (from RS to TS) on the electrostatic contribution of each protein residue. As discussed in our previous studies (e.g., ref (6)), the electrostatic contributions of all the protein residues to the activation barrier can be estimated by the following expression:[3a,24]Here the 332 factor is the conversion to kcal/mol, q are the residual charges of the protein atoms in atomic units (j runs over the protein residues and k runs over the atoms of the jth residues and i over the substrate atoms), r is the distance in A between the kth atom of the jth group and the ith atom of the substrate, ε is the effective dielectric constant for the specific interaction, and ΔQ are the changes in the substrate charges upon going from the RS to TS. Decomposing this expression to the individual group contributions[3a,24] allows one to explore the approximated effect of mutating ionized or polar residues.

Results and Discussion

Accurate estimation of the catalytic effects of the different enzyme construct/mutants can be considered as the most basic requirement for the effective enzyme design or understanding to evolutionary mechanism. Therefore, we started with systematic evaluations of the activation barriers for our systems. Our typical procedure of obtaining activation barrier involved average over 5 free energy profiles, for each enzyme variant (mutant). The details of the calculations are summarized in Table S1 (Supporting Information) and the estimated barriers are summarized in Table 1 and Figure 6).

Table 1

Calculated and Observed Activation Free Energies for the Systems Studied in This Work

systems	Δgobs‡, kcal/mol	Δgcalc‡, kcal/mol
1A4L	27.48	26.42
PT3	22.55	20.97
PT3.1	20.77	20.64
PT3.2	19.31	19.92
PT3.3	18.11	18.59

Figure 6

Correlation between the calculated and observed activation free energies. for the hydrolysis of DECP in the enzymes studied.

The correlation between the calculated and observed activation barriers (Table 1 and Figure 6) suggests that change in activity is driven by the change in transition state binding and not by some other elusive factors (such as substrate binding or dynamics). The successful demonstration of our ability to estimate accurate activation energies also indicates that the binding mode of substrate and the reaction mechanism used are reasonable. It should be noted that this is a designed enzyme, and therefore, no concrete prior information about the binding mode or reaction mechanism is available.

We believe that rational enzyme designing procedure can be improved if we can quantify the contribution of each residue to the transition state binding. Considering the fact that the electrostatic interaction is by far the most important factor in transition state stabilization and therefore enzyme catalysis, we have calculated the electrostatic group contributions of the protein residues. This was done, as discussed in section II.4, by using eq 3 and collecting the contribution of each residue to the overall sum (namely the electrostatic contribution for the energy of moving from the reactant to transition state). Specifically, we have (artificially) changed the charge of protein residues of 1A4L (the “wild type”) from 0 to −1, and then calculated the change in corresponding group contribution upon change of the residual charges of the reacting substrate. As can be seen from Figure 7b, the contributions of residues19 and 296 to the rate limiting C–Olg bond dissociation step,Δg2‡, are positive (note as is clear from the Supporting Information that Figure 7a is for a barrier that does not correspond to the rate limiting step). Thus, changing the charges of the corresponding residues from −1 to 0 should lead to a reduction in Δg2‡. This is consistent with the finding[9] that removing the charges of D19 and D296 (the D19S and D296A mutations) in 1A4L is necessary for effective hydrolysis of DECP. We focus here on these two mutations since they are well-defined experimentally observed electrostatic mutations. In principle we can use the group contributions for further predictions but this is not the purpose of the present work, since these contributions are much less reliable than those obtained from EVB calculations when they involve residues near the substrate.[3a,6a] The group contributions should be, however, very useful for the small contributions of distanced ionized residues, and exploring this point is left to subsequent studies.

Figure 7

Group contributions (in kcal/mol) for (a) the nucleophilic attack and (b) the bond dissociation steps in 1A4L. The group contributions reflect the interactions between the changes in the charge of protein residues from 0 to −1, with the charge change of substrate upon moving from RS to TS1 and TS2. The relatively large positive contributions provide a rough guide for the optimal sites for effective mutations that would enhance the catalytic effect. Since the second step is rate limiting in 1A4L, the corresponding group contributions are those that should be compared to the observed results.

Concluding Remarks

The ability to accurately estimate the activation energy of different variant enzyme of an enzyme can significantly improve the effectiveness of enzyme design efforts. At present, most enzyme design methods rely on directed evolution experiments to refine and increase the activity of the designed enzyme. In principle, in silico procedures can help in increasing the activity of designers enzymes by accurately estimating the effect of proposed mutations on the rate determining activation energies. Gas phase calculations or calculations which explicitly focus on the electrostatic interaction between the protein residues and the TS are very unlikely to have success in estimating the activation barriers as they do not consider the surrounding environment and its reorganization during the reaction. In principle, QM(MO)/MM[25] treatments can account for the enzyme environment. However, the difficulties of obtaining converging free energy calculations make it hard to use such methods in accurately estimating mutational effects.

On the other hand, the EVB has been shown to be capable of estimating the effect of mutational change on activation as early as 1986,[5a] where computer-aided mutations were proposed for rat trypsin. As far as enzyme design is concerned, we like to point out that EVB has been shown to be capable of reproducing the effect of mutations observed in directed evolution of kemp eliminases.[6] However, more studies are clearly needed and thus we have extended here the validation of the EVB to a study of the effects of multiple mutations on the activity of a designed Zn metalloenzyme. In doing so we note that the relatively high reactivity of metalloenzyme, coupled with the wide range of reactions carried out by them, makes them very attractive starting points for introducing new activities. At any rate, in the present study, we have successfully estimated the activities of different variants of the designed metalloenzyme and have reproduced the evolutionary trajectory leading to a new catalytic function (hydrolysis of DECP).

While determining the effect of different mutations on activation energies is the key to effective rational design, it would be useful to have a qualitative guide to propose mutations which can decrease the activation energy and therefore can increase the catalytic activity. Here we provide indications that the electrostatic group contributions can provide an important lead for mutations, which can improve the activity of an enzyme. In particular the group contributions in 1A4L reproduced the experimental trend that mutations that remove the negative charges at position Asp19 and Asp296 increase the activity.

Directed evolution has emerged as a powerful approach that can offers an effective way of optimizing enzyme activity. However, at present such approach has not achieved the same impressive catalytic power on enzymes that evolved by natural evolution. Overcoming this limitation will require exploration of mutational trajectories beyond what has been suggested by directed evolution. The EVB can be very useful in advancing such studies.

Despite the encouraging results of the present study it is important to mention that we did not performed a sufficiently careful study of the reference solution reaction or the effect of the Zn ion and its ligands and used relatively tentative estimates in estimating the reference surface in 1A4L. To further advance in this direction it would be crucial to preform ab initio QM/MM (QM(ai)/MM) free energy calculations for the solution reaction with and without the Zn ion as well as PD calculations of the reaction in the enzyme. It will also be very useful to have LFER experiments for the reaction in the enzyme. This will help in reducing the uncertainties about the rate determining barriers and increases the ability to make quantitative predictions of mutational effects.