Computational methods for the study of carboxylases: The case of crotonyl-CoA carboxylase/reductase

The rising levels of atmospheric CO2 and its impact on climate change call for new methods to transform this greenhouse gas into beneficial compounds. Carboxylases have a significant role in the carbon cycle, converting gigatons of CO2 into biomass annually. One of the most effective and fastest carboxylases is crotonyl-CoA carboxylase/reductase (Ccr). To understand its underlying mechanism, we have developed computational methods and protocols based on all-atom molecular dynamics simulations. These methods provide the CO2 binding locations and free energy inside the active site, dependent on different conformations adopted by Ccr and the presence of the crotonyl-CoA substrate. Furthermore, the adaptive string method and quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations outline the CO2 fixation reaction via two different mechanisms. The direct mechanism involves a hydride transfer creating a reactive enolate, which then binds the electrophilic CO2 molecule, resulting in the carboxylated product. Alternatively, another mechanism involves the formation of a covalent adduct. Our simulations suggest that this adduct serves to store the enolate in a much more stable intermediate avoiding its reduction side reaction, explaining the enzyme's efficiency. Overall, this work presents computational methods for studying carboxylation reactions using Ccr as a model, providing general principles that can be applied to modeling other carboxylases.

Conformational Dynamics of the Most Efficient Carboxylase Contributes to Efficient CO2 Fixation

Crotonyl-CoA carboxylase/reductase (Ccr) is one of the fastest CO2 fixing enzymes and has become part of efficient artificial CO2-fixation pathways in vitro, paving the way for future applications. The underlying mechanism of its efficiency, however, is not yet completely understood. X-ray structures of different intermediates in the catalytic cycle reveal tetramers in a dimer of dimers configuration with two open and two closed active sites. Upon binding a substrate, this active site changes its conformation from the open state to the closed state. It is challenging to predict how these coupled conformational changes will alter the CO2 binding affinity to the reaction's active site. To determine whether the open or closed conformations of Ccr affect binding of CO2 to the active site, we performed all-atom molecular simulations of the various conformations of Ccr. The open conformation without a substrate showed the highest binding affinity. The CO2 binding sites are located near the catalytic relevant Asn81 and His365 residues and in an optimal position for CO2 fixation. Furthermore, they are unaffected by substrate binding, and CO2 molecules stay in these binding sites for a longer time. Longer times at these reactive binding sites facilitate CO2 fixation through the nucleophilic attack of the reactive enolate in the closed conformation. We previously demonstrated that the Asn81Leu variant cannot fix CO2. Simulations of the Asn81Leu variant explain the loss of activity through the removal of the Asn81 and His365 binding sites. Overall, our findings show that the conformational dynamics of the enzyme controls CO2 binding. Conformational changes in Ccr increase the level of CO2 in the open subunit before the substrate is bound, the active site closes, and the reaction starts. The full catalytic Ccr cycle alternates among CO2 addition, conformational change, and chemical reaction in the four subunits of the tetramer coordinated by communication between the two dimers.

A simplified charge projection scheme for long-range electrostatics in ab initio QM/MM calculations

In a previous work [Pan et al., Molecules 23, 2500 (2018)], a charge projection scheme was reported, where outer molecular mechanical (MM) charges [>10 Å from the quantum mechanical (QM) region] were projected onto the electrostatic potential (ESP) grid of the QM region to accurately and efficiently capture long-range electrostatics in ab initio QM/MM calculations. Here, a further simplification to the model is proposed, where the outer MM charges are projected onto inner MM atom positions (instead of ESP grid positions). This enables a representation of the long-range MM electrostatic potential via augmentary charges (AC) on inner MM atoms. Combined with the long-range electrostatic correction function from Cisneros et al. [J. Chem. Phys. 143, 044103 (2015)] to smoothly switch between inner and outer MM regions, this new QM/MM-AC electrostatic model yields accurate and continuous ab initio QM/MM electrostatic energies with a 10 Å cutoff between inner and outer MM regions. This model enables efficient QM/MM cluster calculations with a large number of MM atoms as well as QM/MM calculations with periodic boundary conditions.

