Toward an Enhanced Sampling Molecular Dynamics Method for

Subscriber access provided by The University Library | University of Sheffield

Article

Toward an Enhanced Sampling Molecular Dynamics Method for Studying Ligand-Induced Conformational Changes in Proteins Ole Juul Andersen, Julie Grouleff, Perri Needham, Ross C Walker, and Frank Jensen J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.5b07816 • Publication Date (Web): 20 Oct 2015 Downloaded from http://pubs.acs.org on October 24, 2015

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 B 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

Toward an Enhanced Sampling Molecular Dynamics Method for Studying Ligand-Induced Conformational Changes in Proteins Ole Juul Andersen,†,§ Julie Grouleff,† Perri Needham‡, Ross C. Walker‡,¶ and Frank Jensen*,† †

Department of Chemistry, Aarhus University, Denmark

§

Center for Insoluble Protein Structures (inSPIN) and Interdisciplinary Nanoscience Center

(iNANO), Aarhus University, Denmark ‡

San Diego Supercomputer Center, 10100 Hopkins Drive, La Jolla, California, 92093-0505,

USA

¶

Department of Chemistry and Biochemistry, UC San Diego, 9500 Gilman Drive, La Jolla,

California, 92093, USA

*Corresponding Author: Department of Chemistry, Aarhus University, Langelandsgade 140, 8000 Aarhus C, Denmark, Tel.: +45 8715 5908, Fax: +45 8619 6199, E-mail: [email protected]

ACS Paragon Plus Environment

1


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 30

ABSTRACT

Current enhanced sampling molecular dynamics methods for studying large conformational changes in proteins suffer from certain limitations. These include, among others, the need for user defined collective variables, the prerequisite of both start and endpoint structures of the conformational change, and the need for a priori knowledge of the amount by which to boost specific parts of the potential. In this paper, a framework is proposed for a molecular dynamics method for studying ligand-induced conformational changes, in which the non-bonded interactions between the ligand and the protein are used to calculate a biasing force. The method requires only a single input structure, and does not entail the use of collective variables. We provide a proof-of-concept for accelerating conformational changes in three simple test molecules, as well as promising results for two proteins known to undergo domain closure upon ligand binding. For the ribose-binding protein, backbone root-mean-square deviations as low as 0.75 Å compared to the crystal structure of the closed conformation are obtained within 50 ns simulations, whereas no domain closures are observed in unbiased simulations. A skewed closed structure is obtained for the glutamine-binding protein at high bias values, indicating that specific protein-ligand interactions might suppress important protein-protein interactions.

INTRODUCTION


2

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


Protein function is often coupled to large conformational changes, e.g., closure of the protein structure around a substrate, channel opening or rotation of entire subunits as exemplified by, for example, F0F1-ATP synthase.1 Detailed descriptions of these transitions are difficult to obtain from experiments, as measurements are usually performed on an ensemble of structures. Due to the low probability and short lifetimes of high-energy conformers, these experiments are limited to providing information related to highly populated states.1,2 The development of singlemolecule techniques, such as single-molecule FRET3 and magnetic tweezers4, has made it possible to study the transitions directly; however, the information is usually restricted to properties such as interatomic distances and domain rotation, and thus does not provide a fully atomistic description of the transition pathway. Furthermore, the methods typically require some sort of labeling, often involving mutations in the peptide sequence, which can potentially have a significant impact on the results. An alternative approach is the use of molecular dynamics (MD) simulations, which provide a means of obtaining dynamical information with atomistic resolution. A serious drawback, however, is the computational resources needed to calculate the trajectories. Large domain movements typically take place on the microsecond to millisecond time scale, or potentially even longer.1 Until recently it was infeasible to reach these time scales with traditional brute force MD simulations. The millisecond landmark was reached in 2009 with a 1 ms simulation of the 58-residue protein bovine pancreatic trypsin inhibitor (BPTI), performed on the highly specialized supercomputer Anton,5,6 but routine access to millisecond timescales in traditional MD simulations is far from reality for the general research community. While increases in computational power over time are bringing such timescales closer it will still likely be many years, if ever, before such timescales can be routinely accessed in unbiased MD simulations. As such there is a pressing need to develop alternative methods that can provide for


3


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

an enhancement in sampling. Existing enhanced sampling methods can be broadly divided into the following groups: those that need only a single structure as input, such as normal mode analysis7 (NMA), principal component analysis8 (PCA), steered MD9 (SMD), and accelerated MD10,11 (aMD); those that need two or more structures as input, such as elastic network models2,12-15 (ENMs), targeted MD16,17 (TMD), functional mode analysis18 (FMA), and umbrella sampling19,20 (US); and those that can use both single and multiple structures as input, such as metadynamics21-23. Each of these methods have their own benefits and drawbacks. Many of the methods that require both start and endpoint structures, such as ENMs, TMD and US, are well established and widely used methods; however, the requirement of a known endpoint structure limits their use to cases where such structures are available. NMA and PCA need only a single structure to obtain the low-frequency and large-scale motions, respectively, but the directions obtained from these methods are numerous, and determining which ones are related to the protein function is not a straightforward task. For methods such as SMD, metadynamics and FMA, preconceived knowledge such as a steering direction has to be applied, proper collective variables have to be selected, or quantities directly related to the protein function such as the openness of a channel must be provided. In contrast, aMD simulations require only a single structure as input, and no additional information, beyond selecting a relevant boost potential, has to be applied. Recently, a 500 ns aMD simulation of BPTI was shown to provide a comparable degree of sampling of the conformational space as an unbiased 1 ms simulation performed on Anton.24 Although the method provides greatly enhanced sampling, it has the disadvantage of sampling states that may not be directly related to the function of the protein. In the following sections, a framework for an enhanced sampling MD method, denoted force difference bias (FDB), for studying ligand-induced conformational changes is presented.


