Molecular Dynamics and Quantum Chemical Approach for the

Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Molecular Dynamics and Quantum Chemical Approach for the Estimation of an Intramolecular Hydrogen Bond Strength in Okadaic Acid Toru Matsui, Kanako Yamamoto, Takehiro Fujita, and Kenji Morihashi J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b03272 • Publication Date (Web): 29 Jun 2018 Downloaded from http://pubs.acs.org on June 29, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Dynamics and Quantum Chemical Approach for the Estimation of an Intramolecular Hydrogen Bond Strength in Okadaic Acid Toru Matsuia, *, Kanako Yamamotob, Takehiro Fujitaa, and Kenji Morihashia, b

a: Graduate School of Pure and Applied Sciences, University of Tsukuba, 1-1-1 Tennodai Tsukuba, Ibaraki, 305-8577 Japan. b: Master’s program in Education, University of Tsukuba, 1-1-1 Tennodai Tsukuba, Ibaraki, 305-8577 Japan.

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 33

ABSTRACT

We have evaluated the strength of intramolecular hydrogen bond in a protein based on molecular dynamics and quantum chemical calculation. To estimate the intramolecular hydrogen bond strength in okadaic acid (OA), we analyzed the influence of solvent and protonation states on the hydrogen bond and the entire structure. We performed molecular dynamics calculation and analyzed the strength of the hydrogen bond by measuring bond length and bond angle. The stable structure differs depending on the kind of solvent used and the protonation state of OA. Using the mean interaction energy from the quantum chemical calculation, hydrogen bond length and angle were investigated against bond energy. Although dielectric constant slightly depends on bond energy, the estimation of the intramolecular hydrogen bond strength in OA is possible even in a protein environment. The Coulomb interaction between OA and surrounding arginine produced a more negatively charged O1 in OA. The hydrogen bond energy in the deprotonated state is larger than that in the protonated state.


2

Page 3 of 33 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Hydrogen bond (HB) is a ubiquitous non-covalent bond that has been studied in various fields; e.g., molecular designing,1 biochemistry,2,

3

and drug-design.4 From a biochemical aspect,

inter/intramolecular hydrogen bond stabilizes the structure of DNA.5, 6 HB involves many kinds of chemical interactions such as Coulomb interaction, dipole-dipole, dipole-induced dipole, and van der Waals interaction. Despite the prominence of HB in biomolecules, its strength has not been defined clearly. Many theoretical approaches for the estimation of HB strength have been proposed to further understand the origin of stability in biomolecules.7-10 The International Union of Pure and Applied Chemistry defined HB in 2011 and emphasized the importance of computational chemistry.11, 12 Quantum mechanical (QM) calculation is an effective tool for estimating the HB strength. Hobza’s group has employed a very accurate computational method, CCSD(T), and approximated the complete basis set limit for the HB energy between nuclear acid bases.13,

14

Many researchers have also studied HB in protein.

Deshmukh and Gadre calculated the N-H…O=C intramolecular HB by fragmenting molecules1518

, followed by many researchers.19, 20 Takano’s group also investigated molecular interactions

within the secondary structures of protein.21 QM calculation excels in accurately computing energy, but is limited only to a time-frozen structure. Furthermore, normal DFT calculations succeed in reproducing experimental geometries based on the HB observed by NMR or X-ray structure analyses. Because the geometry of a protein varies with temperature, it is necessary to evaluate HB strength by considering dynamics. Conventionally, molecular dynamics (MD) calculation describes intermolecular interaction with Coulomb and van der Waals interactions. Therefore, MD calculation has frequently been used for evaluating the bond length and angle to understand the behavior of HB in protein.22 Our previous study on a nylon oligomer degrading


3


Page 4 of 33

