Modeling the Complete Catalytic Cycle of Aspartoacylase - The

Article pubs.acs.org/JPCB

Modeling the Complete Catalytic Cycle of Aspartoacylase Ekaterina D. Kots,†,‡ Maria G. Khrenova,† Sofya V. Lushchekina,‡ Sergei D. Varfolomeev,†,‡ Bella L. Grigorenko,†,‡ and Alexander V. Nemukhin*,†,‡ †

Chemistry Department, M.V. Lomonosov Moscow State University, Leninskie Gory 1/3, Moscow 119991, Russian Federation N.M. Emanuel Institute of Biochemical Physics, Russian Academy of Sciences, Kosygina 4, Moscow 119334, Russian Federation

‡

S Supporting Information *

ABSTRACT: The complete catalytic cycle of aspartoacylase (ASPA), a zinc-dependent enzyme responsible for cleavage of N-acetyl-L-aspartate, is characterized by the methods of molecular modeling. The reaction energy profile connecting the enzyme−substrate (ES) and the enzyme−product (EP) complexes is constructed by the quantum mechanics/ molecular mechanics (QM/MM) method assisted by the molecular dynamics (MD) simulations with the QM/MM potentials. Starting from the crystal structure of ASPA complexed with the intermediate analogue, the minimumenergy geometry configurations and the corresponding transition states are located. The stages of substrate binding to the enzyme active site and release of the products are modeled by MD calculations with the replica-exchange umbrella sampling technique. It is shown that the first reaction steps, nucleophilic attack of a zinc-bound nucleophilic water molecule at the carbonyl carbon and the amide bond cleavage, are consistent with the glutamate-assisted mechanism hypothesized for the zinc-dependent hydrolases. The stages of formation of the products, acetate and L-aspartate, and regeneration of the enzyme are characterized for the first time. The constructed free energy diagram from the reactants to the products suggests that the enzyme regeneration, but not the nucleophilic attack of the catalytic water molecule, corresponds to the rate-determining stage of the full catalytic cycle of ASPA.

■

INTRODUCTION Aspartoacylase (ASPA) is a zinc-dependent enzyme that can effectively metabolize N-acetyl-L-aspartate (NAA), the major amino acid derivative in the mammalian brain.1,2 This enzyme performs cleavage of the amide bond in NAA producing the acetate and L-aspartate species (Scheme 1).

reaction mechanism including product formation and release are important to suggest the practical ways controlling the enzyme activity. To the best of our knowledge, a limited set of experimental data directly related to the reaction mechanism have been reported. The first crystal structure of human ASPA in the apo form was considered by Bitto et al. (PDB id: 2I3C).8 On the basis of structural data, the authors assumed that the catalytic mechanism of ASPA should be analogous to that of carboxypeptidases. According to this hypothesis, the zinc ion is coordinated to the His21, Glu24, and His116 side chains, and the Glu178 residue acts as a general base that activates the nucleophilic water molecule. Le Coq et al.9 determined the structure of human ASPA in the complex with a stable intermediate analogue, N-phosphonomethyl-L-aspartate (PDB id: 2O4H). Recently, this research group also characterized the crystal structures of the mutated variants (K213E, Y231C, F295S, E285A) of human ASPA with the same intermediate analogue.10 Two papers reported the observed rate constants kcat of the reaction: 14.22 ± 0.48 s−1 11 and 12.7 ± 0.05 s−1.9 The constants of the Michaelis−Menten kinetics were measured in

Scheme 1. Chemical Reaction of the N-Acetyl-L-aspartate Hydrolysis Catalyzed by ASPA

It is shown that NAA is responsible for the formation and maintenance of myelin being a source of the acetyl groups for the brain fatty acids.3,4 Reduction of NAA concentration is related to several disorders including multiple sclerosis and Alzheimer’s disease. An abnormally high level of NAA is associated with the Canavan disease caused by ASPA deficiency which results in disruption of myelin and spongy degeneration of the white matter of brain.5−7 Details of the ASPA catalyzed hydrolysis © XXXX American Chemical Society

Received: March 11, 2016 Revised: April 13, 2016

A

DOI: 10.1021/acs.jpcb.6b02542 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

to the hydrogen bonding network were assigned by using the Reduce package.19 Water molecules resolved in the crystal structure as well as a shell of almost 600 water molecules solvating the protein were included to the model system. Calculations of the energy values and forces along the reaction profile were carried out using the QM/MM approach. In the majority of these calculations, a QM part included the substrate, the side chains of His21, Glu24, Arg63, His116, Glu178, Phe282, and Tyr288, the peptide bonds from Asn117 and Thr118, the Zn ion, and five water molecules (132 atoms in total). Scheme 2 shows the molecular groups assigned to the QM region.

