Tunnel connects lipid bilayer to occluded odorant-binding site of insect olfactory receptor.

MhOR5, an insect olfactory receptor (OR), has an occluded binding site for the odorant eugenol in both the open and closed states of the ion channel. We used atomistic molecular dynamics simulation (MD) and steered molecular dynamics to examine possible tunnels to the odorant binding site from the protein surface. Four high probability tunnels were identified in the MD results. Surprisingly, three of the tunnels connect the ligand binding site to the lipid bilayer. We found sharp 30%-50% increases or decreases in tunnel bottleneck areas over 70 nsec MD trajectories, both in the ligand-bound and unliganded OR structures. Steered MD showed that eugenol follows the tunnels to the protein surface, and the potential of mean force is quantitatively consistent with the known affinity of eugenol for MhOR5. We examined AlphaFold-generated models of 21 other insect ORs, and we found that 19 had odorant binding sites and tunnels in similar positions to MhOR5. The possibility of a tunnel between the odorant binding site and the lipid bilayer in insect ORs suggests new experiments to test molecular mechanisms for insect odorant reception.

Quantitative characterization of the path of glucose diffusion facilitated by human glucose transporter 1.

Glucose transporter GLUT1 is ubiquitously expressed in the human body from the red cells to the blood-brain barrier to the skeletal muscles. It is physiologically relevant to understand how GLUT1 facilitates diffusion of glucose across the cell membrane. It is also pathologically relevant because GLUT1 deficiency causes neurological disorders and anemia and because GLUT1 overexpression fuels the abnormal growth of cancer cells. This article presents a quantitative investigation of GLUT1 based on all-atom molecular-dynamics (MD) simulations of the transporter embedded in lipid bilayers of asymmetric inner-and-outer-leaflet lipid compositions, subject to asymmetric intra-and-extra-cellular environments. This is in contrast with the current literature of MD studies that have not considered both of the aforementioned asymmetries of the cell membrane. The equilibrium (unbiased) dynamics of GLUT1 shows that it can facilitate glucose diffusion across the cell membrane without undergoing large-scale conformational motions. The Gibbs free-energy profile, which is still lacking in the current literature of GLUT1, quantitatively characterizes the diffusion path of glucose from the periplasm, through an extracellular gate of GLUT1, on to the binding site, and off to the cytoplasm. This transport mechanism is validated by the experimental data that GLUT1 has low water-permeability, uptake-efflux symmetry, and 10 kcal/mol Arrhenius activation barrier around 37 °C.

Effects of the N-terminal dynamics on the conformational states of human dopamine transporter.

Dopamine transporter mediates the neurotransmitter dopamine homeostasis in a sodium-dependent manner. The transport process involves an alternating access of a substrate to the extracellular and intracellular spaces, which is associated with different conformational states of the transporter. However, the underlying mechanism of modulation of the state transition remains elusive. Here we present a computational simulation study of human dopamine transporter to explore its two end states (outward-facing open and inward-facing open) that have not been determined experimentally. We show that the full-length transporter may tend to adopt the inward-facing open state in its free state. The binding of an amphetamine may not trap the transporter in the outward-facing open state with increasing length of the N-terminal. Furthermore, we identify distinct patterns in the interaction networks between the N-terminal and the intracellular region that could stabilize the state of the transporter, independent of substrate binding and phosphorylation. Our results reveal the essential role of the N-terminal dynamics in modulating the functional states of the dopamine transporter, providing molecular insights into the coupling of conformational transition and substrate passage in neurotransmitter transporters.

Aquaglyceroporin AQP7's affinity for its substrate glycerol: Have we reached convergence in the computed values of glycerol-aquaglyceroporin affinity?

