In metabolism research, thermodynamics is usually used to determine the directionality of a reaction or the feasibility of a pathway. However, the relationship between thermodynamic potentials and fluxes is not limited to questions of directionality: thermodynamics also affects the kinetics of reactions through the flux-force relationship, which states that the logarithm of the ratio between the forward and reverse fluxes is directly proportional to the change in Gibbs energy due to a reaction (ΔrG'). Accordingly, if an enzyme catalyzes a reaction with a ΔrG' of -5.7 kJ/mol then the forward flux will be roughly ten times the reverse flux. As ΔrG' approaches equilibrium (ΔrG' = 0 kJ/mol), exponentially more enzyme counterproductively catalyzes the reverse reaction, reducing the net rate at which the reaction proceeds. Thus, the enzyme level required to achieve a given flux increases dramatically near equilibrium. Here, we develop a framework for quantifying the degree to which pathways suffer these thermodynamic limitations on flux. For each pathway, we calculate a single thermodynamically-derived metric (the Max-min Driving Force, MDF), which enables objective ranking of pathways by the degree to which their flux is constrained by low thermodynamic driving force. Our framework accounts for the effect of pH, ionic strength and metabolite concentration ranges and allows us to quantify how alterations to the pathway structure affect the pathway's thermodynamics. Applying this methodology to pathways of central metabolism sheds light on some of their features, including metabolic bypasses (e.g., fermentation pathways bypassing substrate-level phosphorylation), substrate channeling (e.g., of oxaloacetate from malate dehydrogenase to citrate synthase), and use of alternative cofactors (e.g., quinone as an electron acceptor instead of NAD). The methods presented here place another arrow in metabolic engineers' quiver, providing a simple means of evaluating the thermodynamic and kinetic quality of different pathway chemistries that produce the same molecules.

MicroRNAs (miRNAs) are endogenously produced approximately 21-nt riboregulators that associate with Argonaute (Ago) proteins to direct mRNA cleavage or repress the translation of complementary RNAs. Capturing the molecular mechanisms of miRNA interacting with its target will not only reinforce the understanding of underlying RNA interference but also fuel the design of more effective small-interfering RNA strands. To address this, in the present work the RNA-bound (Ago-miRNA, Ago-miRNA-target) and RNA-free Ago forms were analyzed by performing both molecular dynamics simulations and thermodynamic analysis. Based on the principal component analysis results of the simulation trajectories as well as the correlation analysis in fluctuations of residues, we discover that: 1) three important (PAZ, Mid and PIWI) domains exist in Argonaute which define the global dynamics of the protein; 2) the interdomain correlated movements are so crucial for the interaction of Ago-RNAs that they not only facilitate the relaxation of the interactions between residues surrounding the RNA binding channel but also induce certain conformational changes; and 3) it is just these conformational changes that expand the cavity of the active site and open putative pathways for both the substrate uptake and product release. In addition, by thermodynamic analysis we also discover that for both the guide RNA 5'-end recognition and the facilitated site-specific cleavage of the target, the presence of two metal ions (of Mg(2+)) plays a predominant role, and this conclusion is consistent with the observed enzyme catalytic cleavage activity in the ternary complex (Ago-miRNA-mRNA). Our results find that it is the set of arginine amino acids concentrated in the nucleotide-binding channel in Ago, instead of the conventionally-deemed seed base-paring, that makes greater contributions in stabilizing the binding of the nucleic acids to Ago.

Metabolic engineering in the post-genomic era is characterised by the development of new methods for metabolomics and fluxomics, supported by the integration of genetic engineering tools and mathematical modelling. Particularly, constraint-based stoichiometric models have been widely studied: (i) flux balance analysis (FBA) (in silico), and (ii) metabolic flux analysis (MFA) (in vivo). Recent studies have enabled the incorporation of thermodynamics and metabolomics data to improve the predictive capabilities of these approaches. However, an in-depth comparison and evaluation of these methods is lacking. This study presents a thorough analysis of two different in silico methods tested against experimental data (metabolomics and 13C-MFA) for the mesophile Escherichia coli. In particular, a modified version of the recently published matTFA toolbox was created, providing a broader range of physicochemical parameters. Validating against experimental data allowed the determination of the best physicochemical parameters to perform the TFA (Thermodynamics-based Flux Analysis). An analysis of flux pattern changes in the central carbon metabolism between 13C-MFA and TFA highlighted the limited capabilities of both approaches for elucidating the anaplerotic fluxes. In addition, a method based on centrality measures was suggested to identify important metabolites that (if quantified) would allow to further constrain the TFA. Finally, this study emphasised the need for standardisation in the fluxomics community: novel approaches are frequently released but a thorough comparison with currently accepted methods is not always performed.