4

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


The only requirement for the method is access to the apo structure of the protein. In contrast to aMD simulations it focuses the search on conformational changes that are directly induced by the ligand. A proof-of-concept for its ability to accelerate conformational changes is reported for three simple test molecules and the method is tested for two proteins known to undergo domain closure upon binding of a ligand, namely ribose-binding protein (RBP) and glutamine-binding protein (GlnBP).

THEORY AND METHODS Method Framework. The FDB method is based on the hypothesis that the apo and holo configurations of a protein experience a difference in force induced by the ligand that can be used to direct the apo conformation toward the holo conformation. Alternatively stated, if the interactions between a protein and a ligand induce a conformational change, it should be possible to accelerate the conformational change by enhancing the protein-ligand interactions. Obtaining a bias force directly from the difference in the apo and holo potential energy surfaces (PESs) alleviates the user from selecting functionally relevant collective variables or steering directions. The bias force can be considered as a collective variable that is defined implicitly by the system and automatically updated during the simulation. The algorithm of the FDB method is presented in Figure 1, where the apo and holo configurations are termed the reference and target systems, respectively. In brief, the difference in force between the target and reference systems is added to the target system immediately before the update of the accelerations, velocities, and positions. The new coordinates are then copied to the reference system before the next iteration begins. In order to control the amount of bias applied to the system, the force difference between the two systems is multiplied by a coefficient α, termed the bias parameter. The effect of the bias


5


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

potential is illustrated in Figure 2A. By increasing the value of the bias parameter, the energy barrier related to the open-to-closed transition is effectively reduced, while the opposite is true for the reverse transition. Although the target and reference PESs in Figure 2A are greatly simplified compared to the more rugged PESs of proteins, they serve as an illustrative example. Figure 2A also carries the assumption that the reaction coordinates of the reference and target systems are completely identical. This is most likely not the case. It is therefore expected that the bias force contains contributions orthogonal to the reaction coordinate of the target system (Figure 2B). Large bias values would thus result in increasingly off-path transitions and should therefore be avoided. As the force difference originates from the non-bonded interactions, a more efficient implementation of the FDB method from a computational perspective is obtained by splitting the force difference calculation into calculations for the individual vdW and electrostatic pairwise interactions between protein and ligand atoms, rather than including all the force field terms in the calculation (Equation 1). In order to minimize the effects of the orthogonal bias contributions on the conformation of the ligand, bias is only applied to the protein atoms.

∗ = + ∙ ∆ , + ∆ , , ϵ protein Eq. 1

Since there is no ligand present in the reference system, the non-bonded interactions between the ligand and the protein can be set to zero for this system, and the FDB method therefore corresponds to enhancing the effect of the ligand by adding the forces from the non-bonded interactions in the target system, multiplied by the bias parameter, to a standard force calculation (Equation 2). The bias force is calculated for all protein-ligand atom pairs within the complex in order to include any difference present beyond the cutoff distance used in the regular force


6

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


calculation. The electrostatic interactions are calculated in real-space without use of the particle mesh Ewald (PME) method25 and therefore without adjustment by erfc. This is computationally feasible due to the low number of atoms usually present in ligands.

∗ = + ∙ )*+,-.* , + *+,-.* , / , ϵ protein Eq. 2

Figure 1. The initial and final conformations of the system are illustrated in the left part of the figure, where the protein is shown in blue and the ligand in green. The initial target system consists of the apo conformation including a docked ligand. The atoms in the reference system have the same coordinates as in the target system, but no ligand is present (i.e., all non-bonded interactions between the ligand and the protein are set to zero). The algorithm behind the FDB method is shown to the right. Using the force difference as a bias, the system is driven into the ligand-induced conformation.


7


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

Figure 2. (A) Increasing the amount of bias applied to the target system results in a reduction of the energy barrier separating the initial conformation from the final conformation of the protein. (B) As the reference force does not point in the exact same direction as the target force, the bias force will contain orthogonal contributions with respect to the target force. Large values of the bias parameter will lead to a noticeable off-path contribution to the perturbed target force, f*.