enzyme, a kind of serine-protease, necessitated the interpretation of data with a statistical method to understand the role of the HB network in the catalytic triad.23 Combining the advantages of QM and MD will lead to the novel understanding of HB. Using the ensemble of conformation obtained by MD, we can compute the average HB energy by QM calculation. This method will give a new insight to the correlation between HB geometry and energy. For applications, we picked up the case of an okadaic acid (OA, C44H68O13, see Scheme 1 for the structure) in a protein environment. The intramolecular HB between O1 and O24 of OA has been taken much attention by many researchers.

Scheme 1: Chemical structure of OA. OA is isolated from black marine sponge Halichondria okadai.24 25 OA is a tumor-promoting compound, which inhibits protein phosphatase (PP) 1 and 2A.26-30 Owing to these properties, OA has been widely used as benchmark to study the role of many enzymes. A marine sponge accumulates OA for chemical defense. H. okadai prevents self-toxicity by lowering the concentration of OA with two types of enzymes: serine/threonine protease that contains Mn2O2 core, and the OA binding protein, categorized to “okadaic acid binding protein 2” (OABP2).31 Ehara et al. presented the X-ray for the new OABP2 protein, named “OABP2.1.” 32 Accordingly, OA is retained in the cavity between two domains which are the part of OABP2.1. The cavity captures the OA by recognizing the cyclic structure formed by the intramolecular hydrogen bond between the C1 carboxyl group and the C24 OH group. The bond between OABP2.1 and OA is


4

Page 5 of 33 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


further stabilized by the electrostatic interaction of OA with surrounding arginine residues (R44, R185). This intramolecular HB also contributes to the accommodation of OA in marine sponge. Surrounding molecules significantly influence HB. Therefore, it is necessary to investigate the effect of these molecules on HB in OA by varying the conditions. Although the C1 carboxyl group is expected to be deprotonated in neutral solution, it could also be protonated in very low pH. The protonation of C1 carboxyl group expectedly affects HB strength. Additionally, the solvent molecule influences intramolecular HB, Hence, solvent dependency is also included in the “surrounding molecule effect.” In this study, we analyzed the “surrounding molecule effect” on intramolecular HB in OA. The proton content of C1 carboxyl group was also investigated to understand the effect of protonation. This work aims to estimate the HB energy of intramolecular HB in OA in a protein environment, contributing to the significant understanding of the HB in a biomolecule.

Computational Details To evaluate the HB energy in protein environment, two models were considered in this study: (1) homogenous solvent model and (2) protein model. Following (1), the geometry on HB can be interpreted as the interaction energy in solvation. For geometry variety, it is necessary to consider the average interaction energy among the snapshots that satisfy certain geometrical conditions. Hence, a kind of “interaction energy map” was generated as a function of bond length and angle. The interaction energy was also analyzed against the dielectric constant (εvalue) because ε-value is unique for a homogeneous solution. After the MD simulation of (2), we obtained the geometry of intramolecular HB and the interaction energy based on the interaction energy map at a specific ε-value. The details are presented in the following subsection.


5


Page 6 of 33

Pre-calculation We investigated the dynamics associated with OA in four types of solvents: water (WAT) and methanol (MOH) as protic solvents, and N-methyl acetoamide (NMA) and chloroform (CHL) as aprotic solvents. For “the surrounding molecular effect” in a protein, we used the structure of OABP 2.1 from the protein data bank (pdb, PDBID: 4WRI,32 by X-ray crystallization). All methyl groups of trimethyl lysine in the pdb file were replaced by normal lysine to reproduce the actual protein environment. To generate the force field of OA, we adopted the structural information of the deprotonated OA in the pdb file and optimized the whole structure at HF/6-31G(d) level. We considered two patterns for the protonated state of OA as shown in Figure 1. The former, “Protonated I” is more stable at equal computation level. For the latter, “Protonated II,” another intramolecular HB was observed between the carboxyl O1 and O2. Moreover, we computed the electrostatic potential (ESP) of these optimized geometries at HF/6-31G(d) level and determined the ESP charge using ANTECHAMBER program. The optimized geometry and ESP charge of each molecule are presented in the Supporting Information. Molecular dynamics simulation for OA and surroundings All MD simulations were performed with PMEMD module in AMBER14 program.33 A solvation box is generated by tLEaP program implemented in AmberTools15. The detail for the solvation box and the number of atoms is present in the Supporting Information. To generate topology and input card files, we adopted the ff12SB34 and selected a general AMBER force field (GAFF). 35 For the simulation in OABP2.1, we added ca. 6650 water molecules to construct the solvation box (TIP3PBOX)36 around the protein. The system was neutralized by sodium ions (Na+). To obtain a structure at thermal equilibrium, geometry optimization was performed until the absolute gradient of the system was reduced to 0.0001 kcal/mol/Å. We then performed 600 ps MD simulations for all of


