Elucidating the molecular basis of spontaneous activation in an engineered mechanosensitive channel.

Mechanosensitive channel of large conductance (MscL) detects and responds to changes in the pressure profile of cellular membranes and transduces the mechanical energy into electrical and/or chemical signals. MscL can be activated using ultrasonic or chemical activation methods to improve the absorption of medicines and bioactive compounds into cells. However, re-engineering chemical signals such as pH change can trigger channel activation in MscL. This study elucidates the activation mechanism of an engineered MscL at an atomic level through a combination of equilibrium and non-equilibrium (NE) molecular dynamics (MD) simulations. Comparing the wild-type (WT) and engineered MscL activation processes suggests that the two systems are likely associated with different active states and different transition pathways. These findings indicate that (1) periplasmic loops play a key role in the activation process of MscL, (2) the loss of various backbone-backbone hydrogen bonds and salt bridge interactions in the engineered MscL channel causes the spontaneous opening of the channel, and (3) the most significant interactions lost during the activation process are between the transmembrane helices 1 and 2 in engineered MscL channel. The orientation-based biasing approach for producing and optimizing an open MscL model used in this work is a promising way to characterize unknown protein functional states and investigate the activation processes in ion channels and transmembrane proteins in general. This work paves the way for a computational framework for engineering more efficient pH-sensing mechanosensitive channels.

Introduction

Mechanosensitive channel of large conductance (MscL) is an 80 kDa homopentameric membrane protein with each protomer consisting of two TM -helices (Fig. 1) [2], [3]. MscL channels identify and react to changes in the pressure profile of the cell membrane and transduce the mechanical energy into electrical or chemical signals. [4], [5], [6], [7]. Ions, water, and even small proteins can travel through these protein pores [8]. In bacteria, MscL channels serve as an emergency release valve for acute osmolarity in the system, preserving osmotic homeostasis [2]. MscL has been a prominent drug target for many antibiotics [9], [10] due to its important role in preserving bacterial homeostasis [11]. MscL is necessary for the growth of bacteria and any mutation could cause changes in channel kinetics, shifts in the pressure sensitivity curve, and depressed growth rate [10], [12], [13]. MscL is specifically found in many bacterial cell members but is absent in mammalian genomes [10], [12].

Fig. 1

Engineered MscL channel. (A) Cartoon representation of an engineered MscL based on the crystal structure of wild-type MscL from Mycobacterium tuberculosis (Tb-MscL) [1]. Blue and red represent transmembrane (TM) helices 1 and 2, respectively. (B) A representative protomer of engineered Tb-MscL. (C) Transmembrane helix-1 (TM1) of engineered Tb-MscL. A positively charged [2-(trimethylammonium) ethyl] methane thiosulfonate bromide (MTSET) label (shown in ball and stick representation) was attached to A20C of TM1 of Tb-MscL. (D) Different interhelical angles calculated in this study: is the angle between two neighboring TM1 helices, is the angle between two neighboring TM2 helices, is the angle between TM1 and TM2 of neighboring protomers, and is the angle between TM1 and TM2 of the same protomer. A and B represent two neighboring protomers.

MscL channels can be expressed genetically in mammalian cells [14]. MscL can be used for improving the absorption of membrane-impermeable drugs and bioactive materials into cells via ultrasonic [15] or chemical [14] activation processes. The introduction of bacterial MscL in mice has shown a reduced metastasis in the lungs [16]. Furthermore, engineered MscL expressed in rat hippocampal neurons is activated by a low-pressure ultrasound pulse, which controls neuronal excitation [15], [17]. The engineered MscL is more effective in reducing metastasis than WT at low mechanical pressure ultrasound pulse [15]. Thus, engineered MscL can be a viable target in this study because of its ability to activate at low mechanical pressure [17].

By redesigning, the activation of MscL can be set off by chemical signals. For example, pH change [18], [19], [20], [21], [22], [23], [24], [25], which is the basis for using an engineered MscL as a pH-sensing nanovalve in a drug delivery liposome (DDL) [19], [26], [3] (Fig. 1). Targeted modification of cysteine-containing subunits with positively charged thiol groups has been determined to be functional when inserted in liposomes [27]. Under physiological circumstances, a G22C mutant of the Escherichia coli MscL (Ec-MscL) with MTSET labels resulted in the spontaneous opening of the channel even in the absence of a natural trigger [23], [25], [3]. G22 in Ec-MscL corresponds to residue A20 in Mycobacterium tuberculosis MscL (Tb-MscL) [28], [29], [30]. The charge-sensitive activation of A20C in Tb-MscL has also been observed experimentally [29]. In the absence of tension, MTSET triggers spontaneous activation of the A20C Tb-MscL mutant [29]. It has been demonstrated that varying p and hydrophobicity of the labels attached to the G22C engineered MscL mutant results in activation of the channel at different desired pH conditions [19], [20], [26].

MD is a powerful computational tool that is used, among other things, to study various therapeutically important targets at an atomic level [31], [32], [33], [34], [35], [36], [37], [38], [39], [40]. Several MD studies using coarse-grained MD [41], [42], [43] or biased MD [44], [45], [46], [47], [48], [49], [50], [51], [52], [53], [54], [55], [56], [57], have been conducted thus far for studying the activation of wild-type (WT) MscL. Other computational approaches such as static molecular modeling [58] and continuum models [59], [60], [61] have also been used to study the WT MscL activation. However, no study has offered atomic-level details of the engineered channel activation. In this study, the activation mechanism (i.e., transitioning from closed to open state) of Tb-MscL (MscL from Mycobacterium tuberculosis) was investigated at an atomic level for both WT and engineered MscL. For this purpose, residue 20 of all five protomers of Tb-MscL was mutated to cysteine and pH-sensing MTSET labels were attached (Fig. 1 B, C). Characterizing large-scale structural transformations of membrane proteins without compromising the full atomic details is a major problem for conventional computational approaches. Our computational protocol comprised of all-atom equilibrium and non-equilibrium (NE) MD simulations were employed to fully describe the activation process of Tb-MscL. Please note that an X-ray crystal structure of open state MscL is not yet available, however, a closed state crystal structure is available [1]. Therefore, the closed state MscL structure was used for this study. As a control, the above-mentioned computational protocol was implemented on the WT MscL as well.

First, equilibrium MD simulations were conducted, but the channels failed to open completely despite the presence of MTSET labels. Therefore, NE MD simulations were conducted to capture the open state. Further, the open state structures resulting from the NE simulations were relaxed without any restraints (we call them follow-up equilibrium (FUE) simulations) to test the stability of their respective open state conformations. The engineered MscL open state structure remained open in the FUE simulations, whereas the WT open state structure collapsed back to its original closed state. The presence of MTSET labels in the engineered MscL channel resisted the closing of the open state structure. Based on the elaborate computational protocol explained above, a detailed mechanism for activation of the engineered MscL channel has been proposed. The key aspects of this proposed mechanism are that: (1) the extracellular loops play a major role, and (2) the loss of several hydrogen bond (H-Bond) s’ and salt bridges (SB) facilitate transitioning of the channel from a closed to an open state. Overall, this computational study can help explain different aspects of the MscL channel that could not be completely understood due to the limitations of experimental techniques. For instance, the higher dynamic activity of the engineered MscL compared to the WT.

In the next section, methodological details are provided, followed by results and discussion, conclusions and acknowledgments.

Methods

Crystal structure of Tb-MscL in the closed (inactive) state (PDB2OAR) [1] was downloaded from rcsb.org. Initially, the protein was prepared using the Molecular Operating Environment (MOE) software [62]; crystal waters were removed, appropriate protonation states for the residues were assigned using the protonate3D facility, and hydrogens and other missing atoms were added. Further, CHARMM-GUI web-server [63], [64] was used to prepare the protein for MD simulations. The protein was placed in the 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane bilayer, solvated with TIP3P waters, and 0.15 mM NaCl was added to the system. Overall, the total size of the system was . The final system contains: 603 lipids (298 lipids in the upper leaflet and 305 lipids in the lower leaflet), 60,027 TIP3P [65] waters, 179 sodium (Na+), and 184 chloride (Cl−) ions. The total number of atoms in the system were 270,524.

NAMD 2.10/2.13 [66] was used to simulate the system with periodic boundary conditions in the NPT ensemble at 310 K using the Langevin integrator with a damping coefficient of and 1 atm pressure was maintained using the Nosé-Hoover Langevin piston method [67], [68]. A 2 fs time step used was used and the system was equilibrated for 1000 ns. CHARMM36 all-atom additive force field parameters were used to simulate the entire system [69], [70]. Prior to equilibration, each parent system (Set-1) was first energy minimized for 10,000 steps using the conjugate gradient algorithm [71]. Further, the parent structures were relaxed using a sequence of 1 ns [63] multi-step restraining simulations in an NVT ensemble. The non-bonded interactions were cutoff at 10–12 Å and the long-range electrostatic interactions were computed with the particle mesh Ewald (PME) method [72].

A second system (i.e., the first engineered MscL) was generated (5-MTSET) by mutating A20 in each MscL protomer to cysteine. Next, a pH-sensing label, MTSET, was attached to each mutated residue (total five labels). Further, the system was prepared for MD simulations via the procedure explained above for the WT protein. The size of the system and the total number of atoms were approximately similar to the WT, except that engineered MscL has five extra labels attached to the cysteines at position 20. CHARMM General Force Field (CGenFF) [73] parameters were used for the MTSET labels.

A third system (i.e., the second engineered MscL) was also prepared (1-MTSET), which serves as a control. In this system, A20 in only one protomer was mutated to cysteine and one MTSET label was attached to it. The 1-MTSET system was also prepared using the similar protocol explained above. The 5-MTSET and 1-MTSET parent systems were each equilibrated for 1 as well. Two additional simulation repeats (Sets 2 & 3) were also conducted for each of the three systems listed above. Sets 2 & 3 were continued from the parent simulations by using the configurations at 150 ns and 200 ns, respectively, as the input (Table 1). Since the labels attached to the protomers two and four in 5-MTSET were oriented differently from the rest, they were manually reoriented similarly to the other three labels, and we call this system 5-MTSET(*) (Fig. S1), and this system was equilibrated for 500 ns. Sets 2 & 3 were not carried out for this system.

Table 1

Equilibrium MD simulations.

	WT	5-MTSET	1-MTSET	5-MTSET(*)
	(ns)	(ns)	(ns)	(ns)
Set-1	1000	1000	1000	500
Set-2	500	500	500
Set-3	500	500	500

Sets 2 & 3 were started from Set-1. Only one set of simulations was conducted for the 5-MTSET(*) system.

One data point/ns was collected for statistical analysis. Trajectories were visualized using VMD [74]. H-Bond and salt bridge analysis was conducted via VMD plugins [74]. The cut-off distance and angle were 3.5 Å and 30, respectively, for the H-Bond analysis. Only one H-Bond for an interaction pair was counted. For salt bridge analysis, the cut-off distance was 4 Å and the distance between oxygen atoms of the acidic residues and nitrogens of the basic residues was calculated. Principal components analysis (PCA) [75], [76] was carried out using PRODY software [77] on an ensemble of DCD structures and 20 modes were generated. Only  atoms were considered for the PCA calculations.

The NE simulation technique used in this work provides a foundation for better biasing techniques through a series of short simulations. As part of this effort, several mechanically suitable system-specific reaction coordinates for MscL were identified (Please note that the terms reaction coordinates and collective variables (colvars) were used synonymously). The equilibrated structures of WT, 5-MTSET, and 5-MTSET(*) have been opened in a 100 ns NE simulation, each using orientation quaternions of different helices, bundles of helices, or their linear combinations. A set of orientations were used including Q in which Q are orientation quaternions of i TM helix [78], [79]. Quaternion Q=(q0, q1, q2, q3) can be regarded as a composite of a scalar q0 and an ordinary vector q=(q1, q2, q3), or a complex number Q  = q0 + q1i + q2j + q3k with three distinct imaginary parts. A quaternion unit will define the perfect rotation to overlay one set of coordinates, i.e., an orientation quaternion in which the unit is an optimal angle and axis of rotation [78], [79]. A single-molecule FRET (smFRET) based open state model of MscL [80] was used as a target structure to open the pre-equilibrated MscL structures. Each simulation was repeated twice and a force constant of 10,000  was used. The open state models resulting from the NE pulling simulations were further equilibrated for 260 ns, each, without any restraints to validate its conformational stability (Table 2).

Table 2

Protocol for NE and FUE simulations.

Step	Force Constant (K)	Time
	(kcal/mol.rad^2)	(ns)
1	10,000	100
2	10,000	5
3	10,000 to 1000	5
4	1000 to 0	10
5	0	260

Step 1 represents NE pulling simulations; in steps, 2 to 4 structures resulting from step 1 were relaxed by releasing the restraints in a step-wise manner; step 5 represents FUE simulation.

Results and Discussion

MscL spontaneously responds to the attachment of MTSET labels by starting to open

First, to examine the stability during equilibrium simulations, protein RMSD (Fig. S2) was measured. The RMSD of the WT MscL in all three sets remained under 5 Å (Fig. S2A). In the case of 1-MTSET and 5-MTSET, RMSDs plateaued between 6 and 8 Å (Fig. S2B, C) with the latter showing a greater RMSD values than the former on average. This shows that one label may trigger an opening, although five labels trigger a greater response than one and could potentially result in a faster opening (kinetically favored) and/or a more probable opening (thermodynamically favored). The RMSD of the whole protein follows the same pattern as that of the extracellular loops (ECL) (Fig. 2 A-C), highlighting the fact that the ECLs were driving or dominating the changes in the protein. This is interesting, although the labels were attached at the intracellular side, the immediate impact was seen not just on the intracellular side but also on the extracellular side, which was evidenced by the immediate rise in RMSD of the ECLs.

Fig. 2

RMSD () and radius of gyration (R) of extracellular loops in equilibrium MD simulations. (A-C) The root mean squared deviations (RMSD) s of all WT systems (A) were below 11 Å, whereas for the engineered systems (B-C) RMSDs were all greater than 11 Å. (D-F) The R of all WT systems were above 20.5 Å (E-F), whereas for the engineered systems it was fluctuating between 16.5 and 19 Å (D). (E-F) Representative snapshots of WT (G), 1-MTSET (H), and 5-MTSTET(*) (I) systems are shown.

To understand the effect of labels on a specific area of the protein, root mean squared fluctuations (RMSF) of each residue was measured. In all systems, the fluctuations of ECLs were dominated compared to the rest, although the extent of fluctuation varies between the systems. RMSF of the whole protein in the case of WT varied between 0.5–6 Å (Fig. S3A); for 1-MTSET RMSF fluctuated between 0.5–7 Å  except for first the subunit 1 (S1), where it was fluctuating between 0.5–11 Å (Fig. S3B). For the 5-MTSET, RMSF of all subunits were evenly fluctuating between 0.5–7 Å, except for subunit 2 (S2), where it peaked at 9.5 Å (Fig. S3C). RMSF of TM helices in WT (Fig. S3D), 1-MTSET (Fig. S3E), and other systems (Fig. S3F) were varying between 0.5–2 Å, 0.7–3 Å, and 1.0–3.5 Å, respectively. This shows that the impact of the labels in the 5-MTSET systems was symmetrical on the TM helices, although it varies in the ECLs.

MTSET labels trigger a partial spontaneous opening of the channel

Previous computational studies have shown that opening of MscL results in the hydration of the channel pore. [81], [44], [82], [83] To estimate the extent of opening of the channel, water content across the pore was calculated. In all sets of equilibrium simulations, the WT MscL remained completely closed, and there were no waters around the intracellular gate region (Z = to Å, Fig. S2E). In the case of 1-MTSET, the single label failed to open the channel, although there were some waters in the gate region (Fig. S2G). In the case of 5-MTSET (sets 1–3), the channel was more open compared to the WT and 1-MTSET (Fig. S2F), which is consistent with the established phenomenology of pore hydration. In Fig. S2H water content around the intracellular gate region was specifically compared (Figs. S2B). Overall, the extent of the opening of the channel near the intracellular gate region was as follows: 5-MTSET  1-MTSET  WT (Fig. S2H), while the opening on the extracellular side (i.e., near Z  = 15–20 Å) was somewhat in the reverse order (i.e., WT 5-MTSET  1-MTSET (Figs. S2 E–G)). This was due to the collapse of ECLs into the center of the channel pore in the case of 1-MTSET and 5-MTSET and blocking the pore on the extracellular side, which did not happen in WT (Figs. 2 D-F). Overall, based on water content analysis, we conclude that five labels were more effective than 1-MTSET in the opening of the channel. Therefore, 1-MTSET system was not considered for further analysis in the rest of the paper.

MTSET labels disrupt H-Bond network in MscL

A detailed hydrogen bond (H-Bond) interaction analysis was carried out and discovered that the H-Bond interaction pattern was disturbed by MTSET labels in the MscL channel, particularly the backbone-backbone (BB-BB) H-Bonds. The total number of unique BB-BB H-Bonds in the WT was greater than in the 5-MTSET (5%), and this trend was observed in the entire trajectory (Figs. 3A-C). The trend was also reflected in TM1 and TM2 of all protomers, except for TM1s of S3 and S5 and TM2 of S4 (Figs. 3 D, E); TM1 and TM2 refers to the transmembrane helices 1 & 2. This has contributed to our conclusion that the loss of BB-BB H-Bonds in the engineered MscL makes the TM helices more flexible, which in turn might facilitate the conformational changes that were expected to happen during the activation process of the channel.

Fig. 3

WT vs. 5-MTSET hydrogen bond (H-Bond) interaction analysis. (A-C) Time series of backbone-backbone (BB-BB) H-Bonds of WT and 5-MTSET systems. The three MD trials were separately shown in panels A, B, and C. (D-E) BB-BB H-Bonds of TM1 and TM2 helices of WT and 5-MTSET, respectively. Interactions for individual monomers were shown separately. Interactions with 70% interaction frequency were considered. Sets 1–3 were represented as circles, crosses, and triangles, respectively. S1–S5 represents the five protomers, respectively. (F-G) Comparison of N70-N44 and D16-Y94 H-Bonds. D16(TM1, i)-Y94(TM2, i + 4) interaction was on the intracellular side near the labels. i and i + 4 refer to the first and fifth protomers. N44(TM1, i)-N70(TM2, i) interaction was on the extracellular side. (H) The first and fifth protomers were shown in cyan and yellow, respectively. N44, N70, and N94 were colored green, whereas D16 was colored red. WT and 5-MTSET were represented in black and magenta, respectively.

Two inter-helical side chain-side chain (SC-SC) H-Bonds (N70-N44 and D16-Y94) were identified with a greater interaction frequency in WT compared to 5-MTSET. N70-N44 was formed between the TM1 and TM2 of the same protomer (denoted by (i, i); i refers to a protomer) on the extracellular side (Fig. 3H). D16-Y94 was formed between the TM1 (D16) and TM2 (Y94) of neighboring protomers (denoted by (i,i + 4), i.e., TM1 of protomer i interacts with TM2 of protomer i + 4) on the intracellular side (Fig. 3H). The average H-Bond frequency of N70-N44 in the MD simulation sets 1, 2, and 3 in the WT was 50%, 40%, and 52%, and that of the 5-MTSET was 20%, 22%, and 39%, respectively (Fig. 3F); the D16-Y94 H-Bond frequency in WT and 5-MTSET were 49%, 49%, and 40% and 10%, 21%, and 19%, respectively (Fig. 3G). The N70-N44 H-Bond might play a key role in keeping the TM1 and TM2 of the monomer intact on the extracellular side and breaking of this interaction might be crucial for the channel to open, which was also reflected in our analysis. The breaking of D16-Y94 H-Bond in 5-MTSET probably facilitates the formation of R11-E104 (i,i + 1) and E7-R98 (i,i + 2) salt bridges, which are discussed in more detail in the next section.

MTSET labels facilitate the rearrangement of salt bridge interactions

Salt bridge interaction analysis was done in addition to the H-Bond analysis. Instead of considering the entire trajectory, trajectories ranging from 300–500 ns and 800–1000 ns were examined separately (we refer to them as A and B, respectively, in Tables S1, S2). This would help us identify salt bridges that are important for protein stability as well as for conformational differences between WT and 5-MTSET. The salt bridges that were identified in this analysis were categorized into four different classes for the sake of comparison as the following: TM1(i)-TM2(i + 1), TM1(i)-TM2(i + 2), TM2(i)-TM2(i), and TM1(i)-ECL(i + 1). Here i, i + 1, i + 2, i + 3, and i + 4 refer to the five monomers of the channel. These four classes of salt bridges were further categorized as inter and intra-unit salt bridges. Inter-unit salt bridges are interactions between neighboring monomers/units (i.e., TM1(i)-TM2(i + 1), TM1(i)-TM2(i + 2), and TM1(i)-ECL(i + 1); Table S1) and intra-unit salt brides are interactions within the same monomer/unit (i.e., TM2(i)-TM2(i); Table S2).

As shown in Table S1 there was no significant difference between the WT and 5-MTSET in the total number of inter-unit salt bridge s identified ( 7 in either case). However, there was a difference in the distribution of the salt bridges between the WT and 5-MTSET systems. For example, in the case of 5-MTSET, salt bridge s belonging to all three classes were identified, whereas in the WT all identified salt bridge s belong to the classes TM1(i)-TM2(i + 2) and TM1(i)-ECL(i + 1) and no salt bridge s belonging to the class TM1(i)-TM2(i + 1) were found, particularly in part B of the simulation. We believe that the difference in the distribution of these salt bridge s contributed to the extent of opening of the WT and 5-MTSET channels.

Among the salt bridge s that were identified in this study, the salt bridge s we thought to play a role were SB1 and SB7. SB1 (R11-E104) belongs to the class TM1(i)-TM2(i + 1) and is located on the intracellular side near the bottleneck region (Fig. 4A). Two SB1s were identified in 5-MTSET and these were completely absent in the WT. We propose that the MTSET labels facilitated the formation of SB1s in 5-MTSET and that they play a key role in the opening of the channel pore as they are located near the bottle neck. SB7 belongs to the class TM1(i)-ECL(i + 1). SB7 interaction ties the neighboring monomers/units on the extracellular side. Hence, the loss of this interaction on the extracellular side was key for the opening of the channel. In the WT, there were two SB7 (R45-D53) salt bridges keeping the channel intact, whereas these were lost (total zero) in 5-MTSET facilitating the opening of the channel in 5-MTSET. However, two new salt bridge s belonging to the class TM1(i)-ECL(i + 1) (SB5 and SB6; Fig. 4D) were formed in 5-MTSET on the extracellular side; these two were totally absent in the WT. We believe that these two new salt bridge s (SB5 and SB6) were blocking the complete opening of the engineered MscL channel on the extracellular side.

Fig. 4

Salt bridge (SB) interactions that are playing a key role in the activation of the engineered MscL channel in equilibrium MD simulations. Positive and negatively charged residues are colored in blue and red, respectively. Panels A–C represent the WT and panels D-F represent the engineered MscL (5-MTSET). TM1 and TM2 helices are shown as ribbons and cylinders, respectively.

Only one class of intra-unit salt bridges (i.e., TM2(i)-TM2(i)) were identified in our simulations, i.e., D108-R98 (denoted SB8, Fig. 4B). Two SB8 interactions were identified in the WT, whereas only one was found in the 5-MTSET (Table S2). The loss of SB8 interaction in 5-MTSET, which was present at the intersection of TM helices and the intracellular helices (IH), might have facilitated the opening of the channel.

Overall, our detailed salt bridge analysis demonstrates that the two new salt bridge interactions (SB5 and SB6) that were formed in the 5-MTSET in the equilibrium MD analysis were preventing the complete opening of the channel. We believe that these new salt bridge interactions could potentially delay the channel’s opening but it is unclear whether they would play a role once the channel is in the open state. Future computational and experimental studies may shed light on the importance of these salt bridges.

Inter-helical angles confirm the partial opening of engineered MscL

Previous computational studies suggest that the MscL activation mechanism involves expansion of transmembrane complex considerably as it moves radially away from the pore axis, causing the tilt of helices to dramatically rise [82], [83]. The inter-helical angles would be a better indicator for identifying channel activity in response to the introduction of positively charged MTSET labels, therefore various inter-helical angles were measured, which were defined as alpha (), beta (), gamma (), and delta (). is the angle between two neighboring TM1 helices, is the angle between two neighboring TM2 helices, is the angle between TM1 and TM2 of neighboring protomers, and is the angle between TM1 and TM2 of the same protomer (schematically explained in Figs. 1 B and D). The overall behavior of the inter-helical angles (Fig. S4) of all five subunits (we refer to them as S1, S2, S3, S4, and S5) in WT were uniform. Note that we interchangeably use the terms subunits, protomers, and monomers throughout the manuscript. In the 5-MTSET, inter-helical angles show very distinct behavior compared to the WT. All inter-helical angles indicate that the 5-MTSET structure deviated significantly from the WT (Fig. S4 A–D), indicating that the helices of the 5-MTSET system have expanded far more than WT. However, we could not necessarily identify a consistent pattern among different subunits in terms of their inter-helical angles, probably due to the fact that 5-MTSET did not completely open in the equilibrium simulations. Our findings generally support the helical expansion hypothesis proposed in previous computational studies of MscL. [82], [83].

Extracellular loops play a key role in the dynamics of MscL

During channel activation, various computational investigations have shown that the periplasmic loops are disrupted as a result of helical expansion [44], [84], [85], [86]. To uncover the principal variations between the WT and engineered MscL structures, Principal components analysis (PCA) was performed. Projections onto Projections onto principal components (PC) 1 & 2 space clearly separated WT and 5-MTSET systems (Fig. 5A). The contribution of these two PCs to the total variance was 53.5% (Fig. 5C). All triplicate simulations of each system were clustered close to each other along PC1, elucidating the fact that the structural variations along PC1 were reproducible and significant. PCs 2 & 3 did not clearly separate the WT and 5-MTSET systems (Fig. 5 A, B), however, they display how dispersed the 5-MTSET systems were compared to the WT, which is due to the introduction of labels. This supports earlier observations from the inter-helical angle analysis that the inter-helical angles in the 5-MTSET system were dispersed without a coherent trend unlike in the WT.

Fig. 5

Principal components analysis (PCA) of equilibrium trajectories. (A) Projections along PCs’ 1 and 2. (B) Projections along PCs’ 1 and 3. WT and 5-MTSET structures are shown as filled triangles and circles, respectively. (C) The proportion of variance of the first 15 PCs. (D-F) Contribution of individual residues to square displacements in PCs 1 (D), 2 (E), and 3 (F), respectively. The color bar represents the extent of displacement. In all three PCs, maximal displacements were observed in the ECLs. The data shown in these three panels were projected onto the structures in panels G–I. (G-I) Structural variations in PCs 1, 2, and 3. Structures were colored based on the extent of their displacement, which was shown in panels D–F. Bidirectional arrows show the direction of the fluctuation and length of the arrow reflects the magnitude of fluctuation. It is visually evident that the three PCs were overwhelmingly dominated by the fluctuations in the ECLs and that they are playing a key role in discriminating WT and 5-MTSET.

To get a detailed understanding of the structural contribution to the variance in the PCs, the square displacement of each residue was estimated (Fig. 5D–F). The structural variation along PC1 was dominated by the ECLs of S3, S4, and S5 monomers (Fig. 5D), whereas ECLs of S1 and S5 monomers were dominating in PC2 (Fig. 5E), and in PC3 ECLs of monomers S2–S4 contribute the most (Fig. 5F). To demonstrate visually, the square displacements of residues in PCs 1–3 were projected onto the structure as shown in Fig. 5G–I. The square displacements demonstrate that the major structural contributors to the variance in the data were the ECLs. Overall, the key outcome of PC analysis was that the behavior of ECLs of WT and engineered MscL protein was significantly different. PCA supports our earlier hypothesis that the ECLs were playing a key role in the dynamics and functioning of the MscL channel.

Obtaining the open MscL state via non-equilibrium (NE) simulations

The unbiased equilibrium simulations discussed above provide some information on the spontaneous opening of the engineered MscL channel. However, the channel fails to open completely in our unbiased equilibrium simulations. This is due to the fact that the full opening requires various local and global conformational changes that are expected to occur on millisecond or longer timescales. Such timescales are currently beyond the reach of conventional all-atom MD simulations. One may use the all-atom MD to study slow conformational changes with the help of external force. Hence, NE pulling simulations were employed to induce the full opening. In these simulations, a time-dependent biasing potential is used to drive the system from an initial state to a final state by gradually changing the center of a harmonic biasing potential defined in terms of a collective variable. In this study, orientation collective variables (colvars) [78], [79] were used to steer the partially open WT and 5-MTSET structures (resulted from equilibrium simulations) towards the target (i.e., towards a fully open state) [87], [88], in a 100 ns simulation with a force constant (K) of 10,000 (Fig. 6A). A smFRET derived active/open state WT MscL was used as a reference/target structure for these NE simulations [80]. To verify the reproducibility, two sets of NE simulations were conducted for each system (Fig. 6A). To compare and contrast the WT and the 5-MTSET systems, work required to open the MscL structures in NE simulation was calculated. As expected, WT required more work to open compared to the 5-MTSET (Fig. 6A). This supports our earlier observation that labels facilitate opening of the channel. To verify the impact of starting orientation of MTSET labels, NE simulations were also employed to open the 5-MTSET(*) structure, and no significant difference in the work values was observed between the 5-MTSET and 5-MTSET(*) (Fig. 6A) indicating that there was no impact of starting orientation of the labels on the opening/activation of the channel. To evaluate the extent of opening of WT and engineered MscL channels in the NE simulations, RMSDs and water content across the pore were calculated (Figs. 6 B, C, and E). For clarity, only the first set of simulations were compared in Figs. 6B–H. Since a similar NE simulation protocol was employed across all MscL structures, irrespective of the presence/absence of MTSET labels, the extent of the opening was similar in the WT and engineered MscL structures (Fig. 6C), which was also reflected in the RMSDs (Fig. 6B) as well in the water content (Fig. 6E). The  RMSDs at the end of 100 ns NE simulations for the WT and 5-MTSET(*) were 10.5 Å (Fig. 6B) in both cases. Additionally, similar water profiles were observed for the WT and 5-MTSET(*) structures resulting from the NE simulations (Fig. 6E). However, more water molecules were observed in the pore’s bottleneck region (i.e., Z = to Å) in the WT (22) than in the 5-MTSET(*) (7), which we believe was due to the obstruction of the pore by MTSET labels in the 5-MTSET(*) (Fig. 6G).

Fig. 6

NE and FUE simulations result in an open MscL structure. (A) Comparison of NE work values. (B) Protein  RMSD as a function of simulation time. First 100 ns represent the NE pulling simulations and 100–350 ns represents the FUE simulations. Crystal structure (which was in a closed state) was used as a reference for these calculations. (E, F) Water count across the pore resulting from the 100 ns of NE (E) and 260 ns of FUE (F) simulations. The last 2.5 ns of the simulation trajectory was used for the water content calculations. (C, G) Superposition of WT and 5-MTSET(*) structures (we call them ’open(*)’) resulting from NE simulations. WT and 5-MTSET(*) were represented in blue and red, respectively. Side and top views were shown in C and G, respectively. (D,H) Structures (we call them ’Open’) of WT (D) and 5-MTSET(*) (H) resulting from the FUE simulations. Waters within the 3 Å of the protein were shown.

Among mechanistically relevant features of MscL, previously suggested based on both computational and experimental studies, is the presence of a conserved N-terminal helix, which is somewhat perpendicular to the TM1 helix in the inactive state (Fig. 7). During the tension-induced activation process, the pore-lining TM1 helix is shown to align with the N-terminal (NT) helix, resulting in the formation of a straight helix (i.e., straightening of the TM1 helix) [89]. We quantified the angle between the TM1 helix and the NT helix attached to it to determine whether or not the TM1 helix undergoes straightening (Fig. 7). The TM1-NT angle fluctuates roughly between 60 and 90 degrees in equilibrium simulations of the 5-MTSET and WT (Fig. 7 A). The TM1-NT angle behaves somewhat similarly in both systems in equilibrium simulations, not supporting the TM1 straightening hypothesis but not ruling it out either. However, in NE simulations, the TM1-NT angle changes significantly (Fig. 7 B) in both the wild-type system, supporting the TM1 straightening hypothesis. On the other had, we do not observe a clear straightening in the 5-MTSET simulations, even though we employ the same NE strategy for channel opening in both cases. Note that the NT helix is not involved in the biases used in either NE simulation; however, the TM1-NT angle reduces to around 20 degrees in the wild-type simulations and not in the 5-MTSET simulations. The differential behavior of the TM1-NT angle between WT and 5-MTSET is noteworthy because it supports our hypothesis that in the presence of MTSET, MscL is activated in a distinct way.

Fig. 7

Intra-subunit helical angles of N-terminal TM1 and pore-lining TM1 helix. The intra-subunit TM1 Helix angle estimated in this analysis are shown graphically in the left panel. (A-B) N-terminal TM1 and pore-lining TM 1 helix intra-subunit helical angles. Panel A shows the equilibrium simulations results, and Panel B shows NE simulations of the WT and 5-MTSET(*).

The NE pulling simulations used here provide us with yet another way of comparing the activation mechanism in the WT and engineered MscL. While the disadvantage of the unbiased simulation was related to the lack of a complete picture of activation mechanism due to the short simulation times, the main disadvantage of NE pulling method is that it does not necessarily generate a minimum free energy path. In other words, it is generally possible that the pathway generated by the NE method is not the same as that generated in a much longer unbiased simulation. We therefore do not make a claim that the pathways observed for the WT and engineered MscL represent their most relevant transition pathways. Instead, we use this method for a qualitative comparison between the two systems and their pathways based on the fact that both systems use an identical biasing protocol. In addition, since the resulting open states from these NE protocols are out of equilibrium, we have attempted to equilibrate these structures in the next step as discussed below.

Characterizing the open MscL structures resulting from the follow-up equilibrium simulations

We further characterized the open WT(1) and 5-MTSET(*) MscL structure resulting from the follow-up NE simulations. We calculated several variables including RMSF, ECL , number of BB-BB H-Bonds, Average Y94-D16 H-Bond frequency, SC-SC (both inter-/intra-unit) salt bridge interactions and various interhelical angles to compare and contrast between the WT(1) and 5-MTSET(*) open structures. The WT(1) and 5-MTSET(*) MscL structures resulting from the FUE simulations were further characterized. To compare the WT(1) and 5-MTSET(*) open structures, variables such as RMSF, ECL R, number of BB-BB H-Bonds, average Y94-D16 H-Bond frequency, SC-SC (both inter-/intra-unit) salt bridge interactions, and inter-helical angles were analyzed. Significant differences between the open WT(1) and open 5-MTSET(*) structures were observed. The overall fluctuation of the 5-MTSET(*) structure was greater than the WT(1) (Fig. 8A), which probably was due to the presence of labels. The R of ECL of both structures collapsed for the entire length of the trajectory; there was at least 2 Å greater reduction in the WT(1) than 5-MTSET(*) (Fig. 8b) indicating that the extracellular region of 5-MTSET(*) was more open than the WT(1). The number of BB-BB H-Bonds in WT(1) was greater (on average 15) than in the 5-MTSET(*) (Fig. 8C). The formation of these extra H-Bonds facilitated a greater collapse of WT(1) in the FUE simulations than the 5-MTSET(*). This demonstrates that the loss of BB-BB H-Bonds makes the structure more flexible and is essential for its opening/activation. We have also observed a loss of Y94-D16 H-Bond in 5-MTSET(*) (15% vs. 45% H-Bond frequency in 5-MTSET(*) and WT(1), respectively; Fig. 8D). This supports the hypothesis that the loss of this particular H-Bond is key for the activation of engineered MscL.

Fig. 8

Characterizing the open MscL structures resulting from the follow-up equilibrium (FUE) simulations. (A-D) Comparison of RMSF (A), radius of gyration (R) (B), number of BB-BB H-Bonds (C), and average Y94-D16 H-Bond frequency (D) between WT(1) and 5-MTSET(*) systems. WT(1) and 5-MTSET(*) data were represented in blue and red, respectively. The entire 260 ns of each FUE simulation data was considered for H-Bond calculations.

Additionally, we performed salt bridge interaction analysis; SC-SC (both inter-/intra-unit) salt bridges were estimated in both structures resulting from FUE simulations. The total number of inter-unit salt bridge s in the case of 5-MTSET(*) and WT(1) were 4 and 14, respectively. We believe that these extra ten inter-unit salt bridge s in the WT(1) FUE structure, apart from the BB-BB H-Bonds explained above, facilitated its greater collapse compared to the 5-MTSET(*) during FUE simulations. The presence of MTSET labels in 5-MTSET(*) probably resisted the formation of these extra ten inter-unit salt bridge s. Note that only six inter-unit salt bridge s were identified in the WT(1) structure resulting from the earlier equilibrium simulations, whereas seven were identified in the 5-MTSET structure. Therefore, there was a gain of eight new inter-unit salt bridge s in the WT(1) FUE structure, whereas in 5-MTSET(*) there was a loss of three inter-unit salt bridge s, compared to their respective earlier equilibration simulation structures. The 14 inter-unit salt bridge s that were identified in the WT(1) were further classified into seven different classes for the sake of analysis(3). Among the seven classes, only two were identified in the previous equilibrium MD simulations (i.e., TM1(i)-TM2(i + 2) and TM1(i)-ECL(i + 1); Table S1). Therefore, there was a gain of five new classes of salt bridge s in the WT(1) FUE structure. The open WT(1) structure resulting from the FUE simulations have to loose these new five classes of salt bridge s to collapse completely back to its starting structure (i.e., the crystal conformation). The two classes of salt bridge s that were observed in 5-MTSET in the previous equilibrium simulations, i.e., TM1(i)-TM2(i + 1) and TM1(i)-ECL(i + 1), were completely lost in the FUE simulations. This supports our earlier hypothesis that a loss of TM1(i)-ECL(i + 1) salt bridge s are key for the complete opening of the engineered MscL channel. However, there was a gain of three new classes of salt bridge s in 5-MTSET(*) in FUE: TM1(i)-TM1(i + 1), TM2(i)-TM2(i + 1) and IH(i)-IH(i + 1) ((3)). Additionally, there was a gain of two new intra-unit salt bridge s in 5-MTSET(*) (Table 4), which belong to the class ECL(i)-ECL(i). This is another new class that was identified in the FUE simulations. These two salt bridge s were not observed in the WT(1). Overall, our BB-BB H-Bond and salt bridge analysis demonstrate that the loss of these two categories of interactions is key for the activation/opening of the MscL channel.

Table 3

Inter-unit salt bridge interactions in FUE simulations.

Salt bride Class	Salt bride	WT(1)	5-MTSET(*)
TM1(i)-TM2(i + 2)	R11-D108(SB3)	2	1
	E7-K99(SB4)	1	0
TM1(i)-ECL(i + 1)	R45-D53(SB7)	1	0
TM1(i)-TM1(i + 1)	R11-D16(SB9)	2	1
TM1(i)-TM2(i + 1)	K6-E102(SB10)	3	0
TM2(i)-TM2(i + 1)	R98-E104(SB11)	2	1
ECL(i)-ECL(i + 1)	D53-R58(SB12)	1	0
IH(i)-IH(i + 1)	E116-R118(SB13)	2	1
Total		14	4

Last 100 ns of the FUE simulation trajectory was considered. Interactions with 70% were considered. SB is salt bridge. i refer to the first protomer. TM, IH, and ECL refer to transmembrane helices, intracellular helices, and extracellular loops, respectively. Inter-unit interactions refer to interactions between different protomers.

Table 4

Intra-unit salt bridge interactions in FUE simulations.

Salt bride Class	Salt bride	WT(1)	5-MTSET(*)
ECL(i)-ECL(i)	D68-R58(SB14)	0	2
Total		0	2

Last 100 ns of relaxation simulation trajectory was considered. Intra-unit interactions refer to interactions within a protomer.

Proposed mechanism of action

We propose a hypothesis for the opening/activation of the engineered MscL channel based on our extensive analysis. According to this hypothesis, first the channel opens spontaneously due to the repulsion between engineered positively charged labels and also due to steric clashes between the labels and other bottleneck residues, which was also demonstrated through interaction analysis (eight among the top ten residues interacting with the labels were non-polar (Table S3). This sudden jerk-like motion near the bottleneck region leads to the breaking of inter-unit salt bridge s SB4 (TM1(i)-TM2(i + 2); Table S1) and SB7 (TM1(i)-ECL(i + 1); Table S1), and an inter-unit H-Bond (Y94-D16; Fig. 8D). Instead other salt bridges form such as inter-unit salt bridge s SB5 and SB6 (TM1(i)-ECL(i + 1); Table S1) that are absent in the WT MscL but form in the engineered MscL. These salt bridges both contribute to the distortion of the extracellular loops that is missing in the WT MscL. All these changes in the interaction patterns of MscL eventually facilitate transitioning of the structure from a closed state to an open state. The loss of BB-BB H-Bonds (25–30) enables this transition by loosening up the TM helices (Fig. 3 A–C and Fig. 8C). Finally, the gain of inter-unit (SB9, SB11, and SB13; Table 3) and intra-unit salt bridge s (SB14; Table 4) locks the channel in an open state. Furthermore, we propose that the non-polar amino acids L17 and A18 that reside between the labels on the intracellular side were reducing the repulsion between the positively charged labels and were restricting the extent of the opening of the channel.

Conclusions

In this study, the activation mechanism of an engineered MscL channel was investigated at an atomistic level using a combination of all-atom equilibrium and non-equilibrium (NE) MD simulations. In the equilibrium MD simulations, the attached MTSET labels facilitated the opening of the engineered MscL channel to some extent, however, the salt bridge interactions between the extracellular loops and the TM helices prevented its complete opening. Therefore, the NE simulations method were employed to completely open the channel. The major finding of this study is that the open/active conformation of the engineered MscL channel captured via this strategy was not as open as proposed earlier by several other studies. Additionally, we have identified that periplasmic loops were playing a key role in the activation of the channel and that the net loss of several backbone-backbone hydrogen bonds and side chain-side chain salt bridge interactions facilitated the opening of the channel. In this study, we successfully demonstrated that the introduction of labels facilitates the opening of the channel as expected and that the adopted computational protocol is a feasible strategy to elucidate the activation mechanism of ion channels.

CRediT authorship contribution statement

Kalyan Immadisetty: Investigation, Visualization, Formal analysis, Writing - original draft. Adithya Polasa: Visualization, Formal analysis, Writing - review & editing. Reid Shelton: Investigation, Validation. Mahmoud Moradi: Conceptualization, Methodology, Writing - review & editing, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.