Small Molecule Test Systems. As a proof-of-concept, the FDB method was applied to three simple test systems: the dissociation of a positively charged ion from an 18-crown-6 crown ether; the rotation around a single bond in 1,1-diphenylethane-1,2-diamine; and the chair flip of a dihydroxylated hexahydropyridazine. The perturbation in the latter two systems involves a protonation corresponding to the addition of an extra covalent bond. Therefore, for these systems, Equations 1 and 2 are invalid, as they consider conformational changes induced by only non-bonded interactions. For these two systems, the calculation of the bias was augmented with the difference in force obtained from the bonded interactions. To provide an initial proof-ofconcept the FDB method was implemented in the Tinysim MD engine, which is part of the NAMD-Lite package,26 before being fully implemented, as described below, within the AMBER


8

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


PMEMD engine27,28 for production simulations of proteins. Due to the lack of a thermostat in Tinysim, the Berendsen thermostat was implemented based on libraries in the NAMD-Lite package.26 Details concerning the implementation can be found in the Supporting Information.

Crown Ether System. The crown ether system consists of an 18-crown-6 crown ether coordinating to a positively charged ion. In the work describing the synthesis and physicochemical properties of crown ethers, which led to the shared Nobel Prize in Chemistry in 1987, Charles J. Pedersen showed that this particular crown ether has a higher affinity for potassium ions than for sodium ions.29,30 When having a potassium ion in the reference system, and the less favored sodium ion in the target system, the bias force should favor dissociation of the ion, and an increase of the bias parameter should cause a decrease in the dissociation time (Figure 3). The crown ether-ion complex was solvated in a water sphere with a diameter of 36 Å.

Figure 3. Crown ether test system solvated in a 36 Å diameter water sphere. The reference system contains a high-affinity potassium ion while the target system contains a lower-affinity sodium ion. Increasing the bias parameter should cause a decrease in the dissociation time.

Single Bond Rotation System. The single bond rotation system consists of mono- and diprotonated configurations of 1,1-diphenylethane-1,2-diamine (Figure 4A). Torsional energy profiles for the N-C-C-N dihedral angle show that the monoprotonated configuration has two


9


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

equivalent global minima, gauche plus and gauche minus, whereas the global minimum of the diprotonated configuration is the anti conformation (Figure 4B). Starting with either configuration in its more unfavorable minimum, increasing the bias parameter should decrease the time required for the molecule to reach its global minimum conformation through a rotation of the N-C-C-N angle.

Figure 4. (A) Mono- and diprotonated configurations of 1,1-diphenylethane-1,2-diamine. (B) Torsional energy profiles of the N-C-C-N angle for the two configurations. The energy of the global minimum of both configurations has been set to zero. Chair Flip System. A six-membered ring in a chair conformation typically prefers to have its substituents in the equatorial position in order to avoid 1,3-diaxial steric repulsion of axial substituents. However, for certain polyhydroxylated piperidine and hexahydropyridazine ions,


10

Page 11 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 polar substituents favor the axial position to such an extent that the axial conformation becomes the global minimum.31 The chair flip system consists of the trans-4,5dihydroxyhexahydropyridazine molecule, where protonation and deprotonation of the molecule cause it to switch between the two preferred chair conformations (Figure 5). Increasing the bias parameter should cause an increased rate of chair flips.

Figure 5. Chair flip of trans-4,5-dihydroxyhexahydropyridazine upon change in protonation.

Protein-Ligand Test Systems. In order to test the FDB method on protein systems, it was implemented into the PMEMD MD engine of AMBER v12,27,28 as it is not feasible to simulate large systems using Tinysim. The test systems consist of RBP and GlnBP, which are known to undergo domain closure upon binding of D-ribose and glutamine, respectively (Figure 6). As crystal structures of the closed conformations of the proteins exist, it is possible to obtain a measure of how well the method performs by calculating the root-mean-square deviation (RMSD) with respect to the closed states. The PDB IDs of the open and closed structures of RBP used in this study are 1URP32 (chain D) and 2DRI33, respectively. The PDB IDs of the open and closed structures of GlnBP are 1GGG34 (chain A) and 1WDN35, respectively. The binding modes of the two ligands in the open conformations were copied from the closed crystal structures. To determine if the respective ligand should be placed in the upper or the lower domain in the open structure (Figure 6), unbiased simulations were started from both positions. The correct binding site was determined to be in the N-terminal domain for both proteins, as the ligand was seen to


11


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

dissociate from the opposite domain (henceforth termed the lid domain) into the N-terminal pocket within 50 ns, yielding an exact or nearly identical binding mode compared to the closed crystal structure. Both proteins were prepared using the Protein Preparation Wizard available in the Schrödinger Suite version 2014.36-38

Figure 6. Domain closure of RBP and GlnBP upon binding of ribose and glutamine, respectively. The binding modes of the ligands were copied from the closed crystal structures. To determine whether the ligands should be placed in the upper or lower domains in the open structures, simulations were started from both positions.

Simulation Settings.


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


Settings for the small molecule test simulations. For all simulations, a topology was provided for both the target and the reference system. The simulations were run with very short relaxation times and started from a temperature of 5 K, followed by a temperature increase to the desired simulation temperature. This was done in order to avoid a temperature overshoot in the beginning of the simulations, where temperatures could otherwise assume values much higher than the desired temperature, until the Berendsen thermostat equilibrated the temperature of the system. As high temperatures are also able to accelerate conformational changes, this had to be avoided in order to prevent contamination of the results. All simulations were performed using the OPLS-2005 force field39 with a time step of 1 fs. Rotation and translation of the center-ofmass were removed at every step.

