Explicit Solvent Molecular Dynamics Simulations Research Articles

Enhanced Binding of SARS-CoV-2 Spike Protein to Receptor by Distal Polybasic Cleavage Sites.

The receptor-binding domain (RBD) of the SARS-CoV-2 spike protein plays a crucial role in binding the human cell receptor ACE2 that is required for viral entry. Many studies have been conducted to target the structures of RBD–ACE2 binding and to design RBD-targeting vaccines and drugs. Nevertheless, mutations distal from the SARS-CoV-2 RBD also impact its transmissibility and antibody can target non-RBD regions, suggesting the incomplete role of the RBD region in the spike protein–ACE2 binding. Here, in order to elucidate distant binding mechanisms, we analyze complexes of ACE2 with the wild-type spike protein and with key mutants via large-scale all-atom explicit solvent molecular dynamics simulations. We find that though distributed approximately 10 nm away from the RBD, the SARS-CoV-2 polybasic cleavage sites enhance, via electrostatic interactions and hydration, the RBD–ACE2 binding affinity. A negatively charged tetrapeptide (GluGluLeuGlu) is then designed to neutralize the positively charged arginine on the polybasic cleavage sites. We find that the tetrapeptide GluGluLeuGlu binds to one of the three polybasic cleavage sites of the SARS-CoV-2 spike protein lessening by 34% the RBD–ACE2 binding strength. This significant binding energy reduction demonstrates the feasibility to neutralize RBD–ACE2 binding by targeting this specific polybasic cleavage site. Our work enhances understanding of the binding mechanism of SARS-CoV-2 to ACE2, which may aid the design of therapeutics for COVID-19 infection.

Characterizing the water wire in the Gramicidin channel found by Monte Carlo sampling using continuum electrostatics and in molecular dynamics trajectories with conventional or polarizable force fields

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Characterizing the Water Wire in the Gramicidin Channel Found by Monte Carlo Sampling Using Continuum Electrostatics and in Molecular Dynamics Trajectories with Conventional or Polarizable Force Fields

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Conformational distributions of isolated myosin motor domains encode their mechanochemical properties.

Myosin motor domains perform an extraordinary diversity of biological functions despite sharing a common mechanochemical cycle. Motors are adapted to their function, in part, by tuning the thermodynamics and kinetics of steps in this cycle. However, it remains unclear how sequence encodes these differences, since biochemically distinct motors often have nearly indistinguishable crystal structures. We hypothesized that sequences produce distinct biochemical phenotypes by modulating the relative probabilities of an ensemble of conformations primed for different functional roles. To test this hypothesis, we modeled the distribution of conformations for 12 myosin motor domains by building Markov state models (MSMs) from an unprecedented two milliseconds of all-atom, explicit-solvent molecular dynamics simulations. Comparing motors reveals shifts in the balance between nucleotide-favorable and nucleotide-unfavorable P-loop conformations that predict experimentally measured duty ratios and ADP release rates better than sequence or individual structures. This result demonstrates the power of an ensemble perspective for interrogating sequence-function relationships.

Conformational properties of inulin, levan and arabinan studied by molecular dynamics simulations

Inulin, levan and arabinan are the polysaccharides that consist of exclusively furanose units. To date, their conformational features studied at the molecular scale have remained largely unexplained. To tackle this issue, we have performed a series of explicit-solvent molecular dynamics simulations, carried out within the furanose-dedicated force field. None of the polysaccharides exhibits a single, dominating structure type. Instead, they create a large number of separated conformational states originating from the intensive rotation around the φ and ω glycosidic angles. 21-helices are the preferential conformational forms for all compounds but they appear only locally, at the length of several consecutive residues. The flexibility of all three furanose-based polysaccharides is much greater in relation to the (1–4)-linked pyranose polysaccharides and is comparable to that of (1–6)-linked pyranoses. The dynamic geometries of both furanose rings and glycosidic linkages are nearly unchanged independently if considering them at the level of mono-, di- or polysaccharides.

