2.1 Cryo-EM densities help obtain a more reliable initial structure

We used the EmrE protein structure refined by Ovchinnikov et al. [28] as a starting point for further refinement. After placing the structure in a model bacterial membrane, the system was first equilibrated for 1 μs (see S1 Fig). Then, using correlation-based cryo-EM fitting, the protein was restrained to match its electron density map published previously [26]. In this process, a ‘forward’ model of the electron density is created by Gaussian-smoothing the atomistic structure, and forces are applied so as to minimize the cross-correlation between the simulated and actual density. Thanks to this, the spatial arrangement of the protein helices closely matches that of the PDB entry 3B5D (see S2 Fig), increasing the reliability of our initial structure. We then created a relaxed “apo” form by removing the ligand from this system and simulating it under equilibrium conditions for 1.75 μs. Recently, a new NMR-based structure was deposited (PDB 7JK8), and we report its similarity to existing models (PDB 3B5D, our refined structure and the structure obtained by Ovchinnikov V [37]) in S3 Fig.

2.2 Mechanism of ligand translocation is dictated by thermodynamics of binding

In past studies, the ligand translocation pathway of EmrE was broadly characterized in biochemical terms and at timescales accessible to NMR experiments, but these insights into the dynamics of binding were never directly explained on a structural level [30]. Hence, to gain a higher-resolution mechanistic understanding of the transport cycle of EmrE, we employed fully atomistic free energy simulations aimed at reconstruction of thermodynamics of binding at different protonation states of the central glutamates: E14^A and E14^B. The resulting consensus free energy profiles, showed in Fig 2B, are supported by 30 μs of sampling each, and obtained as an average of results from three independent umbrella sampling and metadynamics runs (see Methods and S4 Fig for details). From the depths and positions of free energy minima, the profiles allow for reconstruction of a sequence of events—including protonation changes and the “domain-swap” conformational transition—that lead to ligand translocation driven by proton gradient.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

(A) Probabilities of Grotthuss wire formation. For both glutamates, the fraction of frames in which a continuous wire was found between the carboxyl group and bulk solvent is plotted as a function between the ligand position and the probability calculated as an average from each umbrella sampling window. Note that “open” and “closed” side might refer to either cytoplasm or periplasm, depending on the current conformational state of EmrE. (B) Free energy profiles for ligand entry in the four possible protonation states. The order of affinities is consistent with the intuitive notion that each protonation event decreases affinity for the positively charged ligand. (C) The most thermodynamically plausible mechanism of proton-coupled antiport, postulated based on wire probabilities from panel A and protonation-dependent free energy profiles from panel B (note the dashed profiles are just reflected about Z = 0 and correspond to a state after the conformational change). The mechanism assumes E14^B protonates more eagerly in the ligand-bound state, as shown by experiments and simulations. Note that upon a symmetric transition, the pK_a of E14^B becomes the pK_a of E14^A and vice versa. For a dynamic version of the figure, see S1 Movie. For the discussion of an alternative pathway, cf. S5 Fig.

https://doi.org/10.1371/journal.pcbi.1009454.g002

Upon binding to a deprotonated EmrE dimer (solid orange line in Fig 2B–2D), encounters a series of metastable states, from the entry point at the bilayer surface level (z = 2.2 nm) through two local minima at z = 1.5 and 0.8 nm, to end up stably bound in the free energy basin next to the protein center, between z = -0.1 and 0.5 nm. From here, three different events can follow, with different mechanistic implications: (1) a conformational transition that swaps the domain labels and exposes the binding site to the periplasm; (2) protonation of E14^A; or (3) protonation of E14^B. To narrow down the scope of possibilities, we calculated pK_a values of the two glutamates in a ligand-bound state using equilibrium alchemical free energy simulations. This approach yielded a pK_a for E14^B higher by more than 2 units than that of E14^A (see S1 Table), in agreement with the assignment made in NMR experiments [30]. Even though absolute values of pK_a in proteins cannot be calculated with precision higher than 1.5–2 units [38], and our estimate of absolute pK_as is indeed off by at least 1 unit, this qualitative consistency shows that our model correctly captures the preferential protonation of E14^B.

