Ligand Selectivity Mechanism and Conformational ... - ACS Publications

Subscriber access provided by HACETTEPE UNIVERSITESI KUTUPHANESI

Article

Ligand Selectivity Mechanism and Conformational Changes in Guanine Riboswitch by Molecular Dynamics Simulations and Free Energy Calculations Guodong Hu, Aijing Ma, and Ji Hua Wang J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00139 • Publication Date (Web): 27 Mar 2017 Downloaded from http://pubs.acs.org on March 27, 2017

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.

Journal of Chemical Information and Modeling 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 34

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

Journal of Chemical Information and Modeling

Ligand Selectivity Mechanism and Conformational Changes in Guanine Riboswitch by Molecular Dynamics Simulations and Free Energy Calculations †

†

Guodong Hu,*, Aijing Ma, Jihua Wang*, †

†

Shandong Provincial Key Laboratory of Biophysics, College of Physics and Electronic Information, Dezhou University, Dezhou 253023, China

ABSTRACT Riboswitches regulate gene expression through direct and specific interactions with small metabolite molecules. Binding of a ligand to its RNA target is high selectivity and affinity and induces conformational changes of the RNA’s secondary and tertiary structure. The structural difference of two purine riboswitches aptamers is caused by only one single mutation, where cytosine 74 in the guanine riboswitch is corresponding to a uracil 74 in adenine riboswitch. Here we employed molecular dynamics (MD) simulation, molecular mechanics Poisson Boltzmann

surface

area

(MM-PBSA)

and

thermodynamic

integration

computational

methodologies to evaluate the energetic and conformational changes of ligands binding to purine riboswitches. The snapshots used in MM-PBSA calculation were extracted from ten 50 ns MD simulation trajectories for each complex. These free energy results are in consistent with the experimental data and rationalize the selectivity of the riboswitches for different ligands. In particular, it is found that the loss in binding free energy upon mutation is mainly electrostatic in 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

Page 2 of 34

guanine (GUA) and riboswitch complex. Furthermore, new hydrogen bonds are found in mutated complexes. To reveal the conformational properties of guanine riboswitch, we performed a total of 6 µs MD simulations in both the presence and the absence of the ligand GUA. The MD simulations suggest that the conformation of guanine riboswitch depends on the distance of two groups in the binding pocket of ligand. The conformation is in a close conformation when U51A52 is close to C74-U75.

KEYWORDS: guanine riboswitch; molecular dynamics simulation; binding free energy; ligand selectivity; conformational change

2


Page 3 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION RNA plays a central role in almost every genetic process of the cell. Despite of its polyelectrolyte character, RNA can fold into intricate three-dimensional structures with pockets and cavities that are essential for function and have the potential to bind ligands specifically.1,2 A particularly interesting type of RNA molecules has been named riboswitch, as they have the ability to switch gene expression at special physiological conditions.3 Riboswitches are found with the highest frequency in the 5’ untranslated regions (UTRs) of many bacteria and fungi.4,5 They are composed of two distinct yet mutually interacting domains: an aptamer and an expression platform. The aptamer directly binds the small molecule, and the expression platform undergoes structural changes in response to the changes in the aptamer. Thus riboswitches exercise ligand-dependent control on down stream gene expression.6-8 Riboswitches may assume two different conformations depending on whether the small molecule specific to the riboswitch is bound or unbound.9 Expressions of gene regulated by most of riboswitches are widely found in bacteria; However, only the thiamine pyrophosphate-responsive riboswitch has thus far been employed by eukaryotes. So riboswitches have attracted interest as potential drug targets for antibacterial.10 Ligands to riboswitches are generally complexes of primary metabolism11, including amino acids (glycine and lysine)12,13, nucleotides (adenine and guanine)14,15, cofactors such as Sadenosylmethionine,

S-adenosylhomocysteine,

cobalamine,

adenosylcobalamin,

mononucleotide16, thiamine pyrophosphate17, and metal ions such as Mg2+

18

flavin

. Riboswitches

normally bind their ligands with high affinity and often reject even close chemical relatives of their cognate ligands.19 The purine riboswitches, containing guanine14 and adenine15 classes, were identified in 2003. The crystallographic structures of the guanine and adenine bound aptamer domains revealed the aptamer domain consist of three helices P1, P2 and P3, forming a three-way junctions, with the junction loops J12, J23 and J31 connecting them and with loops L2 and L3 capping P2 and P3, respectively (Figure 1A).14,15 The nucleobase is bound by nucleotides in the center of a three-way 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 34

helical junction with only 5% of its surface still solvent accessible,20 which indicates an induced fit mechanism of binding.3 The nucleotides form a novel intermolecular base triple with the ligand.20,21 Within this binding pocket, it interacts with all functional groups of the ligand and forms a Watson-Crick base pair to C74 for guanine (GUA) (Figure 1B) and the corresponding U74 for the adenine. It is striking that the purine riboswitches discriminate between guanine and adenine by at least 10,000-fold. To understand how the purine riboswitch achieves high degree of specificity for closely related complexes, Batey et al. investigated their interaction with purine analogs with varying functional groups at the 2- and 6-positions (Figure 1C) that have the potential to alter interactions with nucleotide 74 by using a combination of crystallographic and calorimetric approaches. Several crystal structures of the B. subtilis xpt-pbuX guanine riboswitch (GR) and its mutation C74U (GR-C74U) in complexes with different ligands have been solved and published together with binding data.22 Aqvist et al. have calculated the energies of different purine complexes with the guanine and adenine by employing molecular dynamics free energy perturbation calculations and the linear interaction energy method.3 The molecular docking and molecular simulation method were employed to study the interactions and dynamic behavior of GR with ligands by Liu’s group.23,24 On the other hand, a number of experiments have been conducted to investigate the conformational changes of purine riboswitches.8,25 Brenner and co-workers25 studied the conformational dynamics of xpt GR aptamer using single-molecule fluorescence resonance energy transfer spectroscopy. The experiments have demonstrated that ligand binding and dissociation occur on a time scale of seconds,26,27 implying that the conformational rearrangements of riboswitch is a very slow process.19 Molecular dynamics (MD) simulations at the atomic level are at present the most appropriate way to obtain detailed information about the binding of small molecules to RNA.19,28 MD simulations on the aptamer domain of the guanine, adenine, and S-adenosylmethionine sensing riboswitches in both ligand bound and unbound states were reported.19,29-31 Mitra et al. reported 15 ns MD simulations to study the add adenine riboswitch crystal structure, both in the adenine 4


Page 5 of 34

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


