Determination of Base-Flipping Free-Energy Landscapes from Nonequilibrium Stratification

Bubble Relaxation Dynamics in Homopolymer DNA Sequences

Understanding the inherent timescales of large bubbles in DNA is critical to a thorough comprehension of its physicochemical characteristics, as well as their potential role on helix opening and biological function. In this work, we employ the coarse-grained Peyrard-Bishop-Dauxois model of DNA to study relaxation dynamics of large bubbles in homopolymer DNA, using simulations up to the microsecond time scale. By studying energy autocorrelation functions of relatively large bubbles inserted into thermalised DNA molecules, we extract characteristic relaxation times from the equilibration process for both adenine-thymine (AT) and guanine-cytosine (GC) homopolymers. Bubbles of different amplitudes and widths are investigated through extensive statistics and appropriate fittings of their relaxation. Characteristic relaxation times increase with bubble amplitude and width. We show that, within the model, relaxation times are two orders of magnitude longer in GC sequences than in AT sequences. Overall, our results confirm that large bubbles leave a lasting impact on the molecule's dynamics, for times between 0.5-500 ns depending on the homopolymer type and bubble shape, thus clearly affecting long-time evolutions of the molecule.

Energy Landscapes for Base-Flipping in a Model DNA Duplex

We explore the process of base-flipping for four central bases, adenine, guanine, cytosine, and thymine, in a deoxyribonucleic acid (DNA) duplex using the energy landscape perspective. NMR imino-proton exchange and fluorescence correlation spectroscopy studies have been used in previous experiments to obtain lifetimes for bases in paired and extrahelical states. However, the difference of almost 4 orders of magnitude in the base-flipping rates obtained by the two methods implies that they are exploring different pathways and possibly different open states. Our results support the previous suggestion that minor groove opening may be favored by distortions in the DNA backbone and reveal links between sequence effects and the direction of opening, i.e., whether the base flips toward the major or the minor groove side. In particular, base flipping along the minor groove pathway was found to align toward the 5′ side of the backbone. We find that bases align toward the 3′ side of the backbone when flipping along the major groove pathway. However, in some cases for cytosine and thymine, the base flipping along the major groove pathway also aligns toward the 5′ side. The sequence effect may be caused by the polar interactions between the flipping-base and its neighboring bases on either of the strands. For guanine flipping toward the minor groove side, we find that the equilibrium constant for opening is large compared to flipping via the major groove. We find that the estimated rates of base opening, and hence the lifetimes of the closed state, obtained for thymine flipping through small and large angles along the major groove differ by 6 orders of magnitude, whereas for thymine flipping through small angles along the minor groove and large angles along the major groove, the rates differ by 3 orders of magnitude.

Seeding the multi-dimensional nonequilibrium pulling for Hamiltonian variation: indirect nonequilibrium free energy simulations at QM levels

The combination of free energy simulations in the alchemical and configurational spaces provides a feasible route to access the thermodynamic profiles under a computationally demanding target Hamiltonian. Normally, due to the significant differences between the computational cost of ab initio quantum mechanics (QM) calculations and those of semi-empirical quantum mechanics (SQM) and molecular mechanics (MM), this indirect method could be used to obtain the QM thermodynamics by combining the SQM or MM results and the SQM-to-QM or MM-to-QM corrections. In our previous work, a multi-dimensional nonequilibrium pulling framework for Hamiltonian variations was introduced based on bidirectional pulling and bidirectional reweighting. The method performs nonequilibrium free energy simulations in the configurational space to obtain the thermodynamic profile along the conformational change pathway under a selected computationally efficient Hamiltonian, and uses the nonequilibrium alchemical method to correct or perturb the thermodynamic profile to that under the target Hamiltonian. The BAR-based method is designed to achieve the best generality and transferability and thus leads to modest (∼20 fold) speedup. In this work, we explore the possibility of further accelerating the nonequilibrium free energy simulation by employing unidirectional pulling and using the selection criterion to obtain the initial configurations used to initiate nonequilibrium trajectories following the idea of adaptive steered molecular dynamics (ASMD). A single initial condition is used to seed the whole multi-dimensional nonequilibrium free energy simulation and the sampling is performed fully in the nonequilibrium ensemble. Introducing very short ps-length equilibrium sampling to grab more initial seeds could also be helpful. The ASMD scheme estimates the free energy difference with the unidirectional exponential average (EXP), but it does not follow exactly the requirements of the EXP estimator. Another deficiency of the seeding simulation is the inherently sequential or serial pulling due to the inter-segment dependency, which triggers some problems in the parallelizability of the simulation. Numerical tests are performed to grasp some insights and guidelines for using this selection-criterion-based ASMD scheme. The presented selection-criterion-based multi-dimensional ASMD scheme follows the same perturbation network of the BAR-based method, and thus could be used in various Hamiltonian-variation cases.