To distinguish between the remaining two options—conformational change or protonation of E14^B—we calculated the fraction of Grotthuss water wires [39] connecting the glutamate side chains to bulk solvent from each side, cytoplasmic (high-pH) and periplasmic (low-pH), reasoning that E14^B will only be protonated first if it can access protons from the periplasmic side with a ligand bound. As shown in Fig 2A, this is not the case, since E14^B can freely exchange protons with the cytoplasmic side (dashed red line) but not with the periplasm (solid red line); hence the conformational change should occur as the next step, as in the pathway depicted in Fig 2C.

Surprisingly, however, E14^A turned out to be capable of exchanging protons with the periplasm in the presence of a bound ligand (grey solid line), suggesting that under specific conditions (pH of the periplasm lower than the pK_a of E14^A) this residue can still be protonated before the completion of the conformational change, as shown schematically in S5 Fig. As we shall show, both sequences of events lead to slightly distinct but mechanistically viable translocation pathways.

Assuming that the conformational change occurs right after ligand entry, the dimer is now open to the periplasm, and E14^B can easily accept a proton, an event corresponding to a jump from the orange to the lime free energy profile in Fig 2C. Accordingly, the equilibrium position of the ligand shifts towards exit (z = -0.6 nm), where acceptance of another proton is feasible, as shown by the dashed grey line of Fig 2A. Upon second protonation, the system transitions to the blue profile corresponding to a doubly protonated state, where only two minor free energy barriers, of ca. 2 kcal/mol each, separate from full dissociation of from the protein. This complete translocation pathway is illustrated in S1 Movie.

If in turn the proton is instead first accepted by E14^A at sufficiently low periplasmic pH, the ligand will preferentially progress towards more negative values of z (minimum at -0.2 nm, green profile in Fig 2B). Once the conformational change occurs, the symmetrical change of labels “converts” E14^A into E14^B, resulting in a situation identical to the one described above. We note that this dual pathway might represent an evolutionary optimization, allowing the protein to operate at different pH differences, as long as the lower periplasmic pH is capable of protonating both glutamates.

Notably, although the profiles capture the intuitive order of affinities in which each consecutive protonation event decreases the electrostatic attraction between the negatively charged binding site and the aromatic cation, the exact affinities systematically underestimate the experimentally determined values. For the deprotonated and doubly protonated variants, our profiles reach minima of -7 and -4 kcal/mol instead of -11 and -6 kcal/mol [30, 33], likely reflecting the fact that our free energy simulations do not reach the full apo state. However, we point out that it is the structural adaptations of the binding site and its neighborhood that confers efficiency and robustness of the mechanism. In particular, while a watertight dimer interface prohibits Grotthuss wire formation between the closed side and E14^B at all times (Fig 2A), E14^A alternates between connection to the closed (z<0.5 nm) and open (z<0.5 nm) side, thereby enabling adaptability while preventing futile proton transfer. This might address a similar concern regarding the presence of spurious “leaky” states in simulations performed by Vermaas et al. [27].

2.3 Binding site desolvation in the apo form prevents proton leaks

We show above that in the ligand-bound state, E14^A and E14^B are connected to two different proton reservoirs, explaining how EmrE prevents proton leaks in the ligand-bound state. However, a previous simulational study [28] found that the apo form maintains the Grotthuss-like connection between the binding site glutamates and the bulk solvent on a sub-microsecond scale, one that would raise the possibility of “futile” proton transport as the apo form undergoes spontaneous conformational changes.

To verify this finding, we performed a longer simulation of the apo form, extending them to almost 2 μs, and analyzed the resulting trajectory for the presence of continuous hydrogen-bonded water wires. This equilibrium simulation revealed that following the (unphysiologically abrupt) removal of the ligand, EmrE initially remains in a conformation close to the ligand-bound one. The initial reconfiguration of contacts is rather modest, and results in alternating changes in the accessibility of the binding site (Fig 3C). After 1.2 μs, however, a concerted rotameric change of W63^B and M21^B created a hydrophobic barrier around the binding site, so that all Grotthuss wires connected to E14^A and E14^B abruptly disappeared. After additional 0.5 μs, the binding site remained inaccessible to bulk water, with 8 water molecules stably trapped inside the central cavity (Fig 3D). This transition should be easily reversible when the hydrophobic ligand approaches the protein, indicating that the key structural changes between the apo and ligand-bound forms are fast enough to avoid unwanted proton leaks. We need to note, though, that the sole fact of measuring pK_a values for the apo form [40] implies that such spontaneous reopening has to be feasible (so that the protons can reach the central glutamates), albeit likely on a much slower timescale than the one relevant for transfer. Indeed, a comparison of ensemble distributions of interresidue distance measured using spin-label based DEER showed that the apo form has a much less defined structure, hinting at a more complex conformational landscape of the ligand-free form [41].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