AQP7 is one of the four human aquaglyceroporins that facilitate glycerol transport across the cell membrane, a biophysical process that is essential in human physiology. Therefore, it is interesting to compute AQP7's affinity for its substrate (glycerol) with reasonable certainty to compare with the experimental data suggesting high affinity in contrast with most computational studies predicting low affinity. In this study aimed at computing the AQP7-glycerol affinity with high confidence, we implemented a direct computation of the affinity from unbiased equilibrium molecular dynamics (MD) simulations of three all-atom systems constituted with 0.16M, 4.32M, and 10.23M atoms, respectively. These three sets of simulations manifested a fundamental physics law that the intrinsic fluctuations of pressure in a system are inversely proportional to the system size (the number of atoms in it). These simulations showed that the computed values of glycerol-AQP7 affinity are dependent upon the system size (the inverse affinity estimations were, respectively, 47.3 mM, 1.6 mM, and 0.92 mM for the three model systems). In this, we obtained a lower bound for the AQP7-glycerol affinity (an upper bound for the dissociation constant). Namely, the AQP7-glycerol affinity is stronger than 1087/M (the dissociation constant is less than 0.92 mM). Additionally, we conducted hyper steered MD (hSMD) simulations to map out the Gibbs free-energy profile. From the free-energy profile, we produced an independent computation of the AQP7-glycerol dissociation constant being approximately 0.18 mM.

Structural Insights into the Human Mitochondrial Pyruvate Carrier Complexes.

Pyruvate metabolism requires the mitochondrial pyruvate carrier (MPC) proteins to transport pyruvate from the intermembrane space through the inner mitochondrial membrane to the mitochondrial matrix. The lack of the atomic structures of MPC hampers the understanding of the functional states of MPC and molecular interactions with the substrate or inhibitor. Here, we develop the de novo models of human MPC complexes and characterize the conformational dynamics of the MPC heterodimer formed by MPC1 and MPC2 (MPC1/2) by computational simulations. Our results reveal that functional MPC1/2 prefers to adopt an inward-open conformation, with the carrier open to the matrix side, whereas the outward-open states are less populated. The energy barrier for pyruvate transport in MPC1/2 is low enough, and the inhibitor UK5099 blocks the pyruvate transport by stably binding to MPC1/2. Notably, consistent with experimental results, the MPC1 L79H mutation significantly alters the conformations of MPC1/2 and thus fails for substrate transport. However, the MPC1 R97W mutation seems to retain the transport activity. The present de novo models of MPC complexes provide structural insights into the conformational states of MPC complexes and mechanistic understanding of interactions between the substrate/inhibitor and MPC proteins.

Association of sigma-1 receptor with dopamine transporter attenuates the binding of methamphetamine via distinct helix-helix interactions.

Dopamine transporter (DAT) and sigma-1 receptor (σ1R) are potential therapeutic targets to reduce the psychostimulant effects induced by methamphetamine (METH). Interaction of σ1R with DAT could modulate the binding of METH, but the molecular basis of the association of the two transmembrane proteins and how their interactions mediate the binding of METH to DAT or σ1R remain unclear. Here, we characterize the protein-ligand and protein-protein interactions at a molecular level by various theoretical approaches. The present results show that METH adopts a different binding pose in the binding pocket of σ1R and is more likely to act as an agonist. The relatively lower binding affinity of METH to σ1R supports the role of antagonists as inhibitors that protect against METH-induced effects. We demonstrate that σ1R could bind to Drosophila melanogaster DAT (dDAT) through interactions with either the transmembrane helix α12 or α5 of dDAT. Our results showed that the truncated σ1R displays stronger association with dDAT than the full-length σ1R. Although different helix-helix interactions between σ1R and dDAT lead to distinct effects on the dynamics of individual protein, both associations attenuate the binding affinity of METH to dDAT, particularly in the interactions with the helix α5 of dDAT. Together, the present study provides the first computational investigation on the molecular mechanism of coupling METH binding and the association of σ1R with dDAT.

Molecular determinant of substrate binding and specificity of cytochrome P450 2J2.