ref 11 following the release of L-aspartate. The spectroscopy methods were used in this study to quantify the amino group of the released species. In ref 9, the kinetic parameters were determined by a coupled assay wherein the hydrolyzed Laspartate product was converted to fumarate, the production of which was assessed by UV−vis spectrophotometry. From the theoretical side, Zhang et al. reported the results of computational simulations of the initial stages of this reaction.12,13 The authors performed the density functional theory (DFT) calculations for a small cluster composed of a few molecular groups at the active site and computed the potential energy barrier of 34 kcal/mol for the nucleophilic attack of the catalytic water molecule.12 Later,13 the quantum mechanics/ molecular mechanics (QM/MM) based calculations were performed to map the minimum-energy pathway using the semiempirical version of DFT, the self-consistent charge-density functional tight-binding (SCC-DFTB) model, to estimate potentials in the QM region. The reaction pathway was modeled from the ES complex to the first intermediate and then to a system composed of the protonated acetic acid and the doubly charged (2−) L-aspartate with a neutral amine group, but not the true product species as indicated in the right-hand side of Scheme 1. The minimum-energy profile for the pathway computed at the QM(SCC-DFTB)/MM(CHARMM) level again showed a very high potential energy barrier of 24 kcal/mol. Application of the same QM/MM energies to calculate the potential of mean force using the molecular dynamics (MD) umbrella sampling technique resulted in the energy relief with the activation barrier of 16 kcal/mol at the considered reaction segment. The latter energy value is consistent with the measured rate constants for the entire enzyme reaction.9,11 The motivation of this work was to provide the description of all elementary steps of this multistage reaction and to identify the rate-determining states of the entire catalytic cycle of ASPA from the reactants E + S to the products E + P. Suggestion of ref 13 to relate the first reaction step of the nucleophilic attack of water to the experimental rate constant does not seem fully justified. In particular, the results of recent modeling the peptide hydrolysis by the glutamate-assisted mechanism in matrix metalloproteinase 2, another zinc-dependent enzyme, showed14 that the release of the products corresponded to the rate-limiting step. The QM/MM approach15−18 with a large QM subsystem described at the DFT level was applied to compute the energy profile of the chemical transformations between the ES and EP complexes and to locate all reaction intermediates and transition states. By the EP complex, we assumed the product species, acetate (P1) and L-aspartate (P2), trapped at the active site inside the protein. The MD simulations with the QM/MM potentials were carried out for the critical parts of the reaction pathway. The stages of substrate binding to the enzyme active site (E + S → ES) and release of the products (EP → E + P1 + P2) were modeled by the MD calculations with the replica-exchange umbrella sampling technique.

Scheme 2. Molecular Groups Assigned to the QM Regiona

a

The immediate participants of the chemical transformations are shown in red, and the coordination shell of zinc is colored blue.

The critical points on the energy surface were also computed at the QM/MM level with an extended QM subsystem which included in addition to the molecular groups shown in Scheme 2 several charged molecular groups near the active site, Arg71, Asp114, and Arg168, giving rise to 162 atoms in QM. The results obtained with the extended QM model are described in the Supporting Information. Importantly, the basic conclusions of this work were not affected by this extension of the QM region. The structure of the starting point on the computed reaction pathway referred to the first tetrahedral intermediate I1. Then, the minimum-energy reaction pathway was modeled in two directions: to the ES complex and to the EP complex. Location of the minimum-energy stationary points corresponding to the ES and EP configurations and to the reaction intermediates, I1, I2, I3, I4, and I5, was carried out in series of unconstrained QM/MM optimizations and confirmed by computing harmonic vibrational frequencies at the optimized structures. To locate transition state configurations, the Hessian matrices were computed for putative structures followed by optimization along the modes with imaginary frequencies and confirmed by recomputing the frequencies at the located geometry configurations. When the saddle points were located, we verified that the descent in forward and backward directions correctly led to the minimumenergy structures of the corresponding reaction intermediates. Therefore, we believe that the computed minimum-energy profile connects smoothly all the stationary points on reaction pathway from ES to EP. To convert energies at the stationary points to the free energy values, the terms accounting for the

■

MODELS AND METHODS The X-ray structure of the ASPA complex with N-phosphonomethyl-L-aspartate (PDB id: 2O4H)9 presumably mimicking the tetrahedral intermediate analogue was chosen as a source of initial coordinates of heavy atoms. We replaced the phosphorus atom by carbon and added the hydrogen atoms to the substrate and the amino acid side chains. The negatively charged Asp and Glu and positively charged Lys and Arg side chains were assumed. Protonation states of the histidine residues with respect B


Article


Figure 1. Computed energy diagram from the reactants, E + S, to the products, E + P1 + P2. The red asterisk indicates the starting structure for calculations performed either using the QM/MM method (within the limits marked by red arrows) or using MM simulations (within the limits marked by blue arrows).

slightly different Hamiltonians are simulated concurrently. Periodically, exchange between different replicas is performed using the Metropolis Monte Carlo acceptance criterion. By allowing replicas to exchange their configurations frequently, REMD simulations have a better chance to escape from kinetic bottlenecks and efficiently explore all relevant regions of configurational space. Current implementation of the US/HREMD method32 allows one to achieve an excellent performance on large-scale supercomputers. In this work we used the NAMD2.10 package33 and the CHARMM3634−36 force field. First, the structures for each replica were minimized for 2000 steps each with the biasing harmonic potential of 15 kcal/(mol Å2) applied to the distances selected as the reaction coordinates. MD simulations were performed in NPT ensemble at 300 K and 1 atm using Langevin dynamics with a 1 fs time step and Nosé− Hoover Langevin piston method for pressure control. The Coulombic interactions were evaluated with the PME method. In all US/H-REMD simulations the time interval of attempted exchange was set as 100 steps (0.1 ps). Further details of the US/ H-REMD protocol are explained in the Results section and in the Supporting Information. When simulations of product release (EP → E + P1 + P2) were completed, we checked that the obtained structure of the protein without ligands was consistent with the crystal structures of the apoenzyme.

