Nanomechanical properties of MscL α helices: A steered molecular dynamics study

ABSTRACT

Gating of mechanosensitive (MS) channels is driven by a hierarchical cascade of movements and deformations of transmembrane helices in response to bilayer tension. Determining the intrinsic mechanical properties of the individual transmembrane helices is therefore central to understanding the intricacies of the gating mechanism of MS channels. We used a constant-force steered molecular dynamics (SMD) approach to perform unidirectional pulling tests on all the helices of MscL in M. tuberculosis and E. coli homologs. Using this method, we could overcome the issues encountered with the commonly used constant-velocity SMD simulations, such as low mechanical stability of the helix during stretching and high dependency of the elastic properties on the pulling rate. We estimated Young's moduli of the α-helices of MscL to vary between 0.2 and 12.5 GPa with TM2 helix being the stiffest. We also studied the effect of water on the properties of the pore-lining TM1 helix. In the absence of water, this helix exhibited a much stiffer response. By monitoring the number of hydrogen bonds, it appears that water acts like a ‘lubricant’ (softener) during TM1 helix elongation. These data shed light on another physical aspect underlying hydrophobic gating of MS channels, in particular MscL.

KEYWORDS:

• all-atom simulation
• constant velocity
• Escherichia coli
• mechanosensitive channel
• mycobacterium tuberculosis
• nanovalve
• Young's modulus

This article is referred to by:

Channel disassembled: Pick, tweak, and soak parts to soften

Introduction

The mechanosensitive channel of large conductance, MscL acts as a safety valve in bacterial membranes allowing the bacteria to release the turgor pressure under hypoosmotic conditions.^1-3 For this to occur, tension is transmitted directly from the lipid bilayer to MscL resulting in a conformational change that leads to channel gating.^4,5 To date, the structure and function of MscL are well characterized using a plethora of experimental and computational approaches.^1,6-13 MscL is a homopentamer with each monomer consisting of 4 α-helices, 2 transmembrane (TM1 and TM2) and 2 cytoplasmic (N- and C-terminal) helices (Fig. 1).¹⁴ The TM1 helix, which lines the pore, is coupled to the membrane via a juxtaposed horizontal amphipathic N-terminal helix residing at the lipid-solvent interface.^15,16 The TM2 helix faces the bilayer and spans the full bilayer thickness (Fig. 2B) and is connected to a coiled-coil C-terminal helical bundle of the channel pentamer.^14,17 Due to in-plane expansion and thinning of the bilayer, the pore-lining TM1 helix is dragged by the N-terminus and tilts and rotates, resulting in solvation of a hydrophobic gate (Fig. 1C).^6,9,18-20 Therefore, both the TM1 helix and the N-terminus undergo a considerable axial force in the open state.^16,21 MscL activation leads to opening of a large non-selective pore with a diameter of ∼ 3 nm and a unitary conductance in the range of ∼ 3 nS.^6,11,20

Figure 1. Three dimensional structures of the closed (resting) state of (A) MtMscL (PDB code: 2OAR) on the left and a homology model of EcMscL obtained based on 2OAR and 4LKU on the right side.^14,16,17 (B) A subunit of EcMscL has been shown after equilibrated in POPE lipid bilayer. The residues that anchor the protein to the membrane (F7, F10) have been highlighted.¹⁵ (C) TM1 helix becomes aligned with the N-terminus in the open state.¹⁶

Figure 2. Steered molecular dynamics (SMD) simulation using constant-velocity (CV) method. In this method, the α carbons of the first 2 residues (ILE14 and VAL15) in TM1 helix of MtMscL have been fixed on the left hand side. Helical properties such as helix initial length, L₀, pulling rate, v, time t, vector position of the helix end, ${\displaystyle \overrightarrow{r}}$, and the assigned spring constant of the dummy atom, ${\displaystyle k_{dummy}}$, have been indicated. (A) A schematic unidirectional pulling test on the TM1 helix of MtMscL solvated in water. (B) The force versus elongation (ΔL) curve during constant-velocity simulations. The typical spring constant here is 3 kcal/mol/Ų (i.e. ∼210 pN/Å). Four different constant velocities have been assigned to the α carbon of the dummy atom, from 0.1 Å/ps, to 5 Å/ps. The helix elongation trend can be divided into 2 main quasi-linear regions. As shown in the inset multiple regions can be observed in the first regions which are indicative of sequential rupture of the hydrogen bonds (i.e., breakage of the hydrogen bonds during the pulling test). The length of TM1 increases up to a maximum length (>100 % of elongation), when the helix becomes almost unfolded due to the applied stretch. Before this point, helix response is rate dependent, i.e. the higher the pulling rate, the stiffer is the helix response. After this point the helix behavior becomes much stiffer and not rate dependent. It should be noted that the rationale for testing the helix behavior upon application of such high forces and large elongations was to find the reason for the rate dependency observed in the CV simulations. (C) A typical helical behavior under low forces. The spring constant here was 0.6 kcal/mol/Ų (i.e., ∼42 pN/Å) and the pulling rate was 0.1 Å/ps. Similar to previous studies,³⁹ the slope of the best linear fit to the initial part of the diagram (dashed red line) indicates the elasticity modulus, which is ∼0.7 kcal/mol/Ų (i.e., ∼45 pN/Å).

