Computational Study of the Degradation of S-Adenosyl Methionine in

Subscriber access provided by GAZI UNIV

Article

Computational Study of the Degradation of S-Adenosyl Methionine in Water Timm Lankau, Tzu Nung Kuo, and Chin-Hui Yu J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b09639 • Publication Date (Web): 22 Dec 2016 Downloaded from http://pubs.acs.org on December 25, 2016

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry A 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 30

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

Computational Study of the Degradation of S-Adenosyl Methionine in Water Timm Lankau, Tzu Nung Kuo, and Chin Hui Yu∗ Department of Chemistry, National Tsing Hua University, Hsinchu, Taiwan E-mail: [email protected] Phone: +886 (0)3 5162080. Fax: +886 (0)3 5721534

∗

To whom correspondence should be addressed

1 ACS Paragon Plus Environment


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

Abstract The degradation of S -adenosyl methionine (SAM) to homoserine-γ-lactone (HSL) and methyltioadenine (MTA) in water is studied with MD simulations. The AM1 hamiltonian is used for the quantum part and the flexible AMBER force field for the H2 O molecules. The MD simulations predict the free energy barrier for the degradation reaction to be between 109 and 112 kJ mol−1 and an overall gain in free energy of −26 kJ mol−1 . The high barrier and the low energy gain of this reaction can be linked to the interactions among the carboxylate group of the SAM molecule and solvent H2 O molecules which are not observed on the product side. Hence, the H2 O molecules effectively slow down the reaction which otherwise would be much faster.

Introduction S -Adenosyl methionine (variously abbreviated as SAM, SAMe, SAM-e or Adomet) was first reported 1 as active methionine in 1952. Its biochemical usage goes far beyond 2 its most prominent role as a methyl donor 3 which led to the claim of SAM being one of the few molecules required for live. 4 Being such an important molecule, it is hardly surprising that SAM is discussed as a possible cure for many medical conditions such as chronic liver diseases 5 or depressions 6 and marked as a food supplement. 7,8 A wide spread usage of SAM in drug formulations is in great parts hampered by its long known chemical instability. 9–14 Scheme 1 summarises the key steps of the degradation process in water. 9 In neutral and acidic solutions, an oxygen atom from the carboxylate group attacks the γ-carbon atom of the amino acid yielding methylthioadenosine (MTA) and homoserine lactone (HSL, α-amino-γ-butyrolactone) as primary degradation products. The lactone HSL reacts with ambient water in a second step to form the open-chain amino acid homoserine (HS). The hydrogen atoms bound to the carbon atoms adjacent to the sulfonium cation are highly polarised. In alkaline solutions, 10,13 an hydroxide ion attacks one of these hydrogen 2 ACS Paragon Plus Environment

Page 2 of 30

Page 3 of 30


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 ACS Paragon Plus Environment



Page 4 of 30

Page 5 of 30

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


calculations. Truhlar’s M05 25 functional with the 6-311+G(2d,p) basis set was used for the amino acid and the ribose C5 methylene group (high level region), based on a previous evaluation 26 highlighting the superiority of this method for analysis of sulfonium ions in solutions. The low level region (OLYP/6-31G(d,p)) is defined by the remainder of the adenosine nucleoside. The OLYP functional was chosen due to a comprehensive comparison 27 showing its benefits in the analysis of SN 2 reactions.

Figure 1: ONIOM model for the reference calculations, transition state in PCM water. White: H, cyan: C, blue: N, red: O, yellow: S. Table 1: Energies of the degradation process in PCM water obtained from standard quantum calculations. method/PCMa

AM1/H2 O ONIOM/H2 O AM1/none

∆SAM E ∆SAM G ∆PCM E ∆SAM E ∆SAM G ∆SAM E ∆SAM G

SAM 0.0 0.0 −209.739 0.0 0.0 0.0 0.0 TS 74.088 70.240 −163.641 78.086 73.790 27.990 24.892 HSL + MTA −79.415 −128.734 −76.161 −59.197 −110.539 −212.993 −254.085 a Quantum chemical method and the chosen PCM. TS: Transition state. ∆SAM : Difference relative to SAM. ∆PCM E: PCM Stabilisation energy of the solvated species. All energies in kJ mol−1 . Table 1 summarises the energy changes for the degradation reaction calculated with Gaussian 09. 28 The nature of the stationary points was confirmed with frequency calculations. All calculated frequencies for minima were real and transition states featured one and only one imaginary frequency aligned with the reaction path. The obtained harmonic frequencies were used to calculate the changes in Gibbs free energy. The energy changes for the formation of the transition state (Figure 1) agree remarkably well, while the AM1 29 calculations overestimate the overall energy changes associated with 5 ACS Paragon Plus Environment


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 6 of 30

the degradation reaction. Still, the AM1 hamiltonian offers good description of the hypervalent sulfonium ion at very low computational costs. Hence, the AM1 hamiltonian was used for electronic structure calculations with MOPAC 2016 30 for the MD simulations. The QMMM keyword 31,32 extends the one-electron diagonal elements of AM1 Hamiltonian as used in MOPAC by the interactions ϕi of an electron with the Coulomb charges of the surrounding MM water molecules ˜ νν = hνν − ϕi = hνν − h

MM 2 X e qj rij

(1)