A General Picture of Cucurbit[8]uril Host-Guest Binding

Describing, understanding, and designing complex interaction networks within macromolecular systems remain challenging in modern chemical research. Host-guest systems, despite their relative simplicity in both the structural feature and interaction patterns, still pose problems in theoretical modeling. The barrel-shaped supramolecular container cucurbit[8]uril (CB8) shows promising functionalities in various areas, e.g., catalysis and molecular recognition. It can stably coordinate a series of structurally diverse guests with high affinities. In this work, we examine the binding of seven commonly abused drugs to the CB8 host, aiming at providing a general picture of CB8-guest binding. Extensive sampling of the configurational space of these host-guest systems is performed, and the binding pathway and interaction patterns of CB8-guest complexes are investigated. A thorough comparison of widely used fixed-charge models for drug-like molecules is presented. Iterative refitting of the atomic charges suggests significant conformation dependence of charge generation. The initial model generated at the original conformation could be inaccurate for new conformations explored during conformational search, and the newly fitted charge set improves the prediction-experiment correlation significantly. Our investigations of the configurational space of CB8-drug complexes suggest that the host-guest interactions are more complex than expected. Despite the structural simplicities of these molecules, the conformational fluctuations of the host and the guest molecules and orientations of functional groups lead to the existence of an ensemble of binding modes. The insights of the binding thermodynamics, performance of fixed-charge models, and binding patterns of the CB8-guest systems are useful for studying and elucidating the binding mechanism of other host-guest complexes.

Binding thermodynamics and interaction patterns of human purine nucleoside phosphorylase-inhibitor complexes from extensive free energy calculations

Human purine nucleoside phosphorylase (hPNP) plays a significant role in the catabolism of deoxyguanosine. The trimeric protein is an important target in the treatment of T-cell cancers and autoimmune disorders. Experimental studies on the inhibition of the hPNP observe that the first ligand bound to one of three subunits effectively inhibits the protein, while the binding of more ligands to the subsequent sites shows negative cooperativities. In this work, we performed extensive end-point and alchemical free energy calculations to determine the binding thermodynamics of the trimeric protein-ligand system. 13 Immucillin inhibitors with experimental results are under calculation. Two widely accepted charge schemes for small molecules including AM1-BCC and RESP are adopted for ligands. The results of RESP are in better agreement with the experimental reference. Further investigations of the interaction networks in the protein-ligand complexes reveal that several residues play significant roles in stabilizing the complex structure. The most commonly observed ones include PHE200, GLU201, MET219, and ASN243. The conformations of the protein in different protein-ligand complexes are observed to be similar. We expect these insights to aid the development of potent drugs targeting hPNP.

Binding Thermodynamics and Interaction Patterns of Inhibitor-Major Urinary Protein-I Binding from Extensive Free-Energy Calculations: Benchmarking AMBER Force Fields