(A) Presence of the Grotthuss wires during the 1.75 μs equilibrium simulations of the apo state. At 1.2 μs, the wires disappear permanently. Plotted as a function between the time of simulation and the probability calculated as a Gaussian-smoothed wire presence function. (B) A visualization of wires extending from the binding site to the open (orange) and closed (blue) side, with the ligand shown in pink and the two glutamates in yellow. (C,D) Loss of solvation of the ligand binding channel, shown as a comparison between the initial and final structures. (C) Disappearance of the water wires, with only several (8 in this case) water molecules remaining trapped in the binding site cavity. (D) Dominant role of the M21^B “plug” and the W63^B “valve” in the desolvation of the binding channel (see also S2 Movie).

https://doi.org/10.1371/journal.pcbi.1009454.g003

2.4 Solvent-exposed surface of EmrE confers selective affinity for charged aromatic species

EmrE exhibits broad specificity for aromatic cations, compounds whose antibacterial properties often stem from their non-specific interaction with DNA and/or RNA. To investigate the structural basis for this selectivity, as well as provide a clear mechanism of substrate recognition, we studied the behavior of membrane-embedded EmrE in the apo form in the presence of four different ligands, three consisting of four phenyl groups bound to a central atom (cationic , neutral SiPh₄ and anionic ) and one non-aromatic yet hydrophobic cation (). Using pure (protein-free) bilayer systems, we first verified that all charged ligands preferentially reside in the aqueous phase near the bilayer surface, while the neutral tetraphenylsilane shows high affinity for the bilayer interior (see S6 Fig). This indicates that the cationic substrates have to be bound directly from the cytoplasm-bilayer interface—a variant of the popular “molecular vacuum cleaner” mechanism [42], where instead of being captured from inside a membrane leaflet, the charged hydrophobic ligand preferentially adsorbs at the exposed surface of the membrane protein.

Indeed, in the presence of membrane-embedded EmrE, a large fraction of charged aromatic ligands accumulated at the water-protein interface, with a slight preference for the “open” side of the protein (Fig 4A and 4B). In turn, the non-aromatic hydrophobic cation () showed virtually no preference for the protein surface, even though it was found to be an exceptionally good substrate for transport [43]. This indicates that regardless of the ligand’s affinity for the central binding site, aromatic ligands are strongly selected for during initial binding, likely facilitating their expulsion even at low concentrations. Several low-affinity binding sites at the protein surface are directly responsible for this effect, mostly coinciding with surface clusters of aromatic or charged residues. In fact, the observed high ligand densities are almost entirely explainable by the presence of two aromatic residue clusters (Y4/Y6/Y53 and F23/W31/F78/F79/W76, the latter split in two at the “open” side of the dimer; yellow in Fig 4A) and/or charged residues (orange in Fig 4A) on the protein surface.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Binding preference of charged and/or aromatic ligands to the solvent-exposed surface of EmrE.

(A) 3D densities of the central atom near the solvent-exposed surface of EmrE, looking from the “open” (top) and “closed” (bottom) side. Aromatic side chains exposed at the surface are shown in yellow, charged in orange. In each set of simulation, 10 ligands were initially placed in the water phase and allowed to interact with the protein. (B) Ligand density as in (A), but projected on the z-axis. Note that entered into spaces between the protein and the lipid bilayer, occasionally achieving z<0.5 nm but outside of the main ligand-binding channel. (C) A sample structure of ligand-bound EmrE in which the ligand spontaneously passed from bulk solvent to z = 0.8 nm (see extended simulations in S7 Fig).

https://doi.org/10.1371/journal.pcbi.1009454.g004

While it is surprising that both and accumulated at the protein surface to a similar extent, in all replicas only the aromatic cation managed to properly enter the channel leading to the binding site, eventually reaching a distance of 0.8 nm from the center of the protein (S8 Fig). This hints at the presence of another selectivity filter, preventing aromatic anions from accessing the hydrophobic channel leading to the binding site even before they are repelled by the pair of negatively charged E14 residues.

