Structural Investigation of MscL Gating Using Experimental Data and Coarse Grained MD Simulations

The mechanosensitive channel of large conductance (MscL) has become a model system in which to understand mechanosensation, a process involved in osmoregulation and many other physiological functions. While a high resolution closed state structure is available, details of the open structure and the gating mechanism remain unknown. In this study we combine coarse grained simulations with restraints from EPR and FRET experiments to study the structural changes involved in gating with much greater level of conformational sampling than has previously been possible. We generated a set of plausible open pore structures that agree well with existing open pore structures and gating models. Most interestingly, we found that membrane thinning induces a kink in the upper part of TM1 that causes an outward motion of the periplasmic loop away from the pore centre. This previously unobserved structural change might present a new mechanism of tension sensing and might be related to a functional role in osmoregulation.


Effect of restraints
Plots of pore radius vs simulation time ( Figure 3 in main article) demonstrated that the distance restraints are much more effective than the solvent restraints in producing an open pore but also suggested that solvent restraints and tension reinforce the distance restraints and have a stabilising effect on the open pore structure. To further investigate this we analysed the distance restraints as a function of simulation time in the presence and absence of solvent restraints and tension ( Figure S1). The data confirms that both EPR and FRET restraints have a strong effect on the protein structure as seen by the rapid drop off in restraint energy in the first 500 ns of simulation. Similar to the results plots of pore radius vs time the data from the first 1500 ns shows some indication that the distance restraints are reinforced by the solvent restraints and tension as the restraint energy drops off faster and to a slightly lower level for simulations involving tension and solvent restraints. Despite this small difference, the restraint energy clearly remains lower for simulations that combine distance restraints with tension and solvent restraints after the restraints are removed at 1500 ns. This supports the stabilising effect of the solvent restraints and tension on the open pore structure. Furthermore, the data shows that FRET restraints re-enforce the EPR distance restraints and also vice-versa. In all simulations the energy of FRET restraints never decreases to the same low level as the EPR restraint energy and FRET restraints show a larger increase when the restraints are removed. This is most likely a consequence of the different implementation of the EPR and FRET restraint energy. The EPR restraints are implemented as a harmonic half-potential while for the FRET restraints a flat harmonic potential was used. In addition, the EPR restraints are exclusively in the TM helices, a domain of the protein that is less mobile. In contrast, the FRET restraints are distributed over the entire protein. Hence it is harder for the system to reach a structure where all FRET restraints are fulfilled. Figure S1: Energy of FRET and EPR distance restraints as a function of simulation time for different simulation protocols. A/B Effect of combining EPR restraints with FRET restraints and solvent restraints. C/D: Effect of tension on distance restraints. The restraint energy is shown as a % of the maximum value at the start of the simulation. Restraints are introduced over 1000 ns, held constant for the next 500 ns and removed during the last 500 ns of the simulation.  Figure S3: Position of C-terminal in the equilibrated closed pore structure, the final structure from a simulation without restraints (N) and open pore structures from simulations with different combinations of restraints. TM1 is shown in orange, TM2 in red, the C-terminal in cyan, the N-terminal in blue and connecting loops in grey. Figure S4: Top view of structures from the the equilibrated closed pore, the final structure from a simulation without restraints (N) and open pore structures from simulations with different combinations of restraints. TM1 is shown in orange, TM2 in red. Backbone particles of residues V23 (dark grey), G26 (light grey), I92 (blue) and I96 (cyan) were used as an indication of orientation of TM1 and TM2.

