CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
503
Restoring Histone Deacetylase Activity by Waste Product Release. A View from Molecular Mechanics Simulations with Mammalian HDAC8 by Francesco Pietra Accademia Lucchese di Scienze, Lettere e Arti, Classe di Scienze, Palazzo Ducale, Lucca I-55100 (phone/fax: þ 39-0583-417336; e-mail: [email protected])
HDAC8 is a ZnII-based, single-peptide mammalian histone deacetylase that is localized mainly in the cytoskeleton of smooth muscle cells, thus regulating muscle contractility. HDACs are also widely involved in cellular processes, ranging from cell differentiation to proliferation, senescence, and apoptosis; in particular, protecting a telomerase activator from ubiquitin-mediated degradation. How HDACs can eliminate the hydrolytic reaction products, in order that the process of deacetylation of the acetyllysine moiety of histones can take place again, has long been debated in the scientific literature, without reaching any firm conclusion, however. This question is the subject of the present work, carried out along a theoretical line that is capable of describing the whole pathway followed by the acetate product (ACT). A model was built here on the crystal data for the Y306F-mutated HDAC8 complex with a diacetylated peptide of the p53-tumor-suppressor class. That was followed by manually hydrolyzing the acetylated moiety bound to ZnII and discharging the monoacetylated peptide product (MAP). The latter was replaced by a H2O molecule bound to ZnII, while ACT was left free in the reaction cage. This ZnII cluster was DFT-parameterized for the ff99SB force field without any further bias. As the result of random-acceleration molecular dynamics (RAMD) simulations, egress of ACT from the reaction cage toward the aqueous environment can follow three pathways. Two of them utilize the channel for peptide (or histone) uptake and are preferred, if ACT leaves the reaction center before MAP (or the deacetylated histone). The third pathway, developing along the internal channel, is available to ACT even if MAP is still in place.
1. Introduction 1). – 1.1. HDAC Isozymes: Overview. Histone deacetylases (HDACs) are a group of isoenzymes that, along other functions, regulate transcription. They do so by catalyzing deacetylation at the e-N-acetyllysine branch of histones, i.e., by reversing the process of acetylation catalyzed by histone acetyltransferase [1]. As these processes are modulated during malignant transformations of cells, HDACs are widely involved in cellular processes, ranging from cell differentiation to proliferation, senescence, and apoptosis. Various diseases can grow on abnormal HDAC activity, against which the most commonly pursued remedies are to inhibit the enzyme. Therefore, HDAC inhibitors have become most chased targets in cancer research [2], although several HDAC knockouts are known to be lethal during mammalian 1)
Abbreviations: aa, amino acid(s); ACT, acetate product; CPU, central processing unit; CUDA, computer unit device architecture; ff, force field; MD, molecular dynamics; MNA, monoacetylated peptide resulting from the hydrolytic reaction; NPT, number-of-particles pressure temperature (isothermal¢isobaric ensemble); NVT, number-of-particles volume-temperature (canonical ensemble); PDB, protein data bank; QM, quantum mechanics; RAMD, random-acceleration MD; RMSD, root-mean-square deviation; ts, time step; VdW, Van der Waals. Õ 2015 Verlag Helvetica Chimica Acta AG, Zîrich
504
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
embryonic development [3]. Also non-histone HDAC8 substrates have been frequently taken into consideration [1]. HDAC Isoenzymes range in size from 300 to 1000 amino acids, in various classes with different sequences. Despite the differences, the catalytic center of this enzyme family is conserved, pivoting on a divalent metal ion. In vitro, the metal specificity increases in the order NiII < ZnII < FeII < CoII [4], while, in vivo, the metal can only be ZnII or FeII, as far as it is presently known [5]. Studies of mammalian HDACs have been concentrated on HDAC8, as this is a single polypeptide easily obtained in active form [5]. It is localized mainly in the cytoskeleton of smooth muscle cells. The most powerful HDAC8 inhibitor discovered to date is largazole, a cyclodepsipeptide isolated from a cyanobacterium, Symploca sp., which grows on coral reefs of Florida Keys [6]. Largazole acts as a prodrug, becoming bioactive on hydrolysis of the thioester side chain, which unmasks a thiolate function. It is the latter that, by binding to the divalent metal ion, inhibits HDAC8 [7]. The crystal structure of various complexes between HDAC8 and hydroxamate or other thiolate inhibitors have also been elucidated [7], while complexes with acetylated protein substrates in crystalline form have long been elusive. The breakthrough was achieved with a catalytically inactive human Y306F mutant HDAC8 enzyme. This allowed to trap a diacetylated peptide, (Ac)-l-Arg-l-His-l-Lys(e-Ac)-l-Lys(e-Ac), containing a fluorogenic coumarin group at its carboxy terminus, patterned on the p53 tumor suppressor protein, as a model for histones [8]. Present work started from this complex, as illustrated below. 1.2. How HDAC8 Bioactivity Can Be Regenerated Following Substrate Deacetylation. Deacetylation at the e-N-acetyllysine branch of histones was proposed to be carried out by a ZnII-coordinated H2O nucleophile under general acid-base catalysis with either one or two histidine side chains [1]. In order that this process, i.e., acetylation of a fresh histone by histone acetylase followed by deacetylation by histone deacetylase, repeats, the latter enzyme has to be freed from the reaction products at every cycle. How ACT can be eliminated from histone deacetylase was investigated here along MD simulations, starting from the above complex between inactivated HDAC8 and p53-derived diacetylated ligand of PDB code 2 V5W [8]. This represents one of many possible models [8], while ACT is the product of any histone deacetylase. Moreover, while the deacetylated histone, due to its large size, can only retreat along the entrance channel, the small ACT molecule could find alternative exit routes. This is why the fate of the ACT product of the deacetylase reaction was given primary attention here. This work, completed in the year 2012, was waiting for compilation, while new QM/ MM studies have appeared as to the mechanism of hydroxamate inhibition of these enzymes [9], based on the same complex used in our work [8]. Previous hypotheses, outlined above, were mainly consolidated. 2. Results and Discussion. – 2.1. The Model. The above mentioned complex between inactivated HDAC8 and p53-derived diacetylated ligand of PDB ID 2 V5W [8] is constituted of the protein 15 helices and 8 sheets, the catalytic ZnII center, and the diacetylated peptide substrate, as illustrated in Fig. 1. The active center comprises ZnII, coordinated with the D178 and D267 carboxylates, as well the H180 imidazole moiety, the catalytic H2O molecule, and the C¼O O-atom of the lysine-branch e-N-Ac group of
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
505
Fig. 1. Stereoview of the complex between Y306F mutant HDAC8 (PDB ID 2 V5W) and a p53-derived diacetylated peptide, (Ac)-l-Arg-l-His-l-Lys(e-Ac)-l-Lys(e-Ac), containing a fluorogenic coumarin group at its carboxyl terminus [8]. The iron protein comprises 15 helices (helix number (in italics), residue numbers, color: 1, 21 – 29, dark gray; 2, 36 – 48, yellow; 3, 49 – 54, green; 4, 63 – 68, dim gray; 5, 72 – 84, black; 6, 90 – 97, forest green; 7, 107 – 128, hot pink; 8, 156 – 166, deep blue; 9, 182 – 189, sea green; 10, 219 – 223, cyan; 11, 236 – 256, cornflower blue; 12, 280 – 293, magenta; 13, 307 – 324, white; 14, 336 – 341, orange; 15, 358 – 375, purple), other than eight sheets, highlighted in black, and ZnII, shown as a red ball, coordinated with D178, H180 and D267, also highlighted in red. The diacetylated peptide is highlighted in cyan. To aid comparisons, here and in all other figures in this work, the orientation of the complex follows that in [10].
the diacetylated peptide [8]. Mutation of Y306 to F306 prevents stabilization of the incipient amine that would be generated during the hydrolytic reaction, with the result that all the above components are trapped into the reaction cage. A model of the presumed arrangement of this system following the hydrolytic reaction was built here by first manually replacing F306 with Y306, thus restoring the naturally occurring active protein. It was then arbitrarily assumed that the monoacetylated peptide product (MAT) leaves the catalytic center before the ACT. Thus, the peptide was removed and replaced by a H2O molecule, while ACT was left free in the reaction cage (Fig. 2, where residues considered to be important for the functionality of the protein [10] are also highlighted, as described in the figure legend). It was now time to create a ff for the ZnII cluster of HDAC8. So far, ZnII complexes with macromolecules have been treated according to either the nonbonded model, the cationic dummy-atom approach, or the bonded model. The nonbonded model places an integer charge on the metal ion, such as with HDAC complexes with inhibitory hydroxamic acids [11]. This may become computationally expensive, as it requires the inclusion of very-long-range electrostatic interactions, while being also prone to structural deformation on prolonged MD simulations. The cationic dummy approach [12] places dummy cations to mimic valence electrons around the metal ion and is integrated with an electrostatic model. This may present the advantage of avoiding a rigid conformation at the active site, but it has been scarcely validated. The bonded
506
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
Fig. 2. Stereoview of the HDAC8 model used in present work, obtained from the complex of Fig. 1 by replacing F306 with naturally occurring Y306. The diacetylated peptide was replaced by a H2O molecule bound to ZnII and a free molecule of ACT in the hydrolytic reaction cage. By comparison with Fig. 1, it can be recognized that F207 (cyan) lies along the channel for uptake of the diacetylated peptide substrate. Southeastern to F207, the Zn cluster is seen in color by element (C, gray; N, blue; O, red; Zn, dark gray). Also illustrated are R37 (green), G303 (hot pink), and G139 (red), which are considered to take part in the opening/closing of the internal channel [10]. Finally, G305 (orange) is also illustrated.
model has been recently automated for ZnII protein complexes in the framework of the ff99SB ff [13], however neglecting dihedral terms and Lennard¢Jones parameters, which were considered of minor importance. Either the experimental structure was used or the metal cluster was calculated via DFT, and an electrostatic model was added [14]. Present work makes use of automated generation of the ff developed in the framework of the bonded model [15]. As a basis, this takes the ff99SB ff [13], which is known to treat proteins adequately, while allowing minimization with the coordinates of all the protein atoms considered as free variables. Any bias, other than those intrinsic to the ff99SB ff, was avoided, so that dihedral terms and Lennard¢Jones parameters were included. Like other popular non-polarizable ffs, ff99SB represents the electrostatic field through fictitious partial point charges distributed along the various atoms. To this regard, RESP charges [13] were adopted here. This modeling could be straightforwardly applied to ACT at HF/6-31G* level, similarly to any other standard protein component. The strategy for the ZnII cluster pivoted on the additive properties of the ff99SB ff [13], so as to restrict as much as possible the area to submit to QM calculations. To this end, the amino acid ligands D178 and H180 were cut, changing their atom CB to a Me group. In the ff, the Me group was removed and each remain was connected to the protein through a HN¢CH(Me)¢C linker, parameterized at the same theory level as ACT. In contrast, due to the presence of a transition metal, the HF approach could not be used for the ZnII cluster.
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
507
Recourse was made to the DFT theory with b3lyp hybrid functional and 6-311G* basis set. RESP Charges were calculated accordingly. Non-bonded terms for ZnII were taken from [16]. 2.2. MD of ACT Disposal. Classical MD, carried out with NAMD engine [17] for the ff99SB ff [13], proved satisfactory for the stated purpose. Tight convergence minimization resulted in a slight flattening around ZnII and the three bound amino acids with respect to a tetrahedral arrangement (Fig. 2). That slightly distorted arrangement was conserved during MD equilibration (Fig. 2), and RMSD proved satisfactory. It is difficult to appreciate whether this slight conformational discrepancy adds to the many cases of problems encountered by current ffs in representing the experimental structure along minimization and MD [18]. In fact, the ZnII cluster modeled here only represents a possible situation during the hydrolytic reaction, not an experimental structure. The next decision to be taken was about how to follow the pathway of ACT disposal under dynamic conditions. Experimental and theoretical methods to follow moving ligands through proteins have been recently reviewed [19]. While no experimental methodology proved able to describe whole pathways, various theoretical approaches have been applied. Of these [19], RAMD [20] was chosen here. This is a form of accelerated MD, where a tiny external force is applied to the center of mass of the molecule to move. The direction of the external force is applied randomly, and the force is made to continue only if the ligand moves along a direction for a selected period of time. If not, a new direction is chosen, randomly again. A small integration timestep was used here, so that all chemical bonds, except the solvent H2O, could be left unrestrained. Also, a large number of runs with different initial direction of the external force were carried out to reach statistical significance. In defense of RAMD, with the system myoglobin and a diatomic gas [21], this theory allowed identification of all important pathways that had been previsiously discovered through computationally expensive, massive MD [22]. RAMD proved also useful for larger ligands with a variety of proteins, as recently reviewed [23]. The results of RAMD simulations for present HDAC8-ACT model system are compiled in the Table. An acceleration of 0.28 kcal/è · g proved to be a good compromise between applying a small force that does not distort VdW barriers, while gaining in MD speed. Generated in a VdW cubic cage, and bordered by protein components as described in the Table, ACT was observed to reach the external solvent along three major routes, with gates as represented in Fig. 3. Two of them, with ACT egressing through either F207 and F208 or D101 and F208, utilize the channel for peptide (or histone) uptake, and are preferred if ACT leaves the reaction center before MAP (or the deacetylated histone). Residence time into the ACT cage is very substantial in both cases, while an intermediate binding pocket (BP) was only observed along the route toward D101-F208, and with a short residence time. The route through F207 and F208 is described in Fig. 4, with ACT represented by its centroid as a red ball of radius 0.5 è, and the diacetylated peptide of Fig. 1 manually superimposed after completion of the MD movie. It is seen that ACT first attempts to egress eastward, then takes the empty entrance channel. Should this be occupied by MAP, the lysine branch would inhibit the passage. Therefore, this route of egress, and likewise that through D101 and F208, can only be taken by ACT if it leaves before MNA.
508
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
Table. Statistical Recapitulation of ACT Pathways inside HDAC8, and Egress toward the Solvent, under RAMD a ) Representative residues for gates and BPs b )
Gate/total number of trajectories
Pathway/number of trajectories per pathway
F 207-F 208/7
ACT_cage ! gate/4 gate: F 207, F 208 ACT_cage > SBP ! gate/3 SBP: I34, R37, F142, G303, G305, Y307; gate: F 207, F 208
Residence time [ps] of ACT into BPs and total trajectory length [ps] ACT_cage: 31-450; gate: 38-492 ACT_cage: 115-270; SBP: 10-13; gate: 296-344
D101-F 208/5 ACT_cage ! gate
gate: D101, F 207
ACT_cage: 2-172, gate: 5-240
K33-N307/5
gate: K33-N307
ACT_cage: 150-800; gate: 300-950 ACT_cage: 171-495; SBP: 5-10; gate: 374-700
ACT_cage ! gate/2
ACT_cage ! SBP ! gate/3 SBP: I34, R37, F142, G303, G305, Y307; gate: K33-N307 NO EGRESS/9
ACT_cage
ACT_cage : R37, G140, ACT_cage: 500-1000 H143, C153, Q263, M274, G303, Y306
a ) Based on X-ray diffraction data of HDAC8 in complex with with a diacetylated peptide, PDB ID 2 V5W [8]. Residue F 306 was replaced with Y306. The diacetylated peptide substrate was replaced with a H2O molecule bound to ZnII and a free molecule of acetate in the hydrolytic reaction cage. RAMD Acceleration of 0.28 kcal/è · g and ts ¼ 1fs. b ) At 5 è from ACT; HDAC8 numbering.
Fig. 3. Stereoview of gates used by ACT to egress toward the solvent from the hydrolytic center of the complex of Fig. 2, according to RAMD simulations. As compiled in the Table, the three major gates should be imagined through F207-F208, F208-D101, and K33-N307. The first two gates correspond to ACT egress from the diacetylated peptide channel of Fig. 1, while the latter gate fits for one of the exit channels proposed from the examination of static X-ray diffraction structures [10].
The third pathway, with ACT egress through K33 and N305 (Table), develops along the internal channel, thus being available to ACT even if MAP is still in place. A snapshot of the movie at completion of the trajectory (Fig. 5) is instructive. Thus, in line
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
509
Fig. 4. Representative egress pathway of ACT toward the solvent along the diacetylated peptide channel. It is seen that ACT, leaving the hydrolytic reaction cage, egresses through F207 (cornflower blue) and F208 (orange). However, by superimposing the diacetylated peptide of Fig. 1 (cyan), it is seen that the acetylamino branch at the reaction center would constitute a barrier to ACT escape along this route. Therefore, this egress route is only conceivable, if MNA resulting from the hydrolytic reaction leaves the channel before ACT. Note: ACT is represented by its centroid as a red ball of radius 0.5 è, except at the exit portal, where it is represented in full atomic shape. The time interval between consecutive centroids in this movie is 0.5 ps.
Fig. 5. Representative egress pathway of ACT toward the solvent along the internal channel. It is seen that ACT, leaving the hydrolytic cage, first tries unsuccessfully to egress along the diacetylated peptide channel of Fig. 1, unable, however, to overcome the barrier opposed by F207 and F208 (both white, not labeled). Then, ACT tries unsuccessfully to egress from the opposite side, unable to overcome the barrier opposed by W141 (cyan). Finally, ACT finds a way toward the solvent by taking the internal channel [10], egressing through F306 (hot pink) and R37 (blue). Note: as in Fig. 4, except 1 ps between consecutive centroids.
with data in the Table, it is seen that ACT (represented by the centroid, like in Fig. 5) spends a considerable time into the hydrolytic cage, then tries unsuccessfully to egress along the diacetylated peptide channel of Fig. 1, where, however, it proves unable to
510
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
overcome the barrier opposed by F207 and F208 (in white, not labeled). Also attempts at egressing from the opposite side proved unsuccessful, ACT being unable to overcome the barrier opposed by W141 (cyan). Then, RAMD captured the relatively rare event of ACT egressing through F306 (hot pink) and R37 (blue). This corresponds to an exit internal channel that was hypothesized from the examination of static models [10]. It remains unobstructed, even if MNA resulting from the hydrolytic reaction is still in place. 3. Conclusions. – This work investigated how histone deacetylase HCAC8 can get free from the acetate product in order that the hydrolytic process repeats. A model was built on X-ray diffraction data for a complex between inactivated HCAC8 and an acetylated peptide, parameterizing the active ZnII cluster, along a bonded model. Recourse was made to the DFT theory with b3lyp hybrid functional at 6-311G* basis level, without any further bias than those intrinsic to the ff99SB ff. The system, solvated with TIP3P H2O in a periodic box, proved stable on MD simulations. Pathways for expulsion of the acetate product proved amenable to study under RAMD. The outcome is that the acetate product leaves preferentially along the channel for uptake of the diacetylated peptide substrate, if that is free. Otherwise, the internal channel remains available for ACT expulsion. I warmly thank Prof. FranÅois-Yves Dupradeau for invaluable help with the parameterization of the ZnII cluster on his R.E.D. development server. Computational Details. – All MD and RAMD simulations were carried out with software NAMD [17], CUDA v. 2.8, running on Debian GNU/Linux amd64. System setup and analysis were carried out with both VMD [24], v. 1.9, and CHIMERA [25], v. 1.5.3 software, the latter also serving for all the illustrations in this work. The Protein Model. The HDAC8 model was built on the PDB ID 2 V5W, i.e., the X-ray diffraction derived Cartesian coordinates, at 2.0-è resolution, for the point-mutant Y306F human HDAC8¢substrate complex, expressed on Escherichia coli [8]. In physiological solns., the protein exists as a monomer. Therefore, only chain A, with its 80 molecules of crystallization H2O, and the acetylated peptide substrate, were retained. Also, mutation back F306Y was carried out manually, maintaining the F306 orientation. Lacking residues 1 – 14 and 378 – 388 were not regenerated. The H-atoms were added with software REDUCE [26], under the Ð-buildÏ flag. A H2O molecule was made to replace the Ac group of peptide substrate, conserving the position of ligand O-atom. Manually hydrolyzed acetate was displaced just enough to give room to the H2O molecule. All remaining parts of the peptide substrate were removed. Parameterization and MD. Parameterization of the Zn complex, ACT, and linker was carried out on R.E.D. development servers [15], with software GAUSSIAN [27] for the QM computations. The Zn complex comprised ZnII and its ligands D178, H180, and Y306, in addition to the above mentioned H2O molecule. The amino acid ligands were cut so as to have their atom CB becoming a Me group. In the ff generated on the R.E.D. server, the Me group was removed, so that a HN¢CH(Me)¢CO moiety was used to bind the ligands to the protein. HF QM Computations with basis set 6-31G* were carried out for both ACT and the HN¢CH(Me)¢C linker, used to join the Zn complex to the protein. For the Zn complex, recourse was made to the DFT theory with b3lyp hybrid functional and 6-311G* basis set. In all cases, at tight convergence for geometry optimization. Preparation of the topology/parameter files was carried out by merging the ff derived above with ff99SB ff for the protein 363 standard amino acids and its 80 molecules of crystallization H2O by means of software LEAP, as available in AMBERTOOLS [13], v. 1.5. The ensemble was solvated with 15,356 molecules of TIP3 H2O in a rectangular box of x ¼ 86, y ¼ 76, and z ¼ 86 è, and neutralized with 4 Cl ¢ ions. It should be noted that, because of LEAP features [13],
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
511
chain A begins with residue number 1. Therefore, 14 had to be added to the residue number read from simulations in order to refer to the crystal numbering, which was adopted throughout this report. The ensemble was minimized at constant volume and then equilibrated for 10 ns, at 1 atm (Langevin piston), by gradually raising the temp. to 300K (Langevin thermostat with damping, 1). The Zn¢ACT distance was maintained at 4.1 è by applying harmonic forces (colvars). Flags: 1 – 4 scaling, 1 è; cutoff, 9.0 è; switchdist, 6.0 è; and pairlistdist, 11.0 è. The same settings, without applied colvars, were also used for RAMD simulations. Integration for all these computations was at 1-fs timestep with no rigid bonds for any component as a safe side after that it was noticed that the use of rigid bonds altered the RAMD response with the myoglobin¢O2 system [21]. A RAMD acceleration of 0.28 kcal/è · g turned out, from simulations in the range 0.1 – 0.3 kcal/è · g, to be the minimum value to bring consistently ACT out of the solvent box. The direction of acceleration was set to change randomly, if ACT traveled by less than 0.02 è in the evaluated block.
REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]
N. A. Wolfson, C. A. Pitcairn, C. A. Fierke, Biopolymers 2012, 99, 112. M. Dokmanovic, C. Clarke, P. A. Marks, Mol. Cancer Res. 2007, 5, 981. M. Haberland, R. L. Montgomery, E. N. Olson, Nat. Rev. Genet. 2009, 10, 32. S. L. Gantt, S. G. Gattis, C. A. Fierke, Biochemistry 2006, 45, 6170. D. P. Dowling, S. L. Gantt, S. G. Gattis, C. A. Fierke, D. W. Christianson, Biochemistry 2008, 47, 13554. K. Taori, V. J. Paul, H. Luesch, J. Am. Chem. Soc. 2008, 130, 1806. K. E. Cole, D. P. Dowling, M. A. Boone, A. J. Phillips, D. W. Christianson, J. Am. Chem. Soc. 2011, 133, 12474. A. Vannini, C. Volpari, P. Gallinari, P. Jones, M. Mattu, A. Carfı, R. Defrancesco, C. Steinkîhler, S. Di Marco, EMBO Rep. 2007, 8, 879. K. Chen, X. Zhang, Y.-D. Wu, O. Wiest, J. Am. Chem. Soc. 2014, 136, 11636. S. Haider, C. G. Joseph, S. Neidle, C. A. Fierke, M. J. Fuchter, Biorg. Med. Chem. Lett. 2011, 21, 2129. D.-F. Wang, O. Wiest, P. Helquist, H.-Y. Lan-Hargest, N. L. Wiech, J. Med. Chem. 2004, 47, 3409. P. Oelschlaeger, R. D. Schmidt, J. Pleiss, Protein Eng. 2003, 16, 341. D. A. Case, T. E. Cheatham III, T. Darden, H. Gohlke, R. Luo, K. M. Merz Jr., A. Onufriev, C. Simmerling, B. Wang, R. Woods, J. Comp. Chem. 2005, 26, 1668. M. B. Peters, Y. Yang, B. Wang, L. Fîsti-Molnr, M. N. Weaver, K. M. Merz Jr., J. Chem. Theory Comp. 2010, 6, 2935. E. Vanquelef, S. Simon, G. Marquant, E. Garcia, G. Klimerak, J. C. Delepine, P. Cieplak, F.-Y. Dupradeau, Nucl. Acids Res. 2011, 39, W511. S. C. Hoops, K. W. Anderson, K. M. Merz Jr., J. Am. Chem. Soc. 1991, 113, 8262. J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, K. Schulten, J. Comput. Chem. 1998, 26, 1781. P. Koehl, M. Levitt, Nature Struct. Biol. 1999, 6, 108. F. Pietra, Chem. Biodiversity 2013, 10, 86. S. Lîdemann, V. Lounnas, R. C. Wade, J. Mol. Biol. 2000, 303, 797. F. Pietra, Chem. Biodiversity 2011, 8, 1410. J. Z. Ruscio, D. Kumar, M. Shukla, M. G. Prisant, T. M. Murali, A. V. Onufriev, Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 9204. F. Pietra, Chem. Biodiversity 2012, 9, 606. W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph. 1996, 14, 33. E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng, T. E. Ferrin, J. Comp. Chem. 2004, 25, 1605. J. M. Word, S. C. Lovell, J. S. Richardson, D. C. Richardson, J. Mol. Biol. 1999, 285, 1735.
512
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
[27] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, ©. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, D. J. Fox, Gaussian, Inc., Wallingford CT, 2009. Received August 24, 2014
503
Restoring Histone Deacetylase Activity by Waste Product Release. A View from Molecular Mechanics Simulations with Mammalian HDAC8 by Francesco Pietra Accademia Lucchese di Scienze, Lettere e Arti, Classe di Scienze, Palazzo Ducale, Lucca I-55100 (phone/fax: þ 39-0583-417336; e-mail: [email protected])
HDAC8 is a ZnII-based, single-peptide mammalian histone deacetylase that is localized mainly in the cytoskeleton of smooth muscle cells, thus regulating muscle contractility. HDACs are also widely involved in cellular processes, ranging from cell differentiation to proliferation, senescence, and apoptosis; in particular, protecting a telomerase activator from ubiquitin-mediated degradation. How HDACs can eliminate the hydrolytic reaction products, in order that the process of deacetylation of the acetyllysine moiety of histones can take place again, has long been debated in the scientific literature, without reaching any firm conclusion, however. This question is the subject of the present work, carried out along a theoretical line that is capable of describing the whole pathway followed by the acetate product (ACT). A model was built here on the crystal data for the Y306F-mutated HDAC8 complex with a diacetylated peptide of the p53-tumor-suppressor class. That was followed by manually hydrolyzing the acetylated moiety bound to ZnII and discharging the monoacetylated peptide product (MAP). The latter was replaced by a H2O molecule bound to ZnII, while ACT was left free in the reaction cage. This ZnII cluster was DFT-parameterized for the ff99SB force field without any further bias. As the result of random-acceleration molecular dynamics (RAMD) simulations, egress of ACT from the reaction cage toward the aqueous environment can follow three pathways. Two of them utilize the channel for peptide (or histone) uptake and are preferred, if ACT leaves the reaction center before MAP (or the deacetylated histone). The third pathway, developing along the internal channel, is available to ACT even if MAP is still in place.
1. Introduction 1). – 1.1. HDAC Isozymes: Overview. Histone deacetylases (HDACs) are a group of isoenzymes that, along other functions, regulate transcription. They do so by catalyzing deacetylation at the e-N-acetyllysine branch of histones, i.e., by reversing the process of acetylation catalyzed by histone acetyltransferase [1]. As these processes are modulated during malignant transformations of cells, HDACs are widely involved in cellular processes, ranging from cell differentiation to proliferation, senescence, and apoptosis. Various diseases can grow on abnormal HDAC activity, against which the most commonly pursued remedies are to inhibit the enzyme. Therefore, HDAC inhibitors have become most chased targets in cancer research [2], although several HDAC knockouts are known to be lethal during mammalian 1)
Abbreviations: aa, amino acid(s); ACT, acetate product; CPU, central processing unit; CUDA, computer unit device architecture; ff, force field; MD, molecular dynamics; MNA, monoacetylated peptide resulting from the hydrolytic reaction; NPT, number-of-particles pressure temperature (isothermal¢isobaric ensemble); NVT, number-of-particles volume-temperature (canonical ensemble); PDB, protein data bank; QM, quantum mechanics; RAMD, random-acceleration MD; RMSD, root-mean-square deviation; ts, time step; VdW, Van der Waals. Õ 2015 Verlag Helvetica Chimica Acta AG, Zîrich
504
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
embryonic development [3]. Also non-histone HDAC8 substrates have been frequently taken into consideration [1]. HDAC Isoenzymes range in size from 300 to 1000 amino acids, in various classes with different sequences. Despite the differences, the catalytic center of this enzyme family is conserved, pivoting on a divalent metal ion. In vitro, the metal specificity increases in the order NiII < ZnII < FeII < CoII [4], while, in vivo, the metal can only be ZnII or FeII, as far as it is presently known [5]. Studies of mammalian HDACs have been concentrated on HDAC8, as this is a single polypeptide easily obtained in active form [5]. It is localized mainly in the cytoskeleton of smooth muscle cells. The most powerful HDAC8 inhibitor discovered to date is largazole, a cyclodepsipeptide isolated from a cyanobacterium, Symploca sp., which grows on coral reefs of Florida Keys [6]. Largazole acts as a prodrug, becoming bioactive on hydrolysis of the thioester side chain, which unmasks a thiolate function. It is the latter that, by binding to the divalent metal ion, inhibits HDAC8 [7]. The crystal structure of various complexes between HDAC8 and hydroxamate or other thiolate inhibitors have also been elucidated [7], while complexes with acetylated protein substrates in crystalline form have long been elusive. The breakthrough was achieved with a catalytically inactive human Y306F mutant HDAC8 enzyme. This allowed to trap a diacetylated peptide, (Ac)-l-Arg-l-His-l-Lys(e-Ac)-l-Lys(e-Ac), containing a fluorogenic coumarin group at its carboxy terminus, patterned on the p53 tumor suppressor protein, as a model for histones [8]. Present work started from this complex, as illustrated below. 1.2. How HDAC8 Bioactivity Can Be Regenerated Following Substrate Deacetylation. Deacetylation at the e-N-acetyllysine branch of histones was proposed to be carried out by a ZnII-coordinated H2O nucleophile under general acid-base catalysis with either one or two histidine side chains [1]. In order that this process, i.e., acetylation of a fresh histone by histone acetylase followed by deacetylation by histone deacetylase, repeats, the latter enzyme has to be freed from the reaction products at every cycle. How ACT can be eliminated from histone deacetylase was investigated here along MD simulations, starting from the above complex between inactivated HDAC8 and p53-derived diacetylated ligand of PDB code 2 V5W [8]. This represents one of many possible models [8], while ACT is the product of any histone deacetylase. Moreover, while the deacetylated histone, due to its large size, can only retreat along the entrance channel, the small ACT molecule could find alternative exit routes. This is why the fate of the ACT product of the deacetylase reaction was given primary attention here. This work, completed in the year 2012, was waiting for compilation, while new QM/ MM studies have appeared as to the mechanism of hydroxamate inhibition of these enzymes [9], based on the same complex used in our work [8]. Previous hypotheses, outlined above, were mainly consolidated. 2. Results and Discussion. – 2.1. The Model. The above mentioned complex between inactivated HDAC8 and p53-derived diacetylated ligand of PDB ID 2 V5W [8] is constituted of the protein 15 helices and 8 sheets, the catalytic ZnII center, and the diacetylated peptide substrate, as illustrated in Fig. 1. The active center comprises ZnII, coordinated with the D178 and D267 carboxylates, as well the H180 imidazole moiety, the catalytic H2O molecule, and the C¼O O-atom of the lysine-branch e-N-Ac group of
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
505
Fig. 1. Stereoview of the complex between Y306F mutant HDAC8 (PDB ID 2 V5W) and a p53-derived diacetylated peptide, (Ac)-l-Arg-l-His-l-Lys(e-Ac)-l-Lys(e-Ac), containing a fluorogenic coumarin group at its carboxyl terminus [8]. The iron protein comprises 15 helices (helix number (in italics), residue numbers, color: 1, 21 – 29, dark gray; 2, 36 – 48, yellow; 3, 49 – 54, green; 4, 63 – 68, dim gray; 5, 72 – 84, black; 6, 90 – 97, forest green; 7, 107 – 128, hot pink; 8, 156 – 166, deep blue; 9, 182 – 189, sea green; 10, 219 – 223, cyan; 11, 236 – 256, cornflower blue; 12, 280 – 293, magenta; 13, 307 – 324, white; 14, 336 – 341, orange; 15, 358 – 375, purple), other than eight sheets, highlighted in black, and ZnII, shown as a red ball, coordinated with D178, H180 and D267, also highlighted in red. The diacetylated peptide is highlighted in cyan. To aid comparisons, here and in all other figures in this work, the orientation of the complex follows that in [10].
the diacetylated peptide [8]. Mutation of Y306 to F306 prevents stabilization of the incipient amine that would be generated during the hydrolytic reaction, with the result that all the above components are trapped into the reaction cage. A model of the presumed arrangement of this system following the hydrolytic reaction was built here by first manually replacing F306 with Y306, thus restoring the naturally occurring active protein. It was then arbitrarily assumed that the monoacetylated peptide product (MAT) leaves the catalytic center before the ACT. Thus, the peptide was removed and replaced by a H2O molecule, while ACT was left free in the reaction cage (Fig. 2, where residues considered to be important for the functionality of the protein [10] are also highlighted, as described in the figure legend). It was now time to create a ff for the ZnII cluster of HDAC8. So far, ZnII complexes with macromolecules have been treated according to either the nonbonded model, the cationic dummy-atom approach, or the bonded model. The nonbonded model places an integer charge on the metal ion, such as with HDAC complexes with inhibitory hydroxamic acids [11]. This may become computationally expensive, as it requires the inclusion of very-long-range electrostatic interactions, while being also prone to structural deformation on prolonged MD simulations. The cationic dummy approach [12] places dummy cations to mimic valence electrons around the metal ion and is integrated with an electrostatic model. This may present the advantage of avoiding a rigid conformation at the active site, but it has been scarcely validated. The bonded
506
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
Fig. 2. Stereoview of the HDAC8 model used in present work, obtained from the complex of Fig. 1 by replacing F306 with naturally occurring Y306. The diacetylated peptide was replaced by a H2O molecule bound to ZnII and a free molecule of ACT in the hydrolytic reaction cage. By comparison with Fig. 1, it can be recognized that F207 (cyan) lies along the channel for uptake of the diacetylated peptide substrate. Southeastern to F207, the Zn cluster is seen in color by element (C, gray; N, blue; O, red; Zn, dark gray). Also illustrated are R37 (green), G303 (hot pink), and G139 (red), which are considered to take part in the opening/closing of the internal channel [10]. Finally, G305 (orange) is also illustrated.
model has been recently automated for ZnII protein complexes in the framework of the ff99SB ff [13], however neglecting dihedral terms and Lennard¢Jones parameters, which were considered of minor importance. Either the experimental structure was used or the metal cluster was calculated via DFT, and an electrostatic model was added [14]. Present work makes use of automated generation of the ff developed in the framework of the bonded model [15]. As a basis, this takes the ff99SB ff [13], which is known to treat proteins adequately, while allowing minimization with the coordinates of all the protein atoms considered as free variables. Any bias, other than those intrinsic to the ff99SB ff, was avoided, so that dihedral terms and Lennard¢Jones parameters were included. Like other popular non-polarizable ffs, ff99SB represents the electrostatic field through fictitious partial point charges distributed along the various atoms. To this regard, RESP charges [13] were adopted here. This modeling could be straightforwardly applied to ACT at HF/6-31G* level, similarly to any other standard protein component. The strategy for the ZnII cluster pivoted on the additive properties of the ff99SB ff [13], so as to restrict as much as possible the area to submit to QM calculations. To this end, the amino acid ligands D178 and H180 were cut, changing their atom CB to a Me group. In the ff, the Me group was removed and each remain was connected to the protein through a HN¢CH(Me)¢C linker, parameterized at the same theory level as ACT. In contrast, due to the presence of a transition metal, the HF approach could not be used for the ZnII cluster.
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
507
Recourse was made to the DFT theory with b3lyp hybrid functional and 6-311G* basis set. RESP Charges were calculated accordingly. Non-bonded terms for ZnII were taken from [16]. 2.2. MD of ACT Disposal. Classical MD, carried out with NAMD engine [17] for the ff99SB ff [13], proved satisfactory for the stated purpose. Tight convergence minimization resulted in a slight flattening around ZnII and the three bound amino acids with respect to a tetrahedral arrangement (Fig. 2). That slightly distorted arrangement was conserved during MD equilibration (Fig. 2), and RMSD proved satisfactory. It is difficult to appreciate whether this slight conformational discrepancy adds to the many cases of problems encountered by current ffs in representing the experimental structure along minimization and MD [18]. In fact, the ZnII cluster modeled here only represents a possible situation during the hydrolytic reaction, not an experimental structure. The next decision to be taken was about how to follow the pathway of ACT disposal under dynamic conditions. Experimental and theoretical methods to follow moving ligands through proteins have been recently reviewed [19]. While no experimental methodology proved able to describe whole pathways, various theoretical approaches have been applied. Of these [19], RAMD [20] was chosen here. This is a form of accelerated MD, where a tiny external force is applied to the center of mass of the molecule to move. The direction of the external force is applied randomly, and the force is made to continue only if the ligand moves along a direction for a selected period of time. If not, a new direction is chosen, randomly again. A small integration timestep was used here, so that all chemical bonds, except the solvent H2O, could be left unrestrained. Also, a large number of runs with different initial direction of the external force were carried out to reach statistical significance. In defense of RAMD, with the system myoglobin and a diatomic gas [21], this theory allowed identification of all important pathways that had been previsiously discovered through computationally expensive, massive MD [22]. RAMD proved also useful for larger ligands with a variety of proteins, as recently reviewed [23]. The results of RAMD simulations for present HDAC8-ACT model system are compiled in the Table. An acceleration of 0.28 kcal/è · g proved to be a good compromise between applying a small force that does not distort VdW barriers, while gaining in MD speed. Generated in a VdW cubic cage, and bordered by protein components as described in the Table, ACT was observed to reach the external solvent along three major routes, with gates as represented in Fig. 3. Two of them, with ACT egressing through either F207 and F208 or D101 and F208, utilize the channel for peptide (or histone) uptake, and are preferred if ACT leaves the reaction center before MAP (or the deacetylated histone). Residence time into the ACT cage is very substantial in both cases, while an intermediate binding pocket (BP) was only observed along the route toward D101-F208, and with a short residence time. The route through F207 and F208 is described in Fig. 4, with ACT represented by its centroid as a red ball of radius 0.5 è, and the diacetylated peptide of Fig. 1 manually superimposed after completion of the MD movie. It is seen that ACT first attempts to egress eastward, then takes the empty entrance channel. Should this be occupied by MAP, the lysine branch would inhibit the passage. Therefore, this route of egress, and likewise that through D101 and F208, can only be taken by ACT if it leaves before MNA.
508
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
Table. Statistical Recapitulation of ACT Pathways inside HDAC8, and Egress toward the Solvent, under RAMD a ) Representative residues for gates and BPs b )
Gate/total number of trajectories
Pathway/number of trajectories per pathway
F 207-F 208/7
ACT_cage ! gate/4 gate: F 207, F 208 ACT_cage > SBP ! gate/3 SBP: I34, R37, F142, G303, G305, Y307; gate: F 207, F 208
Residence time [ps] of ACT into BPs and total trajectory length [ps] ACT_cage: 31-450; gate: 38-492 ACT_cage: 115-270; SBP: 10-13; gate: 296-344
D101-F 208/5 ACT_cage ! gate
gate: D101, F 207
ACT_cage: 2-172, gate: 5-240
K33-N307/5
gate: K33-N307
ACT_cage: 150-800; gate: 300-950 ACT_cage: 171-495; SBP: 5-10; gate: 374-700
ACT_cage ! gate/2
ACT_cage ! SBP ! gate/3 SBP: I34, R37, F142, G303, G305, Y307; gate: K33-N307 NO EGRESS/9
ACT_cage
ACT_cage : R37, G140, ACT_cage: 500-1000 H143, C153, Q263, M274, G303, Y306
a ) Based on X-ray diffraction data of HDAC8 in complex with with a diacetylated peptide, PDB ID 2 V5W [8]. Residue F 306 was replaced with Y306. The diacetylated peptide substrate was replaced with a H2O molecule bound to ZnII and a free molecule of acetate in the hydrolytic reaction cage. RAMD Acceleration of 0.28 kcal/è · g and ts ¼ 1fs. b ) At 5 è from ACT; HDAC8 numbering.
Fig. 3. Stereoview of gates used by ACT to egress toward the solvent from the hydrolytic center of the complex of Fig. 2, according to RAMD simulations. As compiled in the Table, the three major gates should be imagined through F207-F208, F208-D101, and K33-N307. The first two gates correspond to ACT egress from the diacetylated peptide channel of Fig. 1, while the latter gate fits for one of the exit channels proposed from the examination of static X-ray diffraction structures [10].
The third pathway, with ACT egress through K33 and N305 (Table), develops along the internal channel, thus being available to ACT even if MAP is still in place. A snapshot of the movie at completion of the trajectory (Fig. 5) is instructive. Thus, in line
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
509
Fig. 4. Representative egress pathway of ACT toward the solvent along the diacetylated peptide channel. It is seen that ACT, leaving the hydrolytic reaction cage, egresses through F207 (cornflower blue) and F208 (orange). However, by superimposing the diacetylated peptide of Fig. 1 (cyan), it is seen that the acetylamino branch at the reaction center would constitute a barrier to ACT escape along this route. Therefore, this egress route is only conceivable, if MNA resulting from the hydrolytic reaction leaves the channel before ACT. Note: ACT is represented by its centroid as a red ball of radius 0.5 è, except at the exit portal, where it is represented in full atomic shape. The time interval between consecutive centroids in this movie is 0.5 ps.
Fig. 5. Representative egress pathway of ACT toward the solvent along the internal channel. It is seen that ACT, leaving the hydrolytic cage, first tries unsuccessfully to egress along the diacetylated peptide channel of Fig. 1, unable, however, to overcome the barrier opposed by F207 and F208 (both white, not labeled). Then, ACT tries unsuccessfully to egress from the opposite side, unable to overcome the barrier opposed by W141 (cyan). Finally, ACT finds a way toward the solvent by taking the internal channel [10], egressing through F306 (hot pink) and R37 (blue). Note: as in Fig. 4, except 1 ps between consecutive centroids.
with data in the Table, it is seen that ACT (represented by the centroid, like in Fig. 5) spends a considerable time into the hydrolytic cage, then tries unsuccessfully to egress along the diacetylated peptide channel of Fig. 1, where, however, it proves unable to
510
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
overcome the barrier opposed by F207 and F208 (in white, not labeled). Also attempts at egressing from the opposite side proved unsuccessful, ACT being unable to overcome the barrier opposed by W141 (cyan). Then, RAMD captured the relatively rare event of ACT egressing through F306 (hot pink) and R37 (blue). This corresponds to an exit internal channel that was hypothesized from the examination of static models [10]. It remains unobstructed, even if MNA resulting from the hydrolytic reaction is still in place. 3. Conclusions. – This work investigated how histone deacetylase HCAC8 can get free from the acetate product in order that the hydrolytic process repeats. A model was built on X-ray diffraction data for a complex between inactivated HCAC8 and an acetylated peptide, parameterizing the active ZnII cluster, along a bonded model. Recourse was made to the DFT theory with b3lyp hybrid functional at 6-311G* basis level, without any further bias than those intrinsic to the ff99SB ff. The system, solvated with TIP3P H2O in a periodic box, proved stable on MD simulations. Pathways for expulsion of the acetate product proved amenable to study under RAMD. The outcome is that the acetate product leaves preferentially along the channel for uptake of the diacetylated peptide substrate, if that is free. Otherwise, the internal channel remains available for ACT expulsion. I warmly thank Prof. FranÅois-Yves Dupradeau for invaluable help with the parameterization of the ZnII cluster on his R.E.D. development server. Computational Details. – All MD and RAMD simulations were carried out with software NAMD [17], CUDA v. 2.8, running on Debian GNU/Linux amd64. System setup and analysis were carried out with both VMD [24], v. 1.9, and CHIMERA [25], v. 1.5.3 software, the latter also serving for all the illustrations in this work. The Protein Model. The HDAC8 model was built on the PDB ID 2 V5W, i.e., the X-ray diffraction derived Cartesian coordinates, at 2.0-è resolution, for the point-mutant Y306F human HDAC8¢substrate complex, expressed on Escherichia coli [8]. In physiological solns., the protein exists as a monomer. Therefore, only chain A, with its 80 molecules of crystallization H2O, and the acetylated peptide substrate, were retained. Also, mutation back F306Y was carried out manually, maintaining the F306 orientation. Lacking residues 1 – 14 and 378 – 388 were not regenerated. The H-atoms were added with software REDUCE [26], under the Ð-buildÏ flag. A H2O molecule was made to replace the Ac group of peptide substrate, conserving the position of ligand O-atom. Manually hydrolyzed acetate was displaced just enough to give room to the H2O molecule. All remaining parts of the peptide substrate were removed. Parameterization and MD. Parameterization of the Zn complex, ACT, and linker was carried out on R.E.D. development servers [15], with software GAUSSIAN [27] for the QM computations. The Zn complex comprised ZnII and its ligands D178, H180, and Y306, in addition to the above mentioned H2O molecule. The amino acid ligands were cut so as to have their atom CB becoming a Me group. In the ff generated on the R.E.D. server, the Me group was removed, so that a HN¢CH(Me)¢CO moiety was used to bind the ligands to the protein. HF QM Computations with basis set 6-31G* were carried out for both ACT and the HN¢CH(Me)¢C linker, used to join the Zn complex to the protein. For the Zn complex, recourse was made to the DFT theory with b3lyp hybrid functional and 6-311G* basis set. In all cases, at tight convergence for geometry optimization. Preparation of the topology/parameter files was carried out by merging the ff derived above with ff99SB ff for the protein 363 standard amino acids and its 80 molecules of crystallization H2O by means of software LEAP, as available in AMBERTOOLS [13], v. 1.5. The ensemble was solvated with 15,356 molecules of TIP3 H2O in a rectangular box of x ¼ 86, y ¼ 76, and z ¼ 86 è, and neutralized with 4 Cl ¢ ions. It should be noted that, because of LEAP features [13],
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
511
chain A begins with residue number 1. Therefore, 14 had to be added to the residue number read from simulations in order to refer to the crystal numbering, which was adopted throughout this report. The ensemble was minimized at constant volume and then equilibrated for 10 ns, at 1 atm (Langevin piston), by gradually raising the temp. to 300K (Langevin thermostat with damping, 1). The Zn¢ACT distance was maintained at 4.1 è by applying harmonic forces (colvars). Flags: 1 – 4 scaling, 1 è; cutoff, 9.0 è; switchdist, 6.0 è; and pairlistdist, 11.0 è. The same settings, without applied colvars, were also used for RAMD simulations. Integration for all these computations was at 1-fs timestep with no rigid bonds for any component as a safe side after that it was noticed that the use of rigid bonds altered the RAMD response with the myoglobin¢O2 system [21]. A RAMD acceleration of 0.28 kcal/è · g turned out, from simulations in the range 0.1 – 0.3 kcal/è · g, to be the minimum value to bring consistently ACT out of the solvent box. The direction of acceleration was set to change randomly, if ACT traveled by less than 0.02 è in the evaluated block.
REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]
N. A. Wolfson, C. A. Pitcairn, C. A. Fierke, Biopolymers 2012, 99, 112. M. Dokmanovic, C. Clarke, P. A. Marks, Mol. Cancer Res. 2007, 5, 981. M. Haberland, R. L. Montgomery, E. N. Olson, Nat. Rev. Genet. 2009, 10, 32. S. L. Gantt, S. G. Gattis, C. A. Fierke, Biochemistry 2006, 45, 6170. D. P. Dowling, S. L. Gantt, S. G. Gattis, C. A. Fierke, D. W. Christianson, Biochemistry 2008, 47, 13554. K. Taori, V. J. Paul, H. Luesch, J. Am. Chem. Soc. 2008, 130, 1806. K. E. Cole, D. P. Dowling, M. A. Boone, A. J. Phillips, D. W. Christianson, J. Am. Chem. Soc. 2011, 133, 12474. A. Vannini, C. Volpari, P. Gallinari, P. Jones, M. Mattu, A. Carfı, R. Defrancesco, C. Steinkîhler, S. Di Marco, EMBO Rep. 2007, 8, 879. K. Chen, X. Zhang, Y.-D. Wu, O. Wiest, J. Am. Chem. Soc. 2014, 136, 11636. S. Haider, C. G. Joseph, S. Neidle, C. A. Fierke, M. J. Fuchter, Biorg. Med. Chem. Lett. 2011, 21, 2129. D.-F. Wang, O. Wiest, P. Helquist, H.-Y. Lan-Hargest, N. L. Wiech, J. Med. Chem. 2004, 47, 3409. P. Oelschlaeger, R. D. Schmidt, J. Pleiss, Protein Eng. 2003, 16, 341. D. A. Case, T. E. Cheatham III, T. Darden, H. Gohlke, R. Luo, K. M. Merz Jr., A. Onufriev, C. Simmerling, B. Wang, R. Woods, J. Comp. Chem. 2005, 26, 1668. M. B. Peters, Y. Yang, B. Wang, L. Fîsti-Molnr, M. N. Weaver, K. M. Merz Jr., J. Chem. Theory Comp. 2010, 6, 2935. E. Vanquelef, S. Simon, G. Marquant, E. Garcia, G. Klimerak, J. C. Delepine, P. Cieplak, F.-Y. Dupradeau, Nucl. Acids Res. 2011, 39, W511. S. C. Hoops, K. W. Anderson, K. M. Merz Jr., J. Am. Chem. Soc. 1991, 113, 8262. J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, K. Schulten, J. Comput. Chem. 1998, 26, 1781. P. Koehl, M. Levitt, Nature Struct. Biol. 1999, 6, 108. F. Pietra, Chem. Biodiversity 2013, 10, 86. S. Lîdemann, V. Lounnas, R. C. Wade, J. Mol. Biol. 2000, 303, 797. F. Pietra, Chem. Biodiversity 2011, 8, 1410. J. Z. Ruscio, D. Kumar, M. Shukla, M. G. Prisant, T. M. Murali, A. V. Onufriev, Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 9204. F. Pietra, Chem. Biodiversity 2012, 9, 606. W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph. 1996, 14, 33. E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng, T. E. Ferrin, J. Comp. Chem. 2004, 25, 1605. J. M. Word, S. C. Lovell, J. S. Richardson, D. C. Richardson, J. Mol. Biol. 1999, 285, 1735.
512
CHEMISTRY & BIODIVERSITY – Vol. 12 (2015)
[27] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, ©. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, D. J. Fox, Gaussian, Inc., Wallingford CT, 2009. Received August 24, 2014