Mouse major urinary protein (MUP) plays a key role in the pheromone communication system. The one-end-closed β-barrel of MUP-I forms a small, deep, and hydrophobic central cavity, which could accommodate structurally diverse ligands. Previous computational studies employed old protein force fields and short simulation times to determine the binding thermodynamics or investigated only a small number of structurally similar ligands, which resulted in sampled regions far from the experimental structure, nonconverged sampling outcomes, and limited understanding of the possible interaction patterns that the cavity could produce. In this work, extensive end-point and alchemical free-energy calculations with advanced protein force fields were performed to determine the binding thermodynamics of a series of MUP-inhibitor systems and investigate the inter- and intramolecular interaction patterns. Three series of inhibitors with a total of 14 ligands were simulated. We independently simulated the MUP-inhibitor complexes under two advanced AMBER force fields. Our benchmark test showed that the advanced AMBER force fields including AMBER19SB and AMBER14SB provided better descriptions of the system, and the backbone root-mean-square deviation (RMSD) was significantly lowered compared with previous computational studies with old protein force fields. Surprisingly, although the latest AMBER force field AMBER19SB provided better descriptions of various observables, it neither improved the binding thermodynamics nor lowered the backbone RMSD compared with the previously proposed and widely used AMBER14SB. The older but widely used AMBER14SB actually achieved better performance in the prediction of binding affinities from the alchemical and end-point free-energy calculations. We further analyzed the protein-ligand interaction networks to identify important residues stabilizing the bound structure. Six residues including PHE38, LEU40, PHE90, ALA103, LEU105, and TYR120 were found to contribute the most significant part of protein-ligand interactions, and 10 residues were found to provide favorable interactions stabilizing the bound state. The two AMBER force fields gave extremely similar interaction networks, and the secondary structures also showed similar behavior. Thus, the intra- and intermolecular interaction networks described with the two AMBER force fields are similar. Therefore, AMBER14SB could still be the default option in free-energy calculations to achieve highly accurate binding thermodynamics and interaction patterns.

Distributions of bubble lifetimes and bubble lengths in DNA

We investigate the distribution of bubble lifetimes and bubble lengths in DNA at physiological temperature, by performing extensive molecular dynamics simulations with the Peyrard-Bishop-Dauxois (PBD) model, as well as an extended version (ePBD) having a sequence-dependent stacking interaction, emphasizing the effect of the sequences' guanine-cytosine (GC)/adenine-thymine (AT) content on these distributions. For both models we find that base pair-dependent (GC vs AT) thresholds for considering complementary nucleotides to be separated are able to reproduce the observed dependence of the melting temperature on the GC content of the DNA sequence. Using these thresholds for base pair openings, we obtain bubble lifetime distributions for bubbles of lengths up to ten base pairs as the GC content of the sequences is varied, which are accurately fitted with stretched exponential functions. We find that for both models the average bubble lifetime decreases with increasing either the bubble length or the GC content. In addition, the obtained bubble length distributions are also fitted by appropriate stretched exponential functions and our results show that short bubbles have similar likelihoods for any GC content, but longer ones are substantially more likely to occur in AT-rich sequences. We also show that the ePBD model permits more, longer-lived, bubbles than the PBD system.

SAMPL7 TrimerTrip host-guest binding affinities from extensive alchemical and end-point free energy calculations

The prediction of host-guest binding affinities with computational modelling is still a challenging task. In the 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge, a new host named TrimerTrip was synthesized and the thermodynamic parameters of 16 structurally diverse guests binding to the host were characterized. In the TrimerTrip-guest challenge, only structures of the host and the guests are provided, which indicates that the predictions of both the binding poses and the binding affinities are under assessment. In this work, starting from the binding poses obtained from our previous enhanced sampling simulations in the configurational space, we perform extensive alchemical and end-point free energy calculations to calculate the host-guest binding affinities retrospectively. The alchemical predictions with two widely accepted charge schemes (i.e. AM1-BCC and RESP) are in good agreement with the experimental reference, while the end-point estimates perform poorly in reproducing the experimental binding affinities. Aside from the absolute value of the binding affinity, the rank of binding free energies is also crucial in drug design. Surprisingly, the end-point MM/PBSA method seems very powerful in reproducing the experimental rank of binding affinities. Although the length of our simulations is long and the intermediate spacing is dense, the convergence behavior is not very good, which may arise from the flexibility of the host molecule. Enhanced sampling techniques in the configurational space may be required to obtain fully converged sampling. Further, as the length of sampling in alchemical free energy calculations already achieves several hundred ns, performing direct simulations of the binding/unbinding event in the physical space could be more useful and insightful. More details about the binding pathway and mechanism could be obtained in this way. The nonequilibrium method could also be a nice choice if one insists to use the alchemical method, as the intermediate sampling is avoided to some extent.