However, no single charged or hydrophilic residue seems to single-handedly confer this preference, as the channel entrance is lined with hydrophobic side chains. To observe any differences between an anionic and cationic ligand in a more controlled setting, we extracted 3 channel-bound ligand structures (z-distance to binding site between 1.0 and 1.2 nm) from the equilibrium runs and performed separate 1 μs simulations of both and , running it “as it is” or swapping the charges but keeping the geometry (depending on which system the structure came from; see S7 Fig). Out of the 3 systems (yellow/orange lines in S7 Fig), all remained stably bound to the channel entrance, mostly entering deeper (up to 0.8 nm) into the channel (see Fig 4C and free energy minima in Fig 2B); simultaneously, the least strongly bound dissociated in one of the systems, while the other two progressed deeper into the binding channel similarly to their cationic counterparts (blue/purple lines in S7 Fig). We hence conclude that the selectivity at the channel entrance is weak and might instead stem from the initial partial desolvation of the ligand entering the channel, as would be suggested by the different costs of desolvation, known from experimental data to differ by more than 2 kcal/mol between and [44]. Further down the channel, selectivity for cations would be trivially conferred due to the electrostatic effect of deprotonated (and hence negatively charged) E14^A and E14^B; indeed, the effect of electrostatic attraction on the ligand affinity of the binding site is easily observable in the free energy analysis (cf. orange/blue profiles in Fig 2).

2.5 Ligand entrance and exit are associated with well-defined structural changes

Despite the effectively multi-μs timescale explored in our equilibrium spontaneous entry simulations, we did not observe a complete ligand entrance in which would sample a full pathway connecting the apo and ligand-bound states. To provide a realistic description of such complete transitions on the atomistic level, we turned to the Weighted Ensemble approach [45]—implemented in WESTPA—to generate unbiased trajectory ensembles of entrance into the unprotonated binding site, as well as release from the doubly protonated binding site. In doing so, we expected to observe key structural changes associated with ligand binding and release, and to identify possible intermediate steps along both pathways, potentially generalizable to other EmrE substrates.

Using the WESTPA workflow, we were able to generate successful ligand entry and escape trajectories with a length of up to 10 and 4 ns, respectively (see S3 Movie). Since these simulations represent physical but low-probability instances of such transitions, we ran equilibrium simulations from selected end state configurations to verify that the resulting configurations remain stable over time. In fact, in the shortest binding trajectories the protein clearly failed to adjust to the presence of the ligand, as the ligand quickly escaped the binding site in a plain equilibrium simulation, returning to an intermediate position at z = 0.8 nm. In subsequent trajectories, however, the ligand remained stably bound after entry for at least 100 ns. A visual inspection of differences between these two classes of trajectories indicate that reorganization of E14 and W63 is crucial for the stable accommodation of in the binding site, as shown in S8 Fig.

To identify structural changes accompanying changes in ligand occupancy in a systematic manner, we used the longest available entry/escape trajectories, merged with follow-up equilibrium ones, to separately perform principal component analysis (PCA) on protein-protein hydrogen bonds. The resulting projection of the ligand entry trajectory on the two principal components, shown in Fig 5A, allows to distinguish between two kinds of correlated structural alterations: one associated with the transition from the apo to the bound form (PC1), and another connected to a transient and reversible reorganization of the binding channel (PC2). The residues contributing the most to PC1, shown in Fig 5B, cluster in several groups, mostly surrounding the binding site. In particular, the contacts of Y60^A with E14^A and W63^B with E14^B, present in the apo state and disrupted by the incoming ligand, are partially compensated by new H-bonds between E14^A and W63^A and S43^B with E14^B. Simultaneously, backbone H-bonding in a region of transmembrane helix 3 (TM3) stretching from G65^A to I71^A undergoes rearrangement due to formation of a kink in the exceptionally elastic 65GVG67 region. This kinking was indeed deemed important for ligand binding and the rate of conformational change, as concluded from past mutagenesis, NMR, EPR and cryo-EM studies [26, 29, 46, 47], and might provide an easy way of storing excess free energy for the subsequent ligand release. Finally, minor rearrangements occur at the entry site—S75^B switching from an intermolecular H-bond with S105^A to an intramolecular one with the backbone of a neighboring I75^B—reflecting a more open conformation of the binding channel in the bound than in the apo state.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Spontaneous ligand entry trajectory from weighted ensemble simulations, extended with a 100 ns equilibrium trajectory.

