Structures of rhodopsin in complex with G-protein-coupled receptor kinase 1.

G-protein-coupled receptor (GPCR) kinases (GRKs) selectively phosphorylate activated GPCRs, thereby priming them for desensitization1. Although it is unclear how GRKs recognize these receptors2-4, a conserved region at the GRK N terminus is essential for this process5-8. Here we report a series of cryo-electron microscopy single-particle reconstructions of light-activated rhodopsin (Rho*) bound to rhodopsin kinase (GRK1), wherein the N terminus of GRK1 forms a helix that docks into the open cytoplasmic cleft of Rho*. The helix also packs against the GRK1 kinase domain and stabilizes it in an active configuration. The complex is further stabilized by electrostatic interactions between basic residues that are conserved in most GPCRs and acidic residues that are conserved in GRKs. We did not observe any density for the regulator of G-protein signalling homology domain of GRK1 or the C terminus of rhodopsin. Crosslinking with mass spectrometry analysis confirmed these results and revealed dynamic behaviour in receptor-bound GRK1 that would allow the phosphorylation of multiple sites in the receptor tail. We have identified GRK1 residues whose mutation augments kinase activity and crosslinking with Rho*, as well as residues that are involved in activation by acidic phospholipids. From these data, we present a general model for how a small family of protein kinases can recognize and be activated by hundreds of different GPCRs.

Structural and Molecular Mechanisms of Cytokine-Mediated Endocrine Resistance in Human Breast Cancer Cells.

Human breast cancers that exhibit high proportions of immune cells and elevated levels of pro-inflammatory cytokines predict poor prognosis. Here, we demonstrate that treatment of human MCF-7 breast cancer cells with pro-inflammatory cytokines results in ERα-dependent activation of gene expression and proliferation, in the absence of ligand or presence of 4OH-tamoxifen (TOT). Cytokine activation of ERα and endocrine resistance is dependent on phosphorylation of ERα at S305 in the hinge domain. Phosphorylation of S305 by IKKß establishes an ERα cistrome that substantially overlaps with the estradiol (E2)-dependent ERα cistrome. Structural analyses suggest that S305-P forms a charge-linked bridge with the C-terminal F domain of ERα that enables inter-domain communication and constitutive activity from the N-terminal coactivator-binding site, revealing the structural basis of endocrine resistance. ERα therefore functions as a transcriptional effector of cytokine-induced IKKß signaling, suggesting a mechanism through which the tumor microenvironment controls tumor progression and endocrine resistance.

Non-active Site Residue in Loop L4 Alters Substrate Capture and Product Release in d-Arginine Dehydrogenase.

Numerous studies demonstrate that enzymes undergo multiple conformational changes during catalysis. The malleability of enzymes forms the basis for allosteric regulation: residues located far from the active site can exert long-range dynamical effects on the active site residues to modulate catalysis. The structure of Pseudomonas aeruginosa d-arginine dehydrogenase (PaDADH) shows four loops (L1, L2, L3, and L4) that span the substrate and the FAD-binding domains. Loop L4 comprises residues 329-336, spanning over the flavin cofactor. The I335 residue on loop L4 is ∼10 Šaway from the active site and ∼3.8 Šfrom N(1)-C(2)═O atoms of the flavin. In this study, we used molecular dynamics and biochemical techniques to investigate the effect of the mutation of I335 to histidine on the catalytic function of PaDADH. Molecular dynamics showed that the conformational dynamics of PaDADH are shifted to a more closed conformation in the I335H variant. In agreement with an enzyme that samples more in a closed conformation, the kinetic data of the I335H variant showed a 40-fold decrease in the rate constant of substrate association (k1), a 340-fold reduction in the rate constant of substrate dissociation from the enzyme-substrate complex (k2), and a 24-fold decrease in the rate constant of product release (k5), compared to that of the wild-type. Surprisingly, the kinetic data are consistent with the mutation having a negligible effect on the reactivity of the flavin. Altogether, the data indicate that the residue at position 335 has a long-range dynamical effect on the catalytic function in PaDADH.

Evolving Mutational Buildup in HIV-1 Protease Shifts Conformational Dynamics to Gain Drug Resistance.