Exposure to the environmental toxin β-methylamino-L-alanine (BMAA) is linked to amyotrophic lateral sclerosis (ALS), but its disease-promoting mechanism remains unknown. We propose that incorporation of BMAA into the ALS-linked protein Cu,Zn superoxide dismutase (SOD1) upon translation promotes protein misfolding and aggregation, which has been linked to ALS onset and progression. Using molecular simulation and predictive energetic computation, we demonstrate that substituting any serine with BMAA in SOD1 results in structural destabilization and aberrant dynamics, promoting neurotoxic SOD1 aggregation. We propose that translational incorporation of BMAA into SOD1 is directly responsible for its toxicity in neurodegeneration, and BMAA modification of SOD1 may serve as a biomarker of ALS.

Thermodynamics dictates the structure and function of metabolism. Redox reactions drive cellular energy and material flow. Hence, accurately quantifying the thermodynamics of redox reactions should reveal design principles that shape cellular metabolism. However, only few redox potentials have been measured, and mostly with inconsistent experimental setups. Here, we develop a quantum chemistry approach to calculate redox potentials of biochemical reactions and demonstrate our method predicts experimentally measured potentials with unparalleled accuracy. We then calculate the potentials of all redox pairs that can be generated from biochemically relevant compounds and highlight fundamental trends in redox biochemistry. We further address the question of why NAD/NADP are used as primary electron carriers, demonstrating how their physiological potential range fits the reactions of central metabolism and minimizes the concentration of reactive carbonyls. The use of quantum chemistry can revolutionize our understanding of biochemical phenomena by enabling fast and accurate calculation of thermodynamic values.

Methotrexate (MTX) is widely used for the treatment of childhood acute lymphoblastic leukemia (ALL). The accumulation of MTX and its active metabolites, methotrexate polyglutamates (MTXPG), in ALL cells is an important determinant of its antileukemic effects. We studied 194 of 356 patients enrolled on St. Jude Total XV protocol for newly diagnosed ALL with the goal of characterizing the intracellular pharmacokinetics of MTXPG in leukemia cells; relating these pharmacokinetics to ALL lineage, ploidy and molecular subtype; and using a folate pathway model to simulate optimal treatment strategies. Serial MTX concentrations were measured in plasma and intracellular MTXPG concentrations were measured in circulating leukemia cells. A pharmacokinetic model was developed which accounted for the plasma disposition of MTX along with the transport and metabolism of MTXPG. In addition, a folate pathway model was adapted to simulate the effects of treatment strategies on the inhibition of de novo purine synthesis (DNPS). The intracellular MTXPG pharmacokinetic model parameters differed significantly by lineage, ploidy, and molecular subtypes of ALL. Folylpolyglutamate synthetase (FPGS) activity was higher in B vs T lineage ALL (p<0.005), MTX influx and FPGS activity were higher in hyperdiploid vs non-hyperdiploid ALL (p<0.03), MTX influx and FPGS activity were lower in the t(12;21) (ETV6-RUNX1) subtype (p<0.05), and the ratio of FPGS to γ-glutamyl hydrolase (GGH) activity was lower in the t(1;19) (TCF3-PBX1) subtype (p<0.03) than other genetic subtypes. In addition, the folate pathway model showed differential inhibition of DNPS relative to MTXPG accumulation, MTX dose, and schedule. This study has provided new insights into the intracellular disposition of MTX in leukemia cells and how it affects treatment efficacy.