Cytochrome P450 2J2 (CYP2J2) is responsible for the epoxidation of endogenous arachidonic acid, and is involved in the metabolism of exogenous drugs. To date, no crystal structure of CYP2J2 is available, and the proposed structural basis for the substrate recognition and specificity in CYP2J2 varies with the structural models developed using different computational protocols. In this study, we developed a new structural model of CYP2J2, and explored its sensitivity to substrate binding by molecular dynamics simulations of the interactions with chemically similar fluorescent probes. Our results showed that the induced-fit binding of these probes led to the preferred active poses ready for the catalysis by CYP2J2. Divergent conformational dynamics of CYP2J2 due to the binding of each probe were observed. However, a stable hydrophobic clamp composed of residues I127, F310, A311, V380, and I487 was identified to restrict any substrate access to the active site of CYP2J2. Molecular docking of a series of compounds including amiodarone, astemizole, danazol, ebastine, ketoconazole, terfenadine, terfenadone, and arachidonic acid to CYP2J2 confirmed the role of those residues in determining substrate binding and specificity of CYP2J2. In addition to the flexibility of CYP2J2, the present work also identified other factors such as electrostatic potential in the vicinity of the active site, and substrate strain energy and property that have implications for the interpretation of CYP2J2 metabolism.

Quantitative study of unsaturated transport of glycerol through aquaglyceroporin that has high affinity for glycerol.

The structures of several aquaglyceroporins have been resolved to atomic resolution showing two or more glycerols bound inside a channel and confirming a glycerol-facilitator's affinity for its substrate glycerol. However, the kinetics data of glycerol transport experiments all point to unsaturated transport that is characteristic of low substrate affinity in terms of the Michaelis-Menten kinetics. In this article, we present an in silico-in vitro research focused on AQP3, one of the human aquaglyceroporins that is natively expressed in the abundantly available erythrocytes. We conducted 2.1 µs in silico simulations of AQP3 embedded in a model erythrocyte membrane with intracellular-extracellular asymmetries in leaflet lipid compositions and compartment salt ions. From the equilibrium molecular dynamics (MD) simulations, we elucidated the mechanism of glycerol transport at high substrate concentrations. From the steered MD simulations, we computed the Gibbs free-energy profile throughout the AQP3 channel. From the free-energy profile, we quantified the kinetics of glycerol transport that is unsaturated due to glycerol-glycerol interactions mediated by AQP3 resulting in the concerted movement of two glycerol molecules for the transport of one glycerol molecule across the cell membrane. We conducted in vitro experiments on glycerol uptake into human erythrocytes for a wide range of substrate concentrations and various temperatures. The experimental data quantitatively validated our theoretical-computational conclusions on the unsaturated glycerol transport through AQP3 that has high affinity for glycerol.

Identification of a New Allosteric Binding Site for Cocaine in Dopamine Transporter.

Dopamine (DA) transporter (DAT) is a major target for psychostimulant drugs of abuse such as cocaine that competitively binds to DAT, inhibits DA reuptake, and consequently increases synaptic DA levels. In addition to the central binding site inside DAT, the available experimental evidence suggests the existence of alternative binding sites on DAT, but detection and characterization of these sites are challenging by experiments alone. Here, we integrate multiple computational approaches to probe the potential binding sites on the wild-type Drosophila melanogaster DAT and identify a new allosteric site that displays high affinity for cocaine. This site is located on the surface of DAT, and binding of cocaine is primarily dominated by interactions with hydrophobic residues surrounding the site. We show that cocaine binding to this new site allosterically reduces the binding of DA/cocaine to the central binding pocket, and simultaneous binding of two cocaine molecules to a single DAT seems infeasible. Furthermore, we find that binding of cocaine to this site stabilizes the conformation of DAT but alters the conformational population and thereby reduces the accessibility by DA, providing molecular insights into the inhibitory mechanism of cocaine. In addition, our results indicate that the conformations induced by cocaine binding to this site may be relevant to the oligomerization of DAT, highlighting a potential role of this new site in modulating the function of DAT.