Drug resistance in antiviral treatments is a serious public health problem. Viral proteins mutate very fast, giving them a way to escape drugs by lowering drug binding affinity but with compromised function. Human immunodeficiency virus type I (HIV-1) protease, a critical antiretroviral therapeutic target, represents a model for such viral regulation under inhibition. Drug inhibitors of HIV-1 protease lose effectiveness as the protein evolves through several variants to become more resistant. However, the detailed mechanism of drug resistance in HIV-1 protease is still unclear. Here, we test the hypothesis that mutations throughout the protease alter the protein conformational ensemble to weaken protein-inhibitor binding, resulting in an inefficient protease but still viable virus. Comparing conformational ensembles between variants and the wild type helps detect these function-related dynamical changes. All analyses of over 30 µs simulations converge to the conclusion that conformational dynamics of more drug-resistant variants are more different from that of the wild type. Distinct roles of mutations during viral evolution are discussed, including a mutation predominantly contributing to the increase of drug resistance and a mutation that is responsible (synergistically) for restoring catalytic efficiency. Drug resistance is mainly due to altered flap dynamics that hinder the access to the active site. The mutant variant showing the highest drug resistance has the most ″collapsed″ active-site pocket and hence the largest magnitude of hindrance of drug binding. An enhanced difference contact network community analysis is applied to understand allosteric communications. The method summarizes multiple conformational ensembles in one community network and can be used in future studies to detect function-related dynamics in proteins.

p53 Is Potentially Regulated by Cyclophilin D in the Triple-Proline Loop of the DNA Binding Domain.

The multifunctional protein p53 is the central molecular sensor of cellular stresses. The canonical function of p53 is to transcriptionally activate target genes in response to, for example, DNA damage that may trigger apoptosis. Recently, p53 was also found to play a role in the regulation of necrosis, another type of cell death featured by the mitochondrial permeability transition (mPT). In this process, p53 directly interacts with the mPT regulator cyclophilin D, the detailed mechanism of which however remains poorly understood. Here, we report a comprehensive computational investigation of the p53-cyclophilin D interaction using molecular dynamics simulations and associated analyses. We have identified the specific cyclophilin D binding site on p53 that is located at proline 151 in the DNA binding domain. As a peptidyl-prolyl isomerase, cyclophilin D binds p53 and catalyzes the cis-trans isomerization of the peptide bond preceding proline 151. We have also characterized the effect of such an isomerization and found that the p53 domain in the cis state is overall more rigid than the trans state except for the local region around proline 151. Dynamical changes upon isomerization occur in both local and distal regions, indicating an allosteric effect elicited by the isomerization. We present potential allosteric communication pathways between proline 151 and distal sites, including the DNA binding surface. Our work provides, for the first time, a model for how cyclophilin D binds p53 and regulates its activity by switching the configuration of a specific site.

Tuning Protein Dynamics to Sense Rapid Endoplasmic-Reticulum Calcium Dynamics.

Multi-scale calcium (Ca2+ ) dynamics, exhibiting wide-ranging temporal kinetics, constitutes a ubiquitous mode of signal transduction. We report a novel endoplasmic-reticulum (ER)-targeted Ca2+ indicator, R-CatchER, which showed superior kinetics in vitro (koff ≥2×103  s-1 , kon ≥7×106  M-1 s-1 ) and in multiple cell types. R-CatchER captured spatiotemporal ER Ca2+ dynamics in neurons and hotspots at dendritic branchpoints, enabled the first report of ER Ca2+ oscillations mediated by calcium sensing receptors (CaSRs), and revealed ER Ca2+ -based functional cooperativity of CaSR. We elucidate the mechanism of R-CatchER and propose a principle to rationally design genetically encoded Ca2+ indicators with a single Ca2+ -binding site and fast kinetics by tuning rapid fluorescent-protein dynamics and the electrostatic potential around the chromophore. The design principle is supported by the development of G-CatchER2, an upgrade of our previous (G-)CatchER with improved dynamic range. Our work may facilitate protein design, visualizing Ca2+ dynamics, and drug discovery.

Synergistic mutations in soluble guanylyl cyclase (sGC) reveal a key role for interfacial regions in the sGC activation mechanism.