Efficient Construction of Atomic-Resolution Models of Non-Sulfated Chondroitin Glycosaminoglycan Using Molecular Dynamics Data.

Glycosaminoglycans (GAGs) are linear, structurally diverse, conformationally complex carbohydrate polymers that may contain up to 200 monosaccharides. These characteristics present a challenge for studying GAG conformational thermodynamics at atomic resolution using existing experimental methods. Molecular dynamics (MD) simulations can overcome this challenge but are only feasible for short GAG polymers. To address this problem, we developed an algorithm that applies all conformational parameters contributing to GAG backbone flexibility (i.e., bond lengths, bond angles, and dihedral angles) from unbiased all-atom explicit-solvent MD simulations of short GAG polymers to rapidly construct models of GAGs of arbitrary length. The algorithm was used to generate non-sulfated chondroitin 10- and 20-mer ensembles which were compared to MD-generated ensembles for internal validation. End-to-end distance distributions in constructed and MD-generated ensembles have minimal differences, suggesting that our algorithm produces conformational ensembles that mimic the backbone flexibility seen in simulation. Non-sulfated chondroitin 100- and 200-mer ensembles were constructed within a day, demonstrating the efficiency of the algorithm and reduction in time and computational cost compared to simulation.

Self-Assembly of Allomelanin Dimers and the Impact of Poly(ethylene glycol) on the Assembly: A Molecular Dynamics Simulation Study.