(A) Projection of the trajectory on the two main principal components describing correlated changes in H-bonding patterns. (B) Residues involved in the definition of PC1. Thick, H-bonds involve the backbone; thin, H-bonds engage the side chain. Residues were color-coded to show if they gain (orange), change (gray) or lose (cyan) H-bonds upon ligand binding. (C) Schematic description of the “tumbling” mechanism of ligand binding and unbinding, with individual colors corresponding to distinct phenyl groups of .

https://doi.org/10.1371/journal.pcbi.1009454.g005

Due to a much slower timescale of the bound→apo transition, the structural changes in the exit trajectories are much less pronounced and mostly map to spatially distinct regions, much less informative than the entry pathway described above (see S9 Fig for the corresponding PCA projection plot). However, both entry and exit trajectories highlight an interesting feature of the ligand progression through the binding channel, during which individual phenyl groups of “tumble” in a discrete way (see Fig 5C and S2 Movie). This step-wise movement explains the presence of several free energy minima in Fig 2B, and is enabled by correlating these “tumbles” with rotations of aromatic side chains lining the binding channel (inset of Fig 1), in line with observations from a previous simulation study [48]. This suggests that apart from increasing the affinity for aromatic ligands and sealing the access of water in the apo form, the presence of aromatic residues in the binding site might also serve to accelerate the kinetics of occupancy changes.

4.1 System preparation

The EmrE protein structure refined by Ovchinnikov V. et al. [28, 37], based on the experimental structure (PDB ID: 3B5D) [50], was used as a initial structure. The construction of the full system has been carried out through CHARMM Membrane Builder [51]. Briefly, the protein with its ligand was embedded in a model lipid bacterial membrane composed of 70% (168) POPE and 30% (72)POPG molecules, solvated with 13167 water molecules and 103 Na⁺ and 34 Cl⁻ ions by using the Monte Carlo method. For the ligand, atomic charge calculations were performed using the Gaussian [52] program. The CHARMM36 forcefield [53] was used for bilayer and TIP3 model was used for water. Equilibration and later simulations were performed with the GROMACS 2018.4 and plumed 2.5.0 software [54, 55]. The system was equilibrated in seven steps with the procedure described in the CHARMM-GUI, where the energy of the system was first minimized, next equilibrated in 5 steps with constantly decreasing restraint forces applied to the system (see S1 Fig for a schematic description of simulation lengths and dependencies). The simulations were conducted with Nose-Hoover thermostat at 310 K and semi-isotropic Parrinello-Rahman barostat was used to kept pressure at 1 Bar. The velocity Verlet algorithm equations of motion were integrated with a time step of 2 fs. For a long-range electrostatic interactions Particle Mesh Ewald (PME) algorithm was used with a real-space cutoff of (10 Å) and Van der Waals and short-range Coulomb interactions were represented by using a smooth cutoff of (12 Å) with a switching radius of 10 Å. The system prepared in this way was then subjected to a 1 μs equilibrium simulation. Then, using the correlation-based cryo-EM fitting module recently implemented in Gromacs package, the protein was optimized for its density map [26] for extra 130 ns. The resulting refitted structure, as well as the parameters for , are available at gitlab.com/KomBioMol/emre-setup.

4.2 Altered protonation states

To simulate specific protonation states of the two glutamic acids (main amino acids E14^A and E14^B responsible for the transport activity of the protein) in free energy simulations, additional hydrogens were added to each seeding frame obtained for the deprotonated variant, so that all four runs (umbrella sampling and metadynamics) were initialized from identical coordinates. Our in-house Gromologist library (gitlab.com/KomBioMol/gromologist) was used to automatize this process.

4.3 Simulations of spontaneous entry

The apo state of EmrE was prepared by removing the ligand from the system and then performing an equilibrium simulation of 1.75 μs. The obtained structure was closed to water on both sides of the protein, as shown in Fig 3A, 3C and 3D. Then, the Z-size of the box was increased by 1Å to randomly place 10 of each ligands in the resulting space. During the simulation, the applied pressure closed this gap, allowing ligands to diffuse freely. This process was repeated 3 times and for each ligand, and 1 μs long trajectories were obtained each time. Then, 3 structures closest to full entry into the binding site were extracted from the and systems and their parameters were swapped to produce another set of 1 μs long trajectories.