Results from kink analysis
A first analysis of the formation of a kink in the upper part of TM1 was carried out using a visual inspection of 17 structures from different simulation protocols to identify the presence of a kink and its position. Figure S5 summarises the results by showing the fraction of structures for which a visible kink is present as well as the position of the kink for simulations in POPC without tension, simulations in POPC with tension and simulations in DMPC without tension.  Protocol for structural analysis of simulations Distance Restraints: The energy of the EPR and FRET distance restraints was calculated as a means of testing if the distance restraints are being fulfilled and to investigate the effect of combining distance restraints with solvent restraints or tension. The energy of the restraints was calculated from the difference between the target distance (d open ) and the instantaneous distance (r ij ) between the restrained residues using the potential defined in Eq. 2 and Eq. 3 and summed up over all restrained residues.
Pore radius: The most characteristic feature of the open pore structure is the widening of the hydrophobic gate formed by the backbone particles residues 20 to 28 in TM1. The pore size was estimated by the radius of the circle inscribing the pentagon formed by the 5 subunits and averaged over the residues 20 to 28. The pore radius of open pore structures was determined by averaging the radius calculated from 200 frames in the last 50 ns of the simulations.
Helix tilt: The helix tilt for TM1 and TM2 was estimated by calculating the angle formed by the helix and the normal to the membrane plane. To calculate the helix tilt the direction of TM1 and TM2 were estimated by an axis passing through the backbone particles in each TM helix. The TM1 and TM2 helix tilt was averaged over the 5 subunits and calculated for each frame in the trajectory. The helix tilt for an open pore structure was calculated by determining the average angle of TM1 and TM2 determined by averaging the angles calculated from 200 frames in the last 50 ns of the simulation.

RMSD of functional domains:
To monitor the effect of the various combinations of restraints on the protein, the RMSD between the equilibrated closed structure and the structure at each time step in the simulation was calculated. Plots of RMSD vs time were prepared for the entire protein and the following functional domains: i) the hydrophobic gate formed by residues 20 to 28, ii) transmembrane helices TM1 and TM2 formed by residues 15 to 43 and 72 to 107, iii) the periplasmic loop consisting of residues 44 to 71 and iv) C-terminal formed by 108 to 136.

RMSD between structures:
The RMSD between the closed equilibrated and selected open pore structures was determined by averaging the RMSD values for the specified sections from a RMSD vs time plot calculated from 200 frames in the last 50 ns in the simulation. To ensure the per-residue RMSD is converged the data from the last 50 ns was compared to data sets obtained from 200 frames between 1100 -1150 ns and 1600 -1650 ns. Radial position vs residue: The radial position was used as measure of how far a residue has moved outwards with respect to the pore centre during the simulation and was calculated as the distance between the backbone particle and the pore centre, averaged over 200 frames from the last 50 ns of the simulation. The radial position for each residue was averaged over the 5 subunits. To compare our open pore structures to models from previous studies we calculated the radial position of the α-Carbon in each residue for structures from Sukharev et al [1,2] , Meyer et al [3] and a model from atomistic restrained MD simulations by Corry et al [4].
Kink analysis The kink analysis was carried out on 17 structures from different simulation protocols. The structures were divided into 3 groups; 7 structures from simulations without tension in POPC (with a variety of distance and solvent restraints), 7 structures from simulations with tension in POPC (wit a variety of distance and solvent restraints) and 3 structures from simulations in DMPC. The first to identify the presence of a kink and its position was based on the visual inspection of the average structure calculated from the last 50 ns of each of the different simulations protocols. To confirm our hypothesis we carried out a quantitative analysis that did not use any average structures to remove any potential artefacts caused by the averaging of multiple structures. In the quantitative analysis the kink was defined a RMS deviation from a straight helix as found in the closed pore structure. For each of the 17 structures 25 frames from the last 50 ns of the simulation were analysed using the TM1 helices from the closed pore as a reference structure. Each subunit was compared to the reference structures using a sliding window and the average RMSD was calculated and plotted against residue number. Using this approach the presence of a kink can be detected by an increase in RMSD around the residue where the kink occurs. Note that this approach however does not allow us to identify the position of the kink but only detect its presence. The analysis resulted in 125 RMSD vs residue data sets (5 subunits * 25 frames) for each structure and the presence of a kink was defined as an RMSD value above a given cutoff of 1.5Å. The % occurrence of kinks in the 3 groups of structures was calculated from the RMSD vs residue data sets.

