proteins STRUCTURE O FUNCTION O BIOINFORMATICS

Multiple conformational states and gate opening of outer membrane protein TolC revealed by molecular dynamics simulations Beibei Wang,1 Jingwei Weng,1* and Wenning Wang1,2* 1 Department of Chemistry, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, Fudan University, Shanghai, People’s Republic of China 2 Institutes of Biomedical Sciences, Fudan University, Shanghai, People’s Republic of China

ABSTRACT Outer membrane protein TolC serves as an exit duct for exporting substances out of cell. The occluded periplasmic entrance of TolC is required to open for substrate transport, although the opening mechanism remains elusive. In this study, systematic molecular dynamics (MD) simulations for wild type TolC and six mutants were performed to explore the conformational dynamics of TolC. The periplasmic gate was shown to sample multiple conformational states with various degrees of gating opening. The gate opening was facilitated by all mutations except Y362F, which adopts an even more closed state than wild type TolC. The interprotomer salt-bridge R367–D153 is turned out to be crucial for periplasmic gate opening. The mutations that disrupt the interactions at the periplasmic tip may affect the stability of the trimeric assembly of TolC. Structural asymmetry of the periplasmic gate was observed to be opening size dependent. Asymmetric conformations are found in moderately opening states, while the most and the least opening states are often more symmetric. Finally, it is shown that lowering pH can remarkably stabilize the closed state of the periplasmic gate. Proteins 2014; 82:2169–2179. C 2014 Wiley Periodicals, Inc. V

Key words: TolC; outer membrane protein; MD simulation; conformational change.

INTRODUCTION The Escherichia coli outer membrane protein TolC plays an important role in protein export and drug efflux across the outer membrane.1,2 It acts with an inner membrane transport system pumping substrates out of the cell with the help of membrane fusion proteins.3–5 X-ray crystallography6 revealed that TolC is a homotrimer which folds into a hollow tapered cylinder, comprising a 40 A˚-long bbarrel spanning the outer membrane (the channel domain), a 100 A˚-long a-helical barrel projecting into the periplasm (the tunnel domain), and a third mixed a/b equatorial domain forming a strap around the tunnel [Fig. 1(A)]. The tunnel domain contains 12 a-helices forming six coiled coils, three of which are inclined by 220 toward the tunnel axis, constricting the periplasmic entrance at a closed state.6,7 The narrowest constriction at the entrance is formed by three D374 residues with an effective pore radius of about 2 A˚ [Fig. 1(B)]. A circular network of intra- and interprotomer hydrogen bonds and salt bridges centered on four residues (Y362, R367, D153, and T152) locks the periplasmic entrance in the closed conformation [Fig. 1(C)].

C 2014 WILEY PERIODICALS, INC. V

In spite of the closed conformation revealed by the crystal structure, the periplasmic gate is required to open to facilitate substrate export in vivo.8 Andersen et al.9 explored the opening mechanism of the periplasmic gate by disrupting the circular interaction network at the aperture region of TolC. By introducing several single or double mutations on R367, Y362, D153 and T152, they captured several opening states of TolC which show increased single-channel ion conductance in vitro. Later, crystal structures of the mutants further illustrated the Additional Supporting Information may be found in the online version of this article. Grant sponsor: National Major Basic Research Program of China; grant numbers: 2011CB808505 and 2014CB910201; Grant sponsor: Specialized Research Fund for the Doctoral Program of Higher Education; grant number: 20130071140004; Grant sponsor: Science & Technology Commission of Shanghai Municipality; grant number: 08DZ2270500. *Correspondence to: Jingwei Weng, Department of Chemistry, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, Fudan University, Shanghai, People’s Republic of China. E-mail: [email protected] and Wenning Wang, Department of Chemistry, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, Fudan University, Shanghai, People’s Republic of China. E-mail: [email protected] Received 16 July 2013; Revised 12 March 2014; Accepted 29 March 2014 Published online 7 April 2014 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/prot.24573

PROTEINS

2169

B. Wang et al.

Figure 1 (A) The simulation system of TolC. (B) The periplasimic constrictions and the definition of TCA374 and TCA365 (triangular cross-sectional area). TCA374 is spanned by the Ca atoms of the three Asp374 residues, while TCA365 is spanned by the Ca atoms of the three Gly365 residues. (C) The circular interaction network at the periplasmic end. (D) Time series of the RMSD for the whole protein in the wild type system and the mutant systems. Each system has five trajectories denoted by different colors.

opening states on the atomic scale and revealed that twisting of the a-helices could open the periplasmic gate9 as proposed in the “twist-to-open” transition model.6 However, controversies still remain as these mutants show very different opening states. In the crystal structures of

2170

PROTEINS

mutants R367S (TolCRS hereafter) and Y362F/R367S (TolCYFRS hereafter), the periplasmic gates adopt symmetric conformations.10 The D374 constriction of TolCRS mutant remains closed, while the double mutation Y362F/R367S caused expansion of the pore constriction from 11.3 to

MD Simulation of TolC

50.3 A˚2.10 In the Y362F/R367E mutant (TolCYFRE hereafter), the periplasmic constriction was shown to be partially opened with slight structural asymmetry.11 What is the dominant opening state of TolC and how the mutations induce the gate opening? These questions remain to be answered. Besides the site-directed mutations, variations in pH could also affect the ion conductance of TolC. Lower pH conditions led to significant drop of TolC ion conductance.7,12 This pH dependence further complicates the gating mechanism of the protein. Molecular dynamics (MD) simulations have also been used to investigate the conformational flexibility and gating mechanism of TolC and its mutants.13–15 Equilibration MD simulations of mutants TolCYFRS13 and TolCYFRE14 at physiological salt concentrations did not observe major conformational changes on the timescale of several tens of nanoseconds. It was demonstrated that the high occupation probability of cations (K1 or Na1) at the periplasmic aperture could inhibit the gate opening.14,15 Therefore, to speed up the simulation it would be favorable to use a minimal ion concentration. In this work, we performed MD simulations for WT TolC and six mutants (TolCYF (Y362F), TolCRE (R367E), TolCRS (R367S), TolCDA (D153A), TolCYFRE (Y362F/R367E), and TolCYFRS (Y362F/R367S)) by minimizing the salt concentration toward a solely neutralized system aiming at exploring the conformational space of TolC to the largest extent. To enhance sampling, multicopy simulations were utilized and each system was simulated for 220–370 ns in total. Starting from the closed state, the protein sampled large conformational spaces revealing various conformational states during the close-to-open transition of the periplasmic gate, and both symmetrical and asymmetrical conformations were observed. The effect of pH on TolC was also studied by performing simulations under pH 4.5 or 5.0. Both pH conditions produced a more stable closed conformation at the periplasmic side. MATERIALS AND METHODS The simulation system

All simulations used the crystal structure of the closed state TolC (PDBID: 1EK9). Various mutants were generated on this structure using the Mutator plug-in of VMD.16 The protein was inserted into a 1-palmitoyl, 2oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer following the “shrinking” method.17 The whole system was then solvated with water molecules and sodium cation were randomly positioned to neutralize the simulation system. The final system contains more than 166,000 atoms including 248 POPC molecules (256 before the insertion of the protein) and over 37,800 water molecules. pKa values of titratable residues in TolC were calculated using H1118 and the protonation states of residues were determined according to the calculated

pKa values. All the simulations were done under pH 7.5 if not specified. The protonation states of the residues under this pH is approximately the same as their default states. The MD simulation

All simulations were carried out using NAMD 2.819 or GROMACS 4.6.520 with the CHARMM2721 force field. TIP3P model22 was used for water molecules. In each system, the solvent were first energy minimized for 1000 steps and equilibrated for 200 ps with the protein and lipid fixed, and then the solvent and lipid were energy minimized for 5000 steps and equilibrated for 2 ns with the protein fixed, followed by 1 ns run with the Ca atoms of protein constrained by harmonic forces. Finally, two 20 ns, and three 60 ns production simulations with different initial velocities were conducted without any constraint on the system. In all the simulations, temperature was maintained at 300 K by Langevin dynamics for non-hydrogen atoms with a damping coefficient of 1 ps21. Constant pressure was maintained at 1 bar by the Nose`-Hoover Langevin piston method.23 The van der Waals (vdW) interaction was smoothed at 10 A˚ and truncated at 12 A˚. Long-range electrostatic interactions were treated using the Particle–Mesh Ewald method24 with a 1 A˚ grid spacing. The integration step was set to 2 fs and all bonds involving hydrogen atom was constrained. One additional 150 ns trajectory was conducted for the TolCWT, TolCYF, TolCDA, and TolCYFRS systems, respectively, using GROMACS20 with the CHARMM27 force field. In these simulations, the temperature was maintained at 300 K and the pressure was maintained at 1 bar using the weak coupling method.25 The program VMD16 was used for visualization and data analysis. Targeted MD simulation

Targeted MD guides a known initial structure to a known target structure using an external potential.26 The potential is designed to decrease the root-meansquare deviation (RMSD) of the system relative to the target structure toward a preset value at each step, and can be described as UTMD 5

1k ½RMSDðtÞ2RMSD  ðtÞ2 2N

where RMSD(t) is the instantaneous best-fit RMSD of the current coordinates to the target coordinates, RMSD*(t) is the preset RMSD value for the current time step, while k is the force constant and N the number of targeted atoms. The external forces were imposed on all the 1284 Ca atoms of TolC. The force constant was set to 128.4 or 642 kcal mol21 A˚22 (i.e., 0.1 or 0.5 kcal mol21 A˚22 for each atom), respectively. The time step was set to 1 fs and the coordinates were saved every 500 PROTEINS

2171

B. Wang et al.

steps. Six replicate trajectories of 0.5, 1.0, and 5.0 ns long with different force constants and initial velocities were conducted. Data analysis

The asymmetric coefficient Casy is defined as: 3  X

Casy 5 n51

374 d 374 2d

2

3TCA374

where d 374 is the distance between the Ca atoms of two 374 Asp374 residues, and d is the average value of the three 374 d in the TolC structure. TCA374 is the triangular crosssectional area (TCA) spanned by the Ca atoms of the Asp374 residues. The structures of the first 10 ns of the replicate trajectories are excluded from the clustering analysis for each system. K-Means clustering algorithm was performed with a radius cutoff of 2 A˚ for Ca RMSD using Multiscale Modeling Tools for Structural Biology Tool Set.27 RESULTS MD simulations were performed for the wild-type protein (TolCWT) and six mutants including TolCYF, TolCRE, TolCRS, TolCDA, TolCYFRE, and TolCYFRS. For each system, at least five independent simulations were conducted, comprising two 20 ns and three 60 ns trajectories. The Ca RMSDs relative to the initial structure were calculated for all trajectories [Fig. 1(D)]. The root-meansquare fluctuation (RMSF) analyses (Supporting Information Fig. S1) illustrate that the protein backbone is generally stable with the values of RMSF around 1 A˚ in most regions, while the most flexible parts are located at the periplasmic and the extracellular ends of the protein. Since the RMSDs changed significantly during the first 10 ns for most trajectories, the corresponding snapshots were excluded from data analysis. The effect of mutations on the conformational state of the periplasmic gate

Although simulations of the mutants started from the closed state structure similar to TolCWT, most of them showed larger structural deviations. At the end of 20 ns, most of the replicas of mutants have RMSDs larger than 4 A˚. For each system, three of the five trajectories were elongated to 60 ns. In some cases, the longer simulation time led to larger structural deviations, such as those of TolCDA and TolCYFRS, while the others did not exhibit RMSD increase during the last 40 ns [Fig. 1(D)]. Generally, the mutants showed large structural deviation from the WT crystal structure and remarkable structural diversity. TolCYF is the only exception among the mutants,

2172

PROTEINS

which has high structural stability with RMSD values less than 3 A˚ in all the replicas [Fig. 1(D)]. In order to compare and quantify the conformational changes in these systems, we calculated the average TCA spanned by the Ca atoms of the three D374s or the three G365s in this region [referred to as TCA374 and TCA365, respectively hereafter; Fig. 1(B)]. The average values of both TCA374 and TCA365 show that all systems sampled more opened conformations in the simulations than the crystal structures [Fig. 2(A,B) and Table I]. For TolCWT, the average value of TCA374 (110.7 6 27.2 A˚2) is much larger than that of the TolCWT crystal structure (58.6 A˚2), even larger than that of the crystal structure of TolCYFRS mutant (99.5 A˚2),10 which is the most open conformation reported in crystal structural studies. Previous MD simulation study of TolCWT performed one salt-free simulation for 150 ns, in which the opening of the D374 constriction was also observed, but to a less extent (average TCA374 of 88.7 6 11.6 A˚2).15 For the mutants, our simulation results show that all but TolCYF have larger average values of TCA374 and TCA365 than WT does, indicating that these mutations facilitate the opening of the periplasmic gate [Fig. 2(A,B)]. Enhanced conformational flexibility of the gate in mutants was also manifested by the larger standard deviations of TCAs [Fig. 2(A,B)]. In the previous MD simulation study with applied voltages, the periplasmic mouth of TolC was also observed to open in TolCYFRD and TolCYFRE mutants.14 The values of TCA365 in their simulations for TolCYFRE are smaller than those in our simulations (200 vs. 318.1 A˚2), while the values of TCA374 were not reported.14 Since all the mutants were designed to disrupt the circular network of intra- and interprotomer interactions of Y362–D153, R367–D153, and R367–T152 at the periplasmic bottleneck in the crystal structure of TolC [Fig. 1(C)], one of the surprising observations in our simulations is that the mutant TolCYF has a more constricted periplasmic gate than the TolCWT and showed higher conformational stability [Figs. 1(D) and 2(A,B)]. To examine the roles of the interactions in the conformational changes during simulation, the appearance probabilities of these interactions in the trajectories were calculated [Fig. 2(C)]. In the TolCWT, the interactions between R367 and D153 were most stable, being maintained around 80% of the simulation time, while the other two interaction pairs were relatively weak [Fig. 2(C)]. In TolCYF, the R367–D153 interactions were extremely stable with an appearance probability of more than 90%. This might be attributed to the removal of the hydrogen bond between Y362 and D153 upon Y362F mutation, which allows D153 to form electrostatic and hydrogen bonding interactions with R367 exclusively. It is most likely that the strengthened R367–D153 interaction is responsible for the higher conformational stability of TolCYF. In agreement with this, the other five mutants without this pair of interaction [Fig. 2(C)] have higher

MD Simulation of TolC

Figure 2 The effect of mutations on the periplasmic gate of the protein. (A) The average TCA374 values of the five replicas for each system. The solid lines denote the TCA values in the crystal structures of TolCWT and mutants. (B) The average TCA365 values of the five replicas for each system. (C) The appearance probabilities of the intra- and interprotomer interactions at the periplasmic entrance in various systems.

conformational flexibility at the periplasmic gate and adopt much more opened conformations. The intraprotomer hydrogen bond interaction Y362– D153 was maintained in TolCRS and TolCRE mutants [Fig. 2(C)], although the periplasmic gate opened widely in these two systems. Similarly, the interprotomer interaction R367–T152 was replaced by E367–T152 with even higher strength in mutants TolCRE and TolCYFRE in spite of the significant dilation of the periplasmic constriction in these two systems. Therefore, these results suggest that only the interprotomer interaction R367–D153 is critical for the conformational stability of the periplasmic gate.

Cluster analysis and structural asymmetry evaluation

To further inspect various conformational states of the systems, an RMSD-based cluster analysis was performed for each system (see Methods). Nine clusters with weights higher than 4% were obtained for TolCWT (Table II). These clusters of conformational states differ both in the magnitudes of TCA374 and in the structural asymmetry at the periplasmic gate. A coefficient Casy was defined to quantify the structural asymmetry (see Methods) and calculated for the center structure of each cluster (Table II). For TolCWT, four clusters (C1, C2, C5, and C6) could

Table I

The Periplasmic Entrance TCAs of TolCWT and Mutants Obtained From MD Simulation Trajectories TCA374a (2)

WT YF RE RS DA YFRE YFRS WT (pH 5.0) WT (pH 4.5) a

TCA365b (2)

Average

Minimum

Maximum

Average

Minimum

Maximum

110.7 6 27.2 76.5 6 6.0 154.1 6 24.1 236.2 6 75.7 210.7 6 78.4 166.4 6 105.4 237.9 6 104.5 73.4 6 4.7 71.5 6 4.0

67.7 57.9 96.9 74.8 79.3 54.1 96.3 58.1 57.6

203.5 104.5 280.4 434.6 534.4 441.1 527.7 105.1 89.6

144.5 6 57.9 104.9 6 9.5 232.0 6 50.1 358.3 6 137.7 316.7 6 149.9 318.1 6 180.8 394.3 6 188.6 102.9 6 14.4 123.4 6 11.1

51.3 70.0 88.7 72.1 65.1 132.4 125.2 67.1 75.1

308.2 150.0 498.0 764.3 981.1 915.3 965.5 153.1 177.0

The triangular cross-sectional area spanned by the Ca atoms of the Asp374 residues. The triangular cross-sectional area spanned by the Ca atoms of the Gly365 residues.

b

PROTEINS

2173

B. Wang et al.

Table II

Clusters With a Relatively Large Size (>4%) of the MD Trajectories of TolCWT and its Mutants WT 374

a

Cluster no.

%

C1 C2 C3 C4 C5 C6 C7 C8 C9

28.4 15.9 13.4 12.5 9.1 5.7 5.6 5.2 4.1

TCA (2)

88.1 100.6 156.3 119.3 86.0 87.1 116.2 138.7 156.2

Casy

%

0.0001 0.029 0.271 0.134 0.028 0.004 0.099 0.142 0.373

85.4 8.9 5.7

DA

YF

RE

RS

374

374

Casy 0.225 0.021 0.250 0.081 0.028 0.001 0.140 0.017 0.065

TCA (2)

76.7 71.9 74.9

Casy

%

0.0002 0.001 0.0004

15.7 15.6 15.3 14.2 12.3 8.9 7.5

YFRE

TCA (2)

Casy

%

TCA374 (2)

149.0 130.5 167.6 152.4 135.2 176.9 171.1

0.153 0.074 0.149 0.144 0.267 0.134 0.155

9.8 7.4 7.0 6.2 6.1 5.9 5.9 4.9 4.5

203.8 246.8 170.2 139.6 309.9 89.7 151.8 278.3 327.4

YFRS

Cluster no.

%a

TCA374 (2)

Casy

%

TCA374 (2)

Casy

%

TCA374 (2)

C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14

8.5 6.9 6.8 6.7 6.6 6.6 6.5 6.5 6.2 5.8 5.7 5.1 4.6 4.2

124.7 218.9 95.7 192.0 228.3 128.5 197.4 244.2 274.6 217.5 178.5 174.9 160.7 433.1

0.090 0.171 0.033 0.006 0.006 0.064 0.001 0.045 0.230 0.198 0.046 0.112 0.023 0.046

29.3 17.4 11.9 5.9 5.9 4.6

71.5 116.4 150.4 105.6 106.8 331.7

0.002 0.123 0.224 0.024 0.001 0.047

29.4 5.9 5.9 5.8 5.6 4.7 4.4 4.4 4.2 4.0

132.5 140.1 150.9 401.7 353.3 423.0 373.1 259.6 197.8 291.0

a

Casy 0.127 0.033 0.133 0.065 0.057 0.027 0.074 0.137 0.056 0.088

Percentage of the conformations in the cluster.

be attributed to partially opened states with their TCA374 values falling within the range of 86–100 A˚2. The Casy of these four clusters were very small (0.0001–0.029), and the structure can be considered as symmetric [Fig. 3(A), Table II]. The other five clusters with wider opening gates are much more asymmetric [Fig. 3(A), Table II]. The TolCYF mutant showed higher conformational stability and more closed pore constriction than TolCWT (typical TCA374 less than 80 A˚2), with all three clusters of symmetric states (Table II). The other five mutants sampled much more opened states than TolCWT and TolCYF did [Fig. 2(A,B), Table II]. However, three of them, that is, the TolCRS, TolCDA, and TolCYFRE, have a conformation at the partial opening state similar with that of TolCWT, with TCA374 less than 100 A˚2 (Table II). The C6 cluster of TolCRS is very similar with the C1 cluster of TolCWT [Fig. 3(A)] with a TCA374 value of 89.7 A˚2 and Casy of 0.001 (Table II). The C1 cluster of TolCYFRE has a small value of TCA374 [71.5 A˚2; Fig. 3(A), Table II]. These conformations are reminiscent of the crystal structures of TolCRS and TolCYFRE, which do not show obvious dilation of the constriction relative to TolCWT.10,11 In the other two mutants (TolCRE and TolCYFRS), however, the protein merely sampled widely opening conformations and none of the clusters resembles

2174

PROTEINS

those of TolCWT (Table II). This is generally consistent with the observation that the crystal structure of TolCYFRS adopts more opened conformation than the WT TolC.11 Relationship between the opening size and the structural asymmetry of the periplasmic gate

By analyzing all the conformational clusters of the WT and the mutants, a general trend can be found between the opening size and the structural asymmetry of the periplasmic gate [Fig. 3(B)]. The most constricted conformations in the simulations, that is, those with TCA374 smaller than 120 A˚2, tend to adopt more symmetric structures at the gate, with a Casy below 0.05 [Fig. 3(B)]. The conformations with large structural asymmetry have intermediate size of periplasmic gate, with the values of TCA374 ranging from 120 to 250 A˚2. As the periplasmic gate opens more widely, the structure becomes more symmetric. For example, the C14 cluster of TolCDA and the C6 cluster of TolCYFRS have the most widely opened gates (TCA374 of 433.1 and 423.0 A˚2, respectively) and small asymmetric coefficients [0.046 and 0.027; Table II, Fig. 3(B)]. In these widely opening states, interprotomer interactions below the equatorial domain are all

MD Simulation of TolC

Figure 3 (A) Some representative center structures in cluster analyses of TolCWT and mutants. (B) All the conformational clusters of TolCWT and its mutants mapped on the 2-dimensional space spanned by TCA374 and the asymmetric coefficient Casy. (C) The projection of three targeted MD trajectories of wild type TolC toward the putative open conformation on the above 2-dimensional space. The three trajectories are of 0.5, 1.0, and 5.0 ns long, respectively, colored in a time-dependent manner.

PROTEINS

2175

B. Wang et al.

Figure 4 The 150 ns simulations of TolCWT, TolCYF, TolCDA, TolCYFRS. (A) Time series of the Ca RMSDs with respect to the initial structures of the four systems. (B) Variations of the TCA374 values along the simulation trajectories. (C) Variations of the TCA365 values along the simulation trajectories. (D) The appearance probabilities of the intra- and interprotomer interactions at the periplasmic entrance during the simulations. (E) Conformational snapshots of the four systems mapped on the 2-dimensional space spanned by TCA374 and the asymmetric coefficient Casy.

disrupted, and the periplasmic end splits into three fourhelix bundles, with the intraprotomer packing partly maintained [Fig. 3(A)]. To sum up, both the least and the most opening conformations are more symmetric, while the moderately opening ones could be asymmetric. To verify the symmetry-dependent transition process deduced from the mutants, we performed targeted MD simulations on TolCWT. The closed state was used as the initial structure and a representative widely opening conformation of TolCYFRE was mutated back to WT and set as the target. Three independent simulations of 0.5, 1.0, and 5.0 ns were conducted, respectively, and all of them successfully captured the close-to-open transition. The simulations showed a common feature that the dilation of the pore constriction accompanies variations in the structural symmetry [Fig. 3(C)]. Beginning with a symmetric closed conformation, the gate was moderately opened in the first half of the simulation. Simultaneously, the asymmetric coefficient rapidly increased to more than 0.2 which is comparable to the value of the most asymmetric cluster of the mutants (Table II). As the dilation continued, Casy began to decrease to less than 0.1 at the end of the simulation. In spite of the symmetric initial and target structures, targeted MD simulations sampled asymmetric conformations of the periplasmic gate. Ranging from 120 to 250 A˚2 in TCA374 and from about 0.1 to 0.25 in Casy, these asymmetric

2176

PROTEINS

conformations share similar features with those observed in the above unbiased simulations. The effect of force constant on the targeted MD simulations was also examined, showing that different force constants gave similar results (Supporting Information Fig. S2). Lowering pH increases the conformational stability of TolC

Previous single-channel ion conductance experiment revealed that the conductance decreases with the lowering of the pH.7 This has been explained by the protonation of D371 and D374 at low pH, which may reduce the cation selectivity and the like charge repulsions that helps to expand the tunnel entrance. Here, MD simulations of TolC at pH 5.0 and 4.5 were carried out. The protonation states of all residues were assigned using H1118 (Supporting Information Fig. S3). Specifically, the protonation states of the residues near the periplasmic gate were changed under lower pH levels. Three D374 residues and one of the three D371 residues which were negatively charged in the former simulations were all protonated. Two 20 ns and one 60 ns MD trajectories under pH 4.5 or pH 5.0, respectively, demonstrated significantly decreased conformational flexibility of TolC. The Ca RMSD values were kept below 2 A˚ [Fig. 1(D)] and the structural fluctuation was also evidently suppressed, especially at the periplasmic side

MD Simulation of TolC

(Supporting Information Fig. S1). The average values of TCA374 and TCA365 at pH 5.0 and 4.5 are much smaller than that under pH 7.5 [Table I, Fig. 2(A,B)], and the gate remains symmetric throughout the simulation trajectories. Therefore, the periplasmic gate remains in the closed state under lower pH conditions and the electrostatic repulsion among the ring of six aspartic acids (D371 and D374) is most likely one of the driving force for the opening of the pore constriction. 150 ns simulations of TolCWT, TolCYF, TolCDA, and TolCYFRS are consistent with the findings in the short multicopy simulations

The RMSD profiles of the above 20–60 ns simulations demonstrate obvious lack of convergence in several trajectories [Fig. 1(D)]. In order to verify the main observations based on these short trajectories, we conducted an additional 150 ns trajectory for each of the four systems, including TolCWT, TolCYF, TolCDA, and TolCYFRS, using GROMACS with the same force filed (CHARMM27). The Ca RMSD of the 150 ns simulations clearly show that the trajectories of TolCDA and TolCYFRS still failed to converge at the end of the simulations, while the structures of TolCWT and TolCYF are much more stable with RMSD values less than 2.6 A˚ throughout the simulations [Fig. 4(A)]. In agreement with the multicopy short simulations, the TolCYF exhibited even higher structural stability than the WT protein. The structure variations are mainly due to the conformational change at the periplasmic gate. The variations of the TCA374 and TCA365 values along the trajectories are consistent with those of the RMSDs [Fig. 4(B,C)]. It is worth noting that the TCA values of the four systems in the 150 ns simulations are all smaller than those of the above short simulations. This discrepancy may be due to the systematic difference between NAMD and GROMACS, or simply due to the sampling deficiency. Comparison of the stability of the three interactions, that is, Y362–D153, R367–D153 and R367–T152 that lock the periplasmic gate gives the similar results as those of the short simulations. All these interactions were disrupted in TolCDA and TolCYFRS mutants during the simulations and the critical role of the interaction between R367 and D153 on periplasmic gate stability is confirmed [Fig. 4(D)]. Calculation of the Casy coefficients shows that the relationship between the opening size and the structural asymmetry of the periplasmic gate is the same as that shown in the short simulations [Fig. 4(E)]. The periplasmic gate also experienced an asymmetry stage during the gate opening process [Fig. 4(E)]. DISCUSSION In this study, we performed systematic MD simulations of TolC and its mutants at minimal salt concentra-

tion, which allow us to explore the conformational space of TolC to the largest extent. The multicopy strategy was used to enhance sampling.28 These trajectories revealed much higher conformational flexibility and larger conformational space of TolC than that observed in the crystal structures of TolCWT and mutants,10,11 as well as in the previous simulation studies. Nevertheless, it should be mentioned that the simulations still suffer from the limited sampling. The conformation ensembles obtained in this study are far from a full description of the conformational space of the systems. The sampling deficiency is manifested in several aspects. First, the average TCA values of the five independent trajectories, which are expected to be very similar in the limit of good sampling, differ significantly, especially in the cases of mutant systems (Supporting Information Table S1). Extension of the trajectories to 60 and 150 ns does not decrease the variances between trajectories. Second, in the cluster analysis, the cluster populations of the five independent trajectories only sometimes overlapped with each other, and in many cases, one replica explored a certain region in the conformational space alone and led to an isolated cluster (see TolCYFRE in Table II for instance), indicating the lack of sampling convergence. In the previous multicopy simulation study of TolCWT with much longer simulation trajectories of 150–300 ns long, it was also shown by principle component analysis that different trajectories explored distinct conformational spaces.15 Third, we further analyzed the statistical uncertainty of the averaged TCA374 values of each trajectory using block-averaging analysis.29,30 Supporting Information Figure S4 shows the block-averaged standard errors (BSEs) of the 20 ns trajectories as a function of block length. It is evident that for many trajectories the BSEs failed to converge as the block length increased. As the trajectories were extended to 60 ns, some of them still failed to converge, especially those of TolC mutants, such as TolCRE, TolCRS, TolCDA, and TolCYFRS (Supporting Information Fig. S5). It is noticeable that convergence has not been reached even for the 150 ns simulations (Supporting Information Fig. S6). This is not surprising since the initial conformation of the mutant systems was chosen to be the crystal structure of WT TolC. It is obvious that the structural deviations of the mutants from the WT protein are remarkable and tens of nanoseconds simulation is not sufficient for them to relax to a stable/metastable state. The convergences of the BSEs of TolCWT simulations are much better than those of the mutants, especially for the low pH systems. However, in those trajectories of TolCWT at neutral pH with large structural deviations, BSE estimation still does not converge (Supporting Information Fig. S5). Therefore, the results of BSE calculations imply that the correlation time of periplasmic gate movement is either comparable or longer than hundreds of nanoseconds timescale. Overall, submicrosecond PROTEINS

2177

B. Wang et al.

simulations are not long enough to describe equilibrium behavior of the conformational movement of the periplasmic gate of TolC. Specifically, the periplasmic gate may adopt multiple distinct, well-populated conformations, rather than fluctuating about a single average conformation similar to the crystal structure. This makes the equilibrium sampling extremely difficult with MD simulations. Previous simulation study of membrane protein rhodopsin demonstrated that even the local conformational movement of a loop region is not well sampled in the simulation consisting of 26 independent trajectories of 100 ns long.31 Considering the size of the system and the high conformational flexibility of TolC, unbiased equilibrated sampling of the conformational space at atomic level is very difficult to achieve. In light of the limited sampling, the results of the simulations should be interpreted more cautiously. For example, the results of the cluster analysis do not represent the equilibrium population of conformational states of TolC or its mutants. Nevertheless, as mentioned above, the current simulations have already revealed large conformational spaces of TolC and its mutants away from the crystal structure, which are not observed in the previous studies. In addition, the parallel simulations of seven systems provide hints in understanding the conformational dynamics and working mechanism of TolC. Previous simulation studies discovered that cations occupied at the periplasmic aperture can inhibit the periplasmic gate opening of TolC.14,15 In our simulations, only 17 sodium ions were included to achieve the neutral environment in the system. The density profile of the sodium ions of our simulations [Supporting Information Fig. S7(A)] is very similar to that of the high salt concentration simulations.15 The relationship between the binding probability of Na1 at D374/D371 and the value of TCA374 was plotted in Supporting Information Figure S7(B) for each trajectory of all systems. Correlation exists in some cases. For example, the trajectory with the most widely opening conformation of TolCYFRE has the lowest Na1 residence probability at D374/D371 in TolCYFRE. Generally, the correlation is more obvious in mutants TolCRS, TolCDA, TolCYFRE, and TolCYFRS, while in WT TolC, TolCYF, and TolCRE mutants, the effect of Na1 residence on the gate opening is obscure [Supporting Information Fig. S7(B)]. These observations suggest that the cation occupation plays an important role in gate opening, but it may not the only factor regulating the conformational change. In agreement with this, the previous study showed that the K1 occupation frequency threshold is not applicable to TolCRD and TolCYF mutants.14 The WT TolC sampled a large conformational space, and the most opening conformation (TCA374 5 203.5 A˚2) expands more than three times the size at the D374 constriction relative to the crystal structure (TCA374 5 58.6 A˚2). These simulation results provide a different scenario with respect to the current AcrAB–TolC complex model.

2178

PROTEINS

The current model suggests that the TolC periplasmic entrance remains closed until the interaction with AcrAB triggers the close-to-open transition.32 However, the simulation results suggest that TolC could also be widely opened without AcrAB. Recently, Krishnamoorthy et al.33 provided in vivo experimental data showing that in the absence of AcrAB, substrate molecules can still access the residues at the periplasmic constriction, such as D374. Their results suggest that TolC periplasmic gate is not tightly occluded in vivo. It is worth noting, however, that the time scale of the conformational transitions at physiological salt concentration could be much longer since the previous simulation studies at higher salt concentration failed to visit the open conformations. Consistent with the experimental evidences, most of the mutations which aim to disrupt the bonds tethering the closely packed coiled coils at the periplasmic tip greatly facilitated the dilation of the pore constriction and yielded various open states at the periplasmic gate [Fig. 2(A,B)]. Five of the six mutants showed larger TCAs than the WT [Fig. 2(A,B)], which can be justified by their increase in the single channel conductance.9,12 The degree of pore dilation, however, is much higher than that seen in the crystal structures of three TolC mutants.10,11 In the widely opened conformational states of TolCYFRS, TolCDA, TolCRS, and TolCYFRE mutants, three monomers dissociate completely under the equatorial domain. This is in agreement with the recent experimental data showing that the trimeric AcrAB–TolC assembly is unstable in the TolCYFRE mutant,33,34 and may be misfolded in vivo.33 The high conformational stability of TolCYF and the close correlation between R367–D153 interactions and gate dilation [Fig. 2(C)] suggests that R367–D153 is a critical link for the closed periplasmic gate. This is in agreement with the experimental results that T152V mutant has the same ion conductance as TolCWT.9 However, in the same experimental study, TolCYF showed slightly higher conductance than TolCWT.9 One possible explanation is that the conductance is not simply proportional to the size of the periplasmic gate. Besides the pore radius, the electrostatic interactions between ions and charged residues would also affect the permeation free energy. Previous studies showed that pH condition affects ion conductance significantly.7 Our simulations indicate that this could result from the change of the conformational flexibility of TolC. Both systems with protonation states at low pH conditions showed much smaller TCAs than the pH 7.5 systems, indicating an occluded periplasmic gate. One of the important features in the low-pH systems are the protonated D371 and D374 residues (Supporting Information Fig. S3). On the other hand, it is worth noting that in the simulation systems of low pH conditions, there are no sodium ions at all and the cation inhibited gate opening mechanism can be excluded. We propose that the charge repulsions between the acidic

MD Simulation of TolC

residues (D371 and D374) would greatly affect the state of the periplasmic gate and regulate the substrate export process. Besides TolC, the conformational flexibility of membrane fusion protein, AcrA was also shown to be regulated by pH conditions.35 This reinforces the model in which the proton flux across the cytoplasmic membrane is the signal to trigger the pH-dependent conformational changes of the tripartite AcrAB–TolC complex, regulating the complex assembly and drug efflux. ACKNOWLEDGMENTS We thank the super computer center of Shenzhen and of Fudan University for their allocation of computer time. REFERENCES 1. Schuldiner S. Structural biology—The ins and outs of drug transport. Nature 2006;443:156–157. 2. Higgins CF. Multiple molecular mechanisms for multidrug resistance transporters. Nature 2007;446:749–757. 3. Koronakis V. TolC—the bacterial exit duct for proteins and drugs. FEBS Lett 2003;555:66–71. 4. Koronakis V, Eswaran J, Hughes C. Structure andfunction of tolC: the bacterial exit duct for proteins and drugs. Annu Rev Biochem 2004;73:467–489. 5. Zgurskaya HI. Multicomponent drug efflux complexes: architecture and mechanism of assembly. Future Microbiol 2009;4:919–932. 6. Koronakis V, Sharff A, Koronakis E, Luisi B, Hughes C. Crystal structure of the bacterial membrane protein TolC central to multidrug efflux and protein export. Nature 2000;405:914–919. 7. Andersen C, Hughes C, Koronakis V. Electrophysiological behavior of the TolC channel-tunnel in planar lipid bilayers. J Membr Biol 2002;185:83–92. 8. Eswaran J, Hughes C, Koronakis V. Locking TolC entrance helices to prevent protein translocation by the bacterial type I export apparatus. J Mol Biol 2003;327:309–315. 9. Andersen C, Koronakis E, Bokma E, Eswaran J, Hymphreys D, Hughes C, Koronakis V. Transition to the open state of the ToIC periplasmic tunnel entrance. Proc Natl Acad Sci USA 2002;99: 11103–11108. 10. Pei XY, Hinchliffe P, Symmons MF, Koronakis E, Benz R, Hughes C, Koronakis V. Structures of sequential open states in a symmetrical opening transition of the TolC exit duct. Proc Natl Acad Sci USA 2011;108:2112–2117. 11. Bavro VN, Pietras Z, Furnham N, Perez-Cano L, Fernandez-Recio J, Pei XY, Misra R, Luisi B. Assembly and channel opening in a bacterial drug efflux machine. Mol Cell 2008;30:114–121. 12. Andersen C, Koronakis E, Hughes C, Koronakis V. An aspartate ring at the TolC tunnel entrance determines ion selectivity and presents a target for blocking by large cations. Mol Microbiol 2002; 44:1131–1139. 13. Vaccaro L, Scott KA, Sansom MSP. Gating at both ends and breathing in the middle: conformational dynamics of TolC. Biophys J 2008;95:5681–5691. 14. Schulz R, Kleinekathofer U. Transitions between closed and open conformations of ToIC: the effects of ions in simulations. Biophys J 2009;96:3116–3125. 15. Raunest M, Kandt C. Locked on one side only: ground state dynamics of the outer membrane efflux duct TolC. Biochemistry 2012;51:1719–1729.

16. Humphrey W, Dalke A, Schulten K. VMD: Visual molecular dynamics. J Mol Graph 1996;14:33–38. 17. Kandt C, Ash WL, Tieleman DP. Setting up and running molecular dynamics simulations of membrane proteins. Methods 2007;41:475– 488. 18. Anandakrishnan R, Aguilar B, Onufriev AV. H113.0: automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations. Nucleic Acids Res 2012;40:W537–W541. 19. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kale L, Schulten K. Scalable molecular dynamics with NAMD. J Comput Chem 2005;26:1781–1802. 20. Pronk S, Pall S, Schulz R, Larsson P, Bjelkmar P, Apostolov R, Shirts MR, Smith JC, Kasson PM, van der Spoel D, Hess B, Lindahl E. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013;29:845– 854. 21. MacKerell AD, Feig M, Brooks CL. Improved treatment of the protein backbone in empirical force fields. J Am Chem Soc 2004;126: 698–699. 22. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. J Chem Phys 1983;79:926–935. 23. Feller SE, Zhang YH, Pastor RW, Brooks BR. Constant-pressure molecular-dynamics simulation—the langevin piston method. J Chem Phys 1995;103:4613–4621. 24. Darden T, York D, Pedersen L. Particle mesh ewald—an n.log(n) method for ewald sums in large systems. J Chem Phys 1993;98: 10089–10092. 25. Berendsen HJC, Postma JPM, Vangunsteren WF, Dinola A, Haak JR. Molecular dynamics with coupling to an external bath. J Chem Phys 1984;81:3684–3690. 26. Schlitter J, Engels M, Kruger P, Jacoby E, Wollmer A. Targeted molecular dynamics simulation of conformational change—application to the T[–]R transition in insulin. Mol Simul 1993;10:291–308. 27. Feig M, Karanicolas J, Brooks CL. MMTSB Tool Set: enhanced sampling and multiscale modeling methods for applications in structural biology. J Mol Graph 2004;22:377–395. 28. Caves LSD, Evanseck JD, Karplus M. Locally accessible conformations of proteins: multiple molecular dynamics simulations of crambin. Protein Sci 1998;7:649–666. 29. Flyvbjerg H, Petersen HG. Error-estimates on averages of correlated data. J Chem Phys 1989;91:461–466. 30. Grossfield A, Zuckerman DM. Quantifying uncertainty and sampling quality in biomolecular simulations. Annu Rep Comput Chem 2009;5:23–48. 31. Grossfield A, Feller SE, Pitman MC. Convergence of molecular dynamics simulations of membrane proteins. Proteins: Struct Funct Bioinform 2007;67:31–40. 32. Balakrishnan L, Hughes C, Koronakis V. Substrate-triggered recruitment of the TolC channel-tunnel during type I export of hemolysin by Escherichia coli. J Mol Biol 2001;313:501–510. 33. Krishnamoorthy G, Tikhonova EB, Dhamdhere G, Zgurskaya HI. On the role of TolC in multidrug efflux: the function and assembly of AcrAB–TolC tolerate significant depletion of intracellular TolC protein. Mol Microbiol 2013;87:982–997. 34. Tikhonova EB, Yamada Y, Zgurskaya HI. Sequential mechanism of assembly of multidrug efflux pump AcrAB-ToIC. Chem Biol 2011; 18:454–463. 35. Wang BB, Weng JW, Fan KN, Wang WN. Interdomain flexibility and pH-induced conformational changes of AcrA revealed by molecular dynamics simulations. J Phys Chem B 2012;116:3411– 3420.

PROTEINS

2179