where the orbital ν is centered on the i-th quantum atom and j stretches over all MM atoms of the solvent. The electrostatic potential ϕi incorporates the polarisation of the solute by the solvent into the quantum calculations. A weighing factor was introduced into the calculation of ϕi as used for the MOPAC calculations to avoid un-physically large polarisation terms for small values of rij . ϕ˜i =

MM X

1/2

+ 1/2 tanh

rij − ai bi

e 2 qj rij

(2)

Where ai is typically the Van der Waals radius of the i-th quantum atom, and bi is set for the weight to be 0.01 for rij being 0.95 ai . Hence, solvent atoms within the Van der Waals radius of the quantum atom do not contribute significantly to the value of ϕ˜i . The AMBER 33 force field was used for 1011 H2 O molecules of the solvent. The AMBER force field for water is essentially a flexible version 34 of the TIP3P 35 force field for bulk water. The potential energy of a single water molecule EH2 O arises from the thermal distortion of the molecule itself, and the interactions with its peers and the sum over all water molecules gives the potential energy Esolv of the MM water mantle surrounding the QM solute

Esolv =

NW X i

N

W Coul 1X Eij + EijVdW [E(ωi ) + E(ri1 ) + E(ri2 )] + 2 i6=j

1 E¯H2 O = Esolv NW

(3) (4)

where NW is the number of MM water molecules, E(ωi ) and E(ri ) the harmonic distortion 6 ACS Paragon Plus Environment

Page 7 of 30

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


energies of the HOH bonding angle (ω) and the two OH bonds (ri1 , ri2 ), EijCoul the Coulomb interaction energy between two H2 O molecules and EijVdW the van der Waals interaction energy between the oxygen atoms. E¯H2 O is the average potential energy of a MM water molecule. The interface between the QM and MM region encloses three types of interactions. They are the Coulomb interactions between the Mulliken charges on the QM solute and the fixed charges on MM solvent atoms, the van der Waals interactions between solute and solvent atoms, and a restraint acting on hydrogen bridges. The first two terms define Elink , the interaction energy between the quantum solute and the MM solvent

Elink =

Coul Elink

+

VdW Elink

QM MM M X X q qj e2 Aij Bij i = + 12 − 6 rij rij rij i j

(5)

where qiM is the Mulliken charge 36 of the quantum atom, qj the static charge of the MM atom, and the pair Aij and Bij the van der Waals parameters of the atom pair obtained form the combination rules of the force field used. Equations 1 and 5 contain both the Coulomb interactions between the quantum and MM atoms. To avoid double counting, the quantum energy 37 EAM1 obtained from the MOPAC calculation needs to be corrected for Coul Elink

Ecore =

   E

AM1

ref − EAM1

mechanical embedding (6)

  ref Coul EAM1 − EAM1 − Elink

electronic embedding

ref with Ecore being the corrected AM1 energy and EAM1 the AM1 energy of the optimised

quantum solute in the vacuum as reference. The total potential energy of the system is dominated by Ecore , Elink and Esolv . Changes in Ecore are dominated by changes in the electronic structure of the reactants, the distortion of the quantum solute by thermal motion and its polarisation by the surrounding MM molecules. It contains further the difference between the weighted Coulomb interaction obtained from the quantum calculation and its approximation using Mulliken charges, which 7 ACS Paragon Plus Environment


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 8 of 30

was found to be negligible as the MM water molecules did not penetrate the quantum core. In a system without any QM solvent molecules, Elink can be used to monitor directly the interaction of the organic molecules with the surrounding H2 O molecules and how the analysed reaction changes these interactions. While changes in E¯H2 O are small as the reaction proceeds, changes in Esolv can be big due to the large number of H2 O molecules monitored. These changes in Esolv reflect effectively how the solvent reacts to changes in the solute as reaction proceeds. The van der Waals parameter and the associated combination rules for the calculation of Elink were taken either from the AMBER 33 or from the OPLS-AA 38–40 force fields in its 2008 version (Sections S.1.1 and S.1.2 in the supplement). Both force fields contain hydrogen atoms without van der Waals parameter. Quantum protons without van der Waals parameter can get very close to water oxygen atoms (OW), which leads to un-physical forces in the MD calculations and ultimately to convergence problems with the QM calculations (Figure S.1). A restraint is therefore implemented to the QM/MM interface to avoid unphysical H+ movements. The degree of proton transfer between an electronegative QM atom (Q) and OW is measured with the collective variable g

g = |~r1 | − |~r2 |

(7)

where ~r1 is the vector from Q to H+ and ~r2 that from H+ to OW. The restraint is defined by a composite potential VR (g)

VR (g) =

   0

for g < g0

  1/2 k (g − g0 )2

for g ≥ g0

(8)

which makes it increasingly difficult for the H+ to move beyond the threshold value g0 . Figure 2 shows the number of H+ positions beyond a fixed value for g0 for NH4 Cl in 8 ACS Paragon Plus Environment

Page 9 of 30




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 10 of 30