We present a method to calculate the propensities of regions within a DNA molecule to transition from B-form to Z-form under negative superhelical stresses. We use statistical mechanics to analyze the competition that occurs among all susceptible Z-forming regions at thermodynamic equilibrium in a superhelically stressed DNA of specified sequence. This method, which we call SIBZ, is similar to the SIDD algorithm that was previously developed to analyze superhelical duplex destabilization. A state of the system is determined by assigning to each base pair either the B- or the Z-conformation, accounting for the dinucleotide repeat unit of Z-DNA. The free energy of a state is comprised of the nucleation energy, the sequence-dependent B-Z transition energy, and the energy associated with the residual superhelicity remaining after the change of twist due to transition. Using this information, SIBZ calculates the equilibrium B-Z transition probability of each base pair in the sequence. This can be done at any physiologically reasonable level of negative superhelicity. We use SIBZ to analyze a variety of representative genomic DNA sequences. We show that the dominant Z-DNA forming regions in a sequence can compete in highly complex ways as the superhelicity level changes. Despite having no tunable parameters, the predictions of SIBZ agree precisely with experimental results, both for the onset of transition in plasmids containing introduced Z-forming sequences and for the locations of Z-forming regions in genomic sequences. We calculate the transition profiles of 5 kb regions taken from each of 12,841 mouse genes and centered on the transcription start site (TSS). We find a substantial increase in the frequency of Z-forming regions immediately upstream from the TSS. The approach developed here has the potential to illuminate the occurrence of Z-form regions in vivo, and the possible roles this transition may play in biological processes.

The differential modulation of agonist and antagonist binding to opioid receptors (ORs) by sodium (Na+) has been known for decades. To shed light on the molecular determinants, thermodynamics, and kinetics of Na+ translocation through the μ-OR (MOR), we used a multi-ensemble Markov model framework combining equilibrium and non-equilibrium atomistic molecular dynamics simulations of Na+ binding to MOR active or inactive crystal structures embedded in an explicit lipid bilayer. We identify an energetically favorable, continuous ion pathway through the MOR active conformation only, and provide, for the first time: i) estimates of the energy differences and required timescales of Na+ translocation in inactive and active MORs, ii) estimates of Na+-induced changes to agonist binding validated by radioligand measurements, and iii) testable hypotheses of molecular determinants and correlated motions involved in this translocation, which are likely to play a key role in MOR signaling.

Extensive experimental information supports the formation of ligand-specific conformations of G protein-coupled receptors (GPCRs) as a possible molecular basis for their functional selectivity for signaling pathways. Taking advantage of the recently published inactive and active crystal structures of GPCRs, we have implemented an all-atom computational strategy that combines different adaptive biasing techniques to identify ligand-specific conformations along pre-determined activation pathways. Using the prototypic GPCR β2-adrenergic receptor as a suitable test case for validation, we show that ligands with different efficacies (either inverse agonists, neutral antagonists, or agonists) modulate the free-energy landscape of the receptor by shifting the conformational equilibrium towards active or inactive conformations depending on their elicited physiological response. Notably, we provide for the first time a quantitative description of the thermodynamics of the receptor in an explicit atomistic environment, which accounts for the receptor basal activity and the stabilization of different active-like states by differently potent agonists. Structural inspection of these metastable states reveals unique conformations of the receptor that may have been difficult to retrieve experimentally.