In silico simulations of erythrocyte aquaporins with quantitative in vitro validation.

Modelling water and membrane lipids is an essential element in the computational research of biophysical/biochemical processes such as water transport across the cell membrane. In this study, we examined the accuracies of two popular water models, TIP3P and TIP4P, in the molecular dynamics simulations of erythrocyte aquaporins (AQP1 and AQP3). We modelled the erythrocyte membrane as an asymmetric lipid bilayer with appropriate lipid compositions of its inner and outer leaflet, in comparison with a symmetric lipid bilayer of a single lipid type. We computed the AQP1/3 permeabilities with the transition state theory with full correction for recrossing events. We also conducted cell swelling assays for water transport across the erythrocyte membrane. The experimental results agree with the TIP3P water-erythrocyte membrane model, in confirmation of the expected accuracy of the erythrocyte membrane model, the TIP3P water model, and the CHARMM parameters for water-protein interactions.

Thermodynamic Integration in 3n Dimensions without Biases or Alchemy for Protein Interactions.

Thermodynamic integration (TI), a powerful formalism for computing Gibbs free energy, has been implemented for many biophysical processes with alchemical schemes that require delicate human efforts to choose/design biasing potentials for sampling the desired biophysical events and to remove their artifactitious consequences afterwards. Theoretically, an alchemical scheme is exact but practically, an unsophisticated implementation of this exact formula can cause error amplifications. Small relative errors in the input parameters can be amplified many times in their propagation into the computed free energy [due to subtraction of similar numbers such as (105 ± 5)‒(100 ± 5) = 5 ± 7]. In this paper, we present an unsophisticated implementation of TI in 3n dimensions (3nD) (n=1,2,3…) for the potential of mean force along a 3nD path connecting one state in the bound state ensemble to one state in the unbound state ensemble. Fluctuations in these 3nD are integrated in the bound and unbound state ensembles but not along the 3nD path. Using TI3nD, we computed the standard binding free energies of three protein complexes: trometamol in Salmonella effector SpvD (n=1), biotin-avidin (n=2), and Colicin E9 endonuclease with cognate immunity protein Im9 (n=3). We employed three different protocols in three independent computations of E9-Im9 to show TI3nD's robustness. We also computed the hydration energies of ten biologically relevant compounds (n=1 for water, acetamide, urea, glycerol, trometamol, ammonium and n=2 for erythritol, 1,3-propanediol, xylitol, biotin). Each of the 15 computations is accomplishable within one (for hydration) to ten (for E9-Im9) days on an inexpensive GPU workstation. The computed results all agree with the available experimental data.

Application of the Brown dynamics fluctuation-dissipation theorem to the study of Plasmodium berghei transporter protein PbAQP.

In this article, the Brownian dynamics fluctuation-dissipation theorem (BD-FDT) is applied to the study of transport of neutral solutes across the cellular membrane of Plasmodium berghei (Pb), a disease-causing parasite. Pb infects rodents and causes symptoms in laboratory mice that are comparable to human malaria caused by Plasmodium falciparum (Pf). Due to the relative ease of its genetic engineering, P. berghei has been exploited as a model organism for the study of human malaria. P. berghei expresses one type of aquaporin (AQP), PbAQP, and, in parallel, P. falciparum expresses PfAQP. Either PbAQP or PfAQP is a multifunctional channel protein in the plasma membrane of the rodent/human malarial parasite for homeostasis of water, uptake of glycerol, and excretion of some metabolic wastes across the cell membrane. This FDT-study of the channel protein PbAQP is to elucidate how and how strongly it interacts with water, glycerol, and erythritol. It is found that erythritol, which binds deep inside the conducting pore of PbAQP/PfAQP, inhibits the channel protein's functions of conducting water, glycerol etc. This points to the possibility that erythritol, a sugar substitute, may inhibit the malarial parasites in rodents and in humans.

Extracellular gating of glucose transport through GLUT 1.