bound state and in the adenine free state.8 MD simulations are particularly true if not only dynamic or structural characteristics of RNA-ligand complex formations are considered, but if a comparison also involves the energetics of binding of RNA–ligand complex.32 There exists many computational methods with various levels of computational accuracy and expense to estimate the inhibitor binding affinities,33-35 such as the free energy perturbation (FEP), molecular mechanics Poisson-Boltzmann surface area (MM-PBSA)36,37, and thermodynamic integration (TI)38,39. These methods have been successfully applied not only for the protein and ligand complexes,34,40 but also for the RNA and ligands complexes.3,41-43 The MM-PBSA can be used as a fast method to rank binding affinities of ligands targeting RNA.32 This method combines molecular mechanics (MM) energies and entropic contributions with continuum solvent models to post-process a series of representative snapshots from MD simulation trajectories.44-46 In this work, we analyzed the binding mechanism of four guanine analogues complex with the GR and GR-C74U. The GR and GUA complex is shown in Figure 1A, as well as the four ligands in Figure 1C. The MM-PBSA calculations with the snapshots extracted from ten 50 ns MD simulations for each complex and the TI methods were combined to perform a comparative study. The estimation of the binding free energy terms and its relative changes between GR and GR-C74U provide valuable information into drug target relationships. The dynamics and conformational changes of GR both in the ligand GUA bound and unbound states were investigated via a total of 6 µs MD simulations. The results of long MD simulations offer helpful information to understand how the ligand binding influences the conformational change. METHODS System preparation. Atomic coordinates of GR in complex with 2BP, XAN, 6GU complexes were obtained from the Protein Data Bank (PDB) with PDB IDs 3G4M, 3GAO and 3GER, respectively.22 The structures of ligands are shown in Figure 1C. The starting model of GR and GUA complex was constructed with the crystal structure 3G4M as the template. We mutated C74 of GR to U74 to obtain the starting structures of GR-C74U complexes. All crystallographic water molecules were retained. The AMBER force field (FF12SB)47 was used to describe the RNA 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 34

parameters, ions and water molecules. Single-point calculations with Gaussion 0348 were performed to obtain the electrostatic potential around each ligand by using Hartree-Fock/6-31G* basis set. The general AMBER force field parameters (GAFF)49 and the partial charges were assigned using the RESP method50 with the Antechamber module of AMBER12 package. All complexes were solvated in a rectangular periodic box of TIP3P51 water molecules with a margin distance of 12 Å. To neutralize the charges of systems, an appropriate number of sodium counterions were added. Molecular dynamics simulation. All energy minimizations and MD simulations were carried out using the AMBER12 package52. To relieve the bad contacts and to direct each system toward energetically favorable conformations, we performed a two-step extensive energy minimization process by the steepest descent method and the conjugate gradient algorithm. After minimization, the systems were heated from 0 K to 300 K over 70 ps applying harmonic restraints with force constants of 2 kcal/(mol·Å2) to all atoms of complex. Subsequent a constant-pressure MD simulation was carried out for 90 ps to adjust the solvent density. Finally, the production constant-pressure MD simulation without any restraints for each system lasted 50 ns at 300 K with the Langevin thermostat and a target pressure 1.0 atm using isotropic positional scaling, as well as 1 µs for GR in GUA bound and unbound states. The long range electrostatic interactions were treated with particle mesh Ewald method.53 SHAKE algorithm was employed to treat the covalent bonds involving hydrogen atoms,54 the time step was set to be 2 fs with a cutoff of 12 Å for all MD simulations. Coordinate trajectories were recorded every 10 ps throughout all equilibration and production runs. MM-PBSA calculations. For each complex，ten stable MD trajectories were obtained to estimate the binding free energies ( ∆Gbind ) by using MM-PBSA method. A number of 600 snapshots were chosen evenly from the last 30 ns MD trajectories with an interval of 50 ps for each trajectories, and a total number of snapshots was 6000 for each complex. All counterions and water molecules of each snapshots were stripped before the MM-PBSA calculation. The

6


Page 7 of 34

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


MM-PBSA calculation used in this work is in the same way as in other works. Briefly, the binding free energies ∆Gbind were computed by the following equations55,56:

∆Gbind = ∆EvdW + ∆Eele + ∆Gpol + ∆Gnonpol - T ∆S

（1）

Each term can be computed by the component of complex minus the sum of components of RNA and ligand. ∆Eele and ∆EvdW are the electrostatic and van der Waals components in the gas phase, respectively. These two components were computed using the same parameter set as those used in the MD simulations. The third and the fourth terms are polar ( ∆Gpol ) and nonpolar (

∆Gnonpol ) components of the solvation free energy. The ∆Gpol was calculated by numerically solving the Poisson Boltzmann equation, and the nonpolar component ( ∆Gnonpol ) was determined using ∆Gnonpol = γ SASA + β , where SASA is the solvent-accessible surface area that was determined with MSMS program57,58 with a probe radius of 1.4 Å, γ and β were set to be 0.005 kcal ∙ molିଵ ∙ Åିଶ and 0.0 kcal/mol, respectively59. The conformational entropic contribution to the binding free energy was estimated for a total number of 3000 snapshots using the normal mode analysis.60 Thermodynamic Integration Calculations. From the thermodynamic cycle in Figure 2, we can calculate the relative binding free energies between B1 and B2 as61,62

∆∆Gbind = ∆Gbind ( B2 ) − ∆Gbind ( B1 ) = ∆Gbound ( B1 → B2 ) − ∆G free ( B1 → B2 )

（2）

where ∆Gbind is the binding free energy and ∆G ( B1 → B2 ) is the free energy difference of mutating B1 into B2 , when B is in bound state with A or in free state. The ∆G ( B1 → B2 ) are calculated using TI method which based on the following equation63: n

∆G = ∑ Wi i =1

∂V ∂λ

（3） λi

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 34

Where λ is the coupling parameter, V is the λ coupling potential function. The sampling space was characterized by n=12. The Gaussian quadratic formula was used to design the values of λi and its weight Wi . In all of the calculations, their values successfully tested by previous studies64,65 were taken from the Amber 12 package. The angle brackets represent a Boltzmannweighted average. The process of TI was divided into three steps for each λ . We ran two independent 2 ns MD simulations for each step. In the first step, the charges on atoms that change in the mutation are zeroed. In the second step, B1 is changed into B2 keeping the charges of mutated atoms zeroed. Finally, the charges are introduced also on the mutated atoms of B2 . A soft-core version of the Lennard-Jones potential,66 as implemented in the Amber simulation package, was used to treat the disappearance and appearance of atoms in the second step. RESULTS AND DISCUSSION Analysis of binding free energy from MM-PBSA. The absolute free energies for all complexes of GR and GR-C74U binding with ligands were calculated by using mm_pbsa program in AMBER 12 and summarized in Table 1. Predicted absolute free energies were larger than the experimental results; However, the ranking order is the same. The MM-PBSA method can help us understand binding mechanism of complex in detail, because this method can decompose the total binding free energy into individual energetic terms.36 The van der Waals interactions and the nonpolar solvation energies ( ∆Gvdw+ nonpol ) are responsible for the burial of inhibitor’s hydrophobic groups upon binding. It’s clear that these two terms are favorable for all ligand binding. The electrostatic energies ( ∆Gele + pol ) are unfavorable in all complexes except in GR and GUA complex. The entropic contributions ( -T ∆S ) are unfavorable for GR complexes with a mean value of 12.26 kcal/mol, which is almost equal to the mean value of GR-C74U complexes of 12.16 kcal/mol. According to the experimental data, the 6GU and 2BP ligands bind more favorable to GR than to GR-C74U, the calculated total free energies of these two ligands also give similar results. The calculated results also indicate that the mutation of nucleotide 74 8


Page 9 of 34

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


from C74 to U74 causes an increase in binding free energies for GR-C74U binding to XAN but a large decrease to GUA. Comparisons the energetic terms between the complexes of GR and GRC74U, we found that the differences of ∆Gele + pol between GR and GR-C74U are similar to those of the binding free energies ( ∆Gbind ) for all ligands; However, the differences of ∆Gvdw+ nonpol differ with those of ∆Gbind as the RNA from GR to GR-C74U. The ∆Gbind of GUA becomes weakly, while ∆Gvdw+ nonpol is more favorable. The ∆Gbind of XAN becomes strongly, the