Gating in MscL occurs on a timescale in the order of milliseconds and is accompanied by series of movements resulting in large helical deformations that are dependent on the mechanical properties of each individual helix. In order to understand the intricacies of the gating mechanisms of these MS channels it is important that we fully determine the intrinsic mechanical properties of the individual TM helices. This information will also provide useful inputs and constraints for future computational analyses including coarse-grain molecular dynamics (MD) and continuum mechanics simulations that are routinely employed for investigating the global structural changes during MS channel activation.^22-29 Furthermore, given the growing interest in using MscL as nanosensor and/or nanovalve in biotechnology,^30-37 there are still unanswered questions regarding the intra-molecular mechanics of this family of MS channels. For example, how much force does each α-helix bear during the gating process? How do the helix properties differ in different environments (e.g., in water)? Also, is there any difference between mechanical properties of the helices in different MscL homologs?

The mechanical properties of different proteins have been studied using various theoretical^38-40 and experimental techniques such as steered molecular dynamics (SMD), atomic force microscopy (AFM), and laser optical tweezers.^38,39,41-45 These methods, in combination with X-ray crystallography and cryo-EM structures, have greatly advanced our knowledge of protein structure, mechanical strength and function. Among these approaches, SMD provides atomistic descriptions of the mechanical behavior of different proteins.^46-48 This can be achieved using 2 different approaches, either by applying a constant-force (CF) or a constant-velocity (CV). The advantages and disadvantages of these approaches are however not entirely clear.

In this study, we compared both CV and CF SMD for performing unidirectional pulling tests on the TM helices of MscL from M. tuberculosis and E. coli homologs. Using large number of simulations and the TM1 helix of MscL as our model α-helix, we demonstrated that the CF method is more reliable and accurate compared to the CV method in determining the mechanical properties of α-helices but more expensive in terms of the computational costs. We also examined the effect of water on the elastic behavior of the pore lining TM1 helix of MscL. We showed that in its hydrated state, the TM1 helix of MtMscL became more than 5 times more flexible, while the elasticity of the TM1 helix in EcMscL did not vary significantly. The framework and mechanical properties estimated here have important consequences for ‘hydrophobic lock-dependent gating’^49-54 of ion channels and may apply to many other prokaryotic and eukaryotic ion channels and proteins in general.

Materials and methods

All simulations were performed with the NAMD 2-10 package, where CHARM36 Force Field was employed. Visual Molecular Dynamics (VMD)⁵⁵ and Pymol were used for all visualizations. The crystal structure of the M. tuberculosis MscL (2OAR) was used from the Protein Data Bank (Fig. 1A). The 3D structure of EcMscL (Fig. 1B) was generated based on the crystal structure of the MscL homolog of M. tuberculosis (PDB ID: 2OAR) and crystal structure of E. coli C-terminus (PDB ID: 4LKU)¹⁷ using Phyre2⁵⁶ and Swiss-Model.⁵⁷

Steered molecular dynamics (SMD) simulation

In SMD, external forces/velocities are applied to biomolecules to manipulate them in order to probe and determine their different mechanical properties.^46,58,59 The application of external forces/velocities accelerates processes that are otherwise too slow to model by using non-steered molecular simulations.^6,28 We used 2 different types of steered molecular dynamics simulation, namely constant-velocity (CV) and constant-force (CF). The CV method is also called “Moving Constraints” method in some versions of NAMD and other studies.⁶⁰ In this method, we have restrained the 2 α carbons of the first 2 residues (e.g., ILE14 and VAL15 in the case of TM1 helix of MtMscL) on each helix using a strong harmonic restrain constant of 12 kcal/mol/Ų. Table S1 and S2 in the Supporting Material contain more details about all the helices studied here. The last α carbon at the end of the helix, which is called the SMD atom, is attached to another ‘dummy’ atom by a harmonic spring (e.g., ILE46 in the TM1). The dummy atom is then moved with a given constant-velocity vector (Fig. 2A). We used a wide range of velocities for the TM1 helix (0.1 to 5 Å/ps) to fully investigate the effect of this parameter on estimation of the elastic moduli of our helices. In the CV method the force vector, ${\displaystyle \overrightarrow{F}}$, that has been applied on the SMD atoms can be calculated using Eq. 1(1) ${\displaystyle \overrightarrow{F}=-\nabla \left(\frac{1}{2}k{[vt-(\overrightarrow{r}-{\overrightarrow{r}}_0)\cdot \overrightarrow{n}]}^2\right)}$(1) ,³⁹(1) ${\displaystyle \overrightarrow{F}=-\nabla \left(\frac{1}{2}k{[vt-(\overrightarrow{r}-{\overrightarrow{r}}_0)\cdot \overrightarrow{n}]}^2\right)}$(1)