Crown ether simulations. Simulations were performed in which a potassium ion was present in the reference system, while a sodium ion was present in the target system. The bias force was applied to both the crown ether atoms and the ion. The crown ether simulations were performed using a thermostat relaxation time of 100 fs. The target temperature was set to 300 K. The cutoff value for non-bonded interactions was set to 12.0 Å. A shifting function was used on the electrostatic forces, while a switching function was applied to the vdW forces using a switching distance of 10 Å. The simulations were terminated when the distance from the ion to the centerof-mass of the crown ether exceeded 5 Å. The distance was checked every 1000 steps. A total of 40 simulations were run for each α-value, which ranged from 0.00 to 0.20 with a stepwise increase of 0.04 resulting in an aggregate total of 240 simulations.


13


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

Single bond rotation simulations. Simulations were performed for two setups: one in which the diprotonated configuration was used as the target system (thus having the monoprotonated configuration as the reference system) started from the less favored gauche conformation, termed the Gauche-to-Anti setup; and one in which the monoprotonated configuration was used as the target system (thus having the diprotonated configuration as the reference system) started from the less favored anti conformation, termed the Anti-to-Gauche setup. The single bond rotation simulations were performed using a thermostat relaxation time of 50 fs. The target temperature was set to 300 K. The simulations were run with no non-bonded interaction cutoff. The gaucheto-anti simulations were terminated when the N-C-C-N angle came within 30 degrees of the anti conformation. The anti-to-gauche simulations were terminated when the N-C-C-N angle came within 30 degrees of either of the two gauche conformations. The angle was checked every 100 steps. A total of 40 simulations were run for each α-value for both the gauche-to-anti and anti-togauche transitions. Bias parameter values ranged from 0.00 to 0.25 with a stepwise increase of 0.05 resulting in an aggregate total of 240 simulations.

Chair flip simulations. Simulations were performed for two setups: one in which the neutral configuration was used as the target system (thus having the charged configuration as the reference system) started from the less favored conformation in which the polar substituents are placed axially; and one in which the charged configuration was used as the target system (thus having the neutral configuration as the reference system) started from the less favored conformation in which the polar substituents are placed equatorially. The chair flip simulations were performed using a thermostat relaxation time of 100 fs. The simulations were run with no non-bonded interaction cutoff. The simulations having the neutral and charged configurations as


14

Page 15 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 target system were run at thermostat target temperatures of 300 K and 400 K, respectively. The higher temperature for the simulations having the charged configuration as the target system was chosen to reduce the simulation time needed to observe chair flips for this setup. Chair flips were detected by monitoring the dihedral angles in the cyclic ring every 5000 steps. A total of 200 simulations were run for each α-value for both the neutral and charged conformations. Bias parameter values ranged from 0.00 to 0.25 with a stepwise increase of 0.05 providing an aggregate total of 1200 simulations.

RBP and GlnBP Simulations: As demonstrated by Equation 2, for protein-ligand simulations it is only necessary to provide a topology for the target system, as the interactions between the protein and the non-existing ligand in the reference system is known a priori to be zero. For this type of system, the target-reference terminology thus becomes redundant. The protein was solvated with TIP3P40 water in a truncated octahedron box using a minimum spacing of 12 Å between any solute molecule and the edge of the box. The system was neutralized by addition of one and two sodium ions for the RBP and GlnBP systems, respectively, using the monovalent ion parameters developed by Joung and Cheatham.41 The simulations were performed using the ff12SB42 protein force field. The GLYCAM_06j43 force field was used for the ligand in the RBP system. Antechamber44 was used to obtain parameters for the ligand in the GlnBP system, using the general Amber force field (GAFF)45 and the AM1-BCC charge fitting method46,47. For all calculations, a non-bonded cutoff of 8 Å was applied with long-range electrostatics included via use of the PME method25 with a maximum grid spacing of 1.0 Å. The system was first energy minimized for 1000 steps using a combination of steepest descent (10 steps) followed by conjugate gradient to alleviate strain in the starting configuration. The system was then linearly


15


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

heated from 1 K to 300 K over the course of a 1 ns NVT simulation using Langevin dynamics with a collision frequency of 2 ps-1. This was followed by a 1 ns NPT simulation, in which the pressure was maintained at 1 atm using default settings. This equilibrated system was used for all subsequent production runs. All simulations were performed using a time step of 1 fs. The mdin files used in the individual steps are provided in the Supporting Information.

Coordinate Scan. The coordinate scans of the N-C-C-N angle in the single bond rotation system (Figure 4B) were performed with MacroModel 9.948 using the OPLS-2005 force field39. The scans were run in vacuo with no cutoff and with an increment of 1 degree per step.