zero-point energy, entropic, and thermal contributions were evaluated at each stationary point. The QM/MM calculations of the minimum-energy reaction profile were carried out with the NWChem program.20 The density functional theory with the PBE0 functional21,22 and 631G** basis set was used. The PBE0 functional was shown to perform well in quantum chemical and QM/MM calculations.23−27 All QM atoms, the link hydrogen atoms, and all atoms assigned to the MM region were included to the optimization protocol. Forces and energies in the MM part were computed with the AMBER all atom force field parameters.28 Nonbonded interactions were cut off at 100 Å in the MM part; the same radius was taken for the QM/MM bare charges around the QM part, thus covering all electrostatic contributions in this system. The electronic embedding scheme29 implemented in NWChem was utilized. Molecular dynamics calculations with the QM/MM potentials (MD-QM/MM) were performed with the CP2K program.30 The main objective of the MD-QM/MM simulations carried out at the vicinity of the ES complex was to analyze the flexibility of the zinc coordination sphere. Another part of the reaction pathway where the MD-QM/MM simulations were important referred to the vicinity of the I3 and I4 intermediates. This was a reaction step when the peptide bond of the substrate was already cleaved but the amino group (−NH2) of aspartate remained deprotonated. The MD-QM/MM modeling allowed us to understand the consequence of the proton transfer events and to select proper reaction coordinates. In all MD-QM/MM calculations the QM part consisted of the same set of atoms as used in QM/MM modeling (Scheme 2). The PBE functional with the dispersion correction D3, the mixed Gaussian DZVP and plane wave basis sets, and the Goedecker−Teter−Hutter pseudopotentials were applied. MD trajectories of 15 ps each were executed in the NVT ensemble at 300 K with 1 fs numerical integration time step. We applied the umbrella sampling Hamiltonian replicaexchange molecular dynamics (US/H-REMD) method31 to compute the free energy profiles for substrate binding and product release. In this approach aiming to enhance the sampling efficiency of MD, several copies of the molecular system under

■

RESULTS Computed Energy Diagram. It is convenient to begin consideration of the results showing the computed free energy diagram for the entire enzyme cycle (Figure 1). Every elementary stage along this route is commented in detail in the forthcoming subsections. Here we illustrate only the basic changes in the chemical structures in the corresponding panels in the upper part in Figure 1: from the reactants (NAA + Wat) to the tetrahedral intermediate I1, then to the low-energy intermediate I2, and further to the species appeared after the amide bond cleavage in NAA. The reaction intermediate I1 corresponds to the starting point of our simulation procedure; its position is marked by the asterisk in Figure 1. This model system was constructed computationally by the motifs of the crystal structure PDB id: 2O4H.9 The QM/ C


Article

The Journal of Physical Chemistry B MM-based calculations allowed us to move from the initial conformation toward the ES complex and in another direction toward the EP complex comprising the acetate (P1) and Lasparate (P2) species trapped at the active site inside the protein. Substrate binding to the enzyme active site and release of both products modeled by the MD simulations resulted in the corresponding free energy profiles which were sewed together at the ES and EP1P2 points with the diagram obtained for the stages of chemical transformations. Structures of the First Reaction Intermediate I1 and of the ES Complex. First, we compare the computed equilibrium geometry configuration of the first reaction intermediate I1 of our model system to that of the crystal structure PDB id: 2O4H.9 Such comparison is necessary to validate the computational protocols used in this work. The crystal structure corresponds to the complex of ASPA with N-phosphonomethyl-L-aspartate covalently bound to the protein. Figure 2 demonstrates a

Figure 3. Part of the active site in the ES complex focusing on the position of the substrate. Here and in other figures, carbon atoms are colored in green, oxygen in red, and nitrogen in blue; distances are given in Å.

tightly fixed at the active site by three arginine residues: Arg63, Arg71, and Arg168. The Ow(Wat)−Cs(NAA) distance is 2.66 Å, which is typical for the enzyme catalyzed substrate cleavage by the hydrolases.37 We note that two salt bridges between the positively charged side chains of Arg71 and Arg168 and the negatively charged carboxyl groups of NAA retain at all reaction steps until the EP point, thus keeping the substrate in the productive orientation during the reaction. The configuration of the ligands coordinated to the zinc ion, His21, Glu24, His116, and Wat (see the central part in Figure 4), may be very flexible, as reported in the literature38,39 for other zinc-dependent hydrolases. To characterize this configuration, we carried out MD calculations with the QM/MM potentials near the ES structure. The border panels in Figure 4 show fluctuations of the distances between Zn2+ and respective atoms of the ligands along the MD QM/MM trajectories. Table S2 lists the QM(PBE0/6-31G**)/MM(AMBER) optimized distances, the MD average values obtained with the QM(PBE/DZVP-PW)/MM(CHARMM) potentials, and the energy minimized distances using the QM(PBE/DZVP-PW)/ MM(CHARMM) calculations. These values obtained from the different computational protocols are well consistent. An important conclusion from the MD-QM/MM simulations is that the structure of the metal coordination shell is fairly conservative. The Zn−Oε2(Glu24) distance is the only parameter which fluctuates considerably since the Oε2 atom of Glu24 can be coordinated either by Zn2+ or forms a hydrogen bond with another water molecule present at the active site. However, this particular fluctuation should not be critical for the reaction mechanism. We also show in the Supporting Information the graph of the Ow(Wat)−Cs(NAA) distances along the MD-QM/MM trajectory (Figure S3) demonstrating that the geometry configuration of the reagents favoring the efficient nucleophilic attack of water (Figure 3) remains stable. Reaction Route before the Amide Bond Cleavage: The ES → TS1 → I1 and I1 → TS2 → I2 Steps. We carefully studied the first reaction step between ES and I1 by locating the transition state structure TS1 and analyzing the forward and backward descent from this saddle point to the ES and I1 stationary points. The top left panel in Figure 5 illustrates changes in the reagent structures. The distance between the carbon atom Cs of NAA and the oxygen atom Ow of the catalytic water (2.66 Å in ES, 1.78 in TS1, and 1.48 in I1) was used as a