Synthetic melanin nanoparticles that exhibit properties analogous to naturally found allomelanin can be formed by assembly of dimers/oligomers of the synthetic precursor of allomelanin, 1,8-dihydroxy naphthalene (DHN). To link the nanostructure within these assembled melanin nanoparticles to DHN dimer structure, we use explicit-solvent atomistic molecular dynamics (MD) simulations to study assembly of DHN dimers (2-2', 2-4', and 4-4' and their mixture) into nanoparticles in aqueous solutions. We analyze how the dimer structure and mixture composition impact the molecular interactions that drive assembly, as well as the assembled nanostructure, both internally and on the surface. We find that, prominently, hydrogen-bonding interactions drive the assembly of like-dimers, whereas unlike-dimer stacking interactions play a role in the assembly of dimer mixtures. The aggregate/nanoparticle assembled from 2-2' dimers assumes a spherical morphology as opposed to 4-4' dimers that adopt an anisotropic shape. The surface of the aggregate formed by 2-2' dimers is primarily hydrophobic, while the surface of aggregates formed by 2-4' dimers, 4-4' dimers, and their mixtures is amphiphilic. We also find that the addition of linear poly(ethylene glycol) (PEG) chains to the assembled particle does not alter the aggregate structure formed by a single dimer type. However, the PEG chains prefer to interact with 4-4' dimers more than with 2-2' and 2-4'. In aggregates formed by mixtures of dimers, PEG chains interact preferentially with 2-4' dimers than with 4-4' and 2-2' dimers.

Correlation Length in Concentrated Electrolytes: Insights from All-Atom Molecular Dynamics Simulations.

We study the correlation length of the charge-charge pair correlations in concentrated electrolyte solutions by means of all-atom, explicit-solvent molecular dynamics simulations. We investigate LiCl and NaI in water, which constitute highly soluble, prototypical salts for experiments, as well as two more complex, molecular electrolyte systems of lithium bis(trifluoromethane)sulfonimide (LiTFSI), a salt commonly employed in electrochemical storage systems, in water, and in an organic solvent mixture of dimethoxyethane and dioxolane. Our simulations support the recent experimental observations as well as theoretical predictions of a nonmonotonic behavior of the correlation length with increasing salt concentration. We observe a Debye-Hückel like regime at low concentration, followed by a minimum reached when d/λD ≃ 1, where λD is the Debye correlation length and d is the effective ionic diameter, and an increasing correlation length with salt concentration in very concentrated electrolytes. As in the experiments, we find that the screening length in the concentrated regime follows a universal scaling law as a function d/λD for all studied salts. However, the scaling exponent is significantly lower than the experimentally measured one and lies in the range of the theoretical predictions based on much simpler electrolyte models.

Relaxation Time Spectrum and Dynamics of Stretched Polymer Chain in Dilute θ Solution: Implicit Solvent Model versus Explicit Solvent Model

AbstractRelaxation dynamics of a single stretched polymer chain after releasing the two fixed ends in dilute solutions are investigated by using implicit solvent and explicit solvent molecular dynamic simulations. The phase diagram of polymer in explicit solvent solution is determined, and the θ point is estimated by the form factor of the polymer chain. The relaxation time spectrum is obtained by the inverse Laplace transform of length–time decay curves of the relaxation process. The dependencies of the longest relaxation time in the spectrum on the chain length are compared with the scaling prediction of the Rouse model and the Zimm model of flexible chains. In this paper, it is confirmed that the dynamical scaling analysis of theoretical models are in good agreement with the simulation models in all range of polymerization degree of chain studied. The coarse‐grained explicit solvent simulation model of polymer chain in dilute solution developed in this paper can be a good numerical model for the future study of the dynamics process of macromolecules with hydrodynamic interactions in solutions.

Elongation Factor Tu Switch I Element is a Gate for Aminoacyl-tRNA Selection

Selection of correct aminoacyl (aa)-tRNA at the ribosomal A site is fundamental to maintaining translational fidelity. Aa-tRNA selection is a multistep process facilitated by the guanosine triphosphatase elongation factor (EF)-Tu. EF-Tu delivers aa-tRNA to the ribosomal A site and participates in tRNA selection. The structural mechanism of how EF-Tu is involved in proofreading remains to be fully resolved. Here, we provide evidence that switch I of EF-Tu facilitates EF-Tu's involvement during aa-tRNA selection. Using structure-based and explicit solvent molecular dynamics simulations based on recent cryo-electron microscopy reconstructions, we studied the conformational change of EF-Tu from the guanosine triphosphate to guanine diphosphate conformation during aa-tRNA accommodation. Switch I of EF-Tu rapidly converts from an α-helix into a β-hairpin and moves to interact with the acceptor stem of the aa-tRNA. In doing so, switch I gates the movement of the aa-tRNA during accommodation through steric interactions with the acceptor stem. Pharmacological inhibition of the aa-tRNA accommodation pathway prevents the proper positioning of switch I with the aa-tRNA acceptor stem, suggesting that the observed interactions are specific for cognate aa-tRNA substrates, and thus capable of contributing to the fidelity mechanism.

Going Big: Million Atom Simulations of Ribosomes and Billion Atom Simulations of Chromatin

The rapid increase of computing power has enabled atomistic simulations of larger and larger systems. High performance computing (HPC) platforms with more than one million processor cores are not uncommon among the world's fastest supercomputers. Efficient utilization of these machines requires careful attention to parallel decomposition of the problem. Explicit solvent molecular dynamics (MD) simulations of large macromolecular complexes play an important role in integrating disparate forms of experimental data into a single coherent picture. These simulations help set the stage for reduced model and coarse-grained (CG) simulations capable of simulating large-scale conformational changes important for mechanistic understanding. We used explicit solvent MD simulations to identify the accommodation corridor in the ribosome, critical for tRNA selection during protein synthesis (Sanbonmatsu, et al., PNAS, 2005). Microsecond explicit solvent simulations of the ribosome (2.2 million atoms) also laid the foundations for our energy landscape calculations using all-atom structure-based simulations of spontaneous accommodation events (Whitford, et al., PLoS Comput. Biol., 2013; Whitford, et al., RNA, 2010). We are applying a similar strategy to chromatin architecture, which plays a key role in embryo development, brain function and cancer. As a first step, we have performed the first explicit solvent simulation of an entire gene locus (GATA4), consisting of 427 nucleosomes and over one billion atoms (the first published billion atom biomolecular simulation) (Jung, et al., J. Comp. Chem. 2019). Simulations used the GENESIS molecular dynamics code on the LANL Trinity supercomputer. The multi-disciplinary effort combined computer science, high performance computing, chip design, biophysics, structural biology, and cell biology, including researchers from RIKEN, LANL, NYU, Intel and Cray. Several performance optimizations for the KNL architecture enabled scaling to large numbers of cores.

Comparing Stabilization Strategies Between Engineered and Naturally Thermostable Proteins

Bioengineers strive to create highly stable and functional proteins for use in biomedical and industrial purposes such as biosensors and pharmaceuticals. In order to better design thermostable proteins de novo, we investigated what features and properties contributed to thermostability in both engineered and naturally occurring thermostable proteins. UMC is a de novo-designed protein that is closely related to the thermostable protein UVF (Tm >99°C), which is an engineered thermostable variant of the Engrailed Homeodomain (EnHD) designed by the Mayo group. The Engrailed Homeodomain (EnHD) is a three-helix bundle transcription factor naturally found in fruit flies. We also studied consensus HD (conHD), a consensus design based on the homeodomain family of proteins that is more thermostable than EnHD. In order to compare stabilization strategies between engineered and naturally thermostable proteins, we compared these proteins to the naturally thermostable protein Cbp2Hb. Cbp2Hb is a six-helix protein found in Hyperthermus butylicus, an archaea that lives in temperatures up to 108 °C. We split the structure of Cbp2Hb into its two three-helix-bundle domains, Cbp2HbN, containing the N-terminal bundle, and Cbp2HbC, containing the C-terminal bundle. We performed all-atom, explicit-solvent molecular dynamic simulations at 25 °C and 100 °C on these four thermostable proteins as well as EnHD and used in lucem molecular mechanics to analyze their dynamics. We analyzed the pairwise residue interactions, solvent accessible surface area, and secondary structures in all five proteins to determine the structural and dynamic basis of the proteins’ thermostability.

Interpreting SAXS/WAXS Data with Explicit-Solvent Simulations: A Practical Guide.

Small- and wide-angle X-ray scattering (SAXS/WAXS/SWAXS) have evolved to be accurate tools used to gain structural information of biomolecules in solution. However, the interpretation of SWAXS data remains challenging owing to the low information content of the data and scattering contributions from the solvent. In recent years, methods for the interpretation of SWAXS data based on explicit-solvent molecular dynamics (MD) simulations have become increasingly popular. The physicochemical information in the MD force fields complements the low-information SWAXS data, thereby greatly reducing the risk of overfitting, and the explicit-solvent models may accurately account for scattering contributions from the solvent. In this chapter, we provide a practical introduction to MD-based methods for the interpretation of SWAXS data. First, we present the back-calculation of a SWAXS curve from an MD trajectory as required to validate an MD simulation against experimental SWAXS data. Second, we present the structure refinement of an atomic model against SWAXS data using SAXS-driven MD simulations. Common technical problems together with appropriate solutions are discussed.

Molecular dynamics simulations suggest stabilizing mutations in a de novo designed α/β protein.

Designing functional proteins that can withstand extreme heat is beneficial for industrial and protein therapeutic applications. Thus, elucidating the atomic-level determinants of thermostability is a major interest for rational protein design. To that end, we compared the structure and dynamics of a set of previously designed, thermostable proteins based on the activation domain of human procarboxypeptidase A2 (AYEwt). The mutations in these designed proteins were intended to increase hydrophobic core packing and inter-secondary-structure interactions. To evaluate whether these design strategies were successfully deployed, we performed all-atom, explicit-solvent molecular dynamics (MD) simulations of AYEwt and three designed variants at both 25 and 100°C. Our MD simulations agreed with the relative experimental stabilities of the designs based on their secondary structure content, Cα root-mean-square deviation/fluctuation, and buried-residue solvent accessible surface area. Using a contact analysis, we found that the designs stabilize inter-secondary structure interactions and buried hydrophobic surface area, as intended. Based on our analysis, we designed three additional variants to test the role of helix stabilization, core packing, and a Phe → Met mutation on thermostability. We performed the additional MD simulations and analysis on these variants, and these data supported our predictions.

Implicit Micelle Model for Membrane Proteins Using Superellipsoid Approximation.

Surfactant micelles are often utilized as membrane mimetics for structure determination and functional analysis of membrane proteins. The curved-surface effects of the micelle can perturb membrane protein structure. However, it is difficult to assess such effects and membrane mimetic artifacts by experimental and theoretical methods. Here, we propose an implicit micelle model (IMIC) to be used in molecular dynamics (MD) simulations of membrane proteins. IMIC is an extension of the IMM1 implicit membrane model and additionally introduces a superellipsoid approximation to represent the curved-surface effects. Most of the IMIC parameters are obtained from all-atom explicit solvent MD simulations of 12 membrane proteins in various micelles. The HIV envelope protein gp41, M13 major coat protein gp8, and amyloid precursor protein (APP) dimer are simulated via MD simulations with IMIC. These simulations clearly show how the micelle influences membrane protein structures compared to the bilayer environments. The MD simulations with IMIC provide reliable membrane protein structures in various micelle environments quickly with smaller computational cost than that for an explicit solvent/micelle model.

Ff19SB: Amino-Acid-Specific Protein Backbone Parameters Trained against Quantum Mechanics Energy Surfaces in Solution.

Molecular dynamics (MD) simulations have become increasingly popular in studying the motions and functions of biomolecules. The accuracy of the simulation, however, is highly determined by the molecular mechanics (MM) force field (FF), a set of functions with adjustable parameters to compute the potential energies from atomic positions. However, the overall quality of the FF, such as our previously published ff99SB and ff14SB, can be limited by assumptions that were made years ago. In the updated model presented here (ff19SB), we have significantly improved the backbone profiles for all 20 amino acids. We fit coupled φ/ψ parameters using 2D φ/ψ conformational scans for multiple amino acids, using as reference data the entire 2D quantum mechanics (QM) energy surface. We address the polarization inconsistency during dihedral parameter fitting by using both QM and MM in aqueous solution. Finally, we examine possible dependency of the backbone fitting on side chain rotamer. To extensively validate ff19SB parameters, and to compare to results using other Amber models, we have performed a total of ∼5 ms MD simulations in explicit solvent. Our results show that after amino-acid-specific training against QM data with solvent polarization, ff19SB not only reproduces the differences in amino-acid-specific Protein Data Bank (PDB) Ramachandran maps better but also shows significantly improved capability to differentiate amino-acid-dependent properties such as helical propensities. We also conclude that an inherent underestimation of helicity is present in ff14SB, which is (inexactly) compensated for by an increase in helical content driven by the TIP3P bias toward overly compact structures. In summary, ff19SB, when combined with a more accurate water model such as OPC, should have better predictive power for modeling sequence-specific behavior, protein mutations, and also rational protein design. Of the explicit water models tested here, we recommend use of OPC with ff19SB.

Exploring with Molecular Dynamics the Structural Fate of PLGA Oligomers in Various Solvents.

This study focuses on the solvent effects that promote preferred solvated structures of polylactic-co-glycolic acid (PLGA) oligomers of molecular weight 278, 668, and 1449 u in ethyl acetate, water, and a mixture of both solvents. Our methodology consists of all-atom, explicit solvent molecular dynamics simulations for inspection of the solvated oligomer structures at ambient conditions. Parameters for the generalized Amber force field are developed in this work for the ethyl acetate liquid and the PLGA oligomers. Energetics, oligomer radius of gyration, end-to-end distance, orientational order parameter, flexibility coefficient, and backbone dihedral angles are reported along with a size scaling property yielding a power law for PLGA oligomers in each of the three solvents considered. It is found that the PLGA oligomer has two characteristic states identified by a set of extended structures and a set of collapsed structures, the former being energetically preferred in ethyl acetate and its mixture with water. The two types of PLGA structures occur in the three solvents and although they flip from one to the other in a sporadic fashion, in ethyl acetate, the extended structures may persist for more than 20 ns. The collapsed structures are significantly more frequent in water, occurring seldom in the mixed ethyl acetate-water solvent. PLGA is a biodegradable polymer approved for use in pharmaceutical and biomedical applications. Insights provided therein are of importance for the polymer aggregation process and its glassy state in condensed phases.

TRNA Dissociation from EF-Tu after GTP Hydrolysis: Primary Steps and Antibiotic Inhibition

In each round of ribosomal translation, the translational GTPase elongation factor Tu (EF-Tu) delivers a transfer RNA (tRNA) to the ribosome. After successful decoding, EF-Tu hydrolyzes GTP, which triggers a conformational change that ultimately results in the release of the tRNA from EF-Tu. To identify the primary steps of these conformational changes and how they are prevented by the antibiotic kirromycin, we employed all-atom explicit-solvent molecular dynamics simulations of the full ribosome-EF-Tu complex. Our results suggest that after GTP hydrolysis and Pi release, the loss of interactions between the nucleotide and the switch 1 loop of EF-Tu allows domain D1 of EF-Tu to rotate relative to domains D2 and D3 and leads to an increased flexibility of the switch 1 loop. This rotation induces a closing of the D1-D3 interface and an opening of the D1-D2 interface. We propose that the opening of the D1-D2 interface, which binds the CCA tail of the tRNA, weakens the crucial EF-Tu-tRNA interactions, which lowers tRNA binding affinity, representing the first step of tRNA release. Kirromycin binds within the D1-D3 interface, sterically blocking its closure, but does not prevent hydrolysis. The resulting increased flexibility of switch 1 explains why it is not resolved in kirromycin-bound structures.

String Method for Protein-Protein Binding Free-Energy Calculations.

A powerful computational strategy to determine the equilibrium association constant of two macromolecules with explicit-solvent molecular dynamics (MD) simulations is the "geometric route", which considers the reversible physical separation of the bound complex in solution. Nonetheless, multiple challenges remain to render this type of methodology reliable and computationally efficient in practice. In particular, in one, formulation of the geometric route relies on the potential of mean force (PMF) for physically separating the two binding partners restrained along a straight axis, which must be selected prior to the calculation. However, practical applications indicate that the calculation of the separation PMF along the predefined rectilinear pathway may be suboptimal and slowly convergent. Recognizing that a rectilinear straight separation pathway is generally not representative of how the protein complex physically separates in solution, we put forth a novel theoretical framework for binding free-energy calculations, leaning on the optimal curvilinear minimum free-energy path (MFEP) determined from the string method. The proposed formalism is validated by comparing the results obtained using both rectilinear and curvilinear pathways for a prototypical host-guest complex formed by cucurbit[7]uril (CB[7]) binding benzene, and for the barnase-barstar protein complex. On the basis of multi-microsecond MD calculations, we find that the calculations following the traditional rectilinear pathway and the string-based curvilinear pathway agree quantitatively, but convergence is faster with the latter.

Molecular dynamics simulation of the binding process of ligands to the add adenine riboswitch aptamer.

Riboswitches are RNA-structured elements that modulate gene expression through changing their conformation in response to specific metabolite binding. However, the regulation mechanisms of riboswitches by ligand binding are still not well understood. At present two possible ligand-binding mechanisms have been proposed: conformational selection and induced fit. Based on explicit-solvent molecular dynamics (MD) simulations, we have studied the process of the binding of ligands (adenines) to add adenine riboswitch aptamer (AARA) in detail. Our results show that the relative high flexibility of the junction J23 of AARA allows the ligands to be captured by the binding pocket of AARA in a near-native state, which may be driven by hydrophobic and base-stacking interactions. In addition, the binding of a ligand makes the stem P1 and the junction J23 of AARA more stable than in the absence of the ligand. Generally, our analyses show that the ligand-binding process of the add adenine riboswitch may be partially embodied by a conformational selection mechanism.