Structure of the Brain N-Acetylaspartate Biosynthetic Enzyme NAT8L Revealed by Computer Modeling

We report the results of computational modeling of a three-dimensional all-atom structure of the membrane-associated protein encoded by the NAT8L gene, aspartate N-acetyltransferase, which is essential for brain synthesis of N-acetyl-l-aspartate (NAA). The lack of experimentally derived three-dimensional structures of NAT8L poses one of the obstacles in studies of the mechanism of NAA formation and understanding the precise role of NAA in neurological disorders. We apply a computational protocol employing the contact map prediction, ab initio folding, homology modeling, and refinement to obtain a structure of NAT8L with the aspartate and acetyl coenzyme A cofactors in the protein molecule. To verify the computational protocol, we check its predictive power by reproducing the crystal structure of a related N-acetyltransferase domain, specifically, that from the bacterial N-acetylglutamate synthase. We show that the constructed NAT8L model correlates with structural features of the protein revealed in rare experimental studies. The obtained structure of the enzyme active site with the trapped reactants suggests a mechanism of the acetyl transfer upon NAA formation.

DHHC20 Palmitoyl-Transferase Reshapes the Membrane to Foster Catalysis

Cysteine palmitoylation, a form of S-acylation, is a key posttranslational modification in cellular signaling. This type of reversible lipidation occurs in both plasma and organellar membranes, and is catalyzed by a family of integral membrane proteins known as DHHC acyltransferases. The first step in the S-acylation process is the recognition of free acyl coenzyme A (acyl-CoA) from the lipid bilayer. The DHHC enzyme then becomes autoacylated at a site defined by a conserved Asp-His-His-Cys motif. This reaction entails ionization of the catalytic Cys. Intriguingly, in known DHHC structures, this catalytic Cys appears to be exposed to the hydrophobic interior of the lipid membrane, which would be highly unfavorable for a negatively charged nucleophile, thus hindering autoacylation. Here, we use biochemical and computational methods to reconcile these seemingly contradictory facts. First, we experimentally demonstrate that human DHHC20 is active when reconstituted in POPC nanodiscs. Microsecond-long all-atom molecular dynamics simulations are then calculated for human DHHC20 and for different acyl-CoA forms, also in a POPC membrane. Strikingly, we observe that human DHHC20 induces a drastic deformation in the membrane, particularly on the cytoplasmic side, where autoacylation occurs. As a result, the catalytic Cys becomes hydrated and optimally positioned to encounter the cleavage site in acyl-CoA. In summary, we hypothesize that DHHC enzymes locally reshape the membrane to foster a morphology that is specifically adapted for acyl-CoA recognition and autoacylation.

LICHEM 1.1: Recent Improvements and New Capabilities

The QM/MM method has become a useful tool to investigate various properties of complex systems. We previously introduced the layered interacting chemical models (LICHEM) package to enable QM/MM simulations with advanced potentials by combining various (unmodified) QM and MM codes ( J. Comp. Chem., 2016, 37, 1019). LICHEM provides several capabilities such as the ability to use polarizable force fields, such as AMOEBA, for the MM environment. Here, we describe an updated version of LICHEM (v1.1), which includes several new functionalities including a new method to account for long-range electrostatic effects in QM/MM (QM/MM-LREC), a new implementation for QM/MM with the Gaussian electrostatic model (GEM), and new capabilities for path optimizations using the quadratic string model (QSM) coupled with restrained MM environment optimization.

Structural and dynamic basis of substrate permissiveness in hydroxycinnamoyltransferase (HCT)