∆Gvdw+ nonpol is less favorable. This means that the electrostatic energy would be an important term for the selectivity of GR. Revelation from GR to GR-C74U based on TI calculation. We utilized the more costly and rigorous TI methodology to test the reliable of the results from MM-PBSA method. Table 2 shows the relative binding free energies ( ∆∆ GTI ) for GR and GUA complex relative to the other three complexes. Predicted relative rank order of TI calculations is in agreement with the experimentally determined results. Table 3 shows that ∆∆ GTI for these four ligands and GR complexes relative to their C74U mutated GR-C74U complexes. The results suggest that the binding affinity of GUA binding to GR is stronger than to GR-C74U. The opposite is true for the other three ligands. This is not only accordance with the experimental suggestions but also with the MM-PBSA data. The ∆∆ GTI were calculated by three steps. Two energy differences ∆∆ G step1 caused by removing the charges in step1 and ∆∆ G step 3 caused by adding charge in step

3 reflect the changes in polar interactions. The ∆∆ G step 2 involves in the change of the van der Waals interactions in step 2. Free energy plots for each step of the four GR to GR-C74U TI computations are shown in Figure S1 of the supporting information. Large fluctuations were found for the free energy curves of step3 for all four systems. This suggests that the polar interaction plays an important role for the selectivity of GR.

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 34

Identification of the key nucleotides responsible for the binding of inhibitor. We calculated the contribution of each nucleotide to the binding free energy by utilizing the free energy decomposition method to investigate the differences of binding mechanism between GR and GRC74U binding to ligands. The results were shown in Figure 3 for the complexes of GUA, as well as in Figure S2 for the complexes of the other three ligands. The contribution of a nucleotide is divided into the electrostatic interaction energy, van der Waals energy, polar solvation contribution and nonpolar solvation contribution. All energy components were calculated using the same parameters as the total binding free energy, except for the polar solvation contribution with the generalized Born solvation model. This method has been successfully applied to many protein-inhibitor binding systems35,37. From Figure 3 and Figure S2, we noticed that the interaction spectra are very similar for four ligands and GR complexes. There are six nucleotides (A21, U22, U51, A52, C74 and U75) with favorable average energy less than -1.0 kcal/mol. It is well known that two types of non-covalent interactions affected the structure and dynamics of nucleic acids, one is hydrogen bond (H-bond) mostly governed by electrostatics, the other one is aromatic stacking mostly governed by van der Waals.67,68 Among the six key nucleotides, two nucleotides (U37 and C60) contribute large electrostatic energy, the other four nucleotides contribute large van der Waals energy. Comparisons between the spectra of GR and GR-C74U, only GUA’s complexes show large differences (Figure 3). The crystallographic structure of GUA and GR complex suggests that seven H-bonds were formed between GUA and GR (Figure 1B). The H-bonds analysis were carried out between GR and ligands for four complexes of GR, as well as for four complexes of GR-C74U (Table 4). The results implies that U51 forms three H-bonds with the ligand, which has an amino at 2-position, during all ten MD simulations of the GR complexes. For the complexes of GR, we also found that nucleotide C74 forms three, two and one H-bonds with GUA, 2BP/6GU and XAN, respectively. This means that the substitution of a hydrogen/chlorine atom at 6-position of 2BP/6GU for an oxygen atom at 6-position and a hydrogen atom at 1-position of GUA results in a H-bond lost and the electrostatic contribution of C74 decrease obvious. In fact, the two H-bonds between C74 and 2BP/6GU are reformed, the side chain of nucleotide of C74 is arranged to suit 10


Page 11 of 34

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


for forming these two H-bonds (see Figure 4A). In XAN complex, the substitution of an oxygen atom for an amino at 2-position results in the H-bonds involved the amino lost, and a new stable H-bond between N3 atoms of C74 and ligand formed. We have calculated the average RMSDs of all heavy atoms within 5Å of ligand over a sum 30000 snapshots extracted from ten last 30 ns MD trajectories relative to their crystallographic structures. The average RMSDs of XAN complex (1.13 Å) is larger than that of GUA complex (0.94Å). And the volume of XAN is smaller than that of GUA. Hence, XAN is less stable compared to GUA in the binding pocket. The H-bond analysis also indicates that the U22 forms a H-bond with ligand, however, U22 and the other three nucleotides (A21, A38 and U75) mainly contribute large van der Waals energies. As shown in Figure 4B, the interaction between these four nucleotides and GUA are mainly the aromatic stacking. The nucleobases of U22 and A38 are on one side of the ligand, and the nucleobases of A21 and U75 are on the other side. Comparison among complexes of 2BP, 6GU and GUA. Three nucleotides (A21 U22 and U51) with different contribution large than 0.5 kcal/mol are found by comparing the interaction spectra between 2BP_GR and 6GU_GR (Figure S2). The nucleotide A21 binding to 6GU (-2.82 kcal/mol) is more favorable than to 2BP (-1.30 kcal/mol). As shown in Figure 5, the nucleobase of U22 is located at the other side of ligand relative to that of A21, its contribution in 6GU is less favorable than in 2BP. We calculated the average distances between the ligands and the nucleobase of A21 or U22 over 30000 snapshots, indicating that the 6GU is close to nucleobase of A21 and far from nucleobase of U22 relative to 2BP (Figure 5). As a result, the van der Waals energy of A21 is larger in 6GU complex than in 2BP complex. The H-bond analysis shows U22 forms stable H-bond with 2BP/6GU. Because the relative position of 2BP in the ligand binding pocket is different from that of 6GU. Thus the difference of electrostatic energy is the main factor for U22 in complexes of 2BP and 6GU (Table 4). The H-bond analysis show that U51 forms three stable H-bonds in 6GU and 2BP complexes. However, two H-bonds are found in XAN because the 2-position is an oxygen atom (Table 4). Thus the contribution of U51 is mainly from the electrostatic contribution in the three complexes. The H-bond analysis also suggests that U51 bonds weakly to XAN relative to 6GU/2BP, this is in 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 34

accordance with the decomposition of energy to individual residue (Figure S2). Only one nucleotide C74 was found with different contribution between GUA_GR and 6GU_GR complexes larger than 0.5 kcal/mol by comparing their interaction spectra. The large difference is mainly from the electrostatic contribution. This can be attribute to the H-bond between C74 and the oxygen atom of GUA which is replaced by a chlorine atom in 6GU or a hydrogen atom in 2BP (Figure 1C). The influence of mutation from C74 to U74. The ligands GUA, 6GU and 2BP have the same chemical scaffold. However, the mutation from C74 to U74 causes a significant decrease of the binding affinities for GUA complexes and a increase for

6GU and 2BP complexes. The

