Extracellular loops of BtuB facilitate transport of vitamin B12 through the outer membrane of E. coli

Vitamin B12 (or cobalamin) is an enzymatic cofactor essential both for mammals and bacteria. However, cobalamin can be synthesized only by few microorganisms so most bacteria need to take it up from the environment through the TonB-dependent transport system. The first stage of cobalamin import to E. coli cells occurs through the outer-membrane receptor called BtuB. Vitamin B12 binds with high affinity to the extracellular side of the BtuB protein. BtuB forms a β-barrel with inner luminal domain and extracellular loops. To mechanically allow for cobalamin passage, the luminal domain needs to partially unfold with the help of the inner-membrane TonB protein. However, the mechanism of cobalamin permeation is unknown. Using all-atom molecular dynamics, we simulated the transport of cobalamin through the BtuB receptor embedded in an asymmetric and heterogeneous E. coli outer-membrane. To enhance conformational sampling of the BtuB loops, we developed the Gaussian force-simulated annealing method (GF-SA) and coupled it with umbrella sampling. We found that cobalamin needs to rotate in order to permeate through BtuB. We showed that the mobility of BtuB extracellular loops is crucial for cobalamin binding and transport and resembles an induced-fit mechanism. Loop mobility depends not only on the position of cobalamin but also on the extension of luminal domain. We provided atomistic details of cobalamin transport through the BtuB receptor showing the essential role of the mobility of BtuB extracellular loops. A similar TonB-dependent transport system is used also by many other compounds, such as haem and siderophores, and importantly, can be hijacked by natural antibiotics. Our work could have implications for future delivery of antibiotics to bacteria using this transport system.


Introduction
The outer membrane of Gram-negative bacteria permits downhill diffusion of small hydrophilic molecules into the periplasmic space. However, bulky molecules (over *600 Da) or compounds present in the extracellular medium in scarce amounts, such as organometallic species, cannot permeate through porins, and their transport requires different routes. Since these compounds are fundamental for bacterial growth and viability, and largely determine the virulence of pathogens, bacteria developed specific and selective receptor-dependent active transport systems for these molecules. So far, this kind of system has been best characterized for vitamin B 12 , also called cobalamin (Cbl). Vitamin B 12 is an enzymatic cofactor essential both for mammals and bacteria, but only few microorganisms can synthesize Cbl. Most bacteria including Enterobacteriaceae, Bacillus subtilis and group A Streptococci need to take up vitamin B 12 from the environment.
The BtuB architecture is characteristic for TBDTs and organized in two domains-a 22-stranded β-barrel domain and N-terminal globular-like luminal domain (also known as a hatch or plug domain) that occludes the lumen of the BtuB barrel (Fig 1). The three-dimensional structure of BtuB has been resolved in the apo state (using detergent-based [12] and membrane-based crystallization [13]), in the apo state with bound calcium [12], in the holo state [12], and in the holo state in the complex with the C-terminal part of TonB [14].
The Cbl binding site is located on the two apices of the plug domain surrounded by the extracellular loops that partially sequester the ligand. The crystal structures of BtuB suggest that loops L2 through L5 play key roles in high-affinity binding of vitamin B 12 [12,16] (Fig 1). Also, partial deletions in loops L7, L8, L9, and L11 strongly decreased (but did not abrogate) Cbl binding and transport [17].
In proximity to vitamin B 12 binding pocket, there are two calcium binding sites, contained in a cage formed by five aspartates [12,16]. Intriguingly, for saturating calcium levels, the binding affinity of BtuB for Cbl is 5 nM, while at sub-saturating calcium concentrations, the Cbl affinity decreases 1000-fold [18][19][20]. The electron paramagnetic resonance (EPR) spectroscopy [21] and molecular dynamics (MD) simulations [22][23][24] suggested that calcium ions limit the conformational space sampled by loop L2 preventing its occlusion of the Cbl binding site.
The N-terminal periplasm-exposed region of the luminal domain conserved among the TBDTs [25], called Ton box, was shown crucial for Cbl transport [26][27][28] (Fig 1). In the apo state of BtuB, Ton box exists in equilibrium between the folded and unfolded states. Binding of vitamin B 12 to BtuB shifts this equilibrium towards the disordered state [29][30][31][32]. As a result, Ton box extends approximately 20-30 Å into the periplasm, which may initiate its recognition by the TonB protein [33].
The TonB protein is anchored in the inner membrane by the N-terminal single transmembrane helix, which is connected to the extended proline-rich periplasmic domain. The C-terminal part of this TonB domain reaches out to the outer membrane [34], which enables it to form a non-covalent complex with Ton box of BtuB [14] (Fig 1). If not bound to the BtuB receptor, TonB exists mainly as a dimer, but during transport TonB is converted to a monomeric form to yield the high-affinity BtuB-TonB complex [35,36]. TonB can allosterically affect the binding affinity of vitamin B 12 for BtuB, however the directionality of this activity is still a matter of debate [21,36].
Since the luminal domain of BtuB physically precludes free diffusion and passage of Cbl through BtuB, some conformational rearrangement of this plug, exposing a permeation pathway, has to be induced before TonB disengages from BtuB. The most probable mechanism of luminal domain rearrangement, upon its interaction with TonB, is the mechanical pulling hypothesis [25]. According to this model, TonB mediates a force that changes the plug conformation in a way that Cbl can permeate through BtuB and reach the periplasm. This force is believed to be generated by an electrochemical gradient of protons across the inner membrane utilized by ExbB and ExbD-TonB accessory inner-membrane proteins [37,38]. However, how this proton-motive force is coupled to the outer-membrane BtuB transporter is still  [14]. BtuB extracellular loops are highlighted in green but in terms of topology they belong to the BtuB barrel domain. Right: The Cbl structure (only heavy atoms are shown) and numbering of BtuB extracellular loops in accord with [15]. unknown. Single-molecule force spectroscopy studies have shown that the non-covalent interaction between the Ton box and C-terminal domain of TonB is durable under tension up to about 200 Å extension of the luminal domain [39]. According to steered molecular dynamics (SMD) simulations [39,40], such extension, which corresponds to the unfolding of about 50 amino acids within the luminal domain, is sufficient to create a channel tailored for Cbl passage through the BtuB barrel. Bacterial growth assays suggested that such partial unfolding of the luminal domain precludes free diffusion of undesirable bulky molecules, such as antibiotics, into the cell thus narrowing the set of compounds allowed for passage and transport selectivity to Cbl [39].
One of the crucial issues that remains unanswered is the mechanism of passage of vitamin B 12 through the BtuB itself. Since Cbl binds to the part of the luminal domain that unfolds upon TonB activity, it was speculated that Cbl could be pulled through BtuB. However, SMD simulations of unfolding of the luminal domain have shown that concurrent displacement of vitamin B 12 towards the periplasm was negligible [39]. Thus, in this work, we investigate and describe complete permeation of Cbl through BtuB at atomistic level of detail using enhanced sampling molecular dynamics (MD) simulations, including umbrella sampling (US) simulations. We also identify the roles of BtuB extracellular loops in the transport of vitamin B 12 . To sample the conformational space of the loops in the outer-membrane environment, we developed the Gaussian force-simulated annealing method (GF-SA) and applied it to BtuB in different states of the transport cycle. Overall, we characterize the relation between the vitamin B 12 passage, luminal domain unfolding and geometry of the extracellular loops, providing a whole atomistically-detailed picture of Cbl transport through BtuB.