Soluble guanylyl cyclase (sGC) is the main receptor for nitric oxide (NO) and a central component of the NO-cGMP pathway, critical to cardiovascular function. NO binding to the N-terminal sensor domain in sGC enhances the cyclase activity of the C-terminal catalytic domain. Our understanding of the structural elements regulating this signaling cascade is limited, hindering structure-based drug design efforts that target sGC to improve the management of cardiovascular diseases. Conformational changes are thought to propagate the NO-binding signal throughout the entire sGC heterodimer, via its coiled-coil domain, to reorient the catalytic domain into an active conformation. To identify the structural elements involved in this signal transduction cascade, here we optimized a cGMP-based luciferase assay that reports on heterologous sGC activity in Escherichia coli and identified several mutations that activate sGC. These mutations resided in the dorsal flaps, dimer interface, and GTP-binding regions of the catalytic domain. Combinations of mutations from these different elements synergized, resulting in even greater activity and indicating a complex cross-talk among these regions. Molecular dynamics simulations further revealed conformational changes underlying the functional impact of these mutations. We propose that the interfacial residues play a central role in the sGC activation mechanism by coupling the coiled-coil domain to the active site via a series of hot spots. Our results provide new mechanistic insights not only into the molecular pathway for sGC activation but also for other members of the larger nucleotidyl cyclase family.

Detecting Functional Dynamics in Proteins with Comparative Perturbed-Ensembles Analysis.

Recent advances have made all-atom molecular dynamics (MD) a powerful tool to sample the conformational energy landscape. There are still however three major challenges in the application of MD to biological systems: accuracy of force field, time scale, and the analysis of simulation trajectories. Significant progress in addressing the first two challenges has been made and extensively reviewed previously. This Account focuses on strategies of analyzing simulation data of biomolecules that also covers ways to properly design simulations and validate simulation results. In particular, we examine an approach named comparative perturbed-ensembles analysis, which we developed to efficiently detect dynamics in protein MD simulations that can be linked to biological functions. In our recent studies, we implemented this approach to understand allosteric regulations in several disease-associated human proteins. The central task of a comparative perturbed-ensembles analysis is to compare two or more conformational ensembles of a system generated by MD simulations under distinct perturbation conditions. Perturbations can be different sequence variations, ligand-binding conditions, and other physical/chemical modifications of the system. Each simulation is long enough (e.g., microsecond-long) to ensure sufficient sampling of the local substate. Then, sophisticated bioinformatic and statistical tools are applied to extract function-related information from the simulation data, including principal component analysis, residue-residue contact analysis, difference contact network analysis (dCNA) based on the graph theory, and statistical analysis of side-chain conformations. Computational findings are further validated with experimental data. By comparing distinct conformational ensembles, functional micro- to millisecond dynamics can be inferred. In contrast, such a time scale is difficult to reach in a single simulation; even when reached for a single condition of a system, it is elusive as to what dynamical motions are related to functions without, for example, comparing free and substrate-bound proteins at the minimum. We illustrate our approach with three examples. First, we discuss using the approach to identify allosteric pathways in cyclophilin A (CypA), a member of a ubiquitous class of peptidyl-prolyl cis-trans isomerase enzymes. By comparing side-chain torsion-angle distributions of CypA in wild-type and mutant forms, we identified three pathways: two are consistent with recent nuclear magnetic resonance experiments, whereas the third is a novel pathway. Second, we show how the approach enables a dynamical-evolution analysis of the human cyclophilin family. In the analysis, both conserved and divergent conformational dynamics across three cyclophilin isoforms (CypA, CypD, and CypE) were summarized. The conserved dynamics led to the discovery of allosteric networks resembling those found in CypA. A residue wise determinant underlying the unique dynamics in CypD was also detected and validated with additional mutational MD simulations. In the third example, we applied the approach to elucidate a peptide sequence-dependent allosteric mechanism in human Pin 1, a phosphorylation-dependent peptidyl-prolyl isomerase. We finally present our outlook of future directions. Especially, we envisage how the approach could help open a new avenue in drug discovery.

Comparative structural dynamic analysis of GTPases.

GTPases regulate a multitude of essential cellular processes ranging from movement and division to differentiation and neuronal activity. These ubiquitous enzymes operate by hydrolyzing GTP to GDP with associated conformational changes that modulate affinity for family-specific binding partners. There are three major GTPase superfamilies: Ras-like GTPases, heterotrimeric G proteins and protein-synthesizing GTPases. Although they contain similar nucleotide-binding sites, the detailed mechanisms by which these structurally and functionally diverse superfamilies operate remain unclear. Here we compare and contrast the structural dynamic mechanisms of each superfamily using extensive molecular dynamics (MD) simulations and subsequent network analysis approaches. In particular, dissection of the cross-correlations of atomic displacements in both the GTP and GDP-bound states of Ras, transducin and elongation factor EF-Tu reveals analogous dynamic features. This includes similar dynamic communities and subdomain structures (termed lobes). For all three proteins the GTP-bound state has stronger couplings between equivalent lobes. Network analysis further identifies common and family-specific residues mediating the state-specific coupling of distal functional sites. Mutational simulations demonstrate how disrupting these couplings leads to distal dynamic effects at the nucleotide-binding site of each family. Collectively our studies extend current understanding of GTPase allosteric mechanisms and highlight previously unappreciated similarities across functionally diverse families.