interaction spectra show that there is only one nucleotide 74 with large difference of contribution between GUA_GR and GUA_GR-C74U complexes. Hence, the selectivity of GR is mainly based on the interaction of nucleotide 74, this result is in consistent with Stock’s simulation result.19 As shown in Figure 6A, 6GU and key nucleotides overlapped very well between GR and GR-C74U complexes except the nucleobases of U74 and C74. This conformational difference dues to the reason that 6GU forms two H-bonds with different atoms of nucleobases of C74 and U74 (Figure 6A). In the GUA complexes, the relative position of nucleobases of nucleotide 74 in GR and GR-C74U is reverse compared to in 6GU complexes. GUA forms three H-bonds with C74 in GUA_GR complex, however, only one H-bond with occupancy at 73% was found in GUA_GR-C74U complex (Figure 6B). As shown in Figure 6B, the hydrogen atom at 1-position hampers more H-bond to form. Hence, it would be a good strategy keep this hydrogen atom in the future to design new ligand with high selectivity. The conformational change caused by ligand binding. To evaluate the conformational differences between ligand bound and unbound states, RMSDs of the heavy atoms of GR, heavy atoms of helix P1 and junction loops J23 (P1J23), heavy atoms of GR except the region P1J23 (Not P1J23) and heavy atoms of the ligand GUA were calculated relative to their crystallographic structures versus the MD simulations times (Figure S3). The figure S3 suggests that the RMSDs of the unbound state are larger than those of the bound state. The average RMSDs of GR over 3 MD trajectories of bound state and unbound states are 3.85 and 2.90 Å, respectively. It is clear 12


Page 13 of 34

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


that the conformation of unbound state is different from that of bound state. We noticed that one of the three trajectories (named MD3) of bound state shows a large fluctuation of RMSDs (Figure S3). This may indicate a large conformational change. From the Figure S3, it is clear that the RMSDs of GR are smaller than those of P1J23 and larger than those of Not P1J23 in most of MD simulations times. The region Not P1J23 shows an average value of RMSD of 2.37 Å for the bound state and 2.71 Å for the unbound state, while the region P1J23 shows an average value of RMSD of 3.74 Å for the bound state and 5.52 Å for the unbound state. This suggests that the regions P1J23 are more flexible than the region Not P1J23 in both states, and the major contributions toward increased unbound state RMSD values come from the region P1J23. We calculated the RMSFs from the experimental X-ray temperature factors ( B ) by equation RMSF = 3B / 8π 2 (Figure 7). The RMSFs of nucleotide in region P1J23 are larger than most of the other region. This also implies that the region P1J23 is flexible. The superimposition between two crystallographic structures PDB ID 3GAO and 1Y27 also suggests a large deviation at the P1J23 region, this is in agreement with the study of adenine riboswitch9. Figure 7 also shows the root mean square fluctuations (RMSFs) of C1’ atoms of the nucleotides averaged over the time from 200 ns to 1000 ns of three MD simulations. In general, the RMSFs of the bound state are qualitatively similar to those from the crystallographic structure. It can be seen from Figure 7 that most of nucleotides of P1 and J23 have corresponding lower RMSF values in the unbound state, implying that the region P1J23 would be less flexible in the unbound state than in bound state. This is in accordance with the fact that the GR in the absence of the ligand facilitates transcription via stabilization of the antiterminator stem9,25. To quantitatively evaluate the fluctuation of the region P1J23 during MD simulations, we introduced a parameter which is a distance between two mass center of the backbone atoms of nucleotides 15 and 74 (named D1 in Figure 8A). These distances in crystallographic structures of GUA (PDB ID 1Y27), 2BP(PDB ID 3G4M), 6GU(PDB ID 3GER) and XAN (PDB ID 3GAO) are 10.72, 8.10, 8.11 and 8.67 Å, respectively. The average distances over the trajectories from 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 34

200 ns to 1000 ns for three MD simulations are 14.34 and 17.58 Å for bound state and unbound state, respectively. We defined that the snapshot of MD simulation is in open conformation when its distance D1 greater than the average distance of all six MD simulations (15.96 Å), and it is in close conformation when less than. The percentages, which snapshots are from three MD simulation trajectories for 200 ns to 1000 ns that the specific close conformation exists, are 69.9% and 25.5% for the bound state and unbound state, respectively. This implies that the MD simulations of bound state are mostly in close conformation, and the unbound state in open conformation. The experimental and theoretical information suggested the conformational change is caused by the binding of ligand. Then we carried out a 100 ns MD simulation in which restrained all the atoms of the key nucleotides in binding pocket with a force constants 2 kcal/(mol·Å2) to evaluate the relationship of conformational change between key nucleotides of binding pocket and helix P1. In the restrained MD simulation, the average distance is only 11.80 Å, which is smaller than any one of three MD simulations of bound state for the first 100 ns MD simulation. So the conformational flexibility of the binding pocket would influence the helix P1’s conformation. In order to investigate the conformational differences of binding pocket between open and close conformations, we have calculated the average distances of mass center of backbone atoms of two key nucleotides in the ligand binding pocket and shown in Figure 8A. All the average distances except D1 together can reflect the conformational changes of the binding pocket. The average distances D2 and D3 in bound state are smaller than in unbound state, and the average distances D4 and D5 present opposite case. Little difference between two states for distances D6 and D7 was found. Figure 8B and C show the distances D1, D2 and D3 versus MD simulation times for the bound state of the MD simulations which may show a change of conformation (close or open). It is clear that the line of D1 is the same trend as that of D2 or D3, such as the MD simulation times from 0 ns to 150 ns and from 780 ns to 820 ns in Figure 8B, and from 200 ns to 350 ns in Figure 8C. Thus the conformational change would be caused by the distances D2 and D3. The conformation may be in close conformation when D2 and D3 become small, vice versa. 14


Page 15 of 34

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


Two nucleotides in P3 (U51 and A52) bound tightly to ligand with large free energies, as well as the other two nucleotides in P1 (C74 and U75). This implies that the ligand is a medium making P1 close to P3. The standard deviations (SDs) of D4 and D5 (0.50 and 0.30) for bound state are smaller than those for unbound state (1.29 and 0.81). In the bound state, the nucleotides A21 and U75 form two stable H-bonds with occupancy 99.52% and 99.36% (Figure 8D). In the unbound state, the nucleobases of nucleotides A21 and U75 can move freely for absent ligand. Then the H-bonds become weak with occupancy 55.54% and 55.03% (Figure 8D). The flexible nucleobase of U75 forms H-bond with U51. So the H-bond between U75 and A21 not only strengthen the pocket but also increase the distance D5. We didn’t find the distances D4 and D5 is related to conformational change. CONCLUSIONS The selectivity mechanism of ligands and the mechanism of the GR conformational transition have been studies by molecular dynamics (MD) simulation as well as free energy calculations by molecular mechanics Poisson Boltzmann surface area (MM-PBSA) and thermodynamic integration (TI). The results obtained by the two free-energy methods are consistent with each other. We further showed that the electrostatic energy play an important role for the selectivity of GR. The H-bond analysis indicated that the selectivity of ligand GUA is mainly from the decrease in the number of H-bonds in GR-C74U relative to in GR. The hydrogen atom at the 1position is favorable to form H-bonds with C74 of GR and unfavorable to form H-bonds with U74 of GR-C74U. However, two ligands 2BP and 6GU do not show selectively, because they form H-bonds by rearrangement of the nucleobase of U74. So it would be a good strategy to keep the hydrogen atom for designing better ligands in the future. The MD simulations of the ligand bound and unbound states of GR revealed that the unbound state is in the open conformation for the most MD time, and the bound state in the reverse case. The conformational difference is related with the distances D2 and D3. When the distances become smaller, the structure may be in close conformation.

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 34