Where k is the stiffness of the system, ${\displaystyle v}$ is the pulling velocity, t is time, ${\displaystyle \overrightarrow{r}}$ is the position of the helix end (i.e., ${\displaystyle |\overrightarrow{r}-{\overrightarrow{r}}_0|}$ is the displacement) and ${\displaystyle \overrightarrow{n}}$ is the vector that indicates the pulling direction. It is also noteworthy, that the elongation (or strain) presented in the figures of this study, are the helical elongation which is the distance between the fixed atom and the SMD atom. Therefore, the resulting k from the slope of force versus elongation diagrams only indicates the elasticity constant of the helix. To estimate the Young's modulus, E, from the spring stiffness, k, we used Eq. 2(2) ${\displaystyle E=k \times L_0/A}$(2) as follows.(2) ${\displaystyle E=k \times L_0/A}$(2)

Where ${\displaystyle L_0}$ is the initial length of helix (Table S1 and S2), and A is the average cross-sectional area of helix. For calculating the helical cross-sectional area, we needed to estimate a radius of gyration for each helix as an average radius to be able to calculate their cross-sectional area using Eqs 3 and 4. The radius of gyration here, ${\displaystyle r_{gyr}}$, is the root mean square distance of the atoms from the helix central axis as following.(3) ${\displaystyle r_{gyr}=\sqrt{\frac{\sum {(d-d_{\mathrm{mean}})}^2}{N}}}$(3) (4) ${\displaystyle d=\sqrt{x^2+y^2}\text{\hspace{0.17em}and }{d}_{mean}=\sqrt{x_{\mathrm{mean}}{}^2+y_{\mathrm{mean}}{}^2}}$(4)

Where d and ${\displaystyle d_{\mathrm{mean}}}$ are the in-plane coordinates of each atom and central axis, respectively. N is the total number of atoms. x and y are the coordinates of each atom on the helix except hydrogens when the helix is aligned to the z axis. ${\displaystyle x_{\mathrm{mean}}}$ and ${\displaystyle y_{\mathrm{mean}}}$ are the in-plane coordination of the central axis of the helix when the helix is aligned to the z axis. For example, in the case of the TM1 helix, Eq. 3(3) ${\displaystyle r_{gyr}=\sqrt{\frac{\sum {(d-d_{\mathrm{mean}})}^2}{N}}}$(3) yields the helical radius of 2.98 Å and thus annular cross-sectional area of 27.9 Å², which is similar to previously estimated values for the MscL helices.^11,61

In the CF method, a constant-force vector is applied onto selected atoms on one end of the helix, while the other end has been restrained. In this method, the α carbons of the first 2 residues (e.g., ILE14 and VAL15 in TM1 helix of MtMscL) were restrained, and the force was applied on the α carbon of the last residue (e.g., ILE46). See Table S1 and S2 for the information about the fixed residues and dummy atoms on each helix. In the case of TM1 for example, a constant-force in the range of 0.1-1.5 kcal/mol/Å (∼7 to 104 pN) was applied on the α carbon atoms of the last residue in the direction defined by a vector. This vector links center of the mass of the fixed atoms to the SMD atom, i.e., it is along the helix axis. For the other helices, the force range may alter given that the applied force was gradually increased to reach the helical strain of ∼6 %. We chose to measure the elastic response of each α-helix in this range as it is the maximum range that the MscL helices undergo during gating. This is a common range that has previously been used for other proteins as well.^39,41 After application of each force, the system was allowed to run for 10 ns to stabilize the equilibrium length (Fig. S1). The resultant strain of each helix was calculated by dividing the length of elongation, ${\displaystyle \Delta L}$, by the initial length of the helix, ${\displaystyle L_0}$(i.e., ${\displaystyle \Delta L/L_0 \times 100}$). See Table S1 and S2 for the initial lengths of all the helices. The initial 0% strain corresponds to the non-stretched state. For each applied force, the system was permitted to equilibrate for 10 ns (Fig. S1). For calculating the stress resultant, we calculated force/cross - sectional area of the helix (${\displaystyle \left|\overrightarrow{F}/\overrightarrow{A}\right|}$), where A is calculated using Eq. 3(3) ${\displaystyle r_{gyr}=\sqrt{\frac{\sum {(d-d_{\mathrm{mean}})}^2}{N}}}$(3) . The slope of the stress-strain diagram indicates the Young's modulus of the helix. It should be noted that each point on each stress-strain diagram (e.g., Fig. 3C) is the average of 3 separate simulations. Therefore, with the number of forces we used to create reliable stress-strain curves for each helix, we performed an average of 51 simulations for each individual helix.