Establishing a Framework of Using Residue-Residue Interactions in Protein Difference Network Analysis.

Detailed understanding of interactions between amino acid residues is critical in using promising difference network analysis approaches to map allosteric communication pathways. Using experimental data as benchmarks, we scan values of two essential residue-residue contact parameters: the distance cutoff (dc) and the cutoff of residue separation in sequence (nc). The optimal dc = 4.5 Å is revealed, which defines the upper bound of the first shell of residue-residue packing in proteins, whereas nc is found to have little effects on performance. We also develop a new energy-based contact method for network analyses and find an equivalency between the energy network using the optimal energy cutoff ec = 1.0 kBT and the structure network using dc = 4.5 Å. The simple 4.5-Å contact method is further shown to have comparable prediction accuracy to a contact method using amino acid type-specific distance cutoffs and chemical shift prediction-based methods. This study provides necessary tools in mapping dynamics to functions.

Navigating the conformational landscape of G protein-coupled receptor kinases during allosteric activation.

G protein-coupled receptors (GPCRs) are essential for transferring extracellular signals into carefully choreographed intracellular responses controlling diverse aspects of cell physiology. The duration of GPCR-mediated signaling is primarily regulated via GPCR kinase (GRK)-mediated phosphorylation of activated receptors. Although many GRK structures have been reported, the mechanisms underlying GRK activation are not well-understood, in part because it is unknown how these structures map to the conformational landscape available to this enzyme family. Unlike most other AGC kinases, GRKs rely on their interaction with GPCRs for activation and not phosphorylation. Here, we used principal component analysis of available GRK and protein kinase A crystal structures to identify their dominant domain motions and to provide a framework that helps evaluate how close each GRK structure is to being a catalytically competent state. Our results indicated that disruption of an interface formed between the large lobe of the kinase domain and the regulator of G protein signaling homology domain (RHD) is highly correlated with establishment of the active conformation. By introducing point mutations in the GRK5 RHD-kinase domain interface, we show with both in silico and in vitro experiments that perturbation of this interface leads to higher phosphorylation activity. Navigation of the conformational landscape defined by this bioinformatics-based study is likely common to all GPCR-activated GRKs.

Elucidating Allosteric Communications in Proteins with Difference Contact Network Analysis.

A difference contact network analysis (dCNA) method is developed for delineating allosteric mechanisms in proteins. The new method addresses limitations of conventional network analysis methods and is particularly suitable for allosteric systems undergoing large-amplitude conformational changes during function. Tests show that dCNA works well for proteins of varying sizes and functions. The design of dCNA is general enough to facilitate analyses of diverse dynamic data generated by molecular dynamics, crystallography, or nuclear magnetic resonance.

Structural Determinants Influencing the Potency and Selectivity of Indazole-Paroxetine Hybrid G Protein-Coupled Receptor Kinase 2 Inhibitors.

G protein-coupled receptor kinases (GRKs) phosphorylate activated receptors to promote arrestin binding, decoupling from heterotrimeric G proteins, and internalization. GRK2 and GRK5 are overexpressed in the failing heart and thus have become therapeutic targets. Previously, we discovered two classes of GRK2-selective inhibitors, one stemming from GSK180736A, a Rho-associated coiled-coil containing kinase 1 (ROCK1) inhibitor, the other from paroxetine, a selective serotonin-reuptake inhibitor. These two classes of compounds bind to the GRK2 active site in a similar configuration but contain different hinge-binding "warheads": indazole and benzodioxole, respectively. We surmised from our prior studies that an indazole would be the stronger hinge binder and would impart increased potency when substituted for benzodioxole in paroxetine derivatives. To test this hypothesis, we synthesized a series of hybrid compounds that allowed us to compare the effects of inhibitors that differ only in the identity of the warhead. The indazole-paroxetine analogs were indeed more potent than their respective benzodioxole derivatives but lost selectivity. To investigate how these two warheads dictate selectivity, we determined the crystal structures of three of the indazole hybrid compounds (CCG224061, CCG257284, and CCG258748) in complex with GRK2-Gßγ Comparison of these structures with those of analogous benzodioxole-containing complexes confirmed that the indazole-paroxetine hybrids form stronger interactions with the hinge of the kinase but also stabilize a distinct conformation of the kinase domain of GRK2 compared with previous complexes with paroxetine analogs. This conformation is analogous to one that can be assumed by GRK5, at least partially explaining the loss in selectivity.