System building
To build the BtuB model, the crystal structure of BtuB (PDB ID: 2GSK) [14] in the complex with cyanocobalamin from the OPM (Orientations of Proteins in Membranes) database [41] was used. Thus, in the text, for simplicity, we use the Cbl abbreviation for cyanocobalamin. The first four N-terminal amino acids missing in the BtuB structure were added, and the Cterminal fragment of TonB was removed. Protonation states of BtuB titratable groups were established and hydrogen atoms added with Propka 3.0 using the PDB2PQR server [42]. Cbl, water molecules and calcium ions were preserved as in the crystal structure. The FB15 AMBER force field [43] was used for BtuB. For vitamin B 12 , the AMBER bonded parameters developed by Marques et al [44] were applied. Non-bonded parameters for Cbl were set using General Amber Force Field 2 (GAFF2) and the partial charges were calculated according to the RESP procedure [45] using Gaussian 09 [46] and antechamber (AmberTools17).
The outer membrane of E. coli was built as an asymmetric heterogenous lipid bilayer. Lipids were assembled around BtuB using CellMicrocosmos Membrane Editor 2.2 [47]. The external leaflet of the membrane was made of 29 lipopolysaccharide (LPS) molecules to mimic the K-12 strain of E. coli. This LPS possesses the membrane-forming lipid A, as well as the inner and outer-core oligosaccharides but does not have the O antigen. The interior leaflet of the membrane was composed of six kinds of lipids: POPE, PMPE, PMPG, PSPG, QMPE, and OSPE with the ratio of 8:31:8:8:8:6 [48]. Bonded and non-bonded parameters for lipids were determined using GAFF2 with the partial atomic charges calculated according to the RESP procedure [45] using Gaussian 09 [46] and antechamber (AmberTools17).
Using xleap (AmberTools17) [49] and Visual Molecular Dynamics (VMD) [50], 20 Å and 280 Å layers of water molecules (in the z-coordinate) were added over and below the membrane, respectively. The initial dimensions of system were 95 Å x 95 Å x 380 Å. The negative potential of LPS was neutralized with calcium ions, while the negative charge of the lower-leaflet lipids was balanced by adding sodium ions. In addition, ionic strength of 100 mM of NaCl was applied. The TIP3P-FB water model [51] was used as the one corresponding with the FB15 AMBER force field parameters for the protein. For the sodium and chloride ions, parameters developed by Joung and Cheatham [52], and for the calcium ions, by Li and Merz [53], were applied. The entire membrane-BtuB receptor system with Cbl consisted of about 296,000 atoms (Fig 2).
To confirm the correctness of the prepared asymmetric membrane model, we measured the area per lipid and membrane hydrophobic thickness. The area per lipid was computed using Voronoi tessellation Monte Carlo integration algorithm implemented in the VTMC program [54]. Membrane thickness was calculated as the difference between the average z-coordinate of C2 atoms of the acyl chains 4 and 6 and C4 atoms of the acyl chains 1, 2, 3, and 5 of LPS and average z-coordinate of C2 atoms of acyl tails of all of lower-leaflet lipids [55]. The average area per lipid of LPS is 151 Å 2 , which agrees with the experimental value of 156 Å 2 [56]. The average membrane thickness is 24.05±0.07 Å, which matches the average hydrophobic thickness of all known outer-membrane proteins (*24 Å) [57].

General molecular dynamics simulations protocol
Initially, the systems were energy minimized (under harmonic restraints of 10 kcal/mol/Å 2 set on non-hydrogen atoms) with the steepest descent method (5000 steps) followed by the conjugate gradient method (3000 steps) using the sander program [49]. Subsequent phases were carried out using NAMD 2.12 [58]. After minimization, the solute coordinates were fixed, and the system was gradually heated from 10 to 310 K, for 125 ps in 10 K increments, in the NpT ensemble, with a constant pressure of 1 atm controlled using the Langevin Piston method. The integration time step was set to 0.5 fs. Next, the water molecules and ions were equilibrated at 310 K for 2.5 ns. Further, the solute was restrained with a harmonic constant of 10 kcal/mol/ Å 2 imposed on heavy atoms, and the whole system was gradually thermalized from 10 to 310 K, in 10 K increments of 50 ps length. During equilibration, the restraints were removed over ten rounds of 50 ps each. Next, the integration time step was gradually increased from 0.5 fs to 2 fs, in 0.25 fs increments of 1 ns length, and the whole system was simulated for 400 ns. Periodic boundary conditions and Particle Mesh Ewald method with a grid spacing of 1.0 Å were used. The bonds involving hydrogens within non-water and water molecules were constrained using RATTLE [59] and SETTLE [60] algorithms, respectively. The cutoff for short-range non-bonded interactions was set to 12 Å, with a switching distance of 10 Å. Data were saved every 5 ps.

Steered molecular dynamics
To unfold the BtuB luminal domain in a controlled manner, we used constant-velocity steered molecular dynamics (SMD) [61,62]. The center of mass (COM) of the heavy atoms of the first N-terminal residue was pulled along the z-axis with a constant velocity of 1.0 Å/ns and a force constant of 1 kcal/mol/Å 2 . The force vector was [0, 0, -1]. The COM of the C α atoms of the lower part of the BtuB barrel domain was harmonically restrained, with a force constant of 100 kcal/mol/Å 2 , to prevent the protein-membrane system from drifting along the z-axis upon the BtuB luminal domain pulling.
Our preliminary simulations have shown that the release of vitamin B 12 from the BtuB protein requires additional force applied to vitamin B 12 (see Results and discussion). We observed that unfolding of the luminal domain during constant-velocity SMD does not result in vitamin B 12 transport even though the lumen becomes open and there is enough space for Cbl passage. This finding corroborates the studies of Hickman et al. [39]. However, we found that the release of vitamin B 12 is precluded because it is strongly held by the BtuB extracellular loops. Therefore, to allow for vitamin B 12 passage through BtuB, we applied external constant force on the COM of Cbl with a vector [0, 0, -x], where x was equal to 0, 5, 8, 10, 15, 20 kcal/mol/Å.

Umbrella sampling
To compute the potential of mean force (PMF) of the association and permeation of vitamin B 12 through BtuB, we applied the umbrella sampling (US) methodology [63]. The reaction coordinate was the distance projected onto a normal to the membrane (z-axis) between the COM of Cbl and dummy atom fixed in the center of the Cbl binding pocket. To make the position of the binding pocket comparable between different BtuB systems used in this study, we defined it as the average location of the COM of Cbl in the holo state of BtuB with luminal domain folded (the S2 state described in Results and discussion). The initial association and permeation pathway of vitamin B 12 through BtuB was obtained using SMD. Cbl was pulled from the binding pocket (z = 0) with a constant velocity of 1.0 Å/ns and a force constant of 1 kcal/mol/Å 2 along the z-axis. The force vector was [0, 0, 1] and [0, 0, -1] for association and permeation, respectively. Concurrently, to maintain the plug domain extension of 197 Å, the COM of the N-terminus of the luminal domain was harmonically restrained with a force constant of 100 kcal/mol/Å 2 . The US windows were spaced every 0.5 Å in the range of -30-70 Å of the reaction coordinate, and the force constant was set to 50 kcal/mol/Å 2 . Simulations in each US window lasted 37.5 ns, but first 2.5 ns were treated as an equilibration of the system properties, such as density, total energy and its contributions, as well as RMSD of BtuB-Cbl, thus were discarded from analysis. The PMF was computed using the Weighted Histogram Analysis Method (WHAM) [64,65] with the WHAM program [66]. The convergence of PMF was monitored every 0.2 ns by measuring the root-mean-square deviation (RMSD) of PMF calculated from cumulative data. The reference for PMF RMSD calculations was the last PMF (computed based on all of the collected data).

Gaussian force-simulated annealing method
To optimize the structure and enhance conformational sampling of the extracellular loops in different BtuB states (described in Results and discussion), we developed a Gaussian force-simulated annealing (GF-SA) method. During a typical simulated annealing protocol, a whole system is heated and, subsequently, slowly cooled down to find its global energy minimum. In our approach, the variable in the simulated annealing is not the temperature of the system but external Gaussian forces imposed locally on the BtuB loops.
The idea is to augment the rotation of chemical groups around the rotatable bonds in the protein backbone and side-chains. In the backbone, there are two freely rotatable bonds: C α -C and N-C α . The CO-NH moiety that connects these bonds is a planar rigid structure that prefers the trans position. Therefore, the conformational variety within the backbone arises from the relative positions of the CO-NH groups. In turn, the conformational variability of sidechains dominantly comes from the rotation of the residues around the C α -C β bond. Thus, to enhance protein dynamics within the backbone and side-chains, we applied external forces to the atoms neighboring with the atoms involved in rotatable bonds (mentioned above). In the protein backbone, this applies to the O atom, whereas in side-chains the atom types depend on the residue (e.g., C γ , O γ , etc.).
The forces imposed on the atoms were calculated as normalized vectors. The vectors were normal to the planes formed by the following angles: C α -C-O in the backbone and C α -C β -X in side-chains, where X is C γ , O γ , etc. The forces were multiplied by a pseudo-random number from the Gaussian distribution with the zero average and fixed standard deviation (see below). The C++ 11 Standard Library random class was used to generate the random numbers. In addition, to make random numbers uniformly distributed between the average and standard deviation, we applied the Knuth-b generator. Then, we used the Box-Muller method to transform these numbers to normally distributed. Each Gaussian number was produced using a new seed used by the generator. The seed was assigned nondeterministically by using the random device function. We found this procedure more effective in maintaining the normal distribution of pseudo-random numbers than using only one initial seed followed by the update of the generator state.
The external forces were added every second step of MD simulation and were imposed on the following residues In the GF-SA procedure, the standard deviation of the normally distributed pseudo-numbers, which determines how strongly the dynamics of the loops is enhanced, is changed over the simulation time. The standard deviation s was calculated according to the function: where a is a maximal standard deviation, b is an integer determining the number of stages during heating or cooling, n is an even integer controlling the gradient of the standard deviation. We set a to 35 kcal/mol/Å, n to 2, and b to 100. Heating was conducted in 100 stages x 2 < −100, 0) for 2000 MD steps each. Then, the standard deviation was maintained in its maximum (x = 0) for 250000 MD steps. Subsequently, cooling was carried out in 100 stages x 2 (0, 100 > for 30000 MD steps each. In the systems simulated with vitamin B 12 , during heating and cooling, the ligand was harmonically restrained with a force constant of 10 kcal/mol/ Å 2 . The restraint on vitamin B 12 was released for subsequent 50000 MD steps followed by conventional MD simulation (with no Gaussian forces). Then, the whole procedure, GF-SA-MD, was repeated for a number of cycles. To prevent the protein-membrane system from artificial drifting upon the external Gaussian forces imposed on the BtuB loops, the COM of the C α atoms of the lower part of the BtuB barrel domain was restrained harmonically with a force constant of 100 kcal/mol/Å 2 . The GF-SA method was implemented in NAMD using the TclForces module. The general code of the method was written in Tcl, however all mathematically robust operations were written in C++, and as an extension wrapped into the Tcl script. The wrapper code was generated using the Swig 3.0 program [67].

Data analysis
Trajectories were analyzed with the cpptraj program (AmberTools17) and VMD 1.9.3 using in-house Tcl and C++ scripts. VMD was also used for visualization. Plots were generated with Gnuplot (v5.2). The criteria for hydrogen bonds were set to 3.5 Å between the donor and acceptor, and no less than 135 degrees angle between donor-hydrogen-acceptor. To detect the contacts, the distance of 3.5 Å was used as the sole criterion. Descriptive statistics calculations were performed using in-house C++ scripts. Data sets underwent the Kolmogorov-Smirnov test to check the normality of the distribution. The data series with a not normal distribution were characterized by median, interquartile range (Q1-Q3), and lower-upper extremes range. In GF-SA-MD simulations, only data collected from the conventional MD were analyzed.

Conventional MD does not properly describe the translocation of vitamin B 12 through BtuB
So far, only Hickman et al. [39] performed SMD simulations of unfolding of the BtuB luminal domain in the presence of Cbl, to find a link with their AFM data. In the beginning, we performed analogous SMD simulations, and, similarly, we did not observe any substantial motion of vitamin B 12 towards the periplasmic space (S1 Fig). Cbl was held by the extracellular loops and did not follow the luminal domain upon its unfolding. Thus, to enhance Cbl movement, simultaneously with unfolding of the luminal domain, we imposed a constant force on the Cbl center of mass, directed toward the periplasm along the normal to the membrane (see S1 Video). Depending on the magnitude of this force (Fig 3B), we observed the release of vitamin B 12 at different extensions of the luminal domain. We noted that the earliest full permeation of Cbl was mechanically feasible once the luminal domain was extended to at least 175 Å and upon the external force of 20 kcal/mol/Å. We also found that in the first stage of passage (z = 0-24 Å), vitamin B 12 rotates along the main axis of BtuB in such a way that the Cbl corrin ring becomes oriented parallel to this axis ( Fig 3C). This indicates that plain translation of Cbl through the BtuB barrel would be mechanically ineffective.
However, SMD provides only a qualitative picture of vitamin B 12 permeation. To obtain quantitative data, namely the potential of mean force (PMF) of Cbl passage through BtuB, we applied US. Since Cbl transport is concomitant with unfolding of the BtuB luminal domain, the complete PMF depends on both the position of vitamin B 12 relative to the BtuB barrel and extension of the luminal domain. Thus, the PMF is at least two-dimensional. However, due to huge computational cost of 2D US simulations, we simplified this problem to 1D US calculations.
Hickman et al. [39] have shown by AFM that the interaction between Ton box and TonB is stable up to about 200 Å extension of the luminal domain. Therefore, the window for releasing vitamin B 12 should span the 175-200 Å extension of this domain. During the conventional SMD simulation of luminal domain unfolding, we observed that the highest steering force peaks occur at 160 Å and 215 Å of the luminal domain extension (S1 Fig). Thus, the partially unfolded BtuB becomes most stable while releasing vitamin B 12 . We also identified three local minima of the SMD force within the permeation window, at the 175 Å, 185 Å, and 197 Å extension. The latter, most stable extension, is the maximal up to which TonB can "pull" the luminal domain without breaking of their interactions [39], and is also sufficient to make space for Cbl passage within BtuB. Therefore, we assumed that Cbl passes through BtuB when the luminal domain fluctuates around its 197 Å extension point.

PLOS COMPUTATIONAL BIOLOGY
Subsequently, we performed 1D US calculations of vitamin B 12 permeation under the restrained position of the luminal domain N-terminus at 197 Å extension, which allows for Cbl passage from the mechanical point of view. However, after 35 ns of US simulations, we obtained a highly unstable PMF (S2 Fig), especially in the first stage of Cbl transport, at 0 < z < 30 Å. Notably, during this phase, vitamin B 12 interacts strongly with the BtuB loops requiring steering force for its release (Fig 3C). The loops tightly sequester Cbl at the binding site, and remain stable throughout the simulations. In addition, such rigid and compact configuration of the loops hampers the rotation of vitamin B 12 that is necessary for its permeation ( Fig  3C). Therefore, we concluded that the time scale of conventional MD is insufficient to extensively sample the conformational space of the BtuB extracellular loops.
Furthermore, due to crystallization conditions, the conformations of the loops resolved based on the BtuB crystal may be different than in physiologic state. Numerous studies on BtuB, such as EPR spectroscopy and site-directed spin labeling, demonstrated that either the solutes or osmolytes used in the crystallization buffer alter the conformational sampling of the solvent-exposed protein segments, such as Ton box and extracellular loops, favoring more compact, ordered, and less hydrated structures [32,[68][69][70][71][72]. A similar energy contribution could arise from the crystal lattice [32]. It was also proven that the detergent-based or membrane-based crystallization method affected the ordering of the BtuB extracellular loops and Ton box conformation [12,13]. Additionally, the most recent EPR spectroscopy measurements have shown that the BtuB loops sample a much broader conformational space than revealed by the crystal structures [21]. According to Sikora et al. [21], in the crystals, BtuB extracellular loops are conformationally trapped, representing one of the more compact structures among those sampled in the outer membrane. Thus, efficient conformational sampling of the loops is a must.

Conformational changes of the BtuB loops during the vitamin B 12 transport
To investigate the mobility of the BtuB extracellular loops and their role in Cbl transport, we developed an enhanced sampling Gaussian force-simulated annealing method, GF-SA, and coupled it with conventional MD (for details see Methods and S2 Video). In the GF-SA-MD simulations, we considered different states of the BtuB protein corresponding to different stages of the transport cycle ( Fig 4A) (see S3 Video). Their initial structures were obtained using SMD simulations described above. These states were: the apo state (S1), the holo state with the luminal domain folded (S2), the holo state with the luminal domain partially unfolded (S3), and the holo states with the luminal domain partially unfolded and Cbl shifted towards the periplasm (S4 and S5). Partial unfolding corresponds to the extension of the luminal domain at 197 Å in the z-axis with respect to the initial position. This extension was maintained by applying a harmonic restraint to the N-terminus of the luminal domain with the force constant of 50 kcal/mol/Å 2 . In S4 and S5, vitamin B 12 was positioned, respectively, 25 Å and 50 Å in the z-axis towards the periplasm (for the axis describing the positioning of vitamin B 12 and luminal domain see Fig 3A). For each system, we performed 10 independent simulations, each consisting of 15 consecutive cycles (see S3 Video). Each cycle comprised 7 ns of GF-SA followed by 7 ns of conventional MD.
To assess the conformational sampling of the loops in the GF-SA-MD simulations compared to conventional MD on the microsecond time scale, we measured the distances between several loop residues in state S2, analogously to the EPR spectroscopy measurements of Sikora et al. [21]. Direct comparisons between the distance distributions obtained by EPR and GF-SA-MD are error-prone because the EPR were measured for spin-labeled residues, while the GF-SA-MD for natural, non-labeled ones. However, in GF-SA-MD for all distances, the distance distributions were substantially more spread, suggesting enhanced sampling of loop conformations as compared to conventional MD (S3 Fig). Therefore, we argue that GF-SA-MD may provide a more realistic picture of mobility of the BtuB loops than conventional MD.
To detect the differences in the conformations of the loops in different BtuB states (Fig 4A), on combined MD trajectories, we performed Principal Component Analysis (PCA) of the BtuB barrel domain C α atoms. The first 50 eigenvectors contributed to over 80% of the total mean square fluctuations. Respective trajectories were projected onto the principal components (PCs) and showed largest variations in the two lowest frequency modes PC 1 and PC 2 ( Fig 4B). We found that these PC modes correspond to the opening (PC 1 < 0 or PC 2 > 0) and closing (PC 1 > 0 or PC 2 < 0) motions of BtuB (see S4 and S5 Videos), which mainly involve loops L2, L3, L7, L8, L10, and L11 (S4 Fig). To estimate the probability of the closed and open loop geometries in different BtuB states, we summed the distances between the center of masses of individual loops and the center of mass of Cbl binding site (Fig 5A). The boundary between the closed and open BtuB geometries was defined as the conformation of the BtuB barrel domain (Fig 1) averaged over all trajectories (for which both PC 1 and PC 2 equal 0 Å). We also analyzed the distance distributions of particular loops to Cbl binding pocket to quantify their individual contributions to the closing-opening motions in the different BtuB states (Fig 5B, S1 Table, S5 Fig).
In the apo state of BtuB (S1), as shown in Fig 4B, both PC 1 and PC 2 are widely and almost symmetrically distributed. This agrees with the ratio of the closed to open loop geometries of about 0.5 (Fig 5A). This also suggests an equilibrium between the closing and opening movements of the loops, with BtuB preference to adopt a neutral conformation (in which both PC 1 and PC 2 are positioned around 0 Å). The largest share in these motions comes from the loops L8, L3, and L7 (Fig 5B).

PLOS COMPUTATIONAL BIOLOGY
Binding of Cbl (S2) restricts the distribution of PC 1 and shifts it towards negative values ( Fig 4B). Thus, in S2, the open conformation of BtuB loops is more occupied (0.86) than the closed one (0.14) (Fig 5A). The reason lies in steric hindrance caused by bound Cbl, which makes the loops, especially L3 and L8, protrude outside the Cbl binding site and lowers their mobility. In contrast, loop L7 prefers to move closer to the binding pocket upon Cbl binding (Fig 5B).
In state S3, when the luminal domain gets partially unfolded, PC 1 and PC 2 are again nearly symmetrically distributed (as in S1). However, the range of PC 1 is narrower than is S1 but similar as in S2, suggesting lower amplitude of loop motions in the direction of PC 1 ( Fig  4B). In state S3, the probability of loop closing motions increases to 0.55 (Fig 5A) because loops L3, L8, L10, and L11 approach the binding site more often (Fig 5B). This implies that

PLOS COMPUTATIONAL BIOLOGY
partial unfolding of the luminal domain may serve as a trigger to initiate the closure of the loops during Cbl transport.
In S4, after partial permeation of Cbl, the range of PC 1 shifts towards positive and of PC 2 towards negative values (Fig 4B), and the closed-open loop occurrence becomes opposite to the one of S2 (Fig 5A). This signifies that in S4 the loops prefer an even more closed geometry (0.86) than in S3. The closure of the BtuB barrel is mainly contributed by loops L3 and L8. Besides, in S4, the PC 1 range widens as compared to S3, indicating increased loop mobility. This is mainly observed for loops L2 and L10 (Fig 5B).
When vitamin B 12 is fully permeated (state S5), the distribution of PC 1 and PC 2 resembles a mix of that in S1 and S4 ( Fig 4B). As in S1, both PCs span a wide range of values suggesting large loop mobility. Indeed, most of the loops in S5, including L3, L7, L8, and L10, sample wide conformational space (Fig 5B). However, while PC 2 is symmetrically distributed (like in S1), PC 1 samples more positive values (like in S4). This is because, in S5, loops L2, L3, L10, and L11 approach the Cbl binding site more frequently than in S1 (Fig 5B). Nevertheless, in S5, the closing motions of the loops are less frequent than in S4, mainly due to moving of loop L8 outside the binding pocket. Overall, in states S3-S5, during Cbl permeation with luminal domain partially unfolded, the loops prefer the closed conformations.

Non-specific interactions govern the mobility of BtuB loops
To explain the conformational changes of the extracellular loops in various states of BtuB, we investigated their contacts and energetic properties. The number of contacts within the loops varies in the BtuB states and is maximal in state S2 upon Cbl binding (Fig 6). Then, while moving from state S2 to S5, the intra-loop contacts gradually decrease, with a change of over 100 in the median (Fig 6 and S2 Table). Since the number of hydrogen bonds formed within the loops is similar in these states (S2 Table, S6 Fig), the changes in the intra-loop contacts can be associated with non-specific interactions.
The number of contacts between the loops also depends on the BtuB state and increases during Cbl permeation (states S3-S5) (Fig 6, S7 Fig). Again, the hydrogen bonds between the loops practically do not change between the BtuB states (S6 Fig). This suggests that most interloop contacts during Cbl transport are non-specific, which is also reflected in the energetic properties of these loops.

PLOS COMPUTATIONAL BIOLOGY
Naturally, with more contacts between the loops, the van der Waals energy of their interactions becomes more negative (S7 Fig). However, despite the fact that for most loop pairs (except pairs including L9 and L10), the electrostatic interaction energy becomes more repulsive from state S1 to S5, it is balanced by the stabilizing van der Waals forces.
Not surprisingly, the number of contacts between the loops is the lowest in state S2 because the loops interact also with Cbl and folded luminal domain. Furthermore, in S2 most loops maximize their intra-loop connections (Fig 6). However, one could ask why the number of contacts between the loops is on average lower in S1 than in S5 (Fig 6) since in both states, the loops lack interactions with vitamin B 12 . This occurs because in S1 the number of contacts between the loops and luminal domain is the highest (Fig 6, S3 Table). Thus, it seems that in the apo state of BtuB (S1), the loops, mainly L2, L3, L7, L8, L9, L10 and L11, prefer to interact with the luminal domain rather than with themselves. On the contrary, in S5, after Cbl translocation, many of the interactions of the loops with luminal domain disappear due to its partial unfolding (S8 Fig), and this is compensated by increasing the inter-loop connections.
Analogous phenomenon occurs while moving from state S2 to S3. Interestingly, the intensified contacts between the loops in S3 are accompanied by the increased interactions with Cbl ( Fig 6, S9 Fig).
The contacts between the loops and luminal domain gradually increase during permeation of vitamin B 12 , i.e., between state S3 and S5. When the number of contacts between the loops and Cbl decreases, the loops compensate it by interacting with themselves and luminal domain. Most of these contacts are non-specific, because the number of direct and watermediated hydrogen bonds of the loops with luminal domain remains almost unaffected (S3 Table). In addition, for most loops their electrostatic interaction energy with the luminal domain becomes unfavorable (S10 Fig, S4 Table). However, the overall van der Waals forces become more stabilizing.
We also examined the interactions of the loops with solvent by analyzing the number of contacts and hydrogen bonds formed between all loops and water molecules (Fig 6, S6 Fig, S5  Table). We found that both numbers slightly increase while moving from S2 to S5. As anticipated, the loops were the least hydrated in states S2 and S3 (with bound Cbl) and the most hydrated in state S5 (with Cbl fully permeated). However, note that the hydration of the loops is comparable in states S1, S3, and S4. This implies that during Cbl transport the loops reduce their contacts with water molecules by interacting with other loops (as in S3-S5), luminal domain (S1, S3-S5) or Cbl (S2-S3), which makes sense because the loops contain many aromatic amino acids, such as Tyr, and overall are hydrophobic.

Enhanced sampling of the BtuB loops affects the energetics of vitamin B 12 permeation
The GF-SA simulations suggested that inadequate probing of the conformational space of BtuB extracellular loops in US could be the reason for the highly unstable PMF of Cbl permeation (S2 Fig). Therefore, we coupled US with the GF-SA approach (see Methods). For each window, we carried out 10 consecutive cycles of GF-SA conjugated with US. Each cycle consisted of 7 ns of GF-SA followed by 3 ns of conventional US. So, in aggregate, 30 ns of US with enhanced sampling of the loops for each window were used to compute the final PMF. With this GF-SA-US approach, the PMF stabilized after 7 cycles (S11 Fig). The PMF of Cbl association with BtuB (z < 0) and its permeation through BtuB (z > 0), together with representative structures corresponding to local energy minima, are presented in Fig 7. On the Cbl association pathway, there are no local energy minima; for z < −20 Å, the PMF stabilizes at about 7.7 kcal/mol. Upon further Cbl movement, up to z = 4 Å, the PMF falls to ca. 1 kcal/mol. This almost linear drop in PMF correlates with the decrease in the interaction energy between Cbl and BtuB (Fig 8A).
In z between 4 and 24 Å, the PMF does not exceed 1.4 kcal/mol with the lowest-energy state at z equal 19 Å (Fig 7A). Other multiple energy minima, observed in this region, separated by low barriers might result from imperfect sampling of the degrees of freedom orthogonal to the reaction coordinate, such as the Cbl orientation or loop conformations. Similarly as in SMD simulations, we observed the rotation of vitamin B 12 relative to the main axis of BtuB by about  PLOS COMPUTATIONAL BIOLOGY 120 degrees (Fig 7B). Overall, in this z-region, vitamin B 12 forms many contacts with the BtuB protein ( Fig 9A). Specifically, Cbl interacts with the luminal domain residues 65-67, 72-74, 90-98, and 101, and Cbl interactions with I66, F67, N72, A73, D95, L96, S97 and Q98 occur in almost entire z-range of 4-24 Å (Fig 7).
The interactions between Cbl and BtuB barrel are detectable only when z > 15 Å. However, they are substantially less abundant and more transient than those with the luminal domain or loops. The most frequent contacts are observed with BtuB residues K504, Q521, L546, K564, E585, and T587. In the z-range of 24-31 Å, the PMF increases by about 4 kcal/mol (Fig 7A). This correlates well with the increase in the interaction energy between vitamin B 12 and BtuB (Fig 8A), which is mainly caused by the loss of contacts with all extracellular loops and weakening of interactions of Cbl with the luminal domain. Meanwhile, Cbl contacts with several BtuB barrel residues, including R562, K564, T587, and S589, increase. Furthermore, the tilt of Cbl with respect to the BtuB main axis lowers from 120 to about 80 degrees. After Cbl moving from the local energy maximum at z = 31 Å, the Cbl rotation angle again starts to rise (Fig 7B).
At z = 40 Å, the PMF reaches a local maximum of 5.7 kcal/mol, which most probably is associated with Cbl losing most of its hitherto existing interactions with the luminal domain. Only contacts with a small fraction of luminal domain (residues 133-137) are maintained. Then, vitamin B 12 reaches its maximal tilt relative to BtuB equal about 140 degrees. Next, until z = 47 Å, the angle of Cbl rotation decreases to ca. 80 degrees (Fig 7B).
When z > 55 Å, Cbl no longer contacts the BtuB barrel domain, but interacts with the extended fragment of the luminal domain (residues 50-60) exposed to the periplasmic space.

PLOS COMPUTATIONAL BIOLOGY
The PMF reaches a plateau of ca. 6.3 kcal/mol at z > 60 Å. Since Cbl after leaving the BtuB barrel interacts with the extended luminal domain, the PMF plateau after Cbl permeation is ca. 1.4 kcal/mol lower than the PMF plateau achieved at z about −30 Å (when vitamin B 12 is in the bulk).
On its transport pathway, vitamin B 12 interacts with BtuB via direct and water-mediated hydrogen bonds (Fig 9B) formed primarily through the phosphate, ribose 5'-hydroxyl, and amide groups (S12B and S12C Fig). The same groups preferentially form hydrogen bonds with bulk water (S12D Fig). The hydrogen bonds are augmented by non-specific interactions involving mainly aliphatic moieties of acetamide and propionamide chains, as well as cyanide group (S12A Fig). This seems reasonable because during vitamin B 12 transport Cbl interacts with many amino acids with a hydrophobic side-chain, e.g., Tyr, Ile, Leu, Val, Ala, and Phe.
We believe that the high number of transient and non-specific interactions and watermediated hydrogen bonds formed between Cbl and BtuB helps lower the friction associated with migration of vitamin B 12 within the BtuB protein and enables its rotation in the first phase of permeation (z = 0-24 Å) ( Fig 7B). Next, the average number of direct and watermediated hydrogen bonds formed by Cbl with BtuB decreases and the average number of Cbl hydrogen bonds with bulk water increases (Fig 9B), which helps in the movement of Cbl from BtuB to the periplasm.
To explore how enhancing the sampling of the BtuB loops influenced vitamin B 12 transport, we compared various system properties calculated from the GF-SA-US and conventional US simulations. We observed that in conventional US, in z between 0 and 3 Å, the Cbl-BtuB interaction energy is more favorable (by ca. 40 kcal/mol) as compared to the energy in GF-SA-US ( Fig 8A). This happens because in US, in this z-region, more contacts (by ca. 140) are formed between Cbl and BtuB loops, particularly L8, L11, and L4 ( Fig 8B). As a consequence, the conformations of extracellular loops obtained from conventional MD simulations, even microsecond long, cause tight binding of Cbl in its initial position (z = 0-3 Å), which hampers vitamin B 12 migration towards the periplasm. This corroborates with a high peak of the SMD force necessary to push Cbl through BtuB (Fig 3C).
Intriguingly, in the z-range of 3-24 Å, the interaction energy between Cbl and BtuB is more negative in the GF-SA-US simulations than in conventional US. Again, the more favorable Cbl-BtuB interaction energy results from the higher number of contacts between Cbl and BtuB loops, especially L3, L2, L9, L10, and L11, in GF-SA-US (Fig 8B). This further implies that loops optimize their conformations to favorably interact with vitamin B 12 at the early stage of its permeation through BtuB, which leads to greater stabilization of Cbl, reflected in its low atomic fluctuations (Fig 8C). This stabilization appears important because during this phase of transport, vitamin B 12 rotates, which is crucial for its permeation.

Conclusions
We determined how the conformational dynamics of the BtuB receptor mediates the transport of vitamin B 12 through the E. coli outer-membrane. We performed a series of conventional and enhanced sampling MD simulations of the full-atomistic model of BtuB embedded in an asymmetric and heterogeneous lipid bilayer mimicking the outer membrane of the E. coli K-12 strain.
Our SMD simulations have shown that the release of Cbl from BtuB is mechanically feasible after the luminal domain unfolds at least 175 Å in the direction perpendicular to the membrane. AFM experiments have shown that the non-covalent connection between TonB and Ton box (the N-terminal fragment of the luminal domain) is viable up to ca. 200 Å of the luminal domain extension [39]. Thus, we assumed that vitamin B 12 permeates through BtuB at the domain extension range of 175-200 Å. Therefore, to simulate Cbl permeation, we extended and constrained the luminal domain at 197 Å, to provide enough room inside the BtuB barrel to mechanically allow for Cbl transport. Further, in the simulations, we found that in the early stage of permeation, Cbl needs to rotate so its corrin ring becomes perpendicular to the surface of the membrane.
Conventional US simulations provided an unstable PMF of vitamin B 12 transport because of insufficient probing of the conformational space of BtuB extracellular loops in conventional MD. This suggested that detailed description of loop mobility is essential to properly characterize Cbl permeation. Therefore, to enhance loop sampling, we developed a method named Gaussian force-simulated annealing (GF-SA). Using this approach, we investigated loop mobility in several BtuB states corresponding to different stages in Cbl transport cycle. We found that several loops, especially L2, L3, L7, L8, L10, and L11, participate in low-frequency closing-opening motions of the BtuB barrel, which may be crucial for Cbl transport. In the apo state, the loops span their closed and open conformations equally well. After Cbl binding, the equilibrium shifts toward the open states, and during luminal domain unfolding and subsequent Cbl permeation, the loops (especially L3, L8 and L10) tend to close over the Cbl binding site.
The behavior of BtuB extracellular loops resembles the induced-fit mechanism occurring during ligand-protein binding, in which the protein adapts its conformation in response to ligand binding to yield a high-affinity complex. Such phenomenon has been already observed for other TBDTs such as FecA [73], FhuA [74], FepA [75], and ShuA [76]. In these proteins, extracellular loops perform large conformational motions to fold over the ligand and sequester it at the binding site. For BtuB, we found that the fitting of loops is triggered not only by Cbl binding but also by unfolding of the luminal domain. The former affects mainly the mobility of loops L3 and L8, while the latter largely determines the behavior of L2, L10, and L11.
It appears that the conformational fitting arises primarily from the propensity of the loops, rich in hydrophobic amino acids, to maximize the internal non-specific contacts within the Cbl-BtuB complex. Such loop action stabilizes the Cbl molecule during the initial phase of permeation when it undergoes rotation. This facilitates vitamin B 12 transport through the outermembrane BtuB protein and may help exclude Cbl dissociation to the external environment. This could also prevent some unwanted solutes present in the environment, such as antibiotics, from entering BtuB, thus ensuring cell protection.
Coupling of conventional US with GF-SA, which enhanced loop sampling, considerably stabilized the PMF of vitamin B 12 permeation. We found that Cbl passage is energetically favorable up to ca. 20 Å towards the periplasm. The energy barrier on the Cbl permeation route is only about 6 kcal/mol, which implies that transport is rapid. Besides, many transient interactions and water-mediated hydrogen bonds between Cbl and BtuB assure for smooth motion of vitamin B 12 within and through the BtuB barrel. Note that our model of Cbl transport is limited by the constrained position of the luminal domain. Additionally, despite the dramatic effect of enhanced sampling of the loops on the PMF of Cbl transport, we expect that some relevant degrees of freedom, such as the orientation of vitamin B 12 , may require even longer sampling to achieve convergence. However, it would be computationally extremely demanding to perform US simulations varying two to three collective variables (Cbl position, Cbl orientation, and extension of the luminal domain) combined with prolonged enhanced sampling of the loops. Nevertheless, our simulations clearly suggest the involvement of the loops "pushing" Cbl in an induced-fit mechanism. Yet, what we could have not captured is that the extending luminal domain could "pull" Cbl towards the periplasm, which would lower the permeation barrier even more.