ASSOCIATED CONTENT Supporting Information Available: Figure S1 shows the dV/dλ curves for transformations of the inhibitor and GR complexes in water from 12 λ points of TI calculations. Figure S2 shows the decomposition free energy of inhibitors on a per-nucleotide basis for six complexes. Figure S3 shows the RMSDs of heavy atoms for all the RNA, the region P1J23, all except the region P1J23 and the ligand GUA relative to their crystallographic structure. AUTHOR INFORMATION Corresponding Authors *E-mail: [email protected] (G.H.). [email protected] (J.W.) Notes The authors declare no competing financial interest. ACKNOWLEDGMENTS This work was partially supported by funding from the National Natural Science Foundation of China (11447004 and 61271378), the Natural Science Foundation of Shandong Province (ZR2014JL006), the Taishan Scholars Program of Shandong Province and a Project of Shandong Province Higher Educational Science and Technology Program (J14LJ05). REFERENCES 1

Thomas, J. R. Hergenrother, P. J. Targeting RNA with small molecules. Chem. Rev. 2008, 108, 1171-1224.

2

Daldrop, P., Reyes, F. E., Robinson, D. A., Hammond, C. M., Lilley, D. M., Batey, R. T.Brenk, R. Novel ligands for a purine riboswitch discovered by RNA-ligand docking. Chem. Biol. 2011, 18, 324-335.

3

Sund, J., Lind, C.Aqvist, J. Binding site preorganization and ligand discrimination in the purine riboswitch. J. Phys. Chem. B 2015, 119, 773-782.

4

Mandal, M. Breaker, R. R. Gene regulation by riboswitches. Nat. Rev. Mol. Cell Biol. 2004, 5, 451-463.

5

Soukup, J. K. Soukup, G. A. Riboswitches exert genetic control through metabolite-induced conformational change. Curr. Opin. Struct. Biol. 2004, 14, 344-349.

6

Montange, R. K. Batey, R. T. Riboswitches: emerging themes in RNA structure and function. Annu. Rev.

16


Page 17 of 34

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


Biophys. 2008, 37, 117-133. 7

Gilbert, S. D., Rambo, R. P., Van Tyne, D.Batey, R. T. Structure of the SAM-II riboswitch bound to Sadenosylmethionine. Nat. Struct. Mol. Biol. 2008, 15, 177-182.

8

Sharma, M., Bulusu, G.Mitra, A. MD simulations of ligand-bound and ligand-free aptamer: molecular level insights into the binding and switching mechanism of the add A-riboswitch. RNA 2009, 15, 1673-1692.

9

Priyakumar, U. D. MacKerell, A. D., Jr. Role of the adenine ligand on the stabilization of the secondary and tertiary interactions in the adenine riboswitch. J. Mol. Biol. 2010, 396, 1422-1438.

10

Deigan, K. E. Ferre-D'Amare, A. R. Riboswitches: discovery of drugs that target bacterial gene-regulatory RNAs. Acc. Chem. Res. 2011, 44, 1329-1338.

11

Nguyen, G. T., Scaife, M. A., Helliwell, K. E.Smith, A. G. Role of riboswitches in gene regulation and their potential for algal biotechnology. J. Phycol. 2016, 52, 320-328.

12

Sudarsan, N., Wickiser, J. K., Nakamura, S., Ebert, M. S.Breaker, R. R. An mRNA structure in bacteria that controls gene expression by binding lysine. Genes Dev. 2003, 17, 2688-2697.

13

Mandal, M., Lee, M., Barrick, J. E., Weinberg, Z., Emilsson, G. M., Ruzzo, W. L.Breaker, R. R. A glycinedependent riboswitch that uses cooperative binding to control gene expression. Science 2004, 306, 275-279.

14

Mandal, M., Boese, B., Barrick, J. E., Winkler, W. C.Breaker, R. R. Riboswitches control fundamental biochemical pathways in Bacillus subtilis and other bacteria. Cell 2003, 113, 577-586.

15

Mandal, M. Breaker, R. R. Adenine riboswitches and gene activation by disruption of a transcription terminator. Nat. Struct. Mol. Biol. 2004, 11, 29-35.

16

Serganov, A., Huang, L.Patel, D. J. Coenzyme recognition and gene regulation by a flavin mononucleotide riboswitch. Nature 2009, 458, 233-237.

17

Barrick, J. E. Breaker, R. R. The distributions, mechanisms, and structures of metabolite-binding riboswitches. Genome Biol. 2007, 8, R239.

18

Dann, C. E., 3rd, Wakeman, C. A., Sieling, C. L., Baker, S. C., Irnov, I.Winkler, W. C. Structure and mechanism of a metal-sensing regulatory RNA. Cell 2007, 130, 878-892.

19

Villa, A., Wohnert, J.Stock, G. Molecular dynamics simulation study of the binding of purine bases to the aptamer domain of the guanine sensing riboswitch. Nucleic Acids Res. 2009, 37, 4774-4786.

20

Batey, R. T., Gilbert, S. D.Montange, R. K. Structure of a natural guanine-responsive riboswitch complexed with the metabolite hypoxanthine. Nature 2004, 432, 411-415.

21

Serganov, A., Yuan, Y. R., Pikovskaya, O., Polonskaia, A., Malinina, L., Phan, A. T., Hobartner, C., Micura, R., Breaker, R. R.Patel, D. J. Structural basis for discriminative regulation of gene expression by adenineand guanine-sensing mRNAs. Chem. Biol. 2004, 11, 1729-1741.

22

Gilbert, S. D., Reyes, F. E., Edwards, A. L.Batey, R. T. Adaptive ligand binding by the purine riboswitch in the recognition of guanine and adenine analogs. Structure 2009, 17, 857-868.

23

Ling, B., Wang, Z., Zhang, R., Meng, X., Liu, Y., Zhang, C.Liu, C. Theoretical studies on the interaction of modified pyrimidines and purines with purine riboswitch. J. Mol. Graph. Model. 2009, 28, 37-45.

24

Ling, B., Zhang, R., Wang, Z., Dong, L., Liu, Y., Zhang, C.Liu, C. Theoretical studies on the interaction of guanine riboswitch with guanine and its closest analogues. Mol. Simulat. 2010, 36, 929-938.

25

Brenner, M. D., Scanlan, M. S., Nahas, M. K., Ha, T.Silverman, S. K. Multivector fluorescence analysis of the xpt guanine riboswitch aptamer domain and the conformational role of guanine. Biochemistry 2010, 49, 1596-1605.

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

26

Page 18 of 34

Gilbert, S. D., Stoddard, C. D., Wise, S. J.Batey, R. T. Thermodynamic and kinetic characterization of ligand binding to the purine riboswitch aptamer domain. J. Mol. Biol. 2006, 359, 754-768.

27

Buck, J., Furtig, B., Noeske, J., Wohnert, J.Schwalbe, H. Time-resolved NMR methods resolving ligandinduced RNA folding at atomic resolution. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 15699-15704.

28

Ogrizek, M., Konc, J., Bren, U., Hodoscek, M.Janezic, D. Role of magnesium ions in the reaction mechanism at the interface between Tm1631 protein and its DNA ligand. Chem. Cent. J. 2016, 10, 41.

29

Di Palma, F., Colizzi, F.Bussi, G. Ligand-induced stabilization of the aptamer terminal helix in the add adenine riboswitch. RNA 2013, 19, 1517-1524.

30

Allner, O., Nilsson, L.Villa, A. Loop-loop interaction in an adenine-sensing riboswitch: a molecular

31