Group contribution (GC) methods are conventionally used in thermodynamics analysis of metabolic pathways to estimate the standard Gibbs energy change (ΔrG'o) of enzymatic reactions from limited experimental measurements. However, these methods are limited by their dependence on manually curated groups and inability to capture stereochemical information, leading to low reaction coverage. Herein, we introduce an automated molecular fingerprint-based thermodynamic analysis tool called dGPredictor that enables the consideration of stereochemistry within metabolite structures and thus increases reaction coverage. dGPredictor has comparable prediction accuracy compared to existing GC methods and can capture Gibbs energy changes for isomerase and transferase reactions, which exhibit no overall group changes. We also demonstrate dGPredictor's ability to predict the Gibbs energy change for novel reactions and seamless integration within de novo metabolic pathway design tools such as novoStoic for safeguarding against the inclusion of reaction steps with infeasible directionalities. To facilitate easy access to dGPredictor, we developed a graphical user interface to predict the standard Gibbs energy change for reactions at various pH and ionic strengths. The tool allows customized user input of known metabolites as KEGG IDs and novel metabolites as InChI strings (https://github.com/maranasgroup/dGPredictor).

Standard Gibbs energies of reactions are increasingly being used in metabolic modeling for applying thermodynamic constraints on reaction rates, metabolite concentrations and kinetic parameters. The increasing scope and diversity of metabolic models has led scientists to look for genome-scale solutions that can estimate the standard Gibbs energy of all the reactions in metabolism. Group contribution methods greatly increase coverage, albeit at the price of decreased precision. We present here a way to combine the estimations of group contribution with the more accurate reactant contributions by decomposing each reaction into two parts and applying one of the methods on each of them. This method gives priority to the reactant contributions over group contributions while guaranteeing that all estimations will be consistent, i.e. will not violate the first law of thermodynamics. We show that there is a significant increase in the accuracy of our estimations compared to standard group contribution. Specifically, our cross-validation results show an 80% reduction in the median absolute residual for reactions that can be derived by reactant contributions only. We provide the full framework and source code for deriving estimates of standard reaction Gibbs energy, as well as confidence intervals, and believe this will facilitate the wide use of thermodynamic data for a better understanding of metabolism.

Energy metabolism is central to cellular biology. Thus, genome-scale models of heterotrophic unicellular species must account appropriately for the utilization of external nutrients to synthesize energy metabolites such as ATP. However, metabolic models designed for flux-balance analysis (FBA) may contain thermodynamically impossible energy-generating cycles: without nutrient consumption, these models are still capable of charging energy metabolites (such as ADP→ATP or NADP+→NADPH). Here, we show that energy-generating cycles occur in over 85% of metabolic models without extensive manual curation, such as those contained in the ModelSEED and MetaNetX databases; in contrast, such cycles are rare in the manually curated models of the BiGG database. Energy generating cycles may represent model errors, e.g., erroneous assumptions on reaction reversibilities. Alternatively, part of the cycle may be thermodynamically feasible in one environment, while the remainder is thermodynamically feasible in another environment; as standard FBA does not account for thermodynamics, combining these into an FBA model allows erroneous energy generation. The presence of energy-generating cycles typically inflates maximal biomass production rates by 25%, and may lead to biases in evolutionary simulations. We present efficient computational methods (i) to identify energy generating cycles, using FBA, and (ii) to identify minimal sets of model changes that eliminate them, using a variant of the GlobalFit algorithm.

Many proteins comprising of complex topologies require molecular chaperones to achieve their unique three-dimensional folded structure. The E.coli chaperone, GroEL binds with a large number of unfolded and partially folded proteins, to facilitate proper folding and prevent misfolding and aggregation. Although the major structural components of GroEL are well defined, scaffolds of the non-native substrates that determine chaperone-mediated folding have been difficult to recognize. Here we performed all-atomistic and replica-exchange molecular dynamics simulations to dissect non-native ensemble of an obligate GroEL folder, DapA. Thermodynamics analyses of unfolding simulations revealed populated intermediates with distinct structural characteristics. We found that surface exposed hydrophobic patches are significantly increased, primarily contributed from native and non-native β-sheet elements. We validate the structural properties of these conformers using experimental data, including circular dichroism (CD), 1-anilinonaphthalene-8-sulfonic acid (ANS) binding measurements and previously reported hydrogen-deutrium exchange coupled to mass spectrometry (HDX-MS). Further, we constructed network graphs to elucidate long-range intra-protein connectivity of native and intermediate topologies, demonstrating regions that serve as central "hubs". Overall, our results implicate that genomic variations (or mutations) in the distinct regions of protein structures might disrupt these topological signatures disabling chaperone-mediated folding, leading to formation of aggregates.

Kinetically stable proteins, those whose stability is derived from their slow unfolding kinetics and not thermodynamics, are examples of evolution's best attempts at suppressing unfolding. Especially in highly proteolytic environments, both partially and fully unfolded proteins face potential inactivation through degradation and/or aggregation, hence, slowing unfolding can greatly extend a protein's functional lifetime. The prokaryotic serine protease alpha-lytic protease (alphaLP) has done just that, as its unfolding is both very slow (t(1/2) approximately 1 year) and so cooperative that partial unfolding is negligible, providing a functional advantage over its thermodynamically stable homologs, such as trypsin. Previous studies have identified regions of the domain interface as critical to alphaLP unfolding, though a complete description of the unfolding pathway is missing. In order to identify the alphaLP unfolding pathway and the mechanism for its extreme cooperativity, we performed high temperature molecular dynamics unfolding simulations of both alphaLP and trypsin. The simulated alphaLP unfolding pathway produces a robust transition state ensemble consistent with prior biochemical experiments and clearly shows that unfolding proceeds through a preferential disruption of the domain interface. Through a novel method of calculating unfolding cooperativity, we show that alphaLP unfolds extremely cooperatively while trypsin unfolds gradually. Finally, by examining the behavior of both domain interfaces, we propose a model for the differential unfolding cooperativity of alphaLP and trypsin involving three key regions that differ between the kinetically stable and thermodynamically stable classes of serine proteases.

The role intrinsic statistical fluctuations play in creating avalanches--patterns of complex bursting activity with scale-free properties--is examined in leaky Markovian networks. Using this broad class of models, we develop a probabilistic approach that employs a potential energy landscape perspective coupled with a macroscopic description based on statistical thermodynamics. We identify six important thermodynamic quantities essential for characterizing system behavior as a function of network size: the internal potential energy, entropy, free potential energy, internal pressure, pressure, and bulk modulus. In agreement with classical phase transitions, these quantities evolve smoothly as a function of the network size until a critical value is reached. At that value, a discontinuity in pressure is observed that leads to a spike in the bulk modulus demarcating loss of thermodynamic robustness. We attribute this novel result to a reallocation of the ground states (global minima) of the system's stationary potential energy landscape caused by a noise-induced deformation of its topographic surface. Further analysis demonstrates that appreciable levels of intrinsic noise can cause avalanching, a complex mode of operation that dominates system dynamics at near-critical or subcritical network sizes. Illustrative examples are provided using an epidemiological model of bacterial infection, where avalanching has not been characterized before, and a previously studied model of computational neuroscience, where avalanching was erroneously attributed to specific neural architectures. The general methods developed here can be used to study the emergence of avalanching (and other complex phenomena) in many biological, physical and man-made interaction networks.

We uncovered the universal statistical laws for the biomolecular recognition/binding process. We quantified the statistical energy landscapes for binding, from which we can characterize the distributions of the binding free energy (affinity), the equilibrium constants, the kinetics and the specificity by exploring the different ligands binding with a particular receptor. The results of the analytical studies are confirmed by the microscopic flexible docking simulations. The distribution of binding affinity is Gaussian around the mean and becomes exponential near the tail. The equilibrium constants of the binding follow a log-normal distribution around the mean and a power law distribution in the tail. The intrinsic specificity for biomolecular recognition measures the degree of discrimination of native versus non-native binding and the optimization of which becomes the maximization of the ratio of the free energy gap between the native state and the average of non-native states versus the roughness measured by the variance of the free energy landscape around its mean. The intrinsic specificity obeys a Gaussian distribution near the mean and an exponential distribution near the tail. Furthermore, the kinetics of binding follows a log-normal distribution near the mean and a power law distribution at the tail. Our study provides new insights into the statistical nature of thermodynamics, kinetics and function from different ligands binding with a specific receptor or equivalently specific ligand binding with different receptors. The elucidation of distributions of the kinetics and free energy has guiding roles in studying biomolecular recognition and function through small-molecule evolution and chemical genetics.

The spontaneous dissociation of six small ligands from the active site of FKBP (the FK506 binding protein) is investigated by explicit water molecular dynamics simulations and network analysis. The ligands have between four (dimethylsulphoxide) and eleven (5-diethylamino-2-pentanone) non-hydrogen atoms, and an affinity for FKBP ranging from 20 to 0.2 mM. The conformations of the FKBP/ligand complex saved along multiple trajectories (50 runs at 310 K for each ligand) are grouped according to a set of intermolecular distances into nodes of a network, and the direct transitions between them are the links. The network analysis reveals that the bound state consists of several subbasins, i.e., binding modes characterized by distinct intermolecular hydrogen bonds and hydrophobic contacts. The dissociation kinetics show a simple (i.e., single-exponential) time dependence because the unbinding barrier is much higher than the barriers between subbasins in the bound state. The unbinding transition state is made up of heterogeneous positions and orientations of the ligand in the FKBP active site, which correspond to multiple pathways of dissociation. For the six small ligands of FKBP, the weaker the binding affinity the closer to the bound state (along the intermolecular distance) are the transition state structures, which is a new manifestation of Hammond behavior. Experimental approaches to the study of fragment binding to proteins have limitations in temporal and spatial resolution. Our network analysis of the unbinding simulations of small inhibitors from an enzyme paints a clear picture of the free energy landscape (both thermodynamics and kinetics) of ligand unbinding.

Bacterial growth depends crucially on metabolic fluxes, which are limited by the cell's capacity to maintain metabolic enzymes. The necessary enzyme amount per unit flux is a major determinant of metabolic strategies both in evolution and bioengineering. It depends on enzyme parameters (such as kcat and KM constants), but also on metabolite concentrations. Moreover, similar amounts of different enzymes might incur different costs for the cell, depending on enzyme-specific properties such as protein size and half-life. Here, we developed enzyme cost minimization (ECM), a scalable method for computing enzyme amounts that support a given metabolic flux at a minimal protein cost. The complex interplay of enzyme and metabolite concentrations, e.g. through thermodynamic driving forces and enzyme saturation, would make it hard to solve this optimization problem directly. By treating enzyme cost as a function of metabolite levels, we formulated ECM as a numerically tractable, convex optimization problem. Its tiered approach allows for building models at different levels of detail, depending on the amount of available data. Validating our method with measured metabolite and protein levels in E. coli central metabolism, we found typical prediction fold errors of 4.1 and 2.6, respectively, for the two kinds of data. This result from the cost-optimized metabolic state is significantly better than randomly sampled metabolite profiles, supporting the hypothesis that enzyme cost is important for the fitness of E. coli. ECM can be used to predict enzyme levels and protein cost in natural and engineered pathways, and could be a valuable computational tool to assist metabolic engineering projects. Furthermore, it establishes a direct connection between protein cost and thermodynamics, and provides a physically plausible and computationally tractable way to include enzyme kinetics into constraint-based metabolic models, where kinetics have usually been ignored or oversimplified.

One of the paramount goals of synthetic biology is to have the ability to tune transcriptional networks to targeted levels of expression at will. As a step in that direction, we have constructed a set of 18 unique binding sites for E. coli RNA Polymerase (RNAP) δ⁷⁰ holoenzyme, designed using a model of sequence-dependent binding energy combined with a thermodynamic model of transcription to produce a targeted level of gene expression. This promoter set allows us to determine the correspondence between the absolute numbers of mRNA molecules or protein products and the predicted promoter binding energies measured in k(B)T energy units. These binding sites adhere on average to the predicted level of gene expression over 3 orders of magnitude in constitutive gene expression, to within a factor of 3 in both protein and mRNA copy number. With these promoters in hand, we then place them under the regulatory control of a bacterial repressor and show that again there is a strict correspondence between the measured and predicted levels of expression, demonstrating the transferability of the promoters to an alternate regulatory context. In particular, our thermodynamic model predicts the expression from our promoters under a range of repressor concentrations between several per cell up to over 100 per cell. After correcting the predicted polymerase binding strength using the data from the unregulated promoter, the thermodynamic model accurately predicts the expression for the simple repression strains to within 30%. Demonstration of modular promoter design, where parts of the circuit (such as RNAP/TF binding strength and transcription factor copy number) can be independently chosen from a stock list and combined to give a predictable result, has important implications as an engineering tool for use in synthetic biology.

Ligand binding sites in proteins are often localized to deeply buried cavities, inaccessible to bulk solvent. Yet, in many cases binding of cognate ligands occurs rapidly. An intriguing system is presented by the L99A cavity mutant of T4 Lysozyme (T4L L99A) that rapidly binds benzene (~106 M-1s-1). Although the protein has long served as a model system for protein thermodynamics and crystal structures of both free and benzene-bound T4L L99A are available, the kinetic pathways by which benzene reaches its solvent-inaccessible binding cavity remain elusive. The current work, using extensive molecular dynamics simulation, achieves this by capturing the complete process of spontaneous recognition of benzene by T4L L99A at atomistic resolution. A series of multi-microsecond unbiased molecular dynamics simulation trajectories unequivocally reveal how benzene, starting in bulk solvent, diffuses to the protein and spontaneously reaches the solvent inaccessible cavity of T4L L99A. The simulated and high-resolution X-ray derived bound structures are in excellent agreement. A robust four-state Markov model, developed using cumulative 60 μs trajectories, identifies and quantifies multiple ligand binding pathways with low activation barriers. Interestingly, none of these identified binding pathways required large conformational changes for ligand access to the buried cavity. Rather, these involve transient but crucial opening of a channel to the cavity via subtle displacements in the positions of key helices (helix4/helix6, helix7/helix9) leading to rapid binding. Free energy simulations further elucidate that these channel-opening events would have been unfavorable in wild type T4L. Taken together and via integrating with results from experiments, these simulations provide unprecedented mechanistic insights into the complete ligand recognition process in a buried cavity. By illustrating the power of subtle helix movements in opening up multiple pathways for ligand access, this work offers an alternate view of ligand recognition in a solvent-inaccessible cavity, contrary to the common perception of a single dominant pathway for ligand binding.