Dynamic Coupling and Allosteric Networks in the α Subunit of Heterotrimeric G Proteins.

G protein α subunits cycle between active and inactive conformations to regulate a multitude of intracellular signaling cascades. Important structural transitions occurring during this cycle have been characterized from extensive crystallographic studies. However, the link between observed conformations and the allosteric regulation of binding events at distal sites critical for signaling through G proteins remain unclear. Here we describe molecular dynamics simulations, bioinformatics analysis, and experimental mutagenesis that identifies residues involved in mediating the allosteric coupling of receptor, nucleotide, and helical domain interfaces of Gαi. Most notably, we predict and characterize novel allosteric decoupling mutants, which display enhanced helical domain opening, increased rates of nucleotide exchange, and constitutive activity in the absence of receptor activation. Collectively, our results provide a framework for explaining how binding events and mutations can alter internal dynamic couplings critical for G protein function.

Online interactive analysis of protein structure ensembles with Bio3D-web.

Bio3D-web is an online application for analyzing the sequence, structure and conformational heterogeneity of protein families. Major functionality is provided for identifying protein structure sets for analysis, their alignment and refined structure superposition, sequence and structure conservation analysis, mapping and clustering of conformations and the quantitative comparison of their predicted structural dynamics. AVAILABILITY: Bio3D-web is based on the Bio3D and Shiny R packages. All major browsers are supported and full source code is available under a GPL2 license from http://thegrantlab.org/bio3d-web CONTACT: bjgrant@umich.edu or lars.skjarven@uib.no.

Mapping the Processivity Determinants of the Kinesin-3 Motor Domain.

Kinesin superfamily members play important roles in many diverse cellular processes, including cell motility, cell division, intracellular transport, and regulation of the microtubule cytoskeleton. How the properties of the family-defining motor domain of distinct kinesins are tailored to their different cellular roles remains largely unknown. Here, we employed molecular-dynamics simulations coupled with energetic calculations to infer the family-specific interactions of kinesin-1 and kinesin-3 motor domains with microtubules in different nucleotide states. We then used experimental mutagenesis and single-molecule motility assays to further assess the predicted residue-wise determinants of distinct kinesin-microtubule binding properties. Collectively, our results identify residues in the L8, L11, and α6 regions that contribute to family-specific microtubule interactions and whose mutation affects motor-microtubule complex stability and processive motility (the ability of an individual motor to take multiple steps along its microtubule filament). In particular, substitutions of prominent kinesin-3 residues with those found in kinesin-1, namely, R167S/H171D, K266D, and R346M, were found to decrease kinesin-3 processivity 10-fold and thus approach kinesin-1 levels.

Integrating protein structural dynamics and evolutionary analysis with Bio3D.

BACKGROUND: Popular bioinformatics approaches for studying protein functional dynamics include comparisons of crystallographic structures, molecular dynamics simulations and normal mode analysis. However, determining how observed displacements and predicted motions from these traditionally separate analyses relate to each other, as well as to the evolution of sequence, structure and function within large protein families, remains a considerable challenge. This is in part due to the general lack of tools that integrate information of molecular structure, dynamics and evolution. RESULTS: Here, we describe the integration of new methodologies for evolutionary sequence, structure and simulation analysis into the Bio3D package. This major update includes unique high-throughput normal mode analysis for examining and contrasting the dynamics of related proteins with non-identical sequences and structures, as well as new methods for quantifying dynamical couplings and their residue-wise dissection from correlation network analysis. These new methodologies are integrated with major biomolecular databases as well as established methods for evolutionary sequence and comparative structural analysis. New functionality for directly comparing results derived from normal modes, molecular dynamics and principal component analysis of heterogeneous experimental structure distributions is also included. We demonstrate these integrated capabilities with example applications to dihydrofolate reductase and heterotrimeric G-protein families along with a discussion of the mechanistic insight provided in each case. CONCLUSIONS: The integration of structural dynamics and evolutionary analysis in Bio3D enables researchers to go beyond a prediction of single protein dynamics to investigate dynamical features across large protein families. The Bio3D package is distributed with full source code and extensive documentation as a platform independent R package under a GPL2 license from http://thegrantlab.org/bio3d/ .