Theoretical understanding of the thermodynamics and interactions in transcriptional regulator TtgR-ligand binding

The transcriptional regulator TtgR belongs to the TetR family of transcriptional repressors. It depresses the transcription of the TtgABC operon and itself and thus regulates the extrusion of noxious chemicals with efflux pumps in bacterial cells. As the ligand-binding domain of TtgR is rather flexible, it can bind with a number of structurally diverse ligands, such as antibiotics, flavonoids and aromatic solvents. In the current work, we perform equilibrium and nonequilibrium alchemical free energy simulation to predict the binding affinities of a series of ligands targeting the TtgR protein and an agreement between the theoretical prediction and the experimental result is observed. End-point methods MM/PBSA and MM/GBSA are also employed for comparison. We further study the interaction maps and contacts between the protein and the ligand and identify important interactions in the protein-ligand binding cases. The dynamics fluctuation and secondary structures are also investigated. The current work sheds light on atomic and thermodynamic understanding of the TtgR-ligand interactions.

BAR-based multi-dimensional nonequilibrium pulling for indirect construction of QM/MM free energy landscapes: from semi-empirical to ab initio

The indirect method for the construction of quantum mechanics (QM)/molecular mechanics (MM) free energy landscapes provides a cheaper alternative for free energy simulations at the QM level. The indirect method features a direct calculation of the free energy profile with a computationally efficient but less accurate Hamiltonian (i.e. low-level Hamiltonian) and a low-level-to-high-level correction. In the thermodynamic cycle, the direct low-level calculation along the physically meaningful reaction coordinate is corrected via the alchemical method, which is often achieved with perturbation-based techniques. In our previous work, a multi-dimensional nonequilibrium pulling framework is proposed for the indirect construction of QM/MM free energy landscapes. Previously, we focused on obtaining semi-empirical QM (SQM) results indirectly from direct MM simulations and MM to SQM corrections. In this work, we apply this method to obtain results under ab initio QM Hamiltonians by combining direct SQM results and SQM to QM corrections. A series of SQM and QM Hamiltonians are benchmarked. It is observed that PM6 achieves the best performance among the low-level Hamiltonians. Therefore, we recommend using PM6 as the low-level theory in the indirect free energy simulation. Considering its higher similarity to the high-level Hamiltonians, PM6 corrected with the bond charge correction could be more accurate than the existing AM1-BCC model. Another central result in the current work is a basic protocol of choosing the strength of restraints and an appropriate time step in nonequilibrium free energy simulation at the stiff spring limit. We provide theoretical derivations to emphasize the importance of using a sufficiently large force constant and choosing an appropriate time step. It is worth noting that a general rule of thumb for choosing the time step, according to our derivation, is that a time step of 1 fs or smaller should be used, as long as the stiff spring approximation is employed, even in simulations with constraints on bonds involving hydrogen atoms.

Sulfur-substitution-induced base flipping in the DNA duplex

Base flipping is widely observed in a number of important biological processes. The genetic codes deposited inside the DNA duplex become accessible to external agents upon base flipping. The sulfur substitution of guanine leads to thioguanine, which alters the thermodynamic stability of the GC base pairs and the GT mismatches. Experimental studies conclude that the sulfur substitution decreases the lifetime of the GC base pair. In this work, under three AMBER force fields for nucleotide systems, we firstly performed equilibrium and nonequilibrium free energy simulations to investigate the variation of the thermodynamic profiles in base flipping upon sulfur substitution. It is found that the bsc0 modification, the bsc1 modification and the OL15 modification of AMBER force fields are able to qualitatively describe the sulfur-substitution dependent behavior of the thermodynamics. However, only the two last-generation AMBER force fields are able to provide quantitatively correct predictions. The second computational study on the sulfur substitutions focused on the relative stability of the S6G-C base pair and the S6G-T mismatch. Two conflicting experimental observations were reported by the same authors. One suggested that the S6G-C base pair was more stable, while the other concludes that the S6G-T mismatch was more stable. We answered this question by constructing the free energy profiles along the base flipping pathway computationally.