6

Page 7 of 33 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the system with gradual heating from 0 to 300 K by 50 K (in 100 ps) intervals. The time step was set to 1 fs throughout the simulation. The temperature and pressure were controlled by a weakcoupling algorithm and maintained at 300 K and 1 atm, respectively. After a 60-ns production run, we sampled 50000 data points from 10-60-ns simulations to analyze MD trajectories. The program CPPTRAJ37 was used for the trajectory analyses. We utilized a visualized molecular dynamics (VMD)38 software with graphical user interface and analyzed the bond length, r, and bond angle, θ. Since there are two patterns for hydrogen bonds in protonated states, we denoted the parameters as rHB and θHB as shown in Figure 1. The bond length and angle of OA satisfy both criteria (i.e., (1) rHB is smaller than 3.40 Å and (2) θHB is larger than 100°), conforming with the intramolecular HB from the reference.39 Data points were also counted as an indicator of the HB strength in each solution or protonation state. Evaluating interaction energy Here, HB energy was considered as the interaction energy between two fragments, as it often is. To translate the information on HB geometry into HB energy, snapshots of the molecular dynamics were considered. In this study, we evaluated the HB energy by creating a model for the intramolecular HB in OA. For the simplicity, the smallest model (Model 1) was considered. In modeling the molecule, a CC bond is capped with hydrogen atom, i.e., a C atom was replaced with H atom so that the length of C-H bond is 1.09 Å. Consequently, we focused on the HB between C24 hydroxyl and C1 carboxyl groups, labeled as Fragments A and B, respectively (see Scheme 2). For the validity of the model, we randomly generated 100 input files for Model 2 which consists of 122 atoms in three states.


7


Page 8 of 33

Scheme 2: Complete structure of OA and the molecular model for the estimation of interaction energy. Among the trajectory data of 50000, we selected 1000 patterns of conformation for each production run in any solvent. Therefore, we have 4000 patterns for the conformation. In the polarizable continuum model (PCM),40, 41 the Boys-Bernardi’s counterpoise scheme was used to correct basis set superposition error (BSSE).42 These calculations were performed at MP2/631+G(d) level. HB energy is expressed by: EHB = {E(A) + E(B)} – E(AB) + EBSSE.

(1)

where E(AB) is the total energy of both fragments, E(A or B) is the energy for either fragment A or B, and EBSSE is computed using ghost atoms on each fragment. Since this scheme with PCM is applicable in Gaussian, we can estimate the HB energy at different ε-value. We computed four patterns of calculation with varying ε-values of 184.0 (NMA), 78.4 (WAT), 32.6 (MOH), and 4.7 (CHL).


8

Page 9 of 33 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


By using these computational schemes, the average HB energy was estimated for assigned geometries. All QM calculations were performed with GAUSSIAN 09.43