decreases as the value of g0 increases until it hits its minimum at g0 = −0.4. This decrease can be rationalised by the observation that too small values of g0 keep the H+ too close to QM atom and thereby allow only weak bonds between the H+ and OW. As g0 becomes too large, the disruptive nature of the artificial H+ proton transfer becomes visible and Elink increases again. Considering the results presented in Figures 2 and 3, the MD simulations were done with a value of −0.4 for g0 and a weak force constant k of 50 kcal mol−1 Å−2 . Simulations with selected QM water molecules in the solvation shell carry the risk that these water molecules diffuse out of the shell and are replaced by MM molecules. This problem can be avoided by replacing all MM water molecules by QM molecules, but this would also increase computational costs exorbitantly. Instead, a constraint similar to that for the movement of protons (Equation 8) was used to control movement of QM water molecules

VW (rXO ) =

   0

X for rXO < rXO

  X 2 1/2 kW (rXO − rXO )

for rXO ≥

(9)

X rXO

where rXO is the distance between a hydrogen bonding centre X on the SAM molecule and the oxygen atom of a QM water molecule. Two ways are possible to pair X and O. Either, a fixed X···O pair is assigned at the beginning of the simulation and does not change in course of the simulation. Or, the X···O pairs are assigned on the fly so that each O is paired with the X closest to it. The first scenario ensures that each X always maintains a fixed number of QM water molecules in its solvation shell even if the solute breaks apart, while the QM water molecules are allowed to move freely in the restraints of the solvation shell made X X from overlapping spheres with the radius rXO . The radius rXO of the solvation shell around

the atom X was determined from the un-normalised radial pair the distribution functions 41 X gXO (r) of the simulations with MM water molecules. rXO was assigned to be halfway between

the first maximum and the first minimum (Figure 4) so that a QM water molecule can move around its most likely position in the shell, but gets gently pushed back into the shell (kW =


Page 11 of 30

25 kJ mol−1 ) when it is about to leave the shell. The strongly coordinating oxygen atoms of X the carboxylate group have the smallest value of 3.0 Å for rXO followed by the hydroxy groups

of the sugar (3.9 Å). The weakly coordinating amino groups of methionine and arginine have the largest sphere radius of 4.3 Å. 4.2 CO2 NH2 OH

3.5 2.8 2.1 0.7 0

2

3.90 4.30

1.4 3.00

1000 * g(r)

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


3

4

5

6 7 r [Å]

8

9

10

Figure 4: The un-normalised radial pair the distribution functions g(r) from a simulation X . with MM water molecules to find the radii of the atomic solvation shells rXO MD simulations of an NpT ensemble were done with a house code using a leapfrog integrator 42 (∆t = 1 fs), an Andersen thermostat 43 (probability of a collision with the heat bath 1:75, 298.15 K) and a Berendsen barostat 44 (coupling constant 100 fs, 1 bar). A damping and shifting technique 45 was used for Coulomb and van der Waals forces, and the pressure 46 was calculated on the basis of pair interactions between molecules with the whole quantum region being treated as a single molecule. The centre of the minimum sphere enclosing all quantum atoms was kept with a harmonic constraint (25 kcal mol−1 Å−2 ) close to the centre of the simulation cell (periodic boundary conditions, minimum image convention). The path of the degradation reaction is defined as

λ = λ2 − λ1

(10)

(λ1 , λ2 : OC and CS bond lengths as marked in Figure 1) so that the reactant SAM with its



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 12 of 30

short CS bond is identified by values of λ smaller than zero. A harmonic bias potential

VB (λ) = 1/2 K (λ − λ0 )2

(11)