The ubiquitous glucose transporter 1 (GLUT1) is physiologically and pathologically relevant in energy metabolism of the CNS, skeletal muscles, cancer cells etc. Extensive experiments on GLUT1 produced thorough understandings of its expressions, functions, and structures which were recently resolved to atomic accuracy. However, theoretical understandings are still controversial about how GLUT1 facilitates glucose diffusion across the cell membrane. Molecular dynamics (MD) simulations of the current literature have GLUT1 embedded in a symmetric bilayer of a single lipid type. They provide atomistic illustrations of the alternating access theory (AAT), but the simulation results are inconsistent with the undisputed experimental data of kinetics showing rapid transport of glucose at near-physiological temperatures, high Arrhenius activation barrier in zero-trans uptake, and large trans-acceleration at sub-physiological temperatures. In this research, we embedded GLUT1 in an asymmetric bilayer of multiple lipids to better mimic the erythrocyte membrane. We ran unbiased MD simulations at 37 °C and at 5 °C and found a new mechanism of glucose transport via GLUT1: The extracellular (EC) gate opened wide for EC glucopyranose at 37 °C and, only in the presence of intracellular (IC) glucose, at 5 °C. In the absence of IC glucose at 5 °C, the EC gate opened narrowly for acyclic glucose, gating out glucopyranose. This EC-gating mechanism is simpler than AAT and yet it well explains for the rapid glucose transport at near-physiological temperatures and large trans-acceleration at sub-physiological temperatures. It also explains why zero-trans uptake (involving the pyranose-to-aldehyde transformation) has an Arrhenius barrier ∼20 kcal/mol higher than the equilibrium exchange transport.

Single-channel permeability and glycerol affinity of human aquaglyceroporin AQP3.

For its fundamental relevance, transport of water and glycerol across the erythrocyte membrane has long been investigated before and after the discovery of aquaporins (AQPs), the membrane proteins responsible for water and glycerol transport. AQP1 is abundantly expressed in the human erythrocyte for maintaining its hydrohomeostasis where AQP3 is also expressed (at a level ~30-folds lower than AQP1) facilitating glycerol transport. This research is focused on two of the remaining questions: How permeable is AQP3 to water? What is the glycerol-AQP3 affinity under near-physiological conditions? Through atomistic modelling and large-scale simulations, we found that AQP3 is two to three times more permeable to water than AQP1 and that the glycerol-AQP3 affinity is approximately 500/M. Using these computed values along with the data from the latest literature on AQP1 and on erythrocyte proteomics, we estimated the water and glycerol transport rates across the membrane of an entire erythrocyte. We used these rates to predict the time courses of erythrocyte swelling-shrinking in response to inward and outward osmotic gradients. Experimentally, we monitored the time course of human erythrocytes when subject to an osmotic or glycerol gradient with light scattering in a stopped-flow spectrometer. We observed close agreement between the experimentally measured and the computationally predicted time courses of erythrocytes, which corroborated our computational conclusions on the AQP3 water-permeability and the glycerol-AQP3 affinity.

Gibbs Free-Energy Gradient along the Path of Glucose Transport through Human Glucose Transporter 3.

Fourteen glucose transporters (GLUTs) play essential roles in human physiology by facilitating glucose diffusion across the cell membrane. Due to its central role in the energy metabolism of the central nervous system, GLUT3 has been thoroughly investigated. However, the Gibbs free-energy gradient (what drives the facilitated diffusion of glucose) has not been mapped out along the transport path. Some fundamental questions remain. Here we present a molecular dynamics study of GLUT3 embedded in a lipid bilayer to quantify the free-energy profile along the entire transport path of attracting a ß-d-glucose from the interstitium to the inside of GLUT3 and, from there, releasing it to the cytoplasm by Arrhenius thermal activation. From the free-energy profile, we elucidate the unique Michaelis-Menten characteristics of GLUT3, low KM and high VMAX, specifically suitable for neurons' high and constant demand of energy from their low-glucose environments. We compute GLUT3's binding free energy for ß-d-glucose to be -4.6 kcal/mol in agreement with the experimental value of -4.4 kcal/mol ( KM = 1.4 mM). We also compute the hydration energy of ß-d-glucose, -18.0 kcal/mol vs the experimental data, -17.8 kcal/mol. In this, we establish a dynamics-based connection from GLUT3's crystal structure to its cellular thermodynamics with quantitative accuracy. We predict equal Arrhenius barriers for glucose uptake and efflux through GLUT3 to be tested in future experiments.