Figure 3. Mechanical behavior of the TM1 helix of MtMscL using constant-force (CF) method (A) A schematic unidirectional pulling of the TM1 helix of MtMscL solvated in water. (B) TM1 helix elongates over time as a result of constant pulling force applied on its end. The TM1 length increases to maximum length of ∼52 Å then fluctuates around this value. (C) Stress-strain curve of the unidirectional traction applied to the TM1 helix in water. The range of axial force applied in these simulations ranges from 0.1 to 1.5 kcal/mol/Å (i.e., from ∼7 to 70 pN). The Young's modulus of the TM1 helix can be calculated from the slope of this curve, which in this case is ∼3.2 ± 0.9 GPa (Mean ± SEM). To obtain each point on the stress-strain curve, 3 simulations were performed. The strain at each force has been averaged over 3 simulations. The Young's modulus has been estimated based on 95 % confidence of both stress and strain axis.

In all of our simulations, a modified Nosé-Hoover Langevin piston pressure control provided in NAMD was utilized to control fluctuations in barostats (at 1 atm). This method is combined with a method of temperature control (at 300 K) (Langevin dynamics) to simulate the NPT ensemble. In order to study the role of water on the mechanical properties of the pore-lining α helix TM1, a series of SMD simulations were carried out in the NVT ensembles for 2 cases: (i) with water and (ii) without water (i.e. in vacuum). The TM1 helix was stretched along its axis in water and in vacuum using CF method. In the cases with water, the α helix was solvated in a 40 × 40 × 90–Å water box using the “SOLVATE” module, with the TIP3P water molecule, in visual molecular dynamics (VMD) software. The helix was held fixed for the first 1 ns by restraining the atomic positions and then, only its backbone was maintained fixed while the rest of the system was allowed to relax for the following 1 ns. Thereafter, the whole system was equilibrated in a NPT ensemble for 10 ns with no restraint. An example of RMSD values during equilibration for the TM1 helix of MtMscL is shown in Figure S1. Care was taken of the charge neutrality in the simulations with water. We used HeliQuest and VMD for our bioinformatics data presented in Table S1 and S2.

Results

The mechanical response of the pore lining TM1 helix

We first compared the 2 SMD methods, constant-velocity (CV) and constant-force (CF), in determining Young's elasticity modulus of the TM1 helix of MtMscL solvated in water. In the CV method, the helix elongates as a result of assigning a velocity to the dummy atom attached to the end of the helix by a virtual spring, while the other end is fixed (Fig. 2A, Movie S1). Due to the helix elongation, force is generated in the helix which can be calculated from Eq. 1(1) ${\displaystyle \overrightarrow{F}=-\nabla \left(\frac{1}{2}k{[vt-(\overrightarrow{r}-{\overrightarrow{r}}_0)\cdot \overrightarrow{n}]}^2\right)}$(1) . The force vs. elongation was monitored over the simulation time for 4 different pulling velocities in the range 0.1 to 5 Å/ps (Fig. 2B). Two main semi-linear regions are evident (Fig. 2B). In the first, the slope of force versus elongation, representing the stiffness of the helix, is shallower and rate dependent. This region relates to hydrogen bonds and non-bonded interactions in general, which determine the secondary structure of the helix. In the second region the slope of force vs. elongation is steeper, corresponding to greater stiffness than the first region, and is independent of the pulling velocity. This region corresponds to the scenario after all the hydrogen bonds have been broken, so represents the elasticity or stiffness of the helical backbone, which is mainly determined by bonded interactions. Importantly, the rationale for examining the helix behavior upon application of such high forces and large elongations was merely to understand the rate dependency observed in the CV simulations, which also occurs even at very low forces (e.g. <10 pN). It should be noted that the elasticity measurements (elastic constants and Young's moduli) and comparison between the CV and CF methods were done at forces <100 pN (∼6% strain for TM1 helix) as described in the methods.

We have to note that the first region is not perfectly linear. This can be more clearly seen from the inset in Fig. 2B. As shown, first there was a steep jump in the force value with little to no change in the helix length, but then the helix length increased in multiple steps as the force increased. Each of these semi-linear domains is shown in the force versus elongation graph (Fig. 2B), which indicates the strain dependent behavior of TM1 helix in response to uniaxial force when the CV method was used. Since there was a slight nonlinearity in the helical response and to be able to compare our results with the results obtained from the CF method, we measured the helical properties corresponding to strain < 6 % region (Fig. 2C). This is done similar to what has been done in previous studies,³⁹ by measuring the slope of the best linear fit to the portion of force-elongation diagram which is in the strain range of <6% (e.g., for TM1 helix, this strain value corresponds to 2.9 Å elongation). This is the range we chose to measure the elastic response of each α-helix throughout this study unless otherwise specified.