4.4 1D/2D metadynamics

To obtain the free energy profiles, we used two variants of well-tempered multiple walker metadynamics [56, 57] as implemented in the PLUMED package [55]. For the 1D variant, only the z-coordinate was used, defined as the z-component of the distance vector from the protein center of mass to the phosphorus atom. The ligand was restrained to a cylinder with a diameter of 1.5 nm and height of 6 nm from the center of the protein. In the 2D variant, the difference between the RMSD with respect to two symmetrically swapped domains (corresponding to the major conformational change) was used as an additional reaction coordinate. Overall, 20 walkers with 500 ns per walker (1D) or 26 walkers with 400 ns per walker (2D) were used for each of the four protonation states, yielding a total of 10 μs per protonation state per one free energy calculation. Initial structures were equally spaced along the Z-coordinate, and were taken either from a steered MD run (up to z = 1.1 nm) or from spontaneous entry simulations (z>1.1 nm). In the 2D variant, 13 walkers were initialized with the open side pointing upward and the other 13 downward; the latter were obtained by simple geometric manipulations, i.e. adequate rotation and relabeling of chains.

4.5 Umbrella sampling

25 frames with different ligand positions along the Z axis were created for umbrella sampling. They were obtained by using steered MD to pull the ligand through the membrane using the Gromacs package patched with Plumed. The ligand was restrained along the z-axis using a harmonic potential with a force constant of 500 kJ/nm², and simulated for 500 ns. After that, we unbiased all histograms separately for each ligand, using the weighted histogram analysis (WHAM) method [58] implemented in an in-house script (gitlab.com/KomBioMol/wham).

For ligand transition, the pathway was divided into 36 windows with a spacing of 0.05 nm (from -0.3 to 1.45 nm) and a force constant of 5000 kJ/nm², and 8 windows with a spacing of 0.2 nm (from 1.6 to 3.0 nm) and a force constant of 200 kJ/nm². For each window, 200 to 250 ns of trajectory was generated, totaling ca. 10 μs per protonation state per one free energy calculation, as in the case of metadynamics. The profiles were analyzed in the same manner as above.

Individual profiles from the above free energy methods were averaged to obtain the results shown in Fig 2, following an approach suggested by Laio and Gervasio [59]. These profiles, as well as their convergence, are shown in S4 Fig.

4.6 Grotthuss wire analysis

To identify continuous wires of h-bond-connected water molecules, we used an in-house script (gitlab.com/KomBioMol/proton_wire) that applies a graph search approach to find connections between “source” residues (here: the central glutamates) and a “target” (here: waters above or below the bilayer phosphate level for wires oriented up- and downward, respectively). Standard distance- and angle-based thresholds (≤0.35 nm for acceptor-donor distance and ≥120deg for acceptor-proton-donor angle) were used for h-bond identification.

4.7 Alchemical calculation of pK_a

Calculations of pK_a of glutamic acids were carried out using the Gromacs package patched with Plumed, with an in-house script for topology preparation (gitlab.com/KomBioMol/proton_alchemist). As a reference, a terminally capped glutamic acid was placed and restrained in the water phase with its alchemical states swapped, so that the system was neutral at all times, and the resulting ΔG was identical to ΔΔG corresponding to the pK_a shift. The ΔpK_as were calculated from the Henderson-Hasselbach equation using ΔG values obtained with Bennett Acceptance Ratio (BAR), and added to the reference value for free glutamic acid in solution (4.25).

4.8 Weighted ensemble simulations

WESTPA (The Weighted Ensemble Simulation Toolkit with Parallelization and Analysis) is a high-performance Python framework for simulating rare events on an extended timescale with rigorous kinetics using the Huber-Kim weighted ensemble algorithm [60]. In order to sample the rare event of ligand entry into or escape from the protein, the reaction interval is divided into small bins that can be reached from neighboring ones almost immediately, during a single 10-ps iteration. Within each bin, a fixed number of simulations is carried out in each iteration, and individual runs are then merged or duplicated in a way that preserves the total probability. Initially, 33 such bins were defined, with higher numbers of simulations per bin in areas corresponding to predicted free energy barriers. These values were then adjusted depending on how far the ligand progressed and where it was trapped kinetically. Movies illustrating the respective pathways were generated with Molywood [61].