Substrate permissiveness has long been regarded as the raw materials for the evolution of new enzymatic functions. In land plants, hydroxycinnamoyltransferase (HCT) is an essential enzyme of the phenylpropanoid metabolism. Although essential enzymes are normally associated with high substrate specificity, HCT can utilize a variety of non-native substrates. To examine the structural and dynamic basis of substrate permissiveness in this enzyme, we report the crystal structure of HCT from Selaginella moellendorffii and molecular dynamics (MD) simulations performed on five orthologous HCTs from several major lineages of land plants. Through altogether 17-μs MD simulations, we demonstrate the prevalent swing motion of an arginine handle on a submicrosecond timescale across all five HCTs, which plays a key role in native substrate recognition by these intrinsically promiscuous enzymes. Our simulations further reveal how a non-native substrate of HCT engages a binding site different from that of the native substrate and diffuses to reach the catalytic center and its co-substrate. By numerically solving the Smoluchowski equation, we show that the presence of such an alternative binding site, even when it is distant from the catalytic center, always increases the reaction rate of a given substrate. However, this increase is only significant for enzyme-substrate reactions heavily influenced by diffusion. In these cases, binding non-native substrates ‘off-center’ provides an effective rationale to develop substrate permissiveness while maintaining the native functions of promiscuous enzymes.

Examples abound of enzymes that can process substrates other than their native ones. However, the structural and dynamic basis of this promiscuity remains to be fully understood. In this work, we examine HCT, an intrinsically promiscuous acyltransferase with conserved function in all land plants. We uncover the sub-microsecond swing motion of a key arginine residue facilitating the recognition of both native and non-native substrates of HCT. We also quantify the impact of an off-center binding site on the non-native reaction rate. Although our calculations were inspired by HCT, the results apply in general, i.e., for enzymes heavily influenced by diffusion, binding non-native substrates ‘off-center’, even with rather weak affinity, can accelerate non-native reactions to appreciable levels.

A two-helix motif positions the lysophosphatidic acid acyltransferase active site for catalysis within the membrane bilayer

Phosphatidic acid (PA), the central intermediate in membrane phospholipid synthesis, is generated by two acyltransferases in a pathway conserved in all life forms. The second step in this pathway is catalyzed by 1-acyl-sn-glycerol-3-phosphate acyltransferase, called PlsC in bacteria. Here we present the crystal structure of PlsC from Thermotoga maritima, revealing an unusual hydrophobic/aromatic N-terminal two-helix motif linked to an acyltransferase αβ-domain that contains the catalytic HX₄D motif. PlsC dictates the acyl chain composition of the 2-position of phospholipids, and the acyl chain selectivity 'ruler' is an appropriately placed and closed hydrophobic tunnel. We confirmed this by site-directed mutagenesis and membrane composition analysis of Escherichia coli cells that expressed mutant PlsC. Molecular dynamics (MD) simulations showed that the two-helix motif represents a novel substructure that firmly anchors the protein to one leaflet of the membrane. This binding mode allows the PlsC active site to acylate lysophospholipids within the membrane bilayer by using soluble acyl donors.

A new smoothing function to introduce long-range electrostatic effects in QM/MM calculations

A new method to account for long range electrostatic contributions is proposed and implemented for quantum mechanics/molecular mechanics long range electrostatic correction (QM/MM-LREC) calculations. This method involves the use of the minimum image convention under periodic boundary conditions and a new smoothing function for energies and forces at the cutoff boundary for the Coulomb interactions. Compared to conventional QM/MM calculations without long-range electrostatic corrections, the new method effectively includes effects on the MM environment in the primary image from its replicas in the neighborhood. QM/MM-LREC offers three useful features including the avoidance of calculations in reciprocal space (k-space), with the concomitant avoidance of having to reproduce (analytically or approximately) the QM charge density in k-space, and the straightforward availability of analytical Hessians. The new method is tested and compared with results from smooth particle mesh Ewald (PME) for three systems including a box of neat water, a double proton transfer reaction, and the geometry optimization of the critical point structures for the rate limiting step of the DNA dealkylase AlkB. As with other smoothing or shifting functions, relatively large cutoffs are necessary to achieve comparable accuracy with PME. For the double-proton transfer reaction, the use of a 22 Å cutoff shows a close reaction energy profile and geometries of stationary structures with QM/MM-LREC compared to conventional QM/MM with no truncation. Geometry optimization of stationary structures for the hydrogen abstraction step by AlkB shows some differences between QM/MM-LREC and the conventional QM/MM. These differences underscore the necessity of the inclusion of the long-range electrostatic contribution.