The rationale for doing it is that this is within the range of the strain that the α-helices of EcMscL experience during gating.^6,11,16 Also, this range has previously been used for other proteins in previous studies.^39,41 As described in the methods, for finding the apparent Young's modulus from the spring constant, we had to define the cross-section for the helix. For example, in the case of the TM1 helix, Eq. 3(3) ${\displaystyle r_{gyr}=\sqrt{\frac{\sum {(d-d_{\mathrm{mean}})}^2}{N}}}$(3) yields the radius of 2.98 Å and thus annular cross-sectional area of 27.9 Å², which is similar to previously estimated values for the MscL helices, ∼2.5 Å.^11,61 There is however, uncertainty in estimating the Young's elastic modulus of the α-helices from their spring constant (k), since they are not perfectly isotropic structures. One way to overcome this problem is to divide each helix into multiple domains with different geometries. But this would require many more computations, a more complex theory and a detailed geometric model of the helix. Hence, we consider an α-helix to behave as a homogenous material so that the apparent Young's modulus, E, can be expressed as a simple scalar from Eq. 2(2) ${\displaystyle E=k \times L_0/A}$(2) . Although there are limitations in applying such a simplified model, it provides insight and describes well the overall mechanical behavior of different α-helices.⁶¹ For future work, the helix cross-sectional area could be modeled as its van der Waals (vdW) shape, allowing stress values within the helix to be more precisely calculated based on force over local cross-section, instead of force over a uniform cross-section. This issue may be important if the helix dynamics are associated with bending (flexural) or torsional stiffness—where calculation of the actual cross-sectional area is more important. These more realistic surfaces will also better reflect the “physicochemical” intricacies of MscL in future continuum models, such as when helical charge density has to be taken into account.⁶²

In the constant-force (CF) method, an external constant-force has been applied on one end of the TM1 helix, while the other end was restrained (See Table S1 for the fixed and stressed residues for each individual helix). As a result of axial pulling force, TM1 helix elongates to a stable length very quickly (e.g., in 0.5 ns of simulation). However, we continued the simulation for at least another 10 ns to equilibrate helical length at each force value (Fig. 3A, B and Fig. S1). This procedure is repeated for the forces up to ∼100 pN which is enough to stretch the helix beyond the length observed in the open state structure of MscL.²⁸ The stress-strain curve describing the helix response to a range of pulling forces is semi–linear (Fig. 3C). The slope of this curve indicates the Young's modulus of TM1 solvated in water, which is estimated to be 3.2 ± 0.9 GPa.

Comparison between constant-velocity (CV) and constant-force (CF) method

Given the differences observed in TM1 helix behavior in response to the axial force applied using the CV and CF methods (Figs 2B and 3C) we further investigated the role of the pulling rate and spring constant of the dummy atoms in our CV simulations in low strain regions (<6%). We performed more than 200 SMD simulations using the CV method with pulling rates between 0.1 Å/ps to 5 Å/ps, as well as a wide range of dummy atom spring constant values from 0.1 kcal/mol/Ų to 2 kcal/mol/Ų. Figure 4A shows that the helix Young's modulus is dependent on both the pulling rate as well as the spring constant assigned to the dummy atoms. In general, for faster pulling rates, the Young's modulus was larger, i.e., the helix was stiffer. Specifically, as the rate of pulling decreased from 5.0 Å/ps to 0.1 Å/ps, the average of E values decreased from ∼42 GPa to ∼11 GPa. For some pulling rates (e.g., 5.0 Å/ps), higher spring constants of the dummy atoms resulted in a stiffer helix response. Also at higher pulling rates (i.e., 5 Å/ps), there was a considerable fluctuation in the Young's modulus of the helix which became larger as the spring constant of the dummy atoms increased. When we decreased the rate of pulling to 5 × 10⁻⁴ Å/ps the resulting spring stiffness of the TM1 reduced to ∼25 pN/Å which corresponds to the Young's modulus ∼4.3 GPa (Eq. 2(2) ${\displaystyle E=k \times L_0/A}$(2) ; Fig. S2). This rate is closer to but still faster than the actual extension rate of this helix during the MscL gating, which is <5 × 10⁻⁷ Å/ps, given that MscL gating occurs in several μs and TM1 helix elongates ∼5 Å during gating.¹⁶

Figure 4. Comparing the constant-velocity (CV) method with the constant-force (CF) method for testing the mechanical behavior of α-helices using TM1 helix of MtMscL in water as an example. (A) The effect of pulling rate and stiffness of spring constant of the dummy atoms using constant-velocity method on the mechanical behavior of the TM1 helix of MtMscL. For higher pulling rates (> 1 Å/ps), there is a considerable fluctuation in the Young's modulus depending on the spring constant assigned to the dummy atoms. This fluctuation becomes larger as the spring constant becomes higher. As the rate of pulling decreases from 5.0 Å/ps to 0.1 Å/ps, the average of E values decreases from 42 to 11 GPa. (B) Change in the number of hydrogen bonds of TM1 α helix of MtMscL during the simulation time. The black trace shows the number of hydrogen bonds in the TM1 helix in the CV method, compared to those when the CF method was used. For this typical CV example, the spring constant is 0.6 kcal/mol/Ų (i.e. ∼42 pN/Å) and the pulling velocity is 0.1 Å/ps. The force used in the CF method was 27 pN. The values of Young's moduli are Mean ± SEM. One-way ANOVA was used for statistical analysis with p-value < 0.05 and it was confirmed in groups with low n number using non-parametric Kruskal-Wallis test.

We next compared the number of hydrogen bonds in the TM helix after 10% strain in CV and CF simulations (Fig. 4B). As indicated in this figure, the number of hydrogen bonds in the TM1 helix only marginally decreased for the strain <10% compared to the CF method.