RESULTS AND DISCUSSION The objective of the FDB method is to enhance the rate of conformational transitions by adding a bias force to the system. By gradually increasing the bias parameter, we thus expect to see an increased rate compared to the unbiased case, or equivalently, a decreased time before a transition from the initial conformation occurs. There will in both the biased and unbiased systems be a distribution of initial conformation lifetimes due to the statistical sampling over initial configurations, which we quantify by running multiple simulations and reporting the results in terms of boxplots showing minimum, mean and maximum values, as well as first and third quartile values.

Small molecule test systems. The time needed for a conformational change to occur as a function of the bias parameter is shown in Figure 7. The FDB method successfully accelerates the conformational changes for all the test systems with modest bias parameters of 0.1-0.2. As


16

Page 17 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


exemplified in Figure 7A, the mean time for a change to occur decreases monotonically, but the time for the slowest change observed for the simulations is not correlated with an increase in bias. This is most likely due to lack of proper statistics by the limited number of simulations. Figure 7 provides a proof-of-concept illustrating that at least for simple systems where the conformational transition can be described by a single geometric variable, the FDB method provides a valid acceleration of the dynamics.

Figure 7. Boxplots showing results for the (A) crown ether, (B) single bond rotation, and (C) chair flip simulations.

Protein-Ligand Test Systems. Although the FDB method provides strong evidence for acceleration of the dynamics for the small molecule test systems, the rearrangements observed in ligand-induced conformational changes of proteins are of a far more complex nature. The ability to accelerate the domain closures observed in RBP and GlnBP thus provides a much more realistic test.


17


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

RBP System. Initially, six 50 ns simulations were run with α-values of 0.0, 0.1, 0.2, 0.3, 0.4 and 0.5, respectively. In all simulations, the backbone RMSD with respect to the closed crystal structure was observed to fluctuate between 4.0-5.5 Å with no signs of incipient closure. A closer look at the open structure revealed that D215 in the binding site is part of a salt bridge with R141 (Figure 8B). In the closed structure, this salt bridge is broken, and D215 instead interacts with the ligand (Figure 8A). This should by itself not pose a problem for the FDB method. However, due to the way the ligand was placed in the active site of the open structure, where the coordinates were copied directly from the closed structure after alignment on the N-terminal pocket residues, the ligand became tucked beneath D215, where it could form a hydrogen bond to the backbone carbonyl of D215 (Figure 8B). In the unbiased simulation, where the ligand was placed in the lid domain, D215 rearranges to a conformation more similar to that observed in the closed structure as the ligand moves into the right binding mode (Figure 8C). This rearrangement is not observed in any of the biased simulations. It thus appears that direct placement of the ligand into the binding pocket is not the optimal approach, as rearrangements of the binding pocket during the binding process might be necessary to obtain a correct orientation of the surrounding side chains with respect to the ligand.


18

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


Figure 8. (A) In the closed structure, D215 interacts with the ligand. (B) When the position of the ligand is copied to the open structure, it becomes tucked beneath D215 that forms a salt bridge with R141. In this position, the ligand can form a hydrogen bond to the backbone carbonyl of D215. (C) In the unbiased simulation with the ligand initially placed in the lid domain, D215 rearranges to accommodate the incoming ligand.

A new set of 50 ns simulations was started from the conformation observed in Figure 8C, using α-values of 0.0, 0.1, 0.3, 0.5, 0.75, and 1.0. Five simulations were run for each α-value. The backbone RMSDs with respect to the closed crystal structure are shown in Figure 9. An RMSD of 1.5 Å was chosen to be the criterion for a successful conformational transition. An example of a snapshot with an RMSD of 1.5 Å is provided in Figure 10. For α-values of 0.0 and 0.1, the protein fluctuates around the open conformation (RMSD ~4.5 Å). For a value of 0.3, two of the five simulations reach the closed conformation. One of the two reaches RMSD values as low as 0.75 Å. A third simulation is observed to move toward the closed state, after which it returns to the open state. The two final simulations remain in the open conformation. For values of 0.5, 0.75 and 1.0, all simulations immediately move closer to the closed conformation (RMSD ~3 Å). A single α = 0.5 simulation reaches the closed conformation within 6.5 ns. An α = 0.75


19


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

simulation briefly reaches the fully closed conformation before returning to the semi-closed conformation. Finally, a single α = 1.0 simulation reaches the closed conformation. A larger number of simulations are required to quantify the dependence of the domain closure rate on the bias parameter, and in particular, longer simulations need to be performed for the unbiased case, but these preliminary simulations do at least show that the bias clearly has an accelerating effect since many of the simulations quickly reach a semi-closed state from which transitions into the fully closed structure can occur.

Figure 9. Backbone RMSD values compared to the closed crystal structure of RBP. The initial structure has an RMSD of 4.5 Å and the horizontal line in each plot denotes an RMSD of 1.5 Å.


20

Page 21 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


Figure 10. Example of a simulation snapshot (orange) with an RMSD of 1.5 Å compared to the closed crystal structure (blue). (A) Entire protein. (B) Active site residues surrounding the ligand.