The mechanism of citryl-coenzyme A formation catalyzed by citrate synthase

The enzyme citrate synthase is used by all living cells to catalyze the first step of the citric acid cycle. In this work, we have investigated the enolization and condensation steps catalyzed by citrate synthase, using ab initio (B3LYP/def2-TZVP and MP2/aug-cc-pVDZ) quantum chemical/molecular mechanical hybrid potentials in conjunction with reaction-path-location algorithms and molecular dynamics free energy simulations. The results of the latter indicate that the catalytic His238 residue is in its neutral form, and also argue strongly for the presence of a water molecule in the enzyme's catalytic center. Such a water is observed in some, but not all, of the experimentally resolved structures of the protein. The mechanism itself starts with an enolization that proceeds via an enolate intermediate rather than the enol form, which is much more unstable. This is in agreement with the results of other workers. For the condensation step, we investigated two mechanisms in which there is a direct nucleophilic attack of the enolate intermediate on the oxaloacetate carbonyl carbon, and found the one in which there is no proton transfer from the neighboring arginine to be preferred. Although this residue, Arg329, is not implicated directly in the reaction, it helps to stabilize the negative citryl-CoA formed during the condensation step.

Mechanism of activation of elongation factor Tu by ribosome: catalytic histidine activates GTP by protonation

Elongation factor Tu (EF-Tu) is central to prokaryotic protein synthesis as it has the role of delivering amino-acylated tRNAs to the ribosome. Release of EF-Tu, after correct binding of the EF-Tu:aa-tRNA complex to the ribosome, is initiated by GTP hydrolysis. This reaction, whose mechanism is uncertain, is catalyzed by EF-Tu, but requires activation by the ribosome. There have been a number of mechanistic proposals, including those spurred by a recent X-ray crystallographic analysis of a ribosome:EF-Tu:aa-tRNA:GTP-analog complex. In this work, we have investigated these and alternative hypotheses, using high-level quantum chemical/molecular mechanical simulations for the wild-type protein and its His85Gln mutant. For both proteins, we find previously unsuggested mechanisms as being preferred, in which residue 85, either His or Gln, directly assists in the reaction. Analysis shows that the RNA has a minor catalytic effect in the wild-type reaction, but plays a significant role in the mutant by greatly stabilizing the reaction's transition state. Given the similarity between EF-Tu and other members of the translational G-protein family, it is likely that these mechanisms of ribosome-activated GTP hydrolysis are pertinent to all of these proteins.

A hybrid elastic band string algorithm for studies of enzymatic reactions

A common challenge in theoretical biophysics is the identification of a minimum energy path (MEP) for the rearrangement of a group of atoms from one stable configuration to another. The structure with maximum energy along the MEP approximates the transition state for the process and the energy profile itself permits estimation of the transition rates. In this work we describe a computationally efficient algorithm for the identification of minimum energy paths in complicated biosystems. The algorithm is a hybrid of the nudged elastic band (NEB) and string methods. It has been implemented in the pDynamo simulation program and tested by examining elementary steps in the reaction mechanisms of three enzymes: citrate synthase, RasGAP, and lactate dehydrogenase. Good agreement is found for the energies and geometries of the species along the reaction profiles calculated using the new algorithm and previous versions of the NEB and string techniques, and also those obtained by the common method of adiabatic exploration of the potential energy surface as a function of predefined reaction coordinates. Precisely refined structures of the saddle points along the paths may be subsequently obtained with the climbing image variant of the NEB algorithm. Directions in which the utility of the methods that we have implemented can be further improved are discussed.