Whitford, P. C., Schug, A., Saunders, J., Hennelly, S. P., Onuchic, J. N.Sanbonmatsu, K. Y. Nonlocal helix

dynamics study. RNA 2013, 19, 916-926. formation is key to understanding S-adenosylmethionine-1 riboswitch function. Biophys. J. 2009, 96, L7-9. 32

Fulle, S. Gohlke, H. Molecular recognition of RNA: challenges for modelling interactions and plasticity. J. Mol. Recognit. 2010, 23, 220-231.

33

Li, W., Wang, W.Takada, S. Energy landscape views for interplays among folding, binding, and allostery of calmodulin domains. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 10550-10555.

34

Shen, J., Tan, C., Zhang, Y., Li, X., Li, W., Huang, J., Shen, X.Tang, Y. Discovery of potent ligands for estrogen receptor beta by structure-based virtual screening. J. Med. Chem. 2010, 53, 5361-5365.

35

Hu, G. Wang, J. Ligand selectivity of estrogen receptors by a molecular dynamics study. Eur J Med Chem 2014, 74, 726-735.

36

Yang, Y., Shen, Y., Liu, H.Yao, X. Molecular dynamics simulation and free energy calculation studies of the binding mechanism of allosteric inhibitors with p38alpha MAP kinase. J. Chem. Inf. Model. 2011, 51, 32353246.

37

Gohlke, H., Kiel, C.Case, D. A. Insights into protein–protein binding by binding free energy calculation and free energy decomposition for the Ras–Raf and Ras–RalGDS complexes. J. Mol. Biol. 2003, 330, 891-913.

38

Deng, Y. Roux, B. Computations of standard binding free energies with molecular dynamics simulations. J. Phys. Chem. B 2009, 113, 2234-2246.

39

Roux, B., Nina, M., Pomès, R.Smith, J. C. Thermodynamic stability of water molecules in the bacteriorhodopsin proton channel: a molecular dynamics free energy perturbation study. Biophys. J. 1996, 71, 670-681.

40

Chen, J., Wang, J., Zhu, W.Li, G. A computational analysis of binding modes and conformation changes of MDM2 induced by p53 and inhibitor bindings. J Comput Aided Mol Des 2013, 27, 965-974.

41

Gouda, H., Kuntz, I. D., Case, D. A.Kollman, P. A. Free energy calculations for theophylline binding to an RNA aptamer: Comparison of MM-PBSA and thermodynamic integration methods. Biopolymers 2003, 68, 16-34.

42

Journal of Molecular BiologyKormos, B. L., Benitex, Y., Baranger, A. M.Beveridge, D. L. Affinity and specificity of protein U1A-RNA complex formation based on an additive component free energy model. J. Mol. Biol. 2007, 371, 1405-1419.

43

Reyes, C. M. Kollman, P. A. Structure and thermodynamics of RNA-protein binding: using molecular dynamics and free energy analyses to calculate the free energies of binding and conformational change. J. Mol. Biol. 2000, 297, 1145-1158.

18


Page 19 of 34

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

44


Kollman, P. A., Massova, I., Reyes, C., Kuhn, B., Huo, S., Chong, L., Lee, M., Lee, T., Duan, Y., Wang, W., Donini, O., Cieplak, P., Srinivasan, J., Case, D. A.Cheatham, T. E., 3rd. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33, 889-897.

45

Hu, G., Ma, A., Dou, X., Zhao, L.Wang, J. Computational studies of a mechanism for binding and drug resistance in the wild type and four mutations of HIV-1 protease with a GRL-0519 inhibitor. Int. J. Mol. Sci. 2016, 17, 819.

46

Stoica, I., Sadiq, S. K.Coveney, P. V. Rapid and accurate prediction of binding free energies for saquinavirbound HIV-1 proteases. J. Am. Chem. Soc. 2008, 130, 2639-2648.

47

Zgarbova, M., Otyepka, M., Sponer, J., Mladek, A., Banas, P., Cheatham, T. E., 3rdJurecka, P. Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput. 2011, 7, 2886-2902.

48

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., Montgomery, J. A., Jr.; Vreven, T., Kudin, K. N., Burant, J. C., Millam, J. M., Iyengar, S. S., Tomasi, J., Barone, V., Mennucci, B., Cossi, M., Scalmani, G., Rega, N., Petersson, G. A., Nakatsuji, H., Hada, M., Ehara, M., Toyota, K., Fukuda, R., Hasegawa, J., Ishida, M., Nakajima, T., Honda, Y., Kitao, O., Nakai, H., Klene, M., Li, X., Knox, J. E., Hratchian, H. P., Cross, J. B., Bakken, V., Adamo, C., Jaramillo, J., Gomperts, R., Stratmann, R. E., Yazyev, O., Austin, A. J., Cammi, R., Pomelli, C., Ochterski, J. W., Ayala, P. Y., Morokuma, K., Voth, G. A., Salvador, P., Dannenberg, J. J., Zakrzewski, V. G., Dapprich, S., Daniels, A. D., Strain, M. C., Farkas, O., Malick, D. K., Rabuck, A. D., Raghavachari, K., Foresman, J. B., Ortiz, J. V., Cui, Q., Baboul, A. G., Clifford, S., Cioslowski, J., Stefanov, B. B., Liu, G., Liashenko, A., Piskorz, P., Komaromi, I., Martin, R. L., Fox, D. J., Keith, T., Al-Laham, M. A., Peng, C. Y., Nanayakkara, A., Challacombe, M., Gill, P. M. W., Johnson, B., Chen, W., Wong, M. W., Gonzalez, C.Pople, J. A. Gaussian 03, Revision C.02, . Gaussian, Inc., Wallingford CT, 2004.

49

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.

50

Bayly, C. I., Cieplak, P., Cornell, W.Kollman, P. A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 1993, 97, 10269-10280.

51

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W.Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Comput. Phys. 1983, 79, 926-935.

52

Case, D. A., Darden, T. A., Cheatham, T. E., III, Simmerling, C. L., Wang, J., Duke, R. E., Luo, R., Walker, R. C., Zhang, W., Merz, K. M., Roberts, S., Hayik, S., Roitberg, A., Seabra, G., Swails, J., Götz, A. W., Kolossváry, I., Wong, K. F., Paesani, F., Vanicek, J., Wolf, R. M., Liu, J. W., X., Brozell, S. R., Steinbrecher, T., Gohlke, H., Cai, Q., Ye, X., Wang, J., Hsieh, M.-J., Cui, G., Roe, D. R., Mathews, D. H., Seetin, M. G., Salomon-Ferrer, R., Sagui, C., Babin, V., Luchko, T., Gusarov, S., Kovalenko, A.Kollman, P. A. Amber 12, University of California: San Francisco, 2012.

53

Darden, T., York, D.Pedersen, L. Particle mesh ewald -and.Log(N) method for ewald sums in large systems. J. Comput. Phys. 1993, 98, 10089-10092.

54

Ryckaert, J. P., Ciccotti, G.Berendsen, H. J. C. Numerical-integration of cartesian equations of motion of a

55

Chen, J., Liang, Z., Wang, W., Yi, C., Zhang, S.Zhang, Q. Revealing origin of decrease in potency of

system with constraints - molecular-dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327-341.

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 34

darunavir and amprenavir against HIV-2 relative to HIV-1 protease by molecular dynamics simulations. Sci. Rep. 2014, 4, 6872. 56