Thermodynamics of Amyloid-ß Fibril Elongation: Atomistic Details of the Transition State.

Amyloid-ß (Aß) fibrils and plaques are one of the hallmarks of Alzheimer's disease. While the kinetics of fibrillar growth of Aß have been extensively studied, several vital questions remain. In particular, the atomistic origins of the Arrhenius barrier observed in experiments have not been elucidated. Employing the familiar thermodynamic integration method, we have directly simulated the dissociation of an Aß(15-40) (D23N mutant) peptide from the surface of a filament along its most probable path (MPP) using all-atom molecular dynamics. This allows for a direct calculation of the free energy profile along the MPP, revealing a multipeak energetic barrier between the free peptide state and the aggregated state. By definition of the MPP, this simulated unbinding process represents the reverse of the physical elongation pathway, allowing us to draw biophysically relevant conclusions from the simulation data. Analyzing the detailed atomistic interactions along the MPP, we identify the atomistic origins of these peaks as resulting from the dock-lock mechanism of filament elongation. Careful analysis of the dynamics of filament elongation could prove key to the development of novel therapeutic strategies for amyloid-related diseases.

Affinity and path of binding xylopyranose unto E. coli xylose permease.

Glucose transporters (GLUTs), expressed in all types of human cells, are responsible for the uptake of sugars as the primary energy source for the normal functions of good cells and for the abnormal growth of cancer cells. The E. coli xylose permease (XylE), a homologue of human GLUTs, has been investigated more thoroughly than other major facilitator proteins in the current literature. In this paper, we present a molecular dynamics (MD) study of an all-atom model system to elucidate the atomistic details and the free-energy landscape along the path of binding a xylopyranose (XYP) from the extracellular space to the inside of the transporter protein XylE. From the MD simulations, the Gibbs free energy of binding was found to be -4.4kcal/mol in agreement with the experimental value of -4.7kcal/mol. The accuracy of our study is further shown in the computed hydration energy of XYP of -14.6kcal/mol in comparison with the experimental data of -15.0kcal/mol. Along the binding path, the Gibbs free energy of the XYP-XylE complex first rises from zero in the dissociated state to approximately 4 kcal/mol in the transition state (when XylE slightly increases its opening toward the extracellular side to accommodate XYP) before dropping down to -9.0 kcal/mol in the bound state. These quantitative insights indicate the fast equilibration between the bound and the unbound states of XylE and XYP. They also serve as an atomistic-dynamic corroboration of the experimental conclusion that XylE is a high-affinity sugar transporter.

Computing osmotic permeabilities of aquaporins AQP4, AQP5, and GlpF from near-equilibrium simulations.

Measuring or computing the single-channel permeability of aquaporins/aquaglyceroporins (AQPs) has long been a challenge. The measured values scatter over an order of magnitude but the corresponding Arrhenius activation energies converge in the current literature. Osmotic flux through an AQP was simulated as water current forced through the channel by kilobar hydraulic pressure or theoretically approximated as single-file diffusion. In this paper, we report large scale simulations of osmotic current under sub M gradient through three AQPs (water channels AQP4 and AQP5 and glycerol-water channel GlpF) using the mature particle mesh Ewald technique (PME) for which the established force fields have been optimized with known accuracy. These simulations were implemented with hybrid periodic boundary conditions devised to avoid the artifactitious mixing across the membrane in a regular PME simulation. The computed single-channel permeabilities at 5°C and 25°C are in agreement with recently refined experiments on GlpF. The Arrhenius activation energies extracted from our simulations for all the three AQPs agree with the in vitro measurements. The single-file diffusion approximations from our large-scale simulations are consistent with the current literature on smaller systems. From these unambiguous agreements among the in vitro and in silico studies, we observe the quantitative accuracy of the all-atom force fields of the current literature for water-channel biology. We also observe that AQP4, that is particularly rich in the central nervous system, is more efficient in water conduction and more temperature-sensitive than other water-only channels (excluding glycerol channels that also conduct water when not inhibited by glycerol).