with K = 400 kcal mol−1 Å−2 is used to drive the reaction. The MD simulations for the calculation of the Helmholtz free energy changes were done in several steps. The transition state of a SN 2 reaction is expected to be linear. Hence, the first set of MD simulations was augmented by a weak restraint aiding the simulations to home in on an approximately straight OCS angle. Each of these fast simulations used a larger value for λ0 than its ¯ − λ0 was predecessor and the series continued until a change in sign of converged value of λ observed. The last two simulations bracket the transition state geometry and the one with ¯ − λ0 closest to zero was selected to seed the next step. In the second step, the starting λ point λS for the free energy calculations was obtained from a 600 ps simulation without the angular restraint, in which the temperature T , all energies, the pressure p, the volume V , and the bias (λ1 , λ2 , λ, VB ) converged. Two simulations branched out from the converged starting point in the third step. Depending on the branch followed, λ0 either increased or decreased by 0.1 Å every 80 ps in each simulation to follow the reaction path either to the reactants or to the products. Every 80 ps prior to the change in λ0 , a snapshot (geometry, velocities, forces) of the simulation was dumped. The resulting dump files were used to seed the final simulations for each window i with λi0 and longer simulation times up to 400 ps in the fourth step. The computational load of the last step was spread over a maximum of 16 computers, because the calculations for the windows are independent from each other. ¯ i of λi and its standard deviation σi for the umbrella Each window i provided the average λ integration. 47–49 The unbiased mean force (∂Ai/∂λ) is calculated directly from the distribution Pi (λ) of λ values in the i-th window ¯i ∂Ai −1 ∂ ln Pi (λ) ∂VB (λ) 1 λ−λ ¯i) = − = − K (λ − λ ∂λ RT ∂λ ∂λ R T σi2


(12)

Page 13 of 30

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


with Ai being the unbiassed free Helmholtz energy. The individual windows are then combined into a single expression for the whole integration path windows X ∂Ai ∂A = pi (λ) ∂λ ∂λ i

(13)

where the weight pi (λ) is obtained from the normal distribution Pi (λ) of λ values in each window i during the MD simulation Pi (λ) pi (λ) = P j Pj (λ)

and

¯ i )2 1 −(λ − λ Pi (λ) = √ exp 2 σi2 σi 2 π

(14)

The final expression for the mean force (Equation 13) is a continuous function by its construction and essentially independent from the position of the windows. Its numerical integration by the trapezoidal rule begins and ends the first and the last λ0 value of the simulation with three equidistant support points between any two values of λi0 .

Results and discussion Initial (N = 2048, p = 0.9975 bar, T = 298.19 K) simulations (equilibration: 700 ps, production 1.05 ns) of pure water with the Amber force field yielded a density ρH2 O of 0.987 g cm−3 with a standard deviation of 0.005 g cm−3 which underestimates the experimental value of 0.997 g cm−3 by 0.0093 g cm−3 or 0.9%. The reported 34,35 ρH2 O values for TIP3P water of 0.982 and 0.986 g cm−3 are slightly smaller than the one reported here, but still covered by the standard deviation. The oxygen-oxygen pair distribution has its first maximum at 2.75 Å and its integration up to the first minimum at 3.43 Å indicates that on average 4.33 H2 O form the first solvation shell. The radius of the first solvation shell is slightly smaller than that of 2.77 Å reported for an NVE ensemble, 50 but the number of molecules in the shell agrees well with those reported of the original TIP3P force field (4.35). 34 This good agreement between the AMBER and the TIP3P force field is not surprising since both use 13 ACS Paragon Plus Environment


the same parametrisation for long-range interactions. 140 ΔRAo [kJ/mol]

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 14 of 30

0.20 0.25 0.30 0.15

120 100 80 60 40 20 0 -1.5

-1

-0.5

0 0.5 λ [Å]

1

1.5

2

Figure 5: Comparing different starting points λS for the umbrella integration. To get an estimate of the error in ∆R A in the analysis of the degradation reaction, steps 2 and 3 were repeated four times with different starting values λS ranging between 0.15 and 0.30 Å using the OPLS-AA parameter set for Elink . The initial optimisation (step 1) suggested λS being 0.2 Å to be the best starting point for the umbrella integration. Figure 5 compiles the results from the free energy calculations using the trajectories from step 3 (∆λi0 = 0.1 Å, 41 windows, 80 ps per window). The ∆R A value at λ equal to 2.25 Å was chosen as a common reference point [∆R A (λ = 2.25 Å) = 0] to simplify the comparison. The average position of the maximum is 0.23 Å with the standard deviation of the sample being 0.01 Å suggesting an error of 4% in the position of the transition state. While the branches of the ∆R A (λ) curves overlap very well on the product side (λ > 0), the curve for λS = 0.2 Å does not fit into the crowd on the reactant side. The amino acid part can either be an open chain or coiled up. The coiled up SAM molecule unwinds in a MD simulation without the bias potential; therefore, the dump files from the simulations with λS = 0.2 Å was chosen to seed the final window calculations (step 4). The average of the four curves in Figure 5 suggests a barrier height of 104 kJ mol−1 with a standard deviation of the sample of 5.58 kJ mol−1 and an overall gain in free energy of -33.2 kJ mol−1 with a standard deviation of 4.40 kJ mol−1 . The branches on the product side with much less conformational diversity overlap much better. The predicted barrier for 14 ACS Paragon Plus Environment

Page 15 of 30

the back reaction is 137 kJ mol−1 with a standard deviation of the sample of 1.47 kJ mol−1 , which suggests an error about 1% for the umbrella integration. The final window calculations run much longer (up to 400 ps) than the 80 ps for the generation of the seeds. Hence, it is possible to assume an upper limit for the error in the free energies of roughly 5 kJ mol−1 . 120 ΔRAo [kJ/mol]

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


OPLS-AA AMBER

90 60 30 0 -30

-2 -1.5 -1 -0.5 0 0.5 1 λ [Å]

1.5 2 2.5

Figure 6: Effect of the force field parameter in the link region. Two simulations were done to analyse the influence of the MM parameter sets for the calculation of Elink (Equation 5). λi0 ran from −1.95 to 2.55 Å (∆λi0 = 0.1 Å, 46 windows) for the simulation using the AMBER parameter set and from −1.7 to 2.3 Å (∆λi0 = 0.1 Å, 41 windows) in the simulation with the OPLS-AA force field parameter set. The averages of the bias quadrupel (λ1 , λ2 , λ, VB ) converged reasonably well during the simulations to generate seeds for the window calculation. The analysis of the two branches from the third step showed ¯ i converges faster than other properties of the system. To estimate the influence of that λ the simulation time for a single window on the barrier of the degradation reaction, the simulations times for the fourth step were increased threefold (AMBER parameter, 240 ps) ¯ i converges quickly (Figure S.2 in and fivefold (OPLS-AA parameter, 400 ps). The value of λ the supplement) and the curves for ∆R A (λ) obtained from the umbrella integration of the seeding calculations overlap very well with those from the window calculations (Figures S.3 and S.4) which suggests that the window calculations are sufficiently long for the calculation of free energy changes. Figure 6 shows the changes in free energy as the reaction proceeds. Both curves for 15 ACS Paragon Plus Environment


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

∆R A (λ) are essentially the same despite a shallow minimum in the curve obtained with the AMBER paramentrisation. The minimum is caused by the brief (40 ps) formation of a 1,2-oxathiane ring (Figure 7) in which a carboxylate oxygen atom approaches the sulfonium cation. The shortest S−O distance observed was 2.28 Å. Despite the rapid decay of the oxathiane, the sulfonium cation proved to be a strong attractor for the carboxylate group. The comparison of both ∆R A (λ) curves suggests further that the OPLS-AA parameter set solvates the carboxylate group stronger than the AMBER parameter set so that the oxathiane cannot be observed in calculations with the OPLS-AA parameter set.

Figure 7: Intermediate 1,2-oxathiane ring with a tetravalent sulphur atom. White: H, cyan: C, blue: N, red: O, yellow: S. The ∆R A (λ) curves have their maxima at 108.754 and 112.445 kJ mol−1 which is about 1.5 times higher than the barriers obtained from the standard quantum calculations (Table 1) and the overall gain in free energy is −26.351 and −26.443 kJ mol−1 respectively. The average barrier height of 111 kJ mol−1 for the degradation reaction can be compared with reported barrier heights for the enzymic methyl transfer 15,16 of roughly 88 kJ mol−1 depending on enzyme, substrate and methylation step. It compares much better to the reported barrier of 118 kJ mol−1 for the similar SN 2 reaction of the transfer of a methyl group from (CH3 )3 S+ to Cl– in pure water. 26 Furthermore, reported rate constants 11,13 in the range between 1.80 · 10−5 and 6.0 · 10−6 s−1 for the first order degradation reaction of SAM in water according to Scheme 1 at pH 7.5 and 37◦ C suggest an average barrier height of 106 kJ mol−1 which matches the calculated barrier heights of 109 and 112 kJ mol−1 reasonably well. The solvent molecules stabilise the dipolar SAM molecule much stronger than the neutral products HSL and MTA. The results from the simulation with OPLS-AA parameter set were 16 ACS Paragon Plus Environment

Page 16 of 30

Page 17 of 30

chosen for a more detailed analysis of the degradation reaction, because no oxathiane was observed during the simulation and the longer simulation times ensure a better convergence of the properties of a single window. Figure 8 shows the average values of Elink (Equation 5) and Ecore (Equation 6) for the majority of steps during the MD simulation with the OPLSAA parametrisation. The average of the first 10 (lead in) and the last 13 points (lead out) of each curve were used as reference for the calculation of energy differences. The quantum energy Ecore decreases by −281 kJ mol−1 as the CO–2 ···S+ zwitterion gets annihilated and the uncharged lactone is formed. This stabilisation is counterweighted by the loss of solvation energy as reflected in the increase of Elink . The carboxylate group forms strong hydrogen bonds to the surrounding H2 O molecules, and these hydrogen bonds result in low values for Elink . As the reaction proceeds the carboxylate group is lost and so is the stabilisation by the surrounding H2 O molecules. The changes of Elink (522 kJ mol−1 ) are approximately twice as large as those in Ecore so that the sum increases by 241 kJ mol−1 during the degradation process. 400 200 E [kJ/mol]

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


Ecore

0 -200 -400

Elink

-600 -800

-1000 -1.5

-1

-0.5

0 0.5 λ [Å]

1

1.5

2

Figure 8: The value of Ecore and Elink during the OPLS-AA simulations and their interpolative splines. The grey lines (length not scale) mark the reference points for the calculation of energy differences. Figure 9 shows Esum being the sum of Ecore and Elink as a function of the reaction path λ. The analysis of the smoothed spline linking the data points suggests that Esum increases by 241 kJ mol−1 , which would be a highly endothermic reaction as the energy gain from the formation of the lactone cannot compensate for the loss of favourable interactions 17 ACS Paragon Plus Environment



Page 18 of 30

Page 19 of 30

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


overall gain in Esolv of 194 kJ mol−1 reduces the energy costs of the degradation reaction to 47 kJ mol−1 . Assuming similar pV terms on the reactant and product side, the estimated change in enthalpy is +47 kJ mol−1 while that in free energy is −26 kJ mol−1 which suggest a sizeable increase in entropy as the reaction proceeds. A possible source for this entropy increase in addition to the cleavage of the C−S bond are the H2 O molecules formerly tied to the CO–2 group, which can move more freely when the reaction is completed. Table 2: Position of the transition state on the reaction coordinate. (a)

AM1/none AM1/H2O (a) Ecore (b) Esum (b) Etot (b) ∆R A , OPLS (a)

standard quantum calculation (Table 1),

(b)

λ −0.29 0.16 −0.60 0.54 0.45 0.22

maximum of the fitted spline. λ in Å.

The influence of the solvent on transition state is also visible in its position λ (Equation 10) on the reaction coordinate (Table 2). The addition of a PCM to the standard AM1 calculations moves the early transition state (λ = −0.29 Å) in the vacuum to a late one (λ = 0.16 Å) in water. A similar shift can be observed with Ecore and Esum . Ecore has a maximum of 27 kJ mol−1 relative to the lead in at −0.60 Å. The stabilising interactions between solute and solvent move the transition state up to 0.54 Å as observed with Esum . The reorganisation of the solvent lessens the shift as Etot and ∆R A have their maxima at 0.45 and 0.22 Å respectively. The good agreement between the ∆R A (λ) curves (Figures S.3 and S.4) obtained from ¯ i converges so quickly that the the seed and the final window calculations suggests that λ fast scans with 80 ps per window can be used to probe the setup for possible error sources. Hence, a fast (−3.2 Å ≤ λi0 ≤ 2.6 Å, ∆λi0 = 0.1 Å, 59 windows, 80 ps per window) repetition of the AMBER simulation was done without the polarisation of the quantum atoms to understand the discrepancy between PCM calculations and MD simulations better. The 19 ACS Paragon Plus Environment


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

electronic embedding of the quantum solute in the MM solvent can be turned off by setting ϕi equal to zero (Equation 1). The quantum atoms and the solvent molecules only interact mechanically via the Coulomb and Van der Waals forces as described in the calculation of Elink (Equation 5). This simulation suggests that the 1,2-oxathiane is the dominating species in solution. A small barrier of approximately 15 kJ mol−1 separates the oxathiane from the degradation products HSL and MTA (Figure S.5). The overall gain in free energy is −231 kJ mol−1 relative to the oxathiane and −288 kJ mol−1 relative to the open chain form of SAM, which is similar to the values obtained from standard quantum calculations for the vacuum (Table 1). The polarisation term pushes the simulation from artificial (1,2-oxathiane) to physical (open chain methionine) while stabilising the dipolar SAM molecule. MD simulations with the TIP3P force field for bulk water and its flexible derivatives tend to yield values for the permittivity larger than 100, 34 which would favour the dipolar reactant. A quick (−2.7 Å ≤ λi0 ≤ 3.2 Å, ∆λi0 = 0.1 Å, 60 windows, 80 ps per window) repetition of the simulation with the OPLS-AA parameter for solute-solvent interactions using the SPC/Fw 34 instead of the AMBER force field for water reproduced the ∆R A (λ) curve perfectly (Figure S.6 in the supplement). Hence, the observed stabilisation of the ionic SAM molecule relative to degradation products can be regarded as independent of the force field for bulk water. The charged carboxylate group forms strong hydrogen bonds with the water molecules surrounding it. However, the MM water molecules can not describe the charge transfer typically associated with these hydrogen bonds. To understand the role of charge transfer between solute and solvent, 18 water molecules closest to hydrogen bonding sites were turned into AM1 quantum molecules. Figure 11 shows the the geometry of the SAM molecule and its solvation shell at the starting point (λS = 0.2 Å) for the free energy calculations. The restraints (green spheres) preventing the quantum water molecules from diffusing into the bulk solvent are centred around the adenine NH2 group, the sugar OH groups, the carboxylate O atoms, and the methionine NH2 group. The quantum water molecules can move free


Page 20 of 30

Page 21 of 30

between overlapping spheres of the restraint. The majority of them form a hydrogen bonded network stretching from the carboxylate group to the sugar hydroxy groups.

Figure 11: SAM with its solvation shell at the starting point for the free energy calculation with 18 quantum H2 O molecules. Balls: white: H, cyan: C, blue: N, red: O, yellow: S. Liquorice: red+white: QM H2 O, orange: MM H2 O. The green spheres mark the restraint for the quantum H2 O molecules. 140 ΔRAo [kJ/mol]

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


no H2O, NH2 18 H2O, NH2 + 23 H2O, NH3

120 100 80 60 40 20 0 -1.5 -1 -0.5

0

0.5 1 λ [Å]

1.5

2

2.5

Figure 12: ∆R A (λ) obtained from OPLS-AA simulations with different amounts of quantum H2 O molecules: No quantum water, 18 quantum H2 O, 23 quantum H2 O and a protonated methionine NH2 group with one Cl– counter ion. ∆R A (2.25 Å) = 0. The simulations with QM water molecules converge slower than those with MM water (Figures S.7 and S.8). Hence, the final window simulations (step 4, −2.0 Å ≤ λi0 ≤ 2.5 Å, ∆λi0 = 0.1 Å, 46 windows) ran for 400 ps a piece. Figure 12 shows a comparison of the ∆R A (λ) curves obtained without (Figure 6) and with quantum water. A common zero point was set at λ = 2.25 Å to simplify the comparison. The branches leading to the 21 ACS Paragon Plus Environment


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

products overlap very well, while those to the reactant differ. The quantum water molecules stabilise the zwitterion stronger than MM water molecules which leads to an increase in activation energy (130 kJ mol−1 ) and to a decrease in the overall free energy gain of the reaction (−9.62 kJ mol−1 ).

Figure 13: Quantum part of the starting point for the simulations with a protonated amino group (SAM+ ) and a Cl– counter ion. White: H, cyan: C, blue: N, red: O, yellow: S, green: Cl. The influence of the protonation of the methionine amino group was evaluated with the last setup, in which a proton was added to the amino group and the resulting positive charge compensated with a chloride anion (Figure 13). The restraint on the adenine amino group for the movement of the quantum H2 O molecules was replaced with one on the Cl– anion. The Cl– ion can move freely, and the quantum water molecules were not allowed to move N Cl O between spheres (rNO = 3.0 Å, rClO = rOO = 3.5 Å) to prevent the accidental replacement of

quantum water molecules in the solvation shell of the chloride ion by MM water molecules. Totally, 23 quantum water molecules were used for the final window simulation (−1.6 Å ≤ λi0 ≤ 2.3 Å, ∆λi0 = 0.1 Å, 40 windows, 400 ps per window). Figure 12 shows the nearly perfect overlap of the ∆R A (λ) curves obtained with and without a proton at amino group, which suggests that the ammonium ion has only a very small influence on the formation of the lactone. The difference between the two curves for small values λ is caused by the formation of ion pairs. The chloride ion binds either to the ammonium or to the sulfonium cation. The formation of ion pairs destabilises the reactant 22 ACS Paragon Plus Environment

Page 22 of 30

Page 23 of 30

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


side and the barrier for degradation reaction decreases to 122 kJ mol−1 and overall gain in free energy increases to −15.81 kJ mol−1 . However, the initial simulations to generate the starting point (Figure 13) and the study 26 of the hydrolysis of (CH3 )3 SCl suggest that the zwitterion dissociates so that the curve for non-protonated amino group appears to be more reliable.

Conclusions MD simulations with the AMBER force field for water and the AM1 method for the SAM solute predict the barrier for the degradation of SAM in neutral solutions to be between 109 and 112 kJ mol−1 and the overall gain in free energy to be −26 kJ mol−1 . The high barrier and the low energy gain from the degradation of the high energy SAM molecule can be explained with the loss of the methionine carboxylate group during the reaction. The associated loss of favourable solute-solvent interactions destabilises the transition state and products relative to reactants. An increase in hydrogen bonding among the solvent molecules and more importantly a large positive entropy change turn the degradation reaction into an exergonic process. Ironically, the water molecules which facilitate the degradation process also slow it down.

Acknowledgement This work was supported by the Ministry of Science and Technology (MOST) of Taiwan under Grant numbers MOST 104-2113-M-007-014, MOST 105-2113-M-007-002, and Prof. I. C. Tsai via grant MOST 105-2811-M-007-082.

Supporting Information Available The electronic supplementary material provides the geometry and van der Waals parameters of SAM for the AMBER and the OPLS-AA model. It also contains figures of a rogue 23 ACS Paragon Plus Environment


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 24 of 30

hydrogen atom crashing into an oxygen atom of the H2 O model and plots showing the convergence of λ and the influence of simulation times on the ∆SAM A (λ) curves.

This

material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Cantoni, G. L. The Nature of the Active Methyl Donor Formed Enzymatically from L-Methionine and Adenosinetriphosphate. J. Am. Chem. Soc. 1952, 74, 2942–2943. (2) Roje, S. S -Adenosyl-l-methionine: Beyond the universal methyl group donor. Phytochemistry 2006, 67, 1686–1698. (3) Schubert, H. L.; Blumenthal, R. M.; Cheng, X. Many paths to methyltransfer: a chronicle of convergence. Trends Biochem. Sci. 2003, 28, 329–335. (4) Markham, G. D. Encyclopedia of Life Sciences; John Wiley & Sons, Ltd: Chichester, UK, 2010. (5) Anstee, Q. M.; Day, C. P. S-adenosylmethionine (SAMe) therapy in liver disease: A review of current evidence and clinical utility. J. Hepatol. 2012, 57, 1097–1109. (6) Mischoulon, D.; Fava, M. Role of S-adenosyl-L-methionine in the treatment of depression: a review of the evidence. Am. J. Clin. Nutr. 2002, 76, 1158S–1161S. (7) Fetrow, C. W.; Avila, J. R. Efficacy of the dietary supplement S-adenosyl-L-methionine. Ann. Pharmacother. 2001, 35, 1414–1425. (8) Sarris, J.; Schoendorfer, N.; Kavanagh, D. J. Major depressive disorder and nutritional medicine: a review of monotherapies and adjuvant treatments. Nutr. Rev. 2009, 67, 125–131. (9) Parks, L. W.; Schlenk, F. The stability and hydrolysis of S-adenosyl-methionine; isolation of S-ribosylmethionine. J. Biol. Chem. 1958, 230, 295–305. 24 ACS Paragon Plus Environment

Page 25 of 30

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


(10) Borchardt, R. T. Mechanism of alkaline hydrolysis of S-adenosyl-L-methionine and related sulfonium nucleosides. J. Am. Chem. Soc. 1979, 101, 458–463. (11) Wu, S. E.; Huskey, W. P.; Borchardt, R. T.; Schowen, R. L. Chiral instability at sulfur of S-adenosylmethionine. Biochemistry 1983, 22, 2828–2832. (12) Hoffman, J. L. Chromatographic analysis of the chiral and covalent instability of Sadenosyl-L-methionine. Biochemistry 1986, 25, 4444–4449. (13) Matos, J. R.; Wong, C.-H. S-Adenosylmethionine: Stability and stabilization. Bioorg. Chem. 1987, 15, 71–80. (14) Desiderio, C.; Cavallaro, R. A.; De Rossi, A.; D’Anselmi, F.; Fuso, A.; Scarpa, S. Evaluation of chemical and diastereoisomeric stability of S-adenosylmethionine in aqueous solution by capillary electrophoresis. J. Pharm. Biomed. Anal. 2005, 38, 449–456. (15) Zhang, X.; Bruice, T. C. Catalytic Mechanism and Product Specificity of Rubisco Large Subunit Methyltransferase: QM/MM and MD Investigations. Biochemistry 2007, 46, 5505–5514. (16) Yue, Y.; Chu, Y.; Guo, H. Computational Study of Symmetric Methylation on Histone Arginine Catalyzed by Protein Arginine Methyltransferase PRMT5 through QM/MM MD and Free Energy Simulations. Molecules 2015, 20, 10032–10046. (17) Markham, G. D.; Norrby, P.-O.; Bock, C. W. S-Adenosylmethionine Conformations in Solution and in Protein Complexes: Conformational Influences of the Sulfonium Group. Biochemistry 2002, 41, 7636–7646. (18) Saez, D. A.; Vöhringer-Martinez, E. A consistent S-Adenosylmethionine force field improved by dynamic Hirshfeld-I atomic charges for biomolecular simulation. J. Comput. Aided Mol. Des. 2015, 29, 951–961.



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

(19) Farooqui, J.; Kim, S.; Paik, W. K. Measurement of isoelectric point of S-adenosylL-methionine and its metabolic products by an isoelectric focusing technique. Electrophoresis 1983, 4, 261–265. (20) Barbato, G.; Calabria, R.; Carteni-Farina, M.; D’Auria, G.; De Rosa, M.; Sartorio, R.; Warzburger, S.; Zappia, V. A physico-chemical approach to the study of the binding interaction between S-adenosyl-l-methionine and polyanions: binding constants and nature of the interaction with sodium poly (styrene sulfonate). Biochim. Biophys. Acta 1989, 991, 324–329. (21) Dapprich, S.; Komáromi, I.; Byun, K.; Morokuma, K.; Frisch, M. J. A new ONIOM implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives. J. Mol. Struct.: THEOCHEM 1999, 461-462, 1–21. (22) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F. et al. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678–5796. (23) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999–3094. (24) Scalmani, G.; Frisch, M. J. Continuous surface charge polarizable continuum models of solvation. I. General formalism. J. Chem. Phys. 2010, 132, 114110. (25) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Exchange-correlation functional with broad accuracy for metallic and nonmetallic compounds, kinetics, and noncovalent interactions. J. Chem. Phys. 2005, 123, 161103. (26) Lankau, T.; Yu, C.-H. Solvent effects on the intramolecular conversion of trimethylsulfonium chloride to dimethyl sulfide and methyl chloride. Phys. Chem. Chem. Phys. 2014, 16, 26658–26671. 26 ACS Paragon Plus Environment

Page 26 of 30

Page 27 of 30

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


(27) Swart, M.; Solà, M.; Bickelhaupt, F. M. Energy landscapes of nucleophilic substitution reactions: A comparison of density functional theory and coupled cluster methods. J. Comput. Chem. 2007, 28, 1551–1560. (28) M. J. Frisch et al., Gaussian 09. Gaussian, Inc.: Wallingford CT, 2010. (29) Dewar, M. J.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. Development and use of quantum mechanical molecular models. 76. AM1: a new general purpose quantum mechanical molecular model. J. Am. Chem. Soc. 1985, 107, 3902–3909. (30) Stewart, J. J. P. MOPAC 2016. Stewart Computational Chemistry, 2016. (31) Plotnikov, N. V.; Warshel, A. Exploring, Refining, and Validating the Paradynamics QM/MM Sampling. J. Phys. Chem. B 2012, 116, 10342–10356. (32) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227–249. (33) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179–5197. (34) Wu, Y.; Tepper, H. L.; Voth, G. A. Flexible simple point-charge water model with improved liquid-state properties. J. Chem. Phys. 2006, 124, 024503. (35) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (36) Mulliken, R. S. Electronic Population Analysis on LCAO-MO Molecular Wave Functions. I. J. Chem. Phys. 1955, 23, 1833–1840. 27 ACS Paragon Plus Environment


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

R (37) MOPAC . http://openmopac.net/manual/, last visit: Nov. 2016.

(38) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225–11236. (39) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. L. OPLS all-atom force field for carbohydrates. J. Comput. Chem. 1997, 18, 1955–1970. (40) Jorgensen, W. L. OPLS All-Atom Parameters for Organic Molecules, Ions, Peptides & Nucleic Acids, July 2008. Yale University, 2009. (41) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33–38. (42) Rapaport, D. C. The Art of Molecular Dynamics Simulations, 2nd ed.; Cambridge University Press: Cambridge, 2004. (43) Frenkel, D.; Smit, B. Understanding Molecular Simulation, 2nd ed.; Academic Press: London, 2002. (44) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684. (45) Fennell, C. J.; Gezelter, J. D. Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics. J. Chem. Phys. 2006, 124, 234104. (46) Louwerse, M. J.; Baerends, E. J. Calculation of pressure in case of periodic boundary conditions. Chem. Phys. Lett. 2006, 421, 138–141. (47) Kästner, J.; Thiel, W. Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: “Umbrella integration”. J. Chem. Phys. 2005, 123, 144104. 28 ACS Paragon Plus Environment

Page 28 of 30

Page 29 of 30

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


(48) Kästner, J.; Thiel, W. Analysis of the statistical error in umbrella sampling simulations by umbrella integration. J. Chem. Phys. 2006, 124, 234106. (49) Kästner, J. Umbrella sampling. WIREs: Comput. Mol. Sci. 2011, 1, 932–942. (50) Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105, 9954–9960.



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

Graphical TOC Entry


Page 30 of 30

Computational Study of the Degradation of S-Adenosyl Methionine in

Recommend Documents