Hu, G., Cao, Z., Xu, S., Wang, W.Wang, J. Revealing the binding modes and the unbinding of 14-3-3sigma proteins and inhibitors by computational methods. Scientific reports 2015, 5, 16481.

57

Weiser, J., Shenkin, P. S.Still, W. C. Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 1999, 20, 217-230.

58

Michel, F. S., Arthur, J. O.Jean-Claude, S. Reduced surface: An efficient way to compute molecular surfaces. Biopolymers 1996, 38, 305-320.

59

Onufriev, A., Bashford, D.Case, D. A. Exploring protein native states and large-scale conformational

60

Case, D. A. Normal mode analysis of protein dynamics. Curr. Opin. Struct. Biol. 1994, 4, 285-290.

61

Genheden, S., Nilsson, I.Ryde, U. Binding affinities of factor Xa inhibitors estimated by thermodynamic

changes with a modified generalized born model. Proteins 2004, 55, 383-394.

integration and MM/GBSA. J. Chem. Inf. Model. 2011, 51, 947-958. 62

Leonis, G., Steinbrecher, T.Papadopoulos, M. G. A contribution to the drug resistance mechanism of darunavir, amprenavir, indinavir, and saquinavir complexes with HIV-1 protease due to flap mutation I50V: a systematic MM-PBSA and thermodynamic integration study. J. Chem. Inf. Model. 2013, 53, 2141-2153.

63

Cai, Y. Schiffer, C. A. Decomposing the energetic impact of drug resistant mutations in HIV-1 protease on binding DRV. J. Chem. Theory Comput. 2010, 6, 1358-1368.

64

Tzoupis, H., Leonis, G., Mavromoustakos, T.Papadopoulos, M. G. A comparative molecular dynamics, MM-PBSA and thermodynamic integration study of saquinavir complexes with Wild-Type HIV-1 PR and L10I, G48V, L63P, A71V, G73S, V82A and I84V single mutants. J. Chem. Theory Comput. 2013, 9, 17541764.

65

Chen, J., Wang, X., Zhu, T., Zhang, Q.Zhang, J. Z. A Comparative Insight into Amprenavir Resistance of Mutations V32I, G48V, I50V, I54V, and I84V in HIV-1 Protease Based on Thermodynamic Integration and MM-PBSA Methods. J. Chem. Inf. Model. 2015, 55, 1903-1913.

66

Beutler, T. C., Mark, A. E., van Schaik, R. C., Gerber, P. R.van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222, 529-539.

67

Bommarito, S., Peyret, N.SantaLucia, J., Jr. Thermodynamic parameters for DNA sequences with dangling

68

Mignon, P., Loverix, S., Steyaert, J.Geerlings, P. Influence of the pi-pi interaction on the hydrogen bonding

ends. Nucleic Acids Res. 2000, 28, 1929-1934. capacity of stacked DNA/RNA bases. Nucleic Acids Res. 2005, 33, 1779-1789.

20


Page 21 of 34

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 captions Figure 1

(A) The structure of the aptamer region of guanine riboswitch and its inhibitor in binding

pocket. Riboswitch is shown in cartoon representation colored with each secondary structure. Ligand GUA is show in stick representation as well as dots representation colored with element of atom; (B) Hbonds formed between Ligand GUA and three nucleic acids; (C) The structures of four ligands.

Figure 2

The thermodynamic cycle used to compute the free energy differences. When A is the GR, B1

and B2 are GUA and a ligand, respectively. When A is a ligand, B1 and B2 are GR and GR-C74U, respectively.

Figure 3

The decomposition of inhibitors on a per-nucleotide basis for GUA and GR/GR-C74U

complexes.

Figure 4

(A) The H-bonds formed between ligands and nucleotide C74. The ligand and nucleotide C74

are shown in stick and ball representation for GUA-GR complex, as well as in stick representation for 6GU-GR complex. The other key nucleotides of binding pockets are shown in line representation. The H-bonds are shown in lines and labeled with average distances. (B) The position of the nucleotides interacted with GUA by the aromatic stacking of nucleobases. The GUA and the nucleobases are shown in ball and stick representation.

Figure 5

The superimposition of binding pockets of 6GU and 2BP. The nucleotides and ligand in 2BP

complex are shown in stick and ball representation, as well as stick representation for 6GU complex. The distances of two mass centers and the H-bonds are labeled.

Figure 6

The superimposition of binding pockets of GR and GR-C74U for (A) 6GU and (B) GUA. The

nucleotide 74 and ligand in GR complexes are shown in stick and ball representation, as well as stick representation in GR-C74U complexes. The other key nucleotides are shown in lines. The distances of H-bonds are labeled.

Figure 7

RMSF for the C1’ atoms of the nucleotides. The fluctuation data from the B-factor are also

given. The regions of guanine riboswitch were labeled.

Figure 8

(A) Schematic view of the average distances between the mass centers of the backbone atoms

of two nucleotides, D1 is for G15 and C74, D2 is for A52 and C74, D3 is for U51 and U75, D4 is for U22 and C74, D5 is for A21 and U75, D6 is for U22 and A52, D7 is for A21 and U51. The distance values are shown for bound state and unbound state (in brackets); The D1, D2 and D3 versus the MD simulation time for MD1 (B) and MD3 (C) of bound state; (D) The superimposition of six key nucleotides of two conformations from the MD simulation of unbound state, only showing three nuleotides and the H-bonds formed by U75. The occupancy of H-bonds is shown in pale blue color for 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 34

one conformation and in orange color for the other conformation. Binding pocket is shown in dots. Atoms are presented in stick and ball, as well as in stick representation. .

22


Page 23 of 34

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


Table of Contents Image

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 34

Tables Table 1.

Binding free energies calculated for GR and GR-C74U complexes obtained by

MM-PBSA method a Complexesb

6GU

GUA

Items

GR

GR-C74U

GR

∆E ele

-39.27±0.93

-25.83±0.97

-29.13±0.5

∆E vdw

-26.84±0.28

-27.85±0.27

∆G nonpol

-2.91±0.01

∆G pol

2BP GR

GR-C74U

-30.96±0.76 -29.54±0.41 -32.46±0.46

-34.6±1.96

-34.42±2.04

-27.96±0.3

-28.81±0.17 -25.48±0.21

-25.7±0.31

-27.28±0.22

-26.6±0.41

-2.97±0.01

-3.02±0.02

-3.02±0.01

-2.86±0.01

-2.85±0.01

-2.88±0

-2.91±0.01

38.33±0.76

33.61±1.06

33.51±0.69

34.24±0.76

33.47±0.43

36.17±0.63

40.59±1.64

37.6±1.49

∆G ele+pol

-0.94±0.44

7.78±1.03

4.38±0.51

3.28±0.4

3.93±0.44

3.71±0.36

5.99±0.55

3.18±1.03

∆G vdw+nonpol

-29.75±0.28

-30.82±0.27

-30.98±0.32 -31.83±0.18 -28.34±0.21 -28.55±0.31

-30.16±0.22

-29.51±0.4

∆H

-30.69±0.52

-23.04±1.06

-26.6±0.77

-28.55±0.48 -24.41±0.45 -24.84±0.51

-24.17±0.41

-26.33±1

-T∆S

12.71±0.4

12.69±0.19

12.02±0.22

12.12±0.24

11.91±0.13

11.59±0.19

12.38±0.08

