|Home | About | Journals | Submit | Contact Us | Français|
Protein Arginine Deiminase 4 (PAD4) catalyzes the citrullination of the peptidylarginine via two successive stages: deimination and hydrolysis. Herein, by employing state-of-the-art Born-Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling method, we have characterized the catalytic mechanism of the hydrolysis reaction: first the nucleophilic attack of a water molecule to the Cζ of the thiouronium intermediate yields a stable tetrahedral intermediate, and then the S-Cζ bond breaks to generate the final product, citrulline. Throughout the hydrolysis reaction, His471 and Asp473 play pivotal catalytic roles by first enhancing the nucleophilic ability of the active water through forming shorter and low-barrier hydrogen bonds, and then by serving as proton-accepting groups to deprotonate the water molecule, which is consistent with experimental findings. At the transition state, the spontaneous proton transfer among the reactive water, His471 and Asp473 has been observed. The determined overall free energy barrier for this hydrolysis stage is 16.5 kcal·mol–1, which is lower than the barrier of 20.9 kcal·mol–1 for the deimination stage determined previously with the same computational approach [J. Phys. Chem. B, 113, 12750-12758 (2009)]. Thus the rate-determining step of the PAD4 catalyzed citrullination is the first step of the deimination. Our current work further demonstrates the strength and applicability of the ab initio QM/MM MD approach in simulating enzyme reactions.
Protein Arginine Deiminase 4 (PAD4) is an important post-translational modifying enzyme which catalyzes the deimination of multiple transcriptional regulators, including p300, and histones H2A, H3, H4 [1–5]. It plays a pivotal role in human gene regulation [1–5], and the dysregulation of PAD4 has been linked to the onset and progression of the autoimmune disorder relevant to rheumatoid arthritis [6–9]. Thus, there is significant interest in elucidating its catalytic mechanism and in developing new mechanism-based inhibitors [10–15].
PAD4 is a member of guanidino-group modifying enzyme superfamily (GMSF) , which includes arginine deiminase (ADI) [17–24] and dimethylarginine dimethylaminohydrolase (DDAH) [25–27]. The PAD4 active site is located at the C-terminal domain, of which the structure is comprised of five ββαβ modules. Site directed mutagenesis experiments indicate four key catalytic residues around the active site cleft, including Asp350, His471, Asp473 and Cys645. Based on extensive structural and biochemical studies [13, 19, 20, 25, 27], it has been generally suggested that PAD4 adopts a two-stage catalytic mechanism similar to other GMSF members [11,14,16,19,20,27]: in the initial deimination stage, the nucleophilic attack of the active site Cys645 to the guanidinium carbon of the arginine residue is followed by the cleavage of the Cζ-Nη1 bond to release an ammonia; then in the second reaction stage, the formed covalent thiouronium intermediate is hydrolyzed to generate the product citrulline. Although the proposed reaction scheme is reasonable in general, many detailed mechanistic questions remain elusive.
Very recently, by employing Born-Oppenheimer ab initio QM/MM molecular dynamics simulations and umbrella sampling, we have characterized the Michaelis reactant complex and determined a substrate-assisted proton transfer mechanism for the initial deimination reaction catalyzed by PAD4 , as shown in Fig. 1. Our simulations indicate that both His471 and Cys645 in the PAD4 active site are neutral prior to the reaction; the deprotonation of Cys645 by the substrate arginine occurs in concert with the nucleophilic addition of the Cys thiolate to Cζ of the substrate, which leads to a covalent tetrahedral intermediate; then the Cζ-Nη1 bond cleaves and the resulted ammonia is displaced by a solvent water molecule. The initial deprotonation and nucleophilic attack step is found to be rate-determining. The computed free energy barrier with B3LYP(6-31G*) QM/MM MD simulations and umbrella sampling is 20.9 kcal·mol–1, consistent with the experimental kinetic data. This previous study  sets the stage for our current investigation, which is to theoretically characterize the second hydrolysis stage of the PAD4 catalyzed reaction.
In this work, our employed computational methods and protocols are the same as in our previous study , which center on the Born-Oppenheimer molecular dynamics (MD) simulation with ab initio quantum mechanical/molecular mechanical (QM/MM) potential [28–39] and umbrella sampling method [40–42]. At each time step, the atomic forces as well as the total energy of the QM/MM system are calculated with the B3LYP(6-31G*) QM/MM approach on-the-fly, and Newton equations of motion are integrated. From a series of biased simulations, the potential of mean force (PMF) along the reaction coordinate is obtained with the weighted histogram analysis method (WHAM). This ab initio QM/MM MD approach represents the state-of-the-art method to simulate enzymes [15, 43–46]: it provides an ab initio quantum description of chemical bond formation and breaking process, properly and explicitly models the heterogeneous enzyme environment including both protein and solvent water molecules, and takes account of dynamics of enzyme reaction center and its environment on an equal footing.
The initial structure of the thiouronium intermediate is obtained from our previous study , in which ab initio QM/MM MD simulations had been carried out to study the reaction mechanism of the deimination step of PAD4. For the thiouronium intermediate to be further hydrolyzed, the ammonia molecule needs to di use into the bulk with a water molecule coming into the active site. In this work, the ammonia molecule was replaced by a water molecule, which was orientated to form hydrogen bonds with His471 and Asp473. This system was then neutralized with 7 Cl– atoms and solvated into a ~72×86×74 Å3 TIP3P water rectangular box with a 10 Å buffer distance between the each wall and the nearest solute atoms. The prepared whole simulation system was first equilibrated by a series of minimizations and restrained MD simulations, and then an unrestrained 2.8 ns classic MD simulation was carried out with Amber9 . The trajectory is quite stable. With the resulted snapshot, the QM/MM model was obtained by chopping off the water outside of 27 Å of Cζ atom of the thiouronium intermediate. The resultant system has 9711 atoms in total. Throughout the MD simulations, the periodical boundary condition was employed. Long-range electrostatic interactions were treated with particle mesh Ewald (PME) method [48, 49], and a 8 Å cutoff was introduced for nonbonding interactions, updating the neighbor pair list every 10 steps. The SHAKE algorithm  was applied to constrain all bonds involving hydrogen. A time step of 1 fs was used. All the calculations were conducted with AMBER molecular dynamics package, employing AMBER99SB force field [51, 52] and TIP3P model  for water molecules. Force field parameters for the covalent thiouronium intermediate were derived from that of cysteine and arginine residues, whereas its charge parameters were fitted with the restrained electrostatic potential (RESP) model in Amber9 .
The QM/MM partition is illustrated in Fig. 2. The QM subsystem consists of the side-chains of His471, Asp473, the thiouronium, and one water molecule, which have been described by B3LYP functional with 6-31G(d) basis set. The QM/MM boundary was treated with the pseudobond approach [54, 55]. All other atoms are described by AMBER99SB force field [51, 52] and TIP3P model  for water molecules. The pseudobond QM/MM interface with electrostatic embedding to calculate the energy and forces has been described previously [46,54,56]. The prepared QM/MM system was first minimized with an iterative minimization procedure , and the MM subsystem was further equilibrated by a 500 ps MD simulation with MM force field. Then a 30 ps ab initio QM/MM MD simulation without any restraints for atoms within 20 Å of Cζ of the substrate was carried out. Subsequently, the reaction coordinate driving (RCD) method  was employed to the optimized structure to map out a minimum energy path. For each determined structure along the path, a 500 ps MD simulation with MM force field was performed to equilibrate the MM subsystem, with the QM subsystem being frozen. Finally, the resulting snapshot was used as the starting structure for ab initio QM/MM MD simulation with umbrella sampling. A total of 30 simulation windows were employed. For each window, the total potential energy was biased with a harmonic potential 40 ~ 300 kcal·mol–1 · Å–2 and a 30 ps QM/MM MD simulation was carried out. Consequently, the total length of the ab initio QM/MM MD simulation for the hydrolysis step of PAD4 is 0.9 ns. In QM/MM MD simulations, the forces on atoms in both QM and MM subsystems as well as the total energy are calculated on-the-fly with a pseudobond ab initio QM/MM method at each time step. Newton equations of motion are integrated with Beeman algorithm . Berendsen thermostat method  has been used to control the system temperature at 300 K. The configurations were collected for 20 ps for data analysis after an equilibration period of 10 ps. The probability distribution along the defined reaction coordinate was determined for each window and pieced together with the Weighted Histogram Analysis Method (WHAM) [60, 61] to calculate the PMF.
All the QM/MM calculations have been conducted with modified versions of QCHEM  and TINKER  programs. Throughout the QM/MM minimization and QM/MM MD simulations, the QM subsystem were treated at B3LYP/6-31G(d) level. The spherical boundary condition has been applied so that only the atom inside of 20 Å of the reaction center were allowed to move. A cut off 12 Å for van der Waals interactions, and a cutoff 18 Å for electrostatic interactions among MM atoms were employed. There was no cutoff for electrostatic interactions between QM and MM atoms. The MM atoms were described with AMBER99SB force field [51, 52] and the TIP3P model  was adopted for water molecules.
To characterize the thiouronium intermediate complex which is the reactant state for the hydrolysis reaction, a 30 ps ab initio QM/MM MD simulation has been carried out. Throughout the simulation including an initial 2.8 ns classical MD equilibration, the overall structure of the thiouronium intermediate as well as the enzyme active site has maintained quite well. The active water molecule is held by the three hydrogen bonds with His471, Asp473 and a water molecule, as shown in Figure 3, and the key geometric parameter are listed in Table 1. The average distance between the O atom of the active water (OW for short) and the Cζ atom is 2.89 ± 0.15 Å. The thiouronium intermediate is hydrogen bonded with Asp350 and Asp473, and the Asp473 also forms a hydrogen bond with His644. Overall, the active site of our simulated thiouronium intermediate complex of PAD4 is very consistent with the crystal structure for the thiouronium intermediate complex of ADI , another enzyme in the GMSF family. In their crystal structure , the Cζ-Nε, Cζ-Nη and S-Cζ bonds of the thiouronium intermediate are 1.33, 1.39 and 1.79 Å, respectively, all of which are consistent with the values of 1.33 ± 0.02, 1.33 ± 0.02 and 1.79 ± 0.03 Å in our simulations. It is also observed in the crystal structure  that there is a heavy atom close to the Cζ atom with the distance being 2.56 Å and the heavy atom is very likely to forms hydrogen bonds with His471 and Asp473 with the distances between it and Nη, Oδ2 being 2.66 and 2.49 Å, respectively, which are comparable to the distances of OW-Nη and OW-Oδ2 in our simulated system, as listed in Table 1. This further indicates that our simulated PAD4 thiouronium intermediate for the hydrolysis reaction is very reasonable.
By employing ab initio QM/MM MD simulations with umbrella sampling method, we have determined the free energy profile for the hydrolysis reaction catalyzed by PAD4, as shown in Fig. 4. The overall reaction barrier for this second hydrolysis stage is of 16.5 kcal·mol–1, which is lower than the free energy barrier of 20.9 kcal·mol–1 for the initial deimination stage . This indicates that for the citrullination reaction catalyzed by PAD4, the deimination step is rate-limiting. Our simulations characterize a two-step mechanism for the hydrolysis reaction: first the nucleophilic attack of a water molecule to the Cζ of the thiouronium intermediate yields a stable tetrahedral intermediate, and then the S-Cζ bond breaks to generate the final product, citrulline. The characterized transition states and intermediate are illustrated in Fig. 5 and the key geometric parameters of the corresponding structures are given in Table 1.
In the initial step of the hydrolysis reaction, with the active water approaching the Cζ atom to form the transition state I (TS I), its hydrogen bonds to His471 and Asp473 become significantly shorter. As shown in Table 1, from the thiouronium intermediate to TS I, the OW(Wat)-Nδ(His471) and OW(Wat)-Oδ2(Asp473) distances change from 3.05±0.19 and 2.77±0.11 to 2.73±0.14 and 2.55±0.10 Å, respectively. Meanwhile, Fig. 6 shows the time dependence of two H-OW bonds during the 30 ps QM/MM MD simulation of TS I, in which the elongation of H1-OW or H2-OW bonds to around 1.60 Å indicates the proton transfer to His471 or Asp473. Thus, three representative configurations of TS I have been observed, as illustrated in Fig. 5. It should be noted that no restraint has been applied on hydrogen atoms during the MD simulation, and the proton transfer occurs spontaneously. These results indicate that shorter and low-barrier hydrogen bonds between the water molecule and two nearby basic residues have been formed in the transition state. This finding is consistent with the experimental hypothesis that both His471 and Asp473 play crucial roles in enhancing the nucleophilic ability of the active water . It should be noted that the characterization of such dynamic nature of the transition state demonstrates the strength and power of ab initio QM/MM MD approach.
With the further shortening of the OW-Cζ bond, a stable covalent intermediate is formed with the proton transferring to His471, as illustrated in Fig. 5. At the intermediate, the bond OW-Cζ becomes 1.43 ± 0.03 Å, and bonds S-Cζ, Cζ-Nε and Cζ-Nη get longer to 1.92 ± 0.05, 1.44 ± 0.03 and 1.44 ± 0.03 Å, respectively. H1 is completely bonded to His471 with an bond length of 1.05 ± 0.03 α. For the hydrogen bond network in the active site of the intermediate, most is maintained except the two hydrogen bonds between Asp350 and the substrate (see Fig. 5). By analyzing the snapshots of QM/MM MD simulation on the window of the intermediate, it is found that Asp350 moves away to form hydrogen bonds with the backbone of Gly408 and solvent water molecules. The movement of Asp350 might be caused by the reduction of the positive charge on the substrate during the reaction process, in which a positive charge has been migrated to His471 due to the proton transfer.
In the second step of the hydrolysis reaction, the S-Cζ bond cleaves to generate the final product, citrulline. The barrier of this step is 8.4 kcal·mol–1, while the overall barrier is 16.5 kcal·mol–1 relative to the thiouronium intermediate. Only the S-Cζ bond is chosen as the reaction coordinate for this step, while the proton transfer occurs spontaneously during the QM/MM MD simulations. In the transition state II (TS II), H2 has been mostly transferred to Asp473 with the H2-OW bond being 1.45 ± 0.20 Å. S-Cζ bond is elongated to 2.34 ± 0.01 Å, while OW-Cζ becomes 1.32 ± 0.03 Å, shorter than the normal single O-C bond, which indicates the OW-Cζ partially becomes double bond. With the S-Cζ further getting longer, the citrulline-enzyme complex is formed. To better characterize the citrulline-enzyme complex, a 30 ps ab initio QM/MM MD simulation has been carried out without any restraint, and the structure is illustrated in Fig. 5. The citrulline residue is stable throughout the QM/MM MD simulation, with the bond-lengths of OW-Cζ, Nε-Cζ and Nη-Cζ being 1.26 ± 0.02, 1.36 ± 0.03 and 1.37 ± 0.03 Å, respectively. The hydrogen bonds between the citrulline and His471, Asp473 are conserved. The thiolate group of Cys645 then is hydrogen bonded with Asn588, His644 and a water molecule. No proton transfer is observed among the charged His471, the protonated Asp473 and the ionized Cys645. Therefore, our results indicate His471 and Asp473 not only act as bases to activate the active water but also serve as proton-accepting groups in the hydrolysis step, which is in agreement with the experimental hypothesis . Asp350 is still away from the citrulline, forming hydrogen bonds with Gly408 and a water molecule. Considering that Asp350 is proposed to play a crucial role in substrate-binding [15, 20], this finding indicates that PAD4 is ready for the product release after the hydrolysis step.
In the current work, we have carried out the on-the-fly ab initio QM/MD molecular dynamics with the umbrella sampling method to determine the free energy profile of the hydrolysis stage in the citrullination process catalyzed by PAD4 and characterize the nature of transition states and intermediates. Our calculations indicate that the hydrolysis reaction starts with a water attacking on Cζ of the thiouronium intermediate to form a stable tetrahedral intermediate, and then the S-Cζ bond breaks to generate the final product, citrulline. It is found that His471 and Asp473 play pivotal catalytic roles by first enhancing the nucleophilic ability of the active water through forming shorter and low-barrier hydrogen bonds, and then by serving as proton-accepting groups to deprotonate the water molecule. At the first transition state, the spontaneous proton transfer among the reactive water, His471 and Asp473 has been observed, which further demonstrates the strength and power of ab initio QM/MM MD approach in investigating enzyme reactions. Thus, with our previous study on the initial deimination stage, we have determined the overall reaction mechanism of the citrullination process catalyzed by PAD4. This sets the stage to revisit other important GMSF enzymes, including ADI.
This work was partly supported by NIH (R01-GM079223), NSF (CHE-CAREER-0448156) and the China Scholarship Council. D. X. was supported by the National Natural Science Foundation of China (Grant Nos. 20725312 and 20533060) and the Ministry of Science and Technology (2007CB815201). We thank Dr. Po Hu and Mr. David Rooklin for helpful discussions, and NYU-ITS for providing computational resources.