Comparing the mechanical properties of the α-helices in the MtMscL and EcMscL

A comparison of the mechanical properties of transmembrane helices between homologues of MscL from M. tuberculosis and E. coli is shown in Fig. 5. The TM1 helix and the N-terminus of MtMscL have similar Young's moduli to those in EcMscL, while the TM2 helix of MtMscL is about 4 times stiffer than in EcMscL. Conversely, the C-terminal helix of EcMscL is about 8 times stiffer than the C-terminal helix of MtMscL. There is no apparent correlation between the length and stiffness of the helix in either homolog (Fig. S3). There is also no correlation between the Young's modulus and hydrophobic moment or number of charges (Fig. S4 and S5 and Table S1 and S2). We also showed that the Young's moduli of the TM2 helix in MtMscL and EcMscL are vastly different, although they have very similar number of hydrogen bonds (Fig. S6).

Figure 5. Mechanical properties of MtMscL and EcMscL α-helices in the presence of water measured by all-atom steered molecular dynamic (SMD) simulation. We used constant-force (CF) method here to estimate (A) the elasticity constant and (B) the Young's modulus. The elasticity constant is between 1.0 to 72.2 pN/Å (Young's moduli are between 0.2 to 12.5 GPa). Overall, the mechanical properties of MscL α-helices of both species are similar to each other except for their TM2 and C-terminal helices. TM2 helix in MtMscL is almost 4 times stiffer than the TM2 helix of EcMscL. But the C-terminal helix of EcMscL is about 8 times stiffer than the C-terminal helix of MtMscL. The values of Young's moduli are Mean ± SEM for n = 3. Each Young's modulus is obtained from the stress-strain graphs exemplified in Fig. 3. Student's t-test was used for statistical analysis between each α-helix in MtMscL with its corresponding α-helix in EcMscL. The differences were considered significant for *p-value < 0.05.

Applying steered forces on the specific parts of the pentameric MscL to gate the channel has been studied before to explore the gating mechanism of MscL.^6,63 Applying force on the MscL does not give us clear information about the material properties of each component because the result of such SMD simulation depends on: 1) the position on the protein that these forces are applied to, 2) the effect of surrounding lipid⁶⁴ and 3) the level of difference between the material properties of different helices and/or loops. Therefore, full atomistic calculations are required to derive physical parameters also useful for mesoscopic and continuum modeling of MscL in “membrane protein lattices” models. This “bottom-up” approach enables modeling multiple proteins (channels) in an implicit membrane to capture larger time and length scale phenomena.

Effect of water on the elasticity of TM1 helix

Figure 6 shows the influence of solvation in water on the properties of the TM1 helix obtained using constant-force method SMD simulation. We specifically chose TM1 helix as it forms the MscL pore and thus we were interested to see how its properties altered as it became hydrated. The intention here was to determine a range for the stiffness variation between these 2 conditions (hydrated and dehydrated). Solvation in water decreased the stiffness of the TM1 helix in both MtMscL and EcMscL (E = 3.2 ± 0.9 GPa in water and 8.8 ± 0.2 in vacuo for MtMscL and E = 2.6 ± 0.5 GPa in water and 3.5 ± 0.1 GPa in vacuo for EcMscL) (Fig. 6A and C). To examine why solvation resulted in a larger change in stiffness for MtMscL than EcMscL we determined the evolution of hydrogen bonds on the TM1 helix during solvation in each case. For MtMscL the number of hydrogen bonds was significantly reduced from ∼41 to ∼30 as the helix becomes solvated (Fig. 6B). However, the number of hydrogen bonds did not vary as much when the TM1 helix of EcMscL was solvated in water (Fig. 6D).

Figure 6. Effect of water on mechanical properties of the TM1 helix in MtMscL (A,B) and in EcMscL (C,D) using constant-force SMD simulation. (A) The Young's modulus of TM1 helix in MtMscL is E = 3.2 ± 0.9 GPa when it is solvated in water, and E = 8.8 ± 0.2 GPa in the absence of water (vacuum). Thus, in the absence of water (vacuum), TM1 helix is more than 3 times stiffer in MtMscL. (B) Change in number of the hydrogen bonds of TM1 α helix of MtMscL during 10 ns of SMD simulation. About 41 hydrogen bonds stabilize the secondary structure of TM1 helix in vacuum, while when it is solvated in water, the number of hydrogen bonds becomes significantly reduced to ∼30 bonds. This indicates that water acts as a ‘lubricant’ (softener) during TM1 helix elongation in water. (C) The Young's modulus of the TM1 helix of EcMscL solvated in water is E = 2.6 ± 0.5 GPa whereas it is E = 3.5 ± 0.1 GPa in the absence of water (vacuum). (D) The number of hydrogen bonds does not change considerably when the TM1 helix of EcMscL is solvated in water, thus its elasticity modulus only increases by ∼30 %. For all the hydrogen bond calculations, we used a donor-acceptor distance of 3.5 Å and angle cut-off of 40°. The values presented in (A) and (C) are Mean ± SEM for n = 3. Student's t-test was used for statistical analysis and differences were considered significant for *p-value < 0.05.