Elongation affinity, activation barrier, and stability of Aß42 oligomers/fibrils in physiological saline.

Amyloid-beta (Aß) peptides, Aß40 and the more neurotoxic Aß42, have been the subject of many research efforts for Alzheimer's disease. In two recent independent investigations, the atomistic structure of Aß42 fibril has been clearly established in the S-shaped conformation consisting of three ß-sheets stabilized by salt bridges formed between the Lys28 sidechain and the C-terminus of Ala42. This structure distinctively differs from the long-known structure of Aß40 in the ß-hairpin shaped conformation consisting of two ß-sheets. Recent in silico investigations based on all-atom models have reached closer agreement with the in vitro measurements of Aß40 thermodynamics. In this study, we present an in silico investigation of Aß42 thermodynamics. Using the established force field parameters in seven sets of all-atom simulations, we examined the stability of small Aß42 oligomers in physiological saline. We computed the elongation affinity of the S-shaped Aß42 fibril, reaching agreement with the experimental data. We also estimated the Arrhenius activation barrier along the elongation pathway (from the disordered conformation of a free Aß42 peptide to its S-shaped conformation on a fibril) that amounts to about 16 kcal/mol, which is consistent with the experimental data. Based on these quantitative agreements, we conclude that aggregation of Aß42 peptides into fibrils is thermodynamically slow without precipitation by extrinsic factors such as heparan sulfate proteoglycan and highlight the possibility to prevent Aß42 aggregation by eliminating some precipitation factors or by increasing competitive agents to capture and transport free Aß42 peptides from the cerebrospinal fluid.

Computing the binding affinity of a ligand buried deep inside a protein with the hybrid steered molecular dynamics.

Computing the ligand-protein binding affinity (or the Gibbs free energy) with chemical accuracy has long been a challenge for which many methods/approaches have been developed and refined with various successful applications. False positives and, even more harmful, false negatives have been and still are a common occurrence in practical applications. Inevitable in all approaches are the errors in the force field parameters we obtain from quantum mechanical computation and/or empirical fittings for the intra- and inter-molecular interactions. These errors propagate to the final results of the computed binding affinities even if we were able to perfectly implement the statistical mechanics of all the processes relevant to a given problem. And they are actually amplified to various degrees even in the mature, sophisticated computational approaches. In particular, the free energy perturbation (alchemical) approaches amplify the errors in the force field parameters because they rely on extracting the small differences between similarly large numbers. In this paper, we develop a hybrid steered molecular dynamics (hSMD) approach to the difficult binding problems of a ligand buried deep inside a protein. Sampling the transition along a physical (not alchemical) dissociation path of opening up the binding cavity---pulling out the ligand---closing back the cavity, we can avoid the problem of error amplifications by not relying on small differences between similar numbers. We tested this new form of hSMD on retinol inside cellular retinol-binding protein 1 and three cases of a ligand (a benzylacetate, a 2-nitrothiophene, and a benzene) inside a T4 lysozyme L99A/M102Q(H) double mutant. In all cases, we obtained binding free energies in close agreement with the experimentally measured values. This indicates that the force field parameters we employed are accurate and that hSMD (a brute force, unsophisticated approach) is free from the problem of error amplification suffered by many sophisticated approaches in the literature.