GlnBP system. Six 50 ns simulations were run with α-values of 0.0, 0.1, 0.2, 0.3, 0.4 and 0.5, respectively. The α = 0.0 and α = 0.1 simulations remained in the open state for the entirety of the 50 ns simulation time. The rest of the simulations all undergo domain closure. In all cases, however, the observed closed structure is skewed compared to the crystal structure (Figure 11A). This appears to be due to an electrostatic interaction between K115 in the lid domain and the carboxylate group in the ligand. This interaction is not present in the crystal structure, where K115 instead forms a salt bridge with D10 in the N-terminal domain (Figure 11B). Increasing the amount of bias increases the strength of the K115-ligand interaction. If the bias is too high, this interaction outcompetes the K115-D10 interaction, and the lid domain is pulled into the skewed conformation. The distance between the K115 side chain nitrogen and the D10 carboxylate group is plotted for the simulations with α-values of 0.1, 0.3, and 0.5 in Figure 11C. The distance in the α = 0.1 simulation fluctuates around the distance observed in the open structure. In the α = 0.3 simulation, the lysine moves toward the correct position before ending up in the skewed conformation. In the α = 0.5 simulation, the lysine moves directly into the skewed conformation. The skewed conformation remains stable when extended with 50 ns of unbiased simulation.


21


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

These results demonstrate that the bias force contains excessive orthogonal, off-path contributions to the target force for α-values of 0.2 and above. In an attempt to alleviate this problem, simulations were run in which the bias was applied in intervals. The role of the bias in this case is to nudge the protein in the right direction, while allowing it to relax onto the original PES in the periods without bias. An α = 0.9 simulation, in which bias was applied for 10 ps followed by 30 ps of unbiased simulation, briefly reached an RMSD of ~2 Å and the correct K115-D10 distance within 60 ns. However, a combination that consistently resulted in stable, closed conformations comparable to the crystal structure was not found.

Figure 11. The closed state obtained in the biased GlnBP simulations is different from the crystal structure. (A) The lid domain moves too far in the direction of the black arrow in the simulations (orange structure). (B) In the crystal structure (gray), K115 forms a salt bridge with


22

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


D10. In the skewed, bias induced closed conformation (green) it instead forms a salt bridge with the ligand. The orange, purple and green structures correspond to the dots in panel C. (C) The distance between the K115 side chain nitrogen and the D10 carboxylate group for α-values of 0.1 (orange), 0.3 (purple), and 0.5 (green). The horizontal line indicates the distance observed in the crystal structure. The α = 0.3 simulation moves in the right direction before it goes off-path after 17 ns and finally reaches the skewed conformation.

The results for the two protein-ligand systems reach different conclusions. The domain closure of RBP is successfully accelerated when a bias force is applied to the system. The simulations reveal that direct placement of the ligand in the binding pocket can in some cases hinder the conformational change, as the rearrangements of the side chains that are supposed to occur during the binding process are bypassed. Instead of attempting to obtain a high-scoring pose in a pre-simulation docking calculation, the best approach is likely to place the ligand loosely in the binding pocket and run an unbiased simulation in which the ligand and the binding pocket can adjust to each other prior to the bias simulations. A second observation is that it is not obvious from the simulations when the target conformation has been reached, unless the structure is already known. Several of the simulations reach semi-closed conformations that are stable for the duration of the simulation. Was it not for the simulations that reach the fully closed conformation, these could mistakenly be identified as the end state. This scenario is difficult to avoid, and it was also encountered for the recently published linear response path following method,49 which bears many similarities to the FDB method. An optimal method would be deterministic and always reach the final closed conformation within a short time; however, the path sampled during an MD simulation is largely dependent on the initial velocities, the applied


23


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

thermostat and barostat, and numerical noise. Running multiple simulations will increase the chances of reaching the final structure and thereby help to identify simulations trapped in stable intermediate states. The GlnBP system behaves qualitatively different. Although the bias force in this case also increases transition rates, a large bias value favoring rapid transitions leads to a risk of ending up in the wrong conformation, as specific protein-ligand interactions might suppress important protein-protein interactions. Unless the structure of the closed conformation is already known, this can be difficult to detect. Further studies are required to probe which, if any, of the RBP or GlnBP test cases are representative for the general case. For the protein-ligand simulations presented in this paper, the temperature is controlled equally well for the unbiased and the biased simulations, with similar average temperatures and temperature fluctuations. Subsequent tests of the RBP system have shown that heating of the protein is only observed for excessively large values of the bias parameter. Thus, when increasing the bias parameter, undesirable temperature effects do not arise until after the applied bias forces become unreliable due to off-path contributions. This result might vary depending on the utilized thermostat.

CONCLUSION We propose a force difference bias method for accelerating ligand-induced conformational changes in proteins by molecular dynamics simulations. The bias force is obtained by increasing the non-bonded protein-ligand interactions and is thus calculated from properties of the system itself. This eliminates the need for user defined collective variables and allows the bias force to change as the structure changes. It furthermore only requires that the apo protein structure is


24

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