Results Since the targeted HB does not always form, we considered two schemes for calculating HB energy depending on the situation. (A) For HB for all conformation, we simply calculate their average (labeled as “ALL”). (B) For HB limited to the snapshots that satisfy the criteria for intramolecular HB (labeled as “HBF,” the criteria are shown in the computational details section). Geometry on HB Tables 1 and 2 present the average and median of bond length and bond angle for the HB among 50000 conformations (listed in “ALL”) and among snapshots where intramolecular HB forms (listed in “HBF”), respectively. In deprotonated state, HB forms in more than 30000 patterns (60% of snapshots). Furthermore, in a protic solvent, a solvent molecule can form HB with O1 in place of O24-H, known as “HB alternation.” Consequently, rHB surpassed the threshold value for HB. However, this broken intramolecular HB occasionally recovers, which can affect the sampling. Since we have checked the reproductivity from another sampling computation, this ON-OFF cycle of HB can be observed in protic solution. In CHL solution, several chloride atoms surround the OA, which can bind to either O24-H or O1-H. These weaken the HB and reduce θHB. Especially in the case of CHL in the deprotonated state, a strong Coulomb interaction between the carboxyl group and sodium atom (a counter cation) stabilizes the whole system because of low ε-value. Therefore, the variation of geometrical parameters in CHL, especially in “HBF,” is greater than that in the other solvents.


9


Page 10 of 33

Analysis of protonated states (“Protonated I” and “Protonated II”) shows that protonation weakens intramolecular HB in terms of the average bond length for all HB. In some cases, the computed rHB surpassed the criteria for intramolecular HB. Furthermore, HB-alternation is more probable in protonated states than in deprotonated state. Since O24 atom forms an intramolecular hydrogen bond with other atoms in OA, HB is weak in Protonated II state. In order to visualize the frequency of geometrical parameter of HB, Figure 2 shows a twodimensional histogram with grids of rHB and θHB at 0.05 Å and 5°, respectively. This histogram also aims to understand the stable structure in the equilibrium state for each solvation. Darker (close to black) areas indicate more snapshots that exhibited HB. For example, when focused on CHL in the deprotonated state, the distribution differs from the other cases because of the Coulomb interaction between OA and a counter ion. This effect can be considered as “surrounding molecule effect”. The broadening of the plot, e.g., plots of CHL in all protonated states, suggests variety of the conformation. Since we sampled only 50000 snapshots, the concentration of the frequency is at most 1500-2000 per grid. Although this amount might be insufficient, the results reveal a few behaviors for the average geometry around the HB. For all solvents in the deprotonated state, the plot is concentrated on the average value for “HBF.” Analyses of HB energy First, it is necessary to discuss the validity of using small molecular model. Figure 3 shows the scattering plot for the interaction energy in two models. The detailed numbers are shown in Table S4 in the Supporting Information. There is a linear correlation between two models. The interaction energy in larger model is overestimated compared with the smaller model, because the larger model includes many hydrogen bonds (e.g. C-H…O) other than the targeted HB. The effect of contamination from the other interactions can be estimated to the intercept of the linear fitting, ca. 11-17 kJ/mol in any states. We also applied the molecular


10

Page 11 of 33 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tailoring approach (MTA) proposed by Deshmukh et al.15-18 in the Supporting Information. The results show that the interaction energy in Model 1 is underestimated in 2-10 kJ/mol compared with MTA scheme. Although further researches are necessary to estimate more accurate HB energy, we only focus on the simple model for the purpose of simplifying the HB energy change. Figure 4 shows (1) the 2D-distribution map of the conformations for the calculations by QM and (2) the interaction energy map on a grid (similar to Figure 2), classified by color. Here, the dielectric constant was set to 4.7 (assumed in the CHL solvation). Another case (assumed in MOH solvation), with dielectric constant at 32.6, is shown in the supporting information. The result of “Hydrogen bond” in Table 3 (low “#” values) shows that the average geometry among 1000 conformations agrees with the original data (8000 conformations). Notably, solvent type had no significant effect on the average HB energy because the maximum variation of HB energy is 5 kJ/mol. Although the color of the grid is slightly sparse in this study, this figure will improve when the sampling number increases. The maxima of average HB energy are 35.6 kJ/mol, 28.8 kJ/mol, and 15.6 kJ/mol in deprotonated, Protonated I, and Protonated II states, respectively. The maximum HB energy for Protonated II (15.6 kJ/mol) is categorized as “weak HB” (which is

Molecular Dynamics and Quantum Chemical Approach for the

Recommend Documents