Dissecting the Allosteric Fine-Tuning of Enzyme Catalysis.

Fully understanding the mechanism of allosteric regulation in biomolecules requires separating and examining all of the involved factors. In enzyme catalysis, allosteric effector binding shifts the structure and dynamics of the active site, leading to modified energetic (e.g., energy barrier) and dynamical (e.g., diffusion coefficient) factors underlying the catalyzed reaction rate. Such modifications can be subtle and dependent on the type of allosteric effector, representing a fine-tuning of protein function. The microscopic description of allosteric regulation at the level of function-dictating factors has prospective applications in fundamental and pharmaceutical sciences, which is, however, largely missing so far. Here, we characterize the allosteric fine-tuning of enzyme catalysis, using human Pin1 as an example, by performing more than half-millisecond all-atom molecular dynamics simulations. Changes of reaction kinetics and the dictating factors, including the free energy surface along the reaction coordinate and the diffusion coefficient of the reaction dynamics, under various enzyme and allosteric effector binding conditions are examined. Our results suggest equal importance of the energetic and dynamical factors, both of which can be modulated allosterically, and the combined effect determines the final allosteric output. We also reveal the potential dynamic basis for allosteric modulation using an advanced statistical technique to detect function-related conformational dynamics. Methods developed in this work can be applied to other allosteric systems.

Domain-opening and dynamic coupling in the α-subunit of heterotrimeric G proteins.

Heterotrimeric G proteins are conformational switches that turn on intracellular signaling cascades in response to the activation of G-protein-coupled receptors. Receptor activation by extracellular stimuli promotes a cycle of GTP binding and hydrolysis on the G protein α-subunit (Gα). Important conformational transitions occurring during this cycle have been characterized from extensive crystallographic studies of Gα. However, the link between the observed conformations and the mechanisms involved in G-protein activation and effector interaction remain unclear. Here we describe a comprehensive principal component analysis of available Gα crystallographic structures supplemented with extensive unbiased conventional and accelerated molecular dynamics simulations that together characterize the response of Gα to GTP binding and hydrolysis. Our studies reveal details of activating conformational changes as well as the intrinsic flexibility of the α-helical domain that includes a large-scale 60° domain opening under nucleotide-free conditions. This result is consistent with the recently reported open crystal structure of Gs, the stimulatory G protein for adenylyl cyclase, in complex with the α2 adrenergic receptor. Sets of unique interactions potentially important for the conformational transition are also identified. Moreover simulations reveal nucleotide-dependent dynamical couplings of distal regions and residues potentially important for the allosteric link between functional sites.

Drug uptake pathways of multidrug transporter AcrB studied by molecular simulations and site-directed mutagenesis experiments.

Multidrug resistance has been a critical issue in current chemotherapy. In Escherichia coli , a major efflux pump responsible for the multidrug resistance contains a transporter AcrB. Crystallographic studies and mutational assays of AcrB provided much of structural and overall functional insights, which led to the functionally rotating mechanism. However, the drug uptake pathways are somewhat controversial because at least two possible pathways, the vestibule and the cleft paths, were suggested. Here, combining molecular simulations and site-directed mutagenesis experiments, we addressed the uptake mechanism finding that the drug uptake pathways can be significantly different depending on the properties of drugs. First, in the computational free energy analysis of drug movements along AcrB tunnels, we found a ligand-dependent drug uptake mechanism. With the same molecular sizes, drugs that are both strongly hydrophobic and lipophilic were preferentially taken in via the vestibule path, while other drugs favored the cleft path. Second, direct simulations realized totally about 3500 events of drug uptake by AcrB for a broad range of drug property. These simulations confirmed the ligand-dependent drug uptake and further suggested that a smaller drug favors the vestibule path, while a larger one is taken in via the cleft path. Moreover, the direct simulations identified an alternative uptake path which is not visible in the crystal structure. Third, site-directed mutagenesis of AcrB in E. coli verified that mutations of residues located along the newly identified path significantly reduced the efflux efficiency, supporting its relevance in in vivo function.