available. Proof-of-concept is presented for three small molecule test systems, and ligand induced domain changes for proteins are probed for two systems where the end-point structures are known. Domain closure is successfully accelerated for the RBP system, while simulations of the GlnBP system reveal a methodological hurdle that has to be overcome before the method can reliably be applied to the study of systems in which only the apo conformation is known.

ACKNOWLEDGEMENTS This work was supported by grants from the Danish Council for Independent Research | Technology and Production Sciences (FTP 11-105010) and the Danish National Research Foundation (DNRF59). RCW and PN were supported by an NSF SI2-SSE grant (NSF-1148276) to RCW. RCW also acknowledges funding through fellowships from Intel Corp. and NVIDIA, Inc. Computations were possible through allocations of time at the Centre for Scientific Computing, Aarhus. OJA was supported by a grant from the Novo Scholarship Program.

ASSOCIATED CONTENT Supporting Information. Details related to the implementation of the method in the Tinysim MD engine and the mdin files used for the RBP and GlnBP simulations are provided as Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org.


25


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

REFERENCES

(1) Henzler-Wildman, K.; Kern, D. Dynamic Personalities of Proteins. Nature 2007, 450, 964972. (2) Zheng, W.; Brooks, B. R.; Hummer, G. Protein Conformational Transitions Explored by Mixed Elastic Network Models. Proteins: Struct., Funct., Bioinf. 2007, 69, 43-57. (3) Ha, T.; Enderle, T.; Ogletree, D. F.; Chemla, D. S.; Selvin, P. R.; Weiss, S. Probing the Interaction between Two Single Molecules: Fluorescence Resonance Energy Transfer between a Single Donor and a Single Acceptor. Proc. Natl. Acad. Sci. U. S. A. 1996, 93, 6264-6268. (4) De Vlaminck, I.; Dekker, C. Recent Advances in Magnetic Tweezers. Annu. Rev. Biophys. 2012, 41, 453-472. (5) Shaw, D. E.; Dror, R. O.; Salmon, J. K.; Grossman, J. P.; Mackenzie, K. M.; Bank, J. A.; Young, C.; Deneroff, M. M.; Batson, B.; Bowers, K. J., et al. In Millisecond-Scale Molecular Dynamics Simulations on Anton, Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, Portland, OR, Nov 14-20, 2009; Association for Computing Machinery: New York, NY, 2009. (6) Shaw, D. E.; Deneroff, M. M.; Dror, R. O.; Kuskin, J. S.; Larson, R. H.; Salmon, J. K.; Young, C.; Batson, B.; Bowers, K. J.; Chao, J. C., et al. Anton, A Special-Purpose Machine for Molecular Dynamics Simulation. Commun. ACM 2008, 51, 91-97. (7) Dykeman, E. C.; Sankey, O. F. Normal Mode Analysis and Applications in Biological Physics. J. Phys.: Condens. Matter 2010, 22, 423202. (8) Stein, S. A. M.; Loccisano, A. E.; Firestine, S. M.; Evanseck, J. D. In Annual Reports in Computational Chemistry; Spellmeyer, D. C., Ed.; Elsevier, 2006; Vol. 2, pp 233-262. (9) Adcock, S. A.; McCammon, J. A. Molecular Dynamics: Survey of Methods for Simulating the Activity of Proteins. Chem. Rev. 2006, 106, 1589-1615. (10) Hamelberg, D.; Mongan, J.; McCammon, J. A. Accelerated Molecular Dynamics: A Promising and Efficient Simulation Method for Biomolecules. J. Chem. Phys. 2004, 120, 11919-11929. (11) Hamelberg, D.; de Oliveira, C. A. F.; McCammon, J. A. Sampling of Slow Diffusive Conformational Transitions with Accelerated Molecular Dynamics. J. Chem. Phys. 2007, 127, 155102. (12) Tirion, M. M. Large Amplitude Elastic Motions in Proteins from a Single-Parameter, Atomic Analysis. Phys. Rev. Lett. 1996, 77, 1905-1908.


26

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


(13) Hinsen, K. Analysis of Domain Motions by Approximate Normal Mode Calculations. Proteins: Struct., Funct., Bioinf. 1998, 33, 417-429. (14) Bahar, I.; Rader, A. Coarse-Grained Normal Mode Analysis in Structural Biology. Curr. Opin. Struct. Biol. 2005, 15, 586-592. (15) Chennubhotla, C.; Rader, A. J.; Yang, L. W.; Ivet Bahar, I. Elastic Network Models for Understanding Biomolecular Machinery: From Enzymes to Supramolecular Assemblies. Phys. Biol. 2005, 2, S173-S180. (16) Schlitter, J.; Engels, M.; Krüger, P.; Jacoby, E.; Wollmer, A. Targeted Molecular Dynamics Simulation of Conformational Change-Application to the T ↔ R Transition in Insulin. Mol. Simul. 1993, 10, 291-308. (17) Schlitter, J.; Engels, M.; Krüger, P. Targeted Molecular Dynamics: A New Approach for Searching Pathways of Conformational Transitions. J. Mol. Graphics 1994, 12, 84-89. (18) Hub, J. S.; de Groot, B. L. Detection of Functional Modes in Protein Dynamics. PLoS Comput. Biol. 2009, 5, e1000480. (19) Kästner, J. Umbrella Sampling. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 932942. (20) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo FreeEnergy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187-199. (21) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562-12566. (22) Laio, A.; Gervasio, F. L. Metadynamics: A Method to Simulate Rare Events and Reconstruct the Free Energy in Biophysics, Chemistry and Material Science. Rep. Prog. Phys. 2008, 71, 126601. (23) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 826-843. (24) Pierce, L. C. T.; Salomon-Ferrer, R.; de Oliveira, C. A. F.; McCammon, J. A.; Walker, R. C. Routine Access to Millisecond Time Scale Events with Accelerated Molecular Dynamics. J. Chem. Theory Comput. 2012, 8, 2997-3002. (25) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092. (26) Hardy, D. J. NAMD-Lite, version 1.0.2; University of Illinois at Urbana-Champaign, IL, 2007.