Discussion

Alpha helices play central roles in various biological processes such as molecular and cellular mechanotransduction, cell mechanics and tissue mechanics and are one of the most important elementary building components of proteins.⁶⁵ Due to a lack of understanding of how these secondary structural elements respond to mechanical force and how they interact with each other, the stiffness and unfolding of proteins remain unclear at atomistic resolution. In this study we have used steered molecular dynamics approaches to investigate the mechanical properties of transmembrane α helices in MscL.

Constant-force SMD is a more reliable measure of helix stiffness than constant-velocity

We compared 2 widely used methods of steered molecular dynamics (SMD) namely, constant-velocity (CV) and constant-force (CF) in a more detailed and systematic way compared to previous studies using a large number of SMD simulations.⁶⁶ We monitored the number of hydrogen bonds during pulling the TM1 helix using both methods. We showed that the number of hydrogen bonds in CV method, which determine/stabilize the secondary structure of the helix, quickly and dramatically decreased compared to the CF method during the stretching (Fig. 4B). Hence, we postulate that the most important advantage of CF method over CV method is that the secondary structure of the helix remains better conserved during the pulling test in the former method compared to the latter one. Although, this is the case for common pulling rates (1-5 pÅ/s), we should note that the CV and CF methods will not be much different, in terms of breaking the H-bonds, if both performed at the same force level (Fig. 4B). The hydrogen bonding state of the helix was also reflected in the rate dependent unfolding during CV SMD which has also been shown in previous studies for different proteins.^39,67 Helix behavior was rate dependent only for the strains <100 % because the helical strain in this region is mainly due to the breakage of hydrogen bonds (Fig. 2B and Movie S2). Whereas in higher strains, the helix behavior is not rate dependent as it manifests the helical backbone elasticity.^40,67 Therefore, unlike in the CF method, the Young's moduli obtained from CV method are dependent on the pulling rate and the assigned spring constant to the dummy atoms. For example, the Young's modulus of TM1 helix measured by CV method can vary from ∼11 GPa to ∼42 GPa depending on the puling rate (Figs 2B and 4A).

Consequently, the absolute Young's modulus obtained from CF method is generally more reliable and reproducible compared to those estimated by using CV method. This information is essential given the great need for employing precise tools for determining the mechanical properties of α-helices of membrane proteins, cytoskeletal elements and collagen-like microfibrils. For instance, it is essential to accurately determine the elasticity of tip links in hair cells given its importance in the force transduction in auditory physiology.⁴⁶ Previous studies have rigorously used the CV method combined with experiments for measuring the elasticity of tip links which is formed by protocadherin 15 and cadherin 23.^46,47 Instead, we suggest more accurate computational elasticity measurements by employing the CF method, which may be used to ensure the accuracy of elasticity constants. Also, provided the essential role of tip links in deafness-related structural defects, this may further uncover the molecular elasticity and physical principles underlying the function of tip links.^68,69 However, the down side of this approach is that the number of simulations required to be done to generate a stress-strain curve using the CF method is much larger than using the CV method (51 CF simulations vs. 3 CV simulations for each helix). Thus, the CF method requires a much larger computational cost. Moreover, it has recently been shown that SMD-like techniques can be vastly improved for relatively fast pulling rates, if one combines Minh-Adib's bidirectional estimator with nonlinear WHAM equations to reconstruct and assess PMFs from trajectories.⁷⁰

Comparison of mechanical properties of T helices of MtMscL and EcMscL

Using the CF method, we determined the Young's moduli of all the helices of MtMscL and EcMscL (Fig. 5). The lowest Young's modulus was found for the N-terminus of MtMscL which is ∼0.2 GPa and the highest determined is for the TM2 helix of MtMscL, which is ∼12.5 GPa. This is well in the range of the elasticity moduli previously measured using various experimental and computational approaches for different proteins ranging from 1 GPa to 30 GPa.^71-75 TM1 helix and N-terminal helices of MtMscL have similar mechanical properties to those in EcMscL, which could be due to their higher level of amino acid sequence conservation and identity in these regions. Among all the helices, the Young's modulus of the TM2 helix and the C-terminal helix appear to be vastly different between M. tuberculosis and E. coli. The TM2 helix in MtMscL is much stiffer than TM2 helix in EcMscL. Hence, overall in the transmembrane region, MtMscL is stiffer than EcMscL. Given that TM2 helix is embedded in the lipid bilayer and determines the bilayer hydrophobic length around the channel,^20,76 we suggest that the stiffer TM2 helix may also correspond to less sensitivity of MtMscL to bilayer thinning during the gating compared to EcMscL, as previously shown by patch-clamp experiments.⁷⁷ Conversely, the C-terminal helix in EcMscL is about 8 times stiffer than in MtMscL. However, removal of up to 20 amino acids from the end of EcMscL has been shown to have no significant effect on gating of MscL, but rather may be involved in pH sensitivity.⁷⁸ At this stage, we are not sure what would be the physiological meaning of having stiffer C-terminal bundle. This will remain to be addressed in the future studies to understand the physiological role of the C-terminus in MscL channels.