12.22±0.28

∆G bind

-17.98±0.87

-10.35±1.13

-14.58±0.67

-16.43±0.5

-12.5±0.47

-13.25±0.65

-11.79±0.41

-14.11±0.93

-11.53

-

-8.31

-8.43

-7.35

-9.06

-6.05

-

∆G exp

c

GR-C74U

GR

XAN GR-C74U

a

All values are given in kcal/mol. Standard errors were labeled by the ± signs and calculated by σ= standard diviation/N1/2, where N is ten because we have done ten independent MD simulations for each complex. b The symbols of the energy terms are described in the section of the binding free energy calculations. c The experimental values from refs 22 and 26.


Page 25 of 34

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


Table 2.

Relative binding free energies for GR and GUA complex relative to the other

ligands calculated by thermodynamic integration a GR

∆∆G step1

GUA→6GU

0.66±5.17

GUA→2BP GUA→XAN a

∆∆G step2

∆∆G step3

∆∆G TI

∆∆G exp

∆∆G MM-PBSA

0.29±4.18

0.38±0.52

1.33

3.21

3.4

0.38±5.21

2.78±4.47

-1.36±1.37

1.80

4.12

5.48

1.90±2.77

-1.98±3.98

5.19±5.25

5.11

5.47

6.19

All values are in kcal/mol. Errors labeled by the ± signs represent standard errors.

Table 3.

Relative binding free energies for GR and four ligands complexes relative to their

C74U mutations calculated by thermodynamic integration a

a

GR→GR-C74U

∆∆G step1

GUA

0.25±2.41

6GU

∆∆G step2

∆∆G step3

∆∆G TI

1.16±3.25

0.8±4.72

2.21

7.63

-0.2±2.47

-1.14±3.24

-0.21±4.35

-1.55

-1.85

2BP

0.56±2.47

-1.29±3.21

-0.84±4.48

-1.57

-0.75

XAN

-0.4±2.44

-1.15±3.20

-0.66±4.39

-2.21

-2.32

∆∆G MM-PBSA

All values are in kcal/mol. Errors labeled by the ± signs represent standard errors.



Page 26 of 34

1 2 3 4 5 6 Table 4. Main H-bonds involved in ligands bindinga 7 8 GUA 6GU 2BP XAN 9 Items GR GR-C74U GR GR-C74U GR GR-C74U GR GR-C74U 10 Occu Dis Occu Dis Occu Dis Occu Dis Occu Dis Occu Dis Occu Dis Occu Dis 11 12 [email protected]@H3-N3 99.98 2.88 99.84 2.90 99.89 2.90 99.8 2.89 99.93 2.90 99.86 2.89 13 U51@O4... Lig@ H2-N1 93.31 3.05 93.76 3.02 93.84 3.04 91.55 3.06 91.82 3.05 91.21 3.06 14 88.22 2.97 95.53 3.02 89.11 2.86 94.31 3.04 15 U51@O2... Lig@ H3-N4 16 C74@N3... Lig@ H5-N3 98.74 2.96 98.34 3.09 99.02 3.01 17 Lig@O1... C74@H41-N4 97.16 3.01 97.28 2.93 18 89.68 2.85 86.47 2.96 19 C74@O2... Lig@ H4-N4 20Lig@ N3... C74@H41-N4 99.65 3.00 98.8 3.00 21 C74@N3... Lig@ H3-N4 94.15 3.01 95.74 3.01 22 99.72 2.83 99.62 2.83 96.88 2.89 98.74 2.89 99.9 2.82 99.61 2.84 98.29 2.84 88.56 2.84 [email protected]@HO2`- O2` 24 U74@O4... Lig@ H5-N3 73.43 2.93 25 [email protected]@H2-N5... 97.92 2.83 99.92 2.86 26 68.54 2.85 27 U47@O4... Lig@ H2-N1 28 [email protected]@ H3-N3 93.53 2.99 85.38 2.92 29 a The H-bonds are determined by an acceptor···donor distance (Dis) less than 3.5 Å and an 30 31 acceptor···H donor angle larger than 120°. The occupancy (Occu) is defined as the percentage of 32 the snapshots from MD simulation trajectory that a specific hydrogen bond exists. 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 26 59 60 ACS Paragon Plus Environment

Page 27 of 34

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 1 (A) The structure of the aptamer region of guanine riboswitch and its inhibitor in binding pocket. Riboswitch is shown in cartoon representation colored with each secondary structure. Ligand GUA is show in stick representation as well as dots representation colored with element of atom; (B) H-bonds formed between Ligand GUA and three nucleic acids; (C) The structures of four ligands. 632x384mm (72 x 72 DPI)



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 2 The thermodynamic cycle used to compute the free energy differences. When A is the GR, B1 and B2 are GUA and a ligand, respectively. When A is a ligand, B1 and B2 are GR and GR-C74U, respectively. 349x120mm (96 x 96 DPI)


Page 28 of 34

Page 29 of 34

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 3 The decomposition of inhibitors on a per-nucleotide basis for GUA and GR\GR-C74U complexes. 233x357mm (150 x 150 DPI)



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 4 (A) The H-bonds formed between ligands and nucleotide C74. The ligand and nucleotide C74 are shown in stick and ball representation for GUA-GR complex, as well as in stick representation for 6GU-GR complex. The other key nucleotides of binding pockets are shown in line representation. The H-bonds are shown in lines and labeled with average distances. (B) The position of the nucleotides interacted with GUA by the aromatic stacking of nucleobases. The GUA and the nucleobases are shown in ball and stick representation. 282x436mm (72 x 72 DPI)


Page 30 of 34

Page 31 of 34

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 5 The superimposition of binding pockets of 6GU and 2BP. The nucleotides and ligand in 2BP complex are shown in stick and ball representation, as well as stick representation for 6GU complex. The distances of two mass centers and the H-bonds are labeled. 379x220mm (72 x 72 DPI)



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 6 The superimposition of binding pockets of GR and GR-C74U for (A) 6GU and (B) GUA. The nucleotide 74 and ligand in GR complexes are shown in stick and ball representation, as well as stick representation in GR-C74U complexes. The other key nucleotides are shown in lines. The distances of Hbonds are labeled. 299x400mm (72 x 72 DPI)


Page 32 of 34

Page 33 of 34

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 7 RMSF for the C1’ atoms of the nucleotides. The fluctuation data from the B-factor are also given. The regions of guanine riboswitch were labeled. 289x202mm (150 x 150 DPI)



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) Schematic view of the average distances between the mass centers of the backbone atoms of two nucleotides, D1 is for G15 and C74, D2 is for A52 and C74, D3 is for U51 and U75, D4 is for U22 and C74, D5 is for A21 and U75, D6 is for U22 and A52, D7 is for A21 and U51. The distance values are shown for bound state and unbound state (in brackets); The D1, D2 and D3 versus the MD simulation time for MD1 (B) and MD3 (C) of bound state; (D) The superimposition of six key nucleotides of two conformations from the MD simulation of unbound state, only showing three nuleotides and the H-bonds formed by U75. The occupancy of H-bonds is shown in pale blue color for one conformation and in orange color for the other conformation. Binding pocket is shown in dots. Atoms are presented in stick and ball, as well as in stick representation. 565x424mm (72 x 72 DPI)


Page 34 of 34

Ligand Selectivity Mechanism and Conformational ... - ACS Publications

Recommend Documents