27


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

(27) Case, D. A.; Darden, T. A.; Cheatham III, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M., et al. AMBER 12; University of California, San Francisco, CA, 2012. (28) Salomon-Ferrer, R.; Case, D. A.; Walker, R. C. An Overview of the Amber Biomolecular Simulation Package. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2013, 3, 198-210. (29) Pedersen, C. J. Cyclic Polyethers and their Complexes with Metal Salts. J. Am. Chem. Soc. 1967, 89, 7017-7036. (30) Pedersen, C. J.; Frensdorff, H. K. Macrocyclic Polyethers and their Complexes. Angew. Chem., Int. Ed. Engl. 1972, 11, 16-25. (31) Jensen, H. H.; Lyngbye, L.; Jensen, A.; Bols, M. Stereoelectronic Substituent Effects in Polyhydroxylated Piperidines and Hexahydropyridazines. Chem. - Eur. J. 2002, 8, 12181226. (32) Björkman, A. J.; Mowbray, S. L. Multiple Open Forms of Ribose-Binding Protein Trace the Path of its Conformational Change. J. Mol. Biol. 1998, 279, 651-664. (33) Björkman, A. J.; Binnie, R. A.; Zhang, H.; Cole, L. B.; Hermodson, M. A.; Mowbray, S. L. Probing Protein-Protein Interactions - the Ribose-Binding Protein in Bacterial Transport and Chemotaxis. J. Biol. Chem. 1994, 269, 30206-30211. (34) Hsiao, C.; Sun, Y.; Rose, J.; Wang, B. The Crystal Structure of Glutamine-Binding Protein from Escherichia Coli. J. Mol. Biol. 1996, 262, 225-242. (35) Sun, Y.; Rose, J.; Wang, B.; Hsiao, C. The Structure of Glutamine-Binding Protein Complexed with Glutamine at 1.94 Å Resolution: Comparisons with Other Amino Acid Binding Proteins. J. Mol. Biol. 1998, 278, 219-229. (36) Protein Preparation Wizard, version 2014-1: Epik, version 2.4; Impact, version 5.9; Prime, version 3.2; Schrödinger, LLC, New York, NY, 2014. (37) Sastry, G. M.; Adzhigirey, M.; Day, T.; Annabhimoju, R.; Sherman, W. Protein and Ligand Preparation: Parameters, Protocols, and Influence on Virtual Screening Enrichments. J. Comput.-Aided Mol. Des. 2013, 27, 221-234. (38) Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J. Chem. Theory Comput. 2011, 7, 525-537. (39) Banks, F. L.; Beard, H. S.; Yixiang, C.; Cho, A. R.; Damm, W.; Farid, R.; Felts, A. K.; Halgren, T. A.; Mainz, D. T.; Maple, J. R., et al. Integrated Modeling Program, Applied Chemical Theory (IMPACT). J. Comput. Chem. 2005, 26, 1752-1780.


28

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


(40) 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, 926935. (41) Joung, I. S.; Cheatham III, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020-9041. (42) AmberTools12 Reference Manual, http://ambermd.org/doc12/AmberTools12.pdf (accessed October 16, 2015). (43) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; González-Outeiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates. J. Comput. Chem. 2008, 29, 622-655. (44) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25, 247260. (45) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157-1174. (46) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of HighQuality Atomic Charges. AM1-BCC Model: I. Method. J. Comput. Chem. 2000, 21, 132146. (47) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: II. Parameterization and Validation. J. Comput. Chem. 2002, 23, 1623-1641. (48) MacroModel, version 9.9; Schrödinger, LLC, New York, NY, 2011. (49) Tamura, K.; Hayashi, S. Linear Response Path Following: A Molecular Dynamics Method to Simulate Global Conformational Changes of Protein upon Ligand Binding. J. Chem. Theory Comput. 2015, 11, 2900-2917.


29


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

For Table of Contents Use Only.

Toward an Enhanced Sampling Molecular Dynamics Method for Studying Ligand-Induced Conformational Changes in Proteins Ole Juul Andersen, Julie Grouleff, Perri Needham, Ross C. Walker and Frank Jensen


30

Toward an Enhanced Sampling Molecular Dynamics Method for

Recommend Documents