Implementation of restraints
Distance restraints Table 2 and Figure 2 in the main article summarise the raw experimental data used for the restraints and depicts the position of all residues for both accessibility and distance restraints. For clarity, the residues are only shown on 1 of the 5 subunits but FRET and EPR experiments are carried out with fully labelled proteins where each subunit contains a fluorescent or magnetic tag. The inter-subunit distance hence represents the distance between equivalent residues in neighbouring subunits or diagonally related subunits. Based on the pentameric structure a single distance measurement gives rise to a total of 10 distance restraints, 5 from neighbouring subunits and 5 from diagonal subunits.
FRET is more accurate for measuring changes in distances rather than absolute distances. Consequently, data from FRET experiments was reported as the difference between the where ∆d F RET is the change in distance calculated from the FRET experiment. The restraint was implemented using a flat bottom harmonic potential given by: where r ij is the distance between particle i and particle j at any given time in the simulation, k max is the force constant and r 0 and r 1 are defined by the uncertainty of the inter-subunit distance in the FRET experiment such that r 0 < d open < r 1 . Considering the geometry of the pentameric structure of MscL each of the 9 reported inter-subunit distance resulted in 10 target distances giving a total of 90 distance restraints from FRET experiments. The target distances used for the simulations are listed in Table 1 in [4].
In EPR experiments inter-subunit distances are reported by the proximity parameter Ω, which is only suitable for semi-quantitative analysis. Distances are categorised as far, intermediate and close based on their Ω values where Ω = 1 means that the distance between the residues is > 15Å. The data reported [5] shows that all 65 residues in TM1 and TM2 have values of Ω = 1 in the open state making all inter-subunit distances in the open state > 15 A. The restraints were implemented using harmonic half potentials given by: As for the distance restraints from FRET data each inter-subunit distance results in 10 restraints giving a total of 650 distance restraints for residues in TM1 and TM2. All distance restraints were incorporated into the GROMACS topology of the CG model of the closed channel using tabulated potentials.

Solvent restraints
Solvent restraints were based on the lipid and water accessibility measurements of residues in the TM helices and the periplasmic loop from EPR experiments. The raw experimental data was converted into restraints by changing the interactions of the side chain particles with water to represent the change in solvent exposure of selected residues upon opening of the channel.
In EPR experiments lipid exposure of specific residues is measured by the collision frequency of the probe with lipid accessible O 2 and is reported by the parameter ∆P 1/2 O 2 , representing the change in lipid exposure upon channel opening. Similarly, water exposure of residues is measured by monitoring the collision with water soluble NiEDDA and is represented by the parameter ∆P 1/2 N iEDDA. Although lipid and water accessibility are reported separately they are conceptually related and they can be considered as two opposing trends of the same phenomena i.e. an increase in water exposure can be seen as a decrease in lipid exposure. To simplify the conversion of the lipid and water accessibility data into restraints we combined the data from the two measurements into a single parameter called ∆hydrophobic. This parameter represents a change in hydrophobic character of a given residue upon opening of the channel. Hence an increase in lipid exposure, reported as a positive ∆P 1/2 O 2 to reflect the increased collision frequency, is equivalent to a increase in hydrophobic character. A decrease in lipid exposure represents a decrease in hydrophobicity. In contrast, an increase in water exposure, reported as a positive ∆P 1/2 N iEDDA reflects a decrease in hydrophobic character and a decrease in water accessibility means an increase in hydrophobic character.
Since we are interested in the change in hydrophobic character upon opening of the channel the following relationship holds: where ∆∆P 1/2 O 2 and ∆∆P 1/2 N iEDDA are the difference in ∆P 1/2 O 2 and ∆P 1/2 N iEDDA between the closed and the open pore structure as reported from the EPR experiments. A positive value of ∆hydrophobic corresponds to an increase in hydrophobic character while a negative ∆hydrophobic means the residue experiences a decrease in hydrophobic character. Using the reported values for ∆P 1/2 N iEDDA and ∆P 1/2 O 2 , ∆hydrophobic was determined for 54 residues in TM1 and TM2 [5] and 28 residues in the periplasmic loop [3] (see Table 2 in the main article). Some residues showed no or little change in hydrophobic character and, as a guideline, residues with |∆hydrophobic| < 3.0 were not used for restraints. A total of 43 residues in TM1 and TM2 and 10 residues in the periplasmic loop were selected resulting in a total of 265 solvent restraints.
The next step was to incorporate ∆hydrophobic into the MARTINI force field. In the MARTINI force field the strength of non-bonded interactions are modelled using a Lennard-Jones potential. The strengths of the attraction or repulsion between particles is described by interaction levels ranging from supra attractive (0) to super repulsive (IX) (see Table 1 in [6]).
To implement the solvent restraints the value of ∆hydrophobic for the selected residues was used to increase or decrease the interaction level of the side chain particles with water particles (type P4). Based on the spread of ∆hydrophobic values among the 43 residues, we divided the residues into 3 groups: i) residues with a ∆hydrophobic between 3 and 6 were assigned a change in interaction level of 1, ii) for residues with ∆hydrophobic between 7 and 18 a change of interaction level of II and iii) a change in interaction level of III was used for residues with ∆hydrophobic larger than 18. This assignment is a compromise between introducing a change in the interaction to model the restraint without introducing a change that is too large for a given residue type. According to the above definition of ∆hydrophobic and the interaction matrix of the MARTINI force field a negative ∆hydrophobic value (a 13 decrease in hydrophobic character) will result in a decrease of interaction level between water and the selected residue. Using these changes in interaction level as a guideline, the ∆hydrophobic for each residue was mapped to a change in interaction level to create a new set of particle types which contained the altered interaction levels with water but unaltered interactions with all other particle types. The residues Gly22, Leu19 are used as examples to illustrate the changes in interaction level. In the MARTINI force field water is modelled as a particle type P4 while in the MscL protein, based on residue type and secondary structure, Gly22 is modelled as particle type. N0 with an interaction level IV between water (P4) and Gly22 (N0) as defined in the Interaction matrix of the MARTINI force field (see Table 1 in [6]). Gly22 shows a ∆hydrophobic value of 7 and hence its interaction level with water was increased from IV to VI. Similarly, the interaction level of Leu19 which is modelled as a particle type C1 was decreased from VIII to V, based on its ∆hydrophobic value of -21. Tables S1 to S3 show the changes in interaction levels and the new particle types for the solvent restraints in TM1, TM2 and the periplasmic loop.
Introducing restraints during simulation The final simulation system contained a total of 740 distance restraints and 265 solvent restraints that describe the changes in inter-subunit distances and solvent accessibly between the closed and the open pore structure. Hence, the restraints represent a gradual change in the MscL structure during the gating of the channel. We therefore aimed to slowly introduce the restraints during the simulation rather than impose a sudden change onto the closed structure. This was achieved using the free energy perturbation (FEP) method in GROMACS. A coupling parameter λ is used to linearly interpolate between two states of the simulation system. During the simulation λ is slowly increased from λ = 0, where no restraints are present to λ = 1 where restraints are at their maximum. In the case of the distance restraints this means that the force constant of the harmonic potentials is slowly increased from k = 0 at λ = 0 to k max at λ = 1. Similarly, for the solvent restraints the side chain particles of the restrained residues is slowly morphed from the normal particle type (λ = 0) into the new particle type (λ = 1). This means the interaction level of the side chain particle with water is slowly increased or decreased. Table S1: Changes to the interactions between selected side chain particles and water particles to implement solvent restraints for TM1. ∆hydrophobic is based on lipid and water accessibility from EPR [5]. Normal Particle type and Interaction levels refers to the standard particle types in the Martini force field and its interaction with water particles (type P4) [6]. New Particle type and Interaction level refers to the newly created particle types that have altered interactions with water particles based on ∆hydrophobic (P4) but unaltered interactions with all other particle types.  Table S2: Changes to the interactions between selected side chain particles and water particles to implement solvent restraints for TM2. ∆hydrophobic is based on lipid and water accessibility from EPR [5]. Normal Particle type and Interaction levels refers to the standard particle types in the Martini force field and its interaction with water particles (type P4) [6]. New Particle type and Interaction level refers to the newly created particle types that have altered interactions with water particles based on ∆hydrophobic (P4) but unaltered interactions with all other particle types.  Table S3: Changes to the interactions between selected side chain particles and water particles to implement solvent restraints for the periplasmic loop. ∆hydrophobic is based on lipid and water accessibility from EPR [5]. Normal Particle type and Interaction levels refers to the standard particle types in the Martini force field and its interaction with water particles (type P4) [3]. New Particle type and Interaction level refers to the newly created particle types that have altered interactions with water particles based on ∆hydrophobic (P4) but unaltered interactions with all other particle types.