Figure 2. Superposition of the computed structure of the reaction intermediate I1 (carbon atoms are shown in green, nitrogen in blue, and oxygen in red) with the crystal structure PDB id: 2O4H.9 The molecular groups of the latter are shown as yellow sticks.

superposition of the molecular groups at the active site of the computed and crystal structures. Table S1 in the Supporting Information lists explicitly the most important geometry parameters. We conclude that the geometry parameters obtained computationally are well consistent with those observed in the crystal structure. The equilibrium geometry configuration of ES was derived from the computationally constructed I1 structure as follows. Regeneration of the lytic water molecule Wat was achieved by a gradual elongation of the carbon−oxygen (Cs−Ow) bond in I1 (see Figure 2) accompanied by deprotonation of the Glu178 side chain. Subsequent unconstrained QM/MM optimization resulted in the ES structure. We checked that this minimumenergy point was characterized by all real vibrational frequencies. Figures 3 and 4 illustrate the most important features of the ES geometry. In Figure 3 we focus on the structural aspects of the first reaction step ES → TS1 → I1, the nucleophilic attack of the water molecule (Wat) on the substrate (NAA) assisted by the Glu178 side chain. Configuration of the reactants is favorable for this reaction: the catalytic water molecule is aligned by the hydrogen bonds with Asn117 and Glu178, while the substrate is D


Article


Figure 4. Center: part of the active site in the ES complex focusing on the zinc coordination sphere (His21, Glu24, His116, and Wat). Panels around the central picture show fluctuations of the distances between Zn2+ and the atoms from its coordination shell observed in the MD-QM/MM trajectories.

Figure 5. Changes in the structures upon the reaction steps ES → TS1 → I1 and I1 → TS2 → I2.

reaction step ES → TS1 → I1 should be transferred to the nitrogen atom Ns of the substrate. Correspondingly, the distance between H of the protonated Glu178 and Ns was chosen as a reaction coordinate at the I1 → TS2 → I2 segment of the reaction pathway. Gradual decrease of this distance from its value at I1 allowed us to locate the saddle point (transition state TS2) and the reaction intermediate I2 (see the top right panel in Figure 5). The proton transfer along the properly aligned hydrogen bonding route Oε1(Glu178)...H+...Ns is characterized by a low

reaction coordinate at this step. The obtained activation barrier at this step estimated as 11.3 kcal/mol is lower than 16 kcal/mol reported in ref 13. On the other hand, a fairly low value of this barrier corresponding to the nucleophilic attack of the catalytic water molecule is consistent with the results of modeling the peptide hydrolysis by the glutamate-assisted mechanism in matrix metalloproteinase 2.14 To facilitate the Cs−Ns bond cleavage in NAA, the proton donated from the lytic water molecule to Glu178 at the first E


Article


Figure 6. Stage of the bond cleavage in NAA.

Figure 7. Reaction step I3 → TS4 → I4 corresponding to reorientation of the Glu178 side chain.

energy barrier, less than 1 kcal/mol. The scissile Ns−Cs bond in NAA shows a tendency to elongation (from 1.44 Å in I1 to 1.57 Å in I2) upon the respective conformational changes and proton transfer. Cs−Ns Bond Cleavage in NAA. The distance between the carbon and nitrogen atoms of the substrate (Cs−Ns) was

selected as a reaction coordinate for the next segment on the reaction route started at I2. The third transition state (TS3) was located in QM/MM calculations by gradually increasing the distance from the initial value 1.57 Å in I2 and optimizing all other geometry parameters. The saddle point TS3 separates the F


Article


Figure 8. Evolution of the system between the I4 intermediate and the EP1P2 complex along the MD QM/MM trajectory. The upper panels correspond to different snapshots.

Starting from the different MD-QM/MM snapshots and using the QM/MM optimization, we were able to locate the stationary points corresponding to the enzyme−product (EP or EP1P2) complex, to the intermediate structure I5 between I4 and EP, and two transition state structures TS5 and TS6. Except for the EP configuration, these points on the shallow energy surface are not of a special interest. Figure 9 illustrates the EP structure showing the zinc coordination shell. It is instructive to compare this picture with

I2 and I3 reaction intermediates. Figure 6 illustrates the structural changes at this step. Conformational Changes and Formation of the Reaction Products. It is principally possible to simulate proton transfer from Glu178 to the acetate moiety when starting from I3 and to arrive to the system composed of the protonated acetic acid and the L-aspartate with a neutral amine group. However, to complete the chemical transformations properly (see Scheme 1) and to arrive to the acetate and aspartate species from the intermediate I3, the proton from Glu178 should be transferred to the amine nitrogen of the aspartic moiety. At this point, it was instructive to carry out MD calculations with the QM/MM potentials starting trajectories in the vicinity of the I3 configuration. Analysis of the computed MD trajectories allowed us to visualize conformational changes of the system primarily affecting the side chain of Glu178 and the acetate moiety as well as the back and forth proton transfer between Glu178 and aspartate. Following such analysis we selected the Oε1−Ns distance as a reaction coordinate for the elementary step I3 → TS4 → I4 (Figure 7) corresponding to reorientation of the Glu178 side chain. Upon this conformation change the distance between the acetate and aspartate species (between the Cs and Ns atoms) increased up to 3 Å. The energy barrier separating the I3 and I4 intermediates is 8 kcal/mol. Figure 8 illustrates the dynamical behavior of the system after the I4 point. The upper panels in Figure 8 show different snapshots along the MD-QM/MM trajectory. We observe that the system fluctuates between few conformations moving along a very shallow energy surface. The acetate moiety can be differently coordinated to Zn2+ depending on which oxygen atom is closer to the metal ion (panels a and b vs panels c and d). In both cases, one of the oxygen atoms from the acetate can form a hydrogen bond with the N-product (aspartate); the proton from Glu178 is transferred to Ns (panels a and c). Alternatively, the acetate stays apart from the N-product, and Glu178 remains protonated (panels b and d).

Figure 9. Zinc coordination shell at the EP configuration according to the QM/MM optimization.

the central panel in Figure 4 corresponding to the ES complex. Clearly, the water ligand in the ES is replaced by the acetate product in the EP; in other respects the zinc environment remains the same. Substrate Binding and Product Release. Inspection of the crystal structures of ASPA8,9 and the structure of the G


Article


Figure 10. Illustration of the substrate binding (the upper left panel), the release of L-aspartate (the upper central panel), and the release of acetate (the upper right panel).

ing Information) gated by the side chains of Asn23, Arg63, Gln184, and Glu290. The pulling vector was chosen along the line connecting the Cα-atom of Leu128 and the carboxyl carbon of acetate. The distances were divided into 16 replicas with an interval of 1 Å and a biasing potential force constant was set to 1.5 kcal/(mol Å2). The US/H-REMD trajectories were executed for 10 ns. We considered the acetate (P1) release from both initial points, EP1P2 and EP2 + P1, resulting in close energy profiles. We show the latter pathway in the diagram in Figure 1. A relatively high energy barrier, 6.5 kcal/mol, can be rationalized since the coordination bond between the positively charged zinc ion and the negatively charged acetate species should be disrupted. It is instructive to compare the structure of the protein obtained computationally after release of the products and the crystal structures PDB id: 2I3C8 and PDB id: 2Q5140 of the apoenzyme which have not been employed in our simulations. Figure 11 shows the superposition of the related molecular groups at the active site. Although a general agreement is satisfactory, we see that the most notable difference refers to the position of the Glu178 side chain (cf. yellow and colored sticks). The deviation exceeds the limits of possible conformational fluctuations in the ensemble of protein molecules as shown in Figure 11 for the side chains of Glu178 and His116 (PDB id: 2Q5140). However, this discrepancy can be easily explained. In the computationally derived structure, the Glu178 residue serves as the fourth ligand of the zinc coordination sphere (the acetate product seen in Figure 9 has been removed from the EP complex). In the crystal structures, one of the phosphate molecules added to the system6 occupies the corresponding spot near Zn2+.

computed EP complex allowed us to suggest the possible entrance channel for substrate binding and the exit channels for release of the product species. We illustrate these transport channels in the Supporting Information (Figure S3). The US/H-REMD protocol was used to compute the corresponding free energy profiles. To insert the substrate, we explored the entrance channel gated by the side chains of Arg71, Tyr164, Lys228, and Glu293. The pulling vector for the US/HREMD simulations was defined along the line between the carbonyl carbon atom of NAA and the zinc ion (see upper left panel in Figure 10). The total distance of 17.4 Å was divided into 16 replicas at a spacing of 1 Å. The force constant of a biasing potential was set to 2 kcal/(mol Å2). The trajectory of the 10 ns length was analyzed. As shown in the energy diagrams the activation barrier of 1.7 kcal/mol to insert the substrate was computed while the binding energy of 6.7 kcal/mol was estimated. Prompted by the structure of the EP complex, we exploited the same channel to extract the L-aspartate from the protein cavity. The pulling vector was aligned along the line between the Cαatom of His116 and Cα-atom of aspartate (see the upper central panel in Figure 10). The distance of 14.2 Å was divided into eight replicas at a spacing of 1.5 Å. The force constant of a biasing potential was 1.2 kcal/(mol Å2). The length of the trajectory was 15 ns. According to these calculations, the release of the aspartate from the enzyme active site to the bulk water requires about 2.5 kcal/mol. We explored several options to extract the acetate product from the active site to the exterior of the protein. In particular, we attempted to use the same transport channel as for release of the aspartate. In this case release of both products should be performed sequentially: first, the aspartate, then the acetate, since the aspartate species blocks spatially the channel. However, this route for the acetate resulted in the high energy expenses. Therefore, we explored another possible channel (see Support-

■

DISCUSSION Comparison of the computationally derived structures of the tetrahedral intermediate I1 and of the protein after release of the H


Article


If we apply such analysis for the energy profile between the points ES and EP1P2 (excluding the segments computed by the MD simulations), then the use of the criteria proposed by Kozuch and Shaik43−45 leads to a pair of the rate-determining states: TS4 with the energy +4.1 kcal/mol and I2 with the energy −9.4 kcal/mol (an assignment is performed using the computer code AUTOF44). The corresponding energy span is 13.5 kcal/ mol, which is better consistent with the “experimental” estimate 16 kcal/mol than the value of 11.3 kcal/mol computed for the step ES → TS1 → I1. If we take into consideration the entire energy profile from E + S to E + P1 + P2 and again use the same algorithm,44 then the rate-determining transition state corresponds to the point denoted EP1* in Figure 1 with the energy +6.8 kcal/mol and the rate-determining reaction intermediate I2 with the energy −9.4 kcal/mol. Under such assignment, the energy span is 16.2 kcal/mol. Although it is tempting to claim that the results of our calculations are in a perfect agreement with the kinetic experimental data we prefer to advocate here for another important conclusion. According to the results of the present simulations, the stages of enzyme recovery, but not the stage of the nucleophilic attack of water, are responsible for the observed reaction rates. The key energy barriers (either TS4 or EP1*) refer to the processes of the product formation and product release. In this respect, a fairly slow release of acetate governed by the barrier at EP1*, which determines the rate of substrate hydrolysis by ASPA, might explain a poor performance of the enzyme and an abnormally high level of NAA associated with the Canavan disease.

Figure 11. Superposition of the molecular groups at the active site in the computed structure of the regenerated enzyme (yellow sticks) with those from the crystal structures of the apoenzyme, PDB id: 2I3C8 (colored sticks). For the side chains of Glu178 and His116 we show a variety of conformations40 suggested by the PDB entry 2Q51.

products to the relevant crystal structures as illustrated in Figures 2 and 11 provides a support to the results of our modeling. It is more difficult to compare the computed and experimental data related to the kinetics of this enzyme reaction. The observed kcat values9,11 can be converted to the activation barriers close to 16 kcal/mol if the simple transitions state theory formula is applied. The calculations reported in this work (Figure 1) resulted in the energy barrier of 11.3 kcal/mol at the elementary step ES → TS1 → I1. Extension of the QM subsystem (described in Supporting Information) resulted in even lower value of 10 kcal/ mol. We point out that the computed activation barriers less than 13 kcal/mol for other zinc-dependent hydrolases have been reported in the literature.14,41 Possible errors in the computational protocol using the conventional quantum chemistry methods should be taken into consideration. In this respect, the paper of Navrátil et al.42 gives the proper estimates. The authors carefully studied quantum chemical calculations of activation energies for the hydrolysis reactions by zinc metalloenzymes and concluded that variations of 3−4 kcal/mol (translating to 2−3 orders in rate constants) were encountered for many approaches. By all means we cannot confirm the conclusion from the papers of Zhang et al.12,13 that the ES → TS1 → I1 stage is responsible for the observed reaction kinetics. Analysis of the energy diagram in Figure 1 shows that other stages beyond the nucleophilic attack of water should be also taken into consideration. In this respect, it seems instructive to apply here the proposal known as the energetic span model by Kozuch and Shaik,43−45 actively used in modern catalysis.46−49 This model allows one to estimate the rate-determining states on the multistage reaction route of the catalyst and to evaluate the activation energy of the catalytic cycle (or the energetic span) δE for use in the Arrhenius-like equation. To this goal, the pairs of the rate-determining transition state with the energy TRDTS and the rate-determining reaction intermediate with the energy IRDI should be identified. For the energy diagram presented in Figure 1 the energetic span is δE = TRDTS − IRDI.

■

CONCLUSION

■

ASSOCIATED CONTENT

There are several important results from our theoretical analysis. For the first time we describe the structures and energetics for the complete catalytic cycle of aspartoacylase, a zinc-dependent enzyme responsible for the hydrolysis of the amino acid derivative, N-acetylaspartate. Using a combination of the quantum-based modeling for the chemical steps of this reaction with the molecular dynamics simulations for the stages of substrate binding to the enzyme active site and release of the products, we were able to construct the energy profile from the reactants E + S to the products E + P. Accordingly, all molecular processes of the catalytic cycle can be visualized at the atomic scale. We also show that the enzyme regeneration, but not the nucleophilic attack of the catalytic water molecule, corresponds to the rate-determining stages of the full catalytic cycle of this enzyme.

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b02542. Details of the QM/MM calculations with an extended QM region (section S1); details of the umbrella sampling Hamiltonian replica-exchange molecular dynamics (US/ H-REMD) method (section S2); characterization of the structures of the ES complex and the tetrahedral intermediate I1 (section S3); view of the transport channels for substrate binding and product release (section S4) (PDF) I


Article


■

(16) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198−229. (17) Merz, K. M., Jr. Using Quantum Mechanical Approaches to Study Biological Systems. Acc. Chem. Res. 2014, 47, 2804−2811. (18) van der Kamp, M. V.; Mulholland, A. J. Combined Quantum Mechanics/Molecular Mechanics (QM/MM) Methods in Computational Enzymology. Biochemistry 2013, 52, 2708−2728. (19) Word, J. M.; Lovell, S. C.; Richardson, J. S.; Richardson, D. C. Asparagine and Glutamine: Using Hydrogen Atom Contacts in the Choice of Side-Chain Amide Orientation. J. Mol. Biol. 1999, 285, 1735− 1747. (20) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; et al. NWChem: A Comprehensive and Scalable Open-Source Solution for Large Scale Molecular Simulations. Comput. Phys. Commun. 2010, 181, 1477−1489. (21) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (22) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0Model. J. Chem. Phys. 1999, 110, 6158. (23) Del Campo, J. M.; Gazquez, J. L.; Trickey, S. B.; Vela, A. NonEmpirical Improvement of PBE and Its Hybrid PBE0 for General Description of Molecular Properties. J. Chem. Phys. 2012, 136, 104108. (24) Grigorenko, B. L.; Khrenova, M. G.; Nilov, D. K.; Nemukhin, A. V.; Švedas, V. K. Catalytic Cycle of Penicillin Acylase from Escherichia coli: QM/MM Modeling of Chemical Transformations in the Enzyme Active Site upon Penicillin G Hydrolysis. ACS Catal. 2014, 4, 2521− 2529. (25) Khrenova, M. G.; Mironov, V. A.; Grigorenko, B. L.; Nemukhin, A. V. Modeling the Role of G12V and G13V Ras Mutations in the RasGAP Catalyzed Hydrolysis Reaction of Guanosine Triphosphate. Biochemistry 2014, 53, 7093−7099. (26) Grigorenko, B. L.; Nemukhin, A. V.; Polyakov, I. V.; Morozov, D. I.; Krylov, A. I. First-Principles Characterization of the Energy Landscape and Optical Spectra of Green Fluorescent Protein along the A→I→B Proton Transfer Route. J. Am. Chem. Soc. 2013, 135, 11541−11549. (27) Grigorenko, B. L.; Nemukhin, A. V.; Polyakov, I. V.; Khrenova, M. G.; Krylov, A. I. A Light-Induced Reaction with Oxygen Leads to Chromophore Decomposition and Irreversible Photobleaching in GFPtype Proteins. J. Phys. Chem. B 2015, 119, 5444−5452. (28) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (29) Bakowies, D.; Thiel, W. Hybrid Models for Combined Quantum Mechanical and Molecular Mechanical Approaches. J. Phys. Chem. 1996, 100, 10580−10594. (30) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (31) Park, S.; Kim, T.; Im, W. Transmembrane Helix Assembly by Window Exchange Umbrella Sampling. Phys. Rev. Lett. 2012, 108, 108102. (32) Jiang, W.; Luo, Y.; Maragliano, L.; Roux, B. Calculation of Free Energy Landscape in Multi-Dimensions with Hamiltonian-Exchange Umbrella Sampling on Petascale Supercomputer. J. Chem. Theory Comput. 2012, 8, 4672−4680. (33) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (34) MacKerell, A. D., Jr.; Feig, M.; Brooks, C. L. Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comput. Chem. 2004, 25, 1400−1415.

AUTHOR INFORMATION

Corresponding Author

*Phone +7-495-939-10-96; e-mail [email protected], [email protected] (A.V.N.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. S. Kozuch for advice and helpful discussions. This study was supported by the Russian Science Foundation (project # 14-13-00124). We acknowledge the use of supercomputer resources of the Lomonosov Moscow State University50 and of the Joint Supercomputer Center of the Russian Academy of Sciences.

■

REFERENCES

(1) Clark, J. F.; Doepke, A.; Filosa, J. A.; Wardle, R. L.; Lu, A.; Meeker, T. J.; Pyne-Geithman, G. J. N-Acetylaspartate as a Reservoir for Glutamate. Med. Hypotheses 2006, 67, 506−512. (2) Moffett, J. R.; Ross, B.; Arun, P.; Madhavarao, C. N.; Namboodiri, A. M. A. N-Acetylaspartate in the CNS: From Neurodiagnostics to Neurobiology. Prog. Neurobiol. 2007, 81, 89−131. (3) Chakraborty, G.; Mekala, P.; Yahya, D.; Wu, G.; Ledeen, R. W. Intraneuronal N-Acetylaspartate Supplies Acetyl Groups for Myelin Lipid Synthesis: Evidence for Myelin-Associated Aspartoacylase. J. Neurochem. 2001, 78, 736−745. (4) Nordengen, K.; Heuser, C. Localization of N-Acetylaspartate in Oligodendrocytes/Myelin. Brain Struct. Funct. 2015, 220, 899−917. (5) Matalon, R.; Michals-Matalon, K.; Sebesta, M.; Deanching, M.; Gashkoff, P.; Casanova, J. Aspartoacylase Deficiency and Nacetylaspartic Aciduria in Patients with Canavan Disease. Am. J. Med. Genet. 1988, 29, 463−471. (6) Zayed, H. Canavan Disease: An Arab Scenario. Gene 2015, 560, 9− 14. (7) Zeng, B. J.; Wang, Z. H.; Ribeiro, L. A.; Leone, P.; De Gasperi, R.; Kim, J.; Raghavan, S.; Ong, E.; Pastores, G. M.; Kolodny, E. H. Identification and Characterization of Novel Mutations of the Aspartoacylase Gene in Non-Jewish Patients with Canavan Disease. J. Inherited Metab. Dis. 2002, 25, 557−570. (8) Bitto, E.; Bingman, C. A.; Wesenberg, G. E.; McCoy, J. G.; Phillips, G. N., Jr. Structure of Aspartoacylase, the Brain Enzyme Impaired in Canavan Disease. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 456−461. (9) Le Coq, J.; Pavlovsky, A.; Malik, R.; Sanishvili, R.; Xu, C.; Viola, R. Examination of the Mechanism of Human Brain Aspartoacylase through the Binding of an Intermediate Analogue. Biochemistry 2008, 47, 3484− 3492. (10) Wijayasinghe, Y. S.; Pavlovsky, A. G.; Viola, R. E. Aspartoacylase Catalytic Deficiency as the Cause of Canavan Disease: A Structural Perspective. Biochemistry 2014, 53, 4970−4978. (11) Herga, S.; Berrin, J.-G.; Perrier, J.; Puigserver, A.; Giardina, T. Identification of the Zinc binding Ligands and the Catalytic Residue in Human Aspartoacylase, an Enzyme Involved in Canavan Disease. FEBS Lett. 2006, 580, 5899−5904. (12) Zhang, C.; Gao, J.-Y.; Chen, Z.-Q.; Xue, Y. Molecular Dynamics and Density Functional Theory Studies of Substrate Binding and Catalysis of Human Brain Aspartoacylase. J. Mol. Graphics Modell. 2010, 28, 799−806. (13) Zhang, C.; Liu, X.; Xue, Y. A General Acid-General Base Reaction Mechanism for Human Brain Aspartoacylase: A QM/MM Study. Comput. Theor. Chem. 2012, 980, 85−91. (14) Vasilevskaya, T.; Khrenova, M. G.; Nemukhin, A. V.; Thiel, W. Mechanism of Proteolysis in Matrix Metalloproteinase-2 Revealed by QM/MM Modeling. J. Comput. Chem. 2015, 36, 1621−1630. (15) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. J


Article

The Journal of Physical Chemistry B (35) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257−3273. (36) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Field. J. Comput. Chem. 2010, 31, 671−690. (37) Nemukhin, A. V.; Grigorenko, B. L.; Rogov, A. V.; Topol, I. A.; Burt, S. K. Modeling of Serine Protease Prototype Reactions with the Flexible Effective Fragment Potential Quantum Mechanical/Molecular Mechanical Method. Theor. Chem. Acc. 2004, 111, 36−48. (38) Blumberger, J.; Lamoureux, G.; Klein, M. L. Peptide Hydrolysis in Thermolysin: Ab Initio QM/MM Investigation of the Glu143-Assisted Water Addition Mechanism. J. Chem. Theory Comput. 2007, 3, 1837− 1850. (39) Wu, R.; Hu, P.; Wang, S.; Cao, Z.; Zhang, Y. Flexibility of Catalytic Zinc Coordination in Thermolysin and HDAC8: A Born−Oppenheimer ab Initio QM/MM Molecular Dynamics Study. J. Chem. Theory Comput. 2010, 6, 337−343. (40) Levin, E. J.; Kondrashov, G. E.; Wesenberg, G. N.; Phillips, N. L., Jr. Ensemble Refinement of Protein Crystal Structures: Validation and Application. Structure 2007, 15, 1040−1052. (41) Ion, B. F.; Kazim, E.; Gauld, J. W. A Multi-Scale Computational Study on the Mechanism of Streptococus pneumonia Nicotinamidase (SpNic). Molecules 2014, 19, 15735−15753. (42) Navrátil, V.; Klusák, V.; Rulišek, L. Theoretical Aspects of Hydrolysis of Peptide Bond by Zinc Metalloenzymes. Chem. - Eur. J. 2013, 19, 16634−16645. (43) Kozuch, S.; Shaik, S. Kinetic-Quantum Chemical Model for Catalytic Cycles: The Heber-Bosch Process and the Effect of Reagent Concentration. J. Phys. Chem. A 2008, 112, 6032−6041. (44) Uhe, A.; Kozuch, S.; Shaik, S. Automatic Analysis of Computed Catalytic Cycles. J. Comput. Chem. 2011, 32, 978−985. (45) Kozuch, S.; Shaik, S. How to Conceptualize Catalytic Cycles? The Energetic Span Model. Acc. Chem. Res. 2011, 44, 101−110. (46) García-Melchor, M.; Braga, A. A. C.; Lledós, A.; Ujaque, G.; Feliu Maseras, F. Computational Perspectives on Pd-Catalyzed C-C CrossCoupling Reaction Mechanism. Acc. Chem. Res. 2013, 46, 2626−2634. (47) Mortensen, P. M.; Grunwaldt, J.-D.; Jensen, P. A.; Jensen, A. D. Screening of Catalysts for Hydrodeoxygenation of Phenol as a Model Compound for Bio-oil. ACS Catal. 2013, 3, 1774−1785. (48) Engelin, C.; Jensen, T.; Rodriguez-Rodriguez, S.; Fristrup, P. Mechanistic Investigation of Palladium-Catalyzed Allylic C−H Activation. ACS Catal. 2013, 3, 294−302. (49) Yang, B.; Burch, R.; Hardacre, C.; Headdock, G.; Hu, P. Correction to Understanding the Optimal Adsorption Energies for Catalyst Screening in Heterogeneous Catalysis. ACS Catal. 2014, 4, 943−943. (50) Voevodin, Vl. V.; Zhumatiy, S. A.; Sobolev, S. I.; Antonov, A. S.; Bryzgalov, P. A.; Nikitenko, D. A.; Stefanov, K. S.; Voevodin; Vad, V. Practice of “Lomonosov” Supercomputer. Open Systems J. (Mosc.) 2012, 7, 36−39.

K


Modeling the Complete Catalytic Cycle of Aspartoacylase - The

Recommend Documents