Previous reports have shown a dependence of the elasticity and strength of some helix types on their length.⁶⁷ For example it has been shown that the stiffness of the stack of α-helical ankyrin repeats increases with a decreasing number of repeats.⁷⁹ In our study, the elasticity modulus did not correlate with the helix length or hydrophobic moment or number of charged residues in this study. In addition, we demonstrated here that the charge and other factors such as hydrophobic moment did not correlate with the helix stiffness of MtMscL and EcMscL (Figs S3, S4, S5 and S6).

Effect of hydration on mechanical properties of TM helices of MscL

MD simulations have previously shown that TM1 helix, which forms the pore, is only partially hydrated in the closed state and it becomes fully hydrated in the open state.^6,16 Therefore, we were interested to see how the TM1 helix elasticity changes when it is hydrated. By pulling the TM1 helix of both MtMscL and EcMscL in vacuum and water, we have found that exposure of the TM1 helix to water greatly affects its stiffness. Based on our results (Fig. 6 B, D), we suggest that water by interfering with the hydrogen bonds and electrostatic interactions within the helix, acts as a “lubricant” (softener) during the helix elongation. Thus, the Young's modulus of a helix like TM1 of MtMscL decreases significantly in the presence of water. This is not surprising given the previous reports on similar effects on different proteins such as the effect of water on collagen-like micro-fibrils, immunoglobulin domains and spectrin-like proteins.^60,80,81 In fact, the physiological function of this class of proteins is regulated by water since they should become softer and thus unfold easier as they become hydrated.⁸⁰ However, there are other studies which suggest that the elasticity of those α-helices that are primarily determined by the backbone hydrogen bonds do not change significantly when the helix is solvated in water.⁸² An interesting point here is that water did not change the stiffness of TM1 helix of EcMscL to the same extent as it did in MtMscL. In fact, the effect of water on the elasticity of the TM1 helix in MtMscL was much larger than on the elasticity of the TM1 helix of EcMscL (Fig. 6A, C). We showed that this is analogous to the change in the number of the hydrogen bonds in vacuum and in water (Fig. 6B, D). The exact physiological meaning of this difference is unclear at this point. Given the recent computational MscL study, which showed that MscL gates in the presence of water easier than in vacuum, our results further support the importance of hydrophobic interactions for the gating of MS channels and nanopores.^49,52,83-85

Conclusion

In summary, we calculated Young's moduli and elastic constants for all the helices in MtMscL and EcMscL, which are 2 important members of the bacterial MS channel family. Given a considerable number of all-atom MD simulations, we demonstrated in this study that the widely used constant-velocity approach is less reliable compared to the constant-force method, which is more accurate for determining Young's modulus of α-helices, although it is computationally much more expensive. Interestingly, we showed not only the helical properties among the helices of MscL were quite different, but also the properties of some of the MscL helices (e.g. TM2 helix) turned out to be different between E. coli and M. tuberculosis homologs. Therefore, in these cases we checked if there was any correlation between the elastic properties and helical length, charge or amino acid sequence. We also demonstrated that water could act as a “lubricant” (softener) during TM1 helix elongation, which lines the channel pore in MscL. Consequently, the data and methods resulting from this study have wide implications for understanding the force-enduring properties of the α-helices of MscL and their mechanistic role in MscL gating. In addition, atomistic calculations to derive parameters are also useful for a mesoscopic bead-spring and “hybrid continuum-atomistic” models of channels using energy and force matching. This “bottom-up” approach will allow us to assemble multiple protein structures by springs (or elastic beam/rods) for capturing larger time and length scale phenomena, such as modeling of a cluster of proteins (channels) in lipid patches or attached to other intra- /extra-cellular elements.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Author contributions

N.B., O.B., B.M. and Y.J. designed the research; N.B. and O.B. performed the simulations. M.V., R.N., A.P.H., B.M. and Y.J. contributed the computational/analytic tools; N.B., O.B., B.M. and Y.J. analyzed the data; and N.B., O.B., A.P.H., Y.J. and B.M. wrote the manuscript. All authors read and approved the final manuscript.

Acknowledgments

We would like to thank Dr Charles D Cox for insightful discussions and for feedback on the manuscript. Dr Yoshitaka Nakayama and Dr Massimo Vassalli are kindly acknowledged for proofreading of the manuscript. We also acknowledge the Institute for Research in Fundamental Sciences (IPM), the Institute for Nanoscience and Nanotechnology (INST), Sharif University of Technology of Iran and the supercomputing facility at CSIRO, Australia for providing access to the supercomputer equipped with the standard software.

Funding

N.B. has been supported by a University International Postgraduate Award (UIPA) from the University of New South Wales. This project was supported by a grant (APP1044628) and a Principal Research Fellowship to B.M. from the National Health and Medical Research Council of Australia. No potential conflicts of interest were disclosed.