Introduction

The transmembrane (TM) proteins Piezo 1 and 2 are mechanosensitive ion channels (Coste et al., 2010, 2012) that mediate numerous physiological processes in mammals, including touch sensation and blood pressure control (Wu et al., 2017). Cryo-electron microscopy (cryo-EM) structures of Piezo 1 (Guo and MacKinnon, 2017; Saotome et al., 2018; Zhao et al., 2018) and Piezo 2 (Wang et al., 2019) in detergent micelles revealed three identical monomeric arms with 38 TM helices that spiral out from the central ion channel, which is lined by the innermost TM helices of the arms. The TM domains of the helices in the three arms do not lie in a plane, which led to the suggestions that the Piezo protein curves the cell membrane into a nanodome, and that this protein-membrane nanodome flattens when external forces induce a tension γ in the membrane, leading to channel opening (Guo and MacKinnon, 2017). The flattening of the nanodome is driven by the energy γΔA, where ΔA is the excess area of the nanodome due to its curved shape, compared to the projected area of the nanodome in the plane of the surrounding membrane (Guo and MacKinnon, 2017; Haselwandter and MacKinnon, 2018). Cryo-electron tomography images of reconstituted membrane vesicles with embedded Piezo 1 proteins indeed show that the proteins induce membrane curvature (Guo and MacKinnon, 2017), but elastic modeling of vesicle shapes with varying diameters indicates that this induced curvature is a on average about four times smaller than the curvature of the proteins in the detergent micelles of the high-resolution cryo-EM structures (Haselwandter et al., 2022b,a). High-resolution fluorescence imaging of Piezo 1 in its native cell membrane environment confirms an expansion in the inactivated state from flattening of the arms compared to structural models in detergent micelles, and indicates further flattening upon activation (Mul-hall et al., 2023). A cryo-EM structure of flattened Piezo 1 with widened ion channel has been obtained by embedding Piezo 1 proteins in small membrane vesicles with a diameter of 20 nm in an outside-out orientation in which the vesicle curvature opposes the intrinsic Piezo curvature (Yang et al., 2022). In patch-clamp experiments with Piezo 1 and Piezo 2, an opening of the ion channel induced by membrane tension has been observed for Piezo 1, but not for Piezo 2, which appears to indicate that Piezo 2 activation requires other factors of its cellular context that are not present in the experiments, besides membrane tension (Moroni et al., 2018).

Molecular dynamics (MD) simulations have the potential to complement recent experimental insights on the structure of Piezo proteins in their native membrane environment by providing high-resolution, dynamic information of membrane-embedded Piezo proteins under different membrane tensions (Botello-Smith et al., 2019; Chong et al., 2021; Buyan et al., 2020; Lin et al., 2022; Jiang et al., 2021; De Vecchis et al., 2021). In atomistic simulation trajectories with a length up to 100 ns at membrane tensions between 14.2 mN/m and 67.8 mN/m starting from a tensionless Piezo 1 protein-membrane nanodome with a simulation-box area of 31.4 × 31.4 nm², De Vecchis et al. (2021) observed a flattening of the nanodome, and a significant expansion of the channel volume at the largest tension of 67.8 mN/m. In atomistic simulations of truncated Piezo 1 with TM helices 17 to 38 in a tensionless membrane, a flattening of the protein-membrane nanodome leading to channel opening has been induced by a small simulation-box area of about 480 to 485 nm² that aims to mimic a crowding of Piezo 1 proteins (Jiang et al., 2021). The membrane flattening in these simulations is a consequence of the standard periodic boundary conditions of the simulation box, which lead to an on average vanishing slope of the membrane at the box boundaries.

In this article, we report results from coarse-grained simulations of the Piezo 1 and Piezo 2 protein-membrane nanodome with a length up to 8 μs at membrane tensions between 0 and 20.8 mN/m, starting from tensionless protein-membrane nanodomes with a simulation-box area of 50 × 50 nm². To corroborate these coarse-grained simulation results, we performed additional atomistic simulations of the Piezo 2 protein-membrane nanodome with a length up to 300 ns and an initial, tensionless simulation-box area of 33 × 33 nm² at membrane tensions between 0 and 18 mN/m. A main focus in our analysis is on determining the excess area ΔA of the protein-membrane nanodomes as a function of the membrane tension γ. In both our coarse-grained simulations with the Martini 2.2 force field (de Jong et al., 2013a) and atomistic simulations with the CHARMM36 force field (Huang and MacKerell Jr, 2013), we obtain an excess area ΔA for the protein-membrane nanodome of about 40 nm² in the tensionless state, a reduction to values of about 5 nm² and below at tension values larger than 10 mN/m at which the nanodome is nearly completely flattened, and a half-maximal response of ΔA at tension values of about 3 to 4 mN/m, which is within the range of experimentally determined values for the half-maximal activation of Piezo 1 (Lewis and Grandl, 2015; Cox et al., 2016). In agreement with experimental observations for membrane-embedded Piezo 1 (Haselwandter et al., 2022b,a; Mulhall et al., 2023), the excess area ΔA as well as the height and curvature of the tensionless protein-membrane nanodomes is significantly reduced compared to corresponding values for the cryo-EM structure of the proteins in detergent micelles. At high membrane tensions, we observe a widening of the Piezo1 ion channel as in the flattened cryo-EM structure of Piezo 1 in small membrane vesicles (Yang et al., 2022). The Piezo2 ion channel, in contrast, does not respond to membrane tension in our simulations, which is in line with patch-clamp experiments in which Piezo 1, but not Piezo 2, was shown to be activated by membrane tension alone (Moroni et al., 2018). In addition to ΔA, we determine the bending energy of the lipid membrane in the nanodomes, which allows to construct an elasticity model that identifies the energetic contributions of the protein and the membrane in the response of the Piezo protein-membrane nanodome to membrane tension.

Results

Piezo 2 protein structure and nanodome shape in tensionless membranes

Our simulations of the Piezo 2 protein-membrane nanodome are based on the cryo-EM structure of the Piezo 2 trimer (PDB ID 6KG7) (Wang et al., 2019), which includes all 38 TM helices of the Piezo 2 monomers. In these simulations, the Piezo 2 trimer is embedded in an asymmetric membrane that mimics essential aspects of the lipid composition of cell membranes (see Methods). We find that the relaxed simulation structure of Piezo 2 in tensionless membranes is significantly flattened compared to the cryo-EM structure of the Piezo 2 trimer in detergent micelles (see Figure 1). To compare these structures, we have performed coarse-grained simulations of membrane-embedded Piezo 2 in which the protein structure is constrained to the cryo-EM structure, in addition to simulations without such constraints in which the protein structure relaxes in the tensionless membrane within the 8 μs of the simulation trajectories. The two exemplary conformations in Figure 1(a) compare the protein-membrane nanodome at the end of simulation trajectories with and without constraints on the protein structure. In these conformations, the TM helices of Piezo 2 are indicated as rods, and the shape of the membrane midplane is shown as a continuous surface determined from the positions of the lipid head beads in the two leaflets of the membranes (see Methods). The relaxed simulation structure of Piezo 2 is clearly less curved than the constrained cryo-EM structure. In addition to these different curvatures of the protein-membrane nanodome, the membrane shape in both exemplary conformations of Figure 1(a) exhibits undulations from thermally excited shape fluctuations.

Figure 1.

To quantify the protein-induced curvature and excess area of the nanodome, and to eliminate the additional curvature and excess area of thermal membrane shape fluctuations present also in protein-free membranes, we average the continuous membrane-midplane shapes of individual conformations over the same time intervals of all trajectories. The average membrane shapes, membrane contours, and lipid densities in Figure 1(a) are averaged over the conformations in the last microseconds of all 10 simulation trajectories of length 8 μs with the relaxed Piezo 2 protein, and over the last microseconds of the 3 simulation trajectories with the constrained protein (see Methods). Because of the translational and rotational diffusion of the Piezo proteins along the simulation trajectories, the averaging of conformations requires an alignment of the protein (see Methods). The average membrane shapes and lipid densities in Figure 1(a) therefore have a circular projected area with a diameter that corresponds to the average x and y dimensions of the simulation boxes of the conformations, because only these ‘inscribed’ circular membrane segments are overlying in all rotationally aligned conformations. From the simulation trajectories in which the protein structure is constrained to the cryo-EM structure of Piezo 2 in detergent micelles, we obtain a nanodome midplane shape with a height of 10 nm, which is significantly larger than the height of about 4 nm for the nanodome shape with the relaxed Piezo 2 protein. The mean curvature profiles of the nanodome shapes in Figure 1(c), which we calculate from the radially averaged nanodome profiles of Figure 1(b), exhibit maxima at a radial distance of about 4 nm from the channel. In these mean curvature profiles, the maximal mean curvature of the nanodome with constrained protein structure is about two times larger than the maximal mean curvature of the nanodome with the relaxed protein. The mean curvature of the nanodome with the relaxed Piezo 2 protein drops to 0 for radial distances larger than about 17 nm from the channel center, which indicates a catenoidal shape of the non-planar surrounding membrane with a principal curvature in the radial direction of the shape profile in Figure 1(b) that is oppositely equal to the principal curvature in the perpendicular, tangential direction (Haselwandter and MacKinnon, 2018).

Tension-induced flattening of the protein-membrane nanodome

With increasing membrane tension, the nanodome shape and coarse-grained simulation structure of the membrane-embedded Piezo 2 protein flattens further (see Figure 2). The average height of the nanodome decreases from 4.4 nm in tensionless membranes to about 3.5 nm, 2.9 nm, and 2.1 nm at the intermediate membrane tensions γ of 1.4 mN/m, 2.8 mN/m, and 5.5 mN/m of our simulations, and to 1.6 nm and 0.9 nm at the largest tensions of 10.8 mN/m and 20.8 nM/m (see radial nanodome profiles in Figure 2(d)). In atomistic simulations of membrane-embedded Piezo 2, we observe a comparable tension-induced flattening of the Piezo 2 protein-membrane nanodome, and a comparable height of the nanodome in tensionless membranes, which corroborates our coarse-grained simulation results (see Figure 3).

Figure 2.

Figure 3.

To obtain values of the excess area ΔA that are not limited by the finite simulation time, in particular for the atomistic simulations, we divide all simulation trajectories into 5 equal intervals and determine the nanodome shape in each interval by averaging over the conformations of all independent simulation runs in this interval. In Figure 4(a) to (c), we plot the excess area ΔA of the nanodome shape in the last time intervals versus the inverse time t⁻¹ of the interval midpoints, for our coarse-grained simulations of membrane-embedded Piezo 1 and Piezo 2 proteins and our atomistic simulations of the membrane-embedded Piezo 2 protein. The excess area ΔA here was calculated as the area difference between the curved nanodome shape and the planar projection of the nanodome (see Methods). As starting points of our coarse-grained simulations of membrane-embedded Piezo 1, we devised a structure of full-length Piezo 1 by adding the unresolved 12 N-terminal TM helices to the Piezo 1 cryo-EM structure with PDB ID 6B3R (Guo and MacKinnon, 2017) via homology modelling (see Methods). In all simulation systems, we observe a decrease of ΔA with increasing simulation time, and therefore extrapolate the ΔA values to t⁻¹ = 0 by linear fitting (dashed lines and shaded prediction error bands). In Figure 3(e), the extrapolated ΔA values are plotted versus the membrane tension γ in the simulations. In all simulation systems, we obtain a maximal excess area ΔA for the protein-membrane nanodome of about 40 nm² in tensionless membrane, and half-maximal ΔA values at tensions of about 3 to 4 mN/m, which is within the range of experimentally determined values for the half-maximal activation of Piezo 1 (Lewis and Grandl, 2015; Cox et al., 2016). The nanodome excess area ΔA of about 40 nm² in tensionless membranes is significantly smaller than the excess area ΔA = 165 ±1 nm² of the average membrane nanodome shape with constrained cryo-EM structure of Piezo 2 shown in Figure 1(a).

Figure 4.

Figure 4—figure supplement 1.

Elasticity model of the Piezo protein-membrane nanodome

We now aim to construct a model for the elastic energies of the Piezo protein and the membrane in the tension-induced nanodome flattening observed in our simulations. The nanodome height in tensionless membranes and the tension-induced flattening of the nanodome are a consequence of opposing protein and membrane energies. The energy of the membrane is the sum E_m = E_b + γΔA of the membrane bending energy E_b and the tension-associated energy γΔA, and ‘prefers’ a planar state in which this energy tends to zero. The energy E_p of the Piezo protein conformation, in contrast, prefers a curved nanodome state. The nanodome height is determined by the sum of these energies (Guo and MacKinnon, 2017; Haselwandter and MacKinnon, 2018), i.e. by the total energy

To describe the tension-dependent Piezo protein conformation, we first note that the vertical displacement z₄ of the fourth 4-TM unit of the Piezo 2 arms relative to the channel center is larger than the vertical displacements z_i of the other 4-TM units in Figure 2(e). We therefore use the vertical displacement z₄ of 4-TM unit 4 to describe the vertical height of the Piezo 2 TM domain, and the shift Δz(γ) = z₄(0) − z₄(γ) of this displacement in comparison to the tensionless nanodome as parameter for the tension-induced flattening of Piezo 2. Figure 5(a) illustrates that the Piezo 2 height change Δz in our coarse-grained simulations is proportional to the membrane tension γ for tension values up to 6 mN/m.

Figure 5.

Figure 5—figure supplement 1.

At the equilibrated Piezo 2 height obtained at a given tension value, the vertical forces resulting from the protein elasticity and from the membrane elasticity are oppositely equal. If this was not the case, the Piezo 2 height would further increase or decrease until such a force balance is achieved at this tension value. We can therefore calculate the vertical force exerted by the flattening-resisting protein from the opposing force exerted by the membrane. At a given membrane tension γ, the vertical force resulting from the membrane energy E_m is

This force is the sum of two terms, and calculating these terms requires to determine the membrane bending energy E_b and the excess area ΔA of the nanodome as functions of the Piezo height change Δz. In Figure 5(b), we plot the ΔA values obtained from the coarse-grained simulations of membrane-embedded Piezo 2 at different tensions versus the Δz values of the protein at these tensions, which indicates a linear relation within the error margins obtained from the simulations. From the linear fit in Figure 5(b), we obtain dΔA/dΔz = −14.6 ± 1.8 nm.

For calculating the membrane bending energy of the averaged nanodome shapes, we use a discretized version of the bending energy and include the lipid densities shown in Figure 2(c) in order to limit the bending energy calculations to the lipid membrane of the nanodome (see Methods). To avoid limitations due to the finite simulation times, we extrapolate the bending energy values obtained for the average nanodome shapes in 5 intervals of the simulation trajectories to long timescales akin to the temporal extrapolations of the excess area ΔA (see Figure 4(d) and (e)). The resulting, extrapolated bending energy values E_b of the lipid membrane in the Piezo 2 protein-membrane nanodome decrease linearly with the Piezo 2 height change Δz (see Figure 5(c)). From linear fitting, we obtain dΔE_b/dΔz = −0.82 ± 0.17 κ/nm where κ is the bending rigidity of the lipid membrane. For the coarse-grained membrane of the lipid composition used in our simulations with main component POPC, we estimate a bending rigidity value κ = 25 ± 5 k_BT based on simulation results for coarse-grained membranes composed of POPC and of various lipid mixtures (Fowler et al., 2016).

From equation 2 and the derivatives dΔA/dΔz and dΔE_b/dΔz determined above via linear modelling, we obtain the four data points in Figure 5(d) for the four Δz values at the membrane tensions γ = 0, 1.4, 2.8, and 5.5 mN/m, which indicate a linear increase of the force F with the Piezo height change Δz. From linear fitting, we obtain the force F (Δz = 0) = 88 ± 15 pN for the tensionless membrane and the force increase dF /dΔz = 60 ± 20 pN/nm. This linear force relation is consistent with a harmonic-spring model for the elasticity of the Piezo protein. Extrapolation to zero force in Figure 5(d) indicates a force-free state of Piezo 2 for Δz = −1.5 ± 0.6 nm, which is – within the statistical accuracy – in accordance with the difference of −2.0 ± 0.2 nm between the vertical displacement z₄(γ = 0) = 3.7±0.2 nm of the membrane-embedded Piezo 2 protein in tensionless membranes and the vertical displacement z₄ = 5.7 nm in the cryo-EM structure of Piezo 2 in detergent micelles (see Figure 1(d)). Our force modelling thus is consistent with the Piezo cryo-EM structure in detergent micelles as force-free conformation of the protein.

From our coarse-grained simulation data of membrane-embedded Piezo 1, we obtain a linear force relation with force F (Δz = 0) = 113 ± 28 pN in the tensionless state and a force increase dF /dΔz = 63 ± 30 pN/nm, which agrees with our results for Piezo 2 within errors (see Figure 5– figure supplement 1). The larger errors in the force relation for Piezo 1 compared to Piezo 2 are also a consequence of fewer tension values considered in our coarse-grained simulations of membrane-embedded Piezo 1 and, thus, fewer data points in the modelling.

Response of the Piezo 1 and 2 ion channel to large membrane tensions

Curved and flattened cryo-EM structures obtained for Piezo 1 embedded in small membrane vesicles with outside-in and outside-out protein orientation indicate a widening of the Piezo 1 ion channel in response to flattening (Yang et al., 2022). In the outside-in orientation, the extracellular side of Piezo 1 points towards the inside of the vesicles with mean diameter of 20 nm, which results in highly curved Piezo structure. In the outside-out orientation, the extracellular side of Piezo 1 points towards the vesicle outside, and the Piezo protein structure is nearly completely flattened because the vesicle curvature opposes the intrinsic Piezo curvature in this orientation. In the flattened Piezo 1 structure, the distance between the residues L2469 of the three TM helices 38 that line the ion channel is increased by 5 Å compared to the curved Piezo 1 structure, which indicates a widening of the outermost constriction site of the Piezo 1 ion channel formed by these residues (Yang et al., 2022).

In our coarse-grained simulations of membrane-embedded Piezo 1, we observe a comparable widening of the ion channel lined by the three TM helices 38 in response to tension-induced flattening of the Piezo 1 protein. At the largest simulated membrane tension γ = 32 mN/m, the Piezo 1 protein is nearly completely flattened, and distance between the channel-constricting residues L2469 of the three TM helices 38 is increased by about 3 ± 1 Å compared to the Piezo 1 channel in our tensionless membranes (see Figure 6). This increased distance between the residues L2469 at large membrane tensions is associated with an increased tilt angle of the TM helices 38 relative to the vertical z direction of our simulation conformations. In our coarse-grained simulations of membrane-embedded Piezo 2, in contrast, we do not observe a response of the ion channel to tension (see Figure 6), apparently in line with patch-clamp experiments of Piezo 1 and Piezo 2, in which an tension-induced opening of the ion channel has only been observed for Piezo 1, but not for Piezo 2 (Moroni et al., 2018).

Figure 6.

Discussion and Conclusions

Based on coarse-grained and atomistic simulations of membrane-embedded Piezo proteins, we have investigated the shape and excess area ΔA of the protein-membrane nanodome as well as the conformation of the Piezo protein in tensionless membranes and for a physiologically relevant range of membrane tensions. Our coarse-grained simulations of membrane-embedded full-length Piezo 1 and 2 proteins allowed to consider relatively large segments of about 50 × 50 nm² on simulation timescales up to 8 μs. To assess and corroborate the validity of these coarse-grained simulations, we performed atomistic simulations of membrane-embedded full-length Piezo 2 in smaller membrane segments of about 33 × 33 nm² on simulation timescales up to 300 ns. To quantify the protein-induced curvature and excess area of the nanodome, and to eliminate the additional curvature and excess area of thermal membrane shape fluctuations, we have determined average nanodome shapes from the simulation conformations (see Figure 1). We have determined the excess area ΔA from extrapolation of ΔA values for average nanodome shapes in different time intervals of the trajectories to reduce sampling errors from the finite length of the simulation trajectories (see Figure 4).

As a main result, we obtained an excess area ΔA of the Piezo protein-membrane nanodome of about 40 nm² in tensionless membranes from both coarse-grained and atomistic simulations. This excess area of the relaxed protein-membrane nanodome in tensionless membranes is significantly smaller than the excess area ΔA = 165 ± 1 nm² obtained from our simulations with the constrained cryo-EM structure of Piezo 2, which in turn is comparable in magnitude to the estimate ΔA = 250 nm² from an analysis of the Piezo 2 dimensions in the cryo-EM structure (Wang et al., 2019). The excess area ΔA is reduced to values of about 5 nm² at tension values larger than 10 mN/m at which the nanodome is nearly completely flattened in our coarse-grained and atomistic simulations, and is half-maximal at tension values of about 3 to 4 mN/m, which are within the range of experimentally determined values for the half-maximal activation of Piezo 1 (Lewis and Grandl, 2015; Cox et al., 2016). Based on our results, the change of the nanodome excess area ΔA required for Piezo activation thus can be estimated to be about 20 nm², which is significantly smaller than previous estimates based on cryo-EM structures, but still large compared to the area difference between the open and closed state of the ion channel and, thus, in line with the suggested flattening-induced Piezo activation (Guo and MacKinnon, 2017).

In addition, our elasticity analysis of the nanodome shapes allowed to construct an elastic model for Piezo activation that distinguishes the different energy components in the tension-induced flattening of the protein-membrane nanodome. According to this model, the Piezo protein elasticity can be approximately described as a harmonic spring, with a force-free conformation that agrees with the curved protein conformation in the cryo-EM structures within the modeling accuracy. The elastic energy of the membrane, which opposes the protein elasticity, is the sum of the membrane bending energy E_b and the tension-associated energy term γΔA that drives the tension-induced flattening. In tensionless membranes, the interplay of the protein elastic energy E_p with the membrane bending energy E_b leads to a flattening of the protein relative to cryo-EM structure that is associated with a force of roughly 100 pN acting on the protein (see Figure 5(d)). With further tension-induced flattening, the force on the protein increases by about 60 pN per nanometer in reduced protein height. With atomic force microscopy (AFM), Lin et al. (2019) measured the force-induced height change of Piezo1 protein-membrane nanodomes adsorbed to substrates. In this system, the force-induced height change is not only affected by the elastic energy of the protein and the membrane, but also by the adhesion energy of the membrane to the substrate, which led to the spreading of membranes vesicles with embedded Piezo 1 on the substrate in the generation of the substrate-adsorbed membranes, and may also be affected by the water pocket between the nanodome and the substrate. The smaller force response of about 7 pN per nanometer in reduced nanodome height measured in the AFM experiments therefore cannot be compared to the force response of the protein obtained from our modeling of simulation data for the tension-induced flattening of the Piezo protein-membrane nanodome.

Materials and methods

Coarse-grained simulations

Protein modeling

Our Piezo 2 simulation system is based on the cryo-EM structure of the full-length Piezo 2 trimer (PDB ID 6KG7) (Wang et al., 2019), which includes all 38 TM helices of the Piezo 2 monomers. We added shorter loops of less than 25 residues not resolved in this cryo-EM structure with the software MODELLER v10.2 (Šali and Blundell, 1993; Webb and Šali, 2016). The remaining unresolved longer loops not included in our Piezo 2 simulation system are indicated by dashed lines in Fig. 7. In cryo-EM structures of Piezo 1 (Guo and MacKinnon, 2017; Saotome et al., 2018; Zhao et al., 2018), several of the outer, N-terminal TM helices are not resolved. We therefore used the cryo-EM structure of full-length Piezo2 as template structure to add the unresolved 12 N-terminal TM helices to the Piezo 1 cryo-EM structure with PDB ID 6B3R (Guo and MacKinnon, 2017) via homology modelling with MODELLER v10.2. Shorter loops of less then 25 residues not resolved in the resulting structure were added as for Piezo 2. To avoid artefacts from missing long protein loops, the distances between the terminal residues of these missing loops were constrained to the distances in the cryo-EM structures using a harmonic potential with spring constant 10 kJ mol⁻¹ nm⁻² in the coarse-grained simulations.

Figure 7.

Force field

We used the Martini 2.2 force field, which exhibits a general 4 to 1 mapping of non-hydrogen atoms into coarse-grained simulation beads (de Jong et al., 2013a; Periole and Marrink, 2013; Marrink et al., 2007; Monticelli et al., 2008), in all coarse-grained simulations performed with the software GROMACS 2018.3 (Abraham et al., 2015). To allow for tension-induced conformational changes of Piezo 2, we employed the standard protein secondary structure constraints of Martini 2.2 with secondary structures identified with DSSP (Kabsch and Sander, 1983), but did not constrain the protein tertiary structure by an additional elastic network (Periole and Marrink, 2013). Instead, we relied on the capabilities of the Martini coarse-grained force field for modeling membrane systems with TM helix assemblies (Sharma and Juffer, 2013; Chavent et al., 2014; Majumder and Straub, 2021).

Membrane embedding

To set up the protein-membrane nanodome of the coarse-grained Piezo 2 simulation system, we first embedded an single monomer of the Piezo 2 trimer into a membrane because the transmembrane region of a Piezo monomer is not strongly curved and can be placed reasonably well in a planar membrane. To this end, we converted the Piezo 2 monomer structure into the coarse-grained Martini 2.2 representation (de Jong et al., 2013a; Marrink et al., 2007; Monticelli et al., 2008) using the martinize script (de Jong et al., 2013b) and generated a coarse-grained planar, asymmetric membrane of area 50 nm × 50 nm with the composition indicated in Table 1 using the INSert membrANE (INSANE) software script (Wassenaar et al., 2015). We then oriented and centered the transmembrane domain of the Piezo 2 monomer along the x-y plane of the membrane and solvated the system with coarse-grained water beads in GROMACS 2018.3 (Abraham et al., 2015). To equilibrate the membrane around the Piezo 2 monomer, we constrained the protein beads using harmonic potentials with force constant 10 000 kJ mol⁻¹ nm⁻², performed an energy minimization with 5000 steps of the steepest descent algorithm, and run a subsequent MD simulation with a length of 50 ns in the NPT ensemble. In this simulation, the pressure was kept at 1 bar using the Berendsen barostat (Berendsen et al., 1984) and the temperature was maintained at 310 K using a velocity-rescale thermostat (Bussi et al., 2007) with individual coupling of protein, membrane, and solvent. The distance cutoff for the non-bonded interactions was 1.2 nm in this simulation, and bonds were constrained with the LINCS algorithm (Hess et al., 1997).

Table 1.

To embed the Piezo 2 trimer, we first generated a Piezo 2 monomer with ‘lipid envelop’ by selecting all lipid molecules of the membrane-embedded protein monomer that contain at least one bead with distance smaller than 3 nm from a protein backbone bead. We then superimposed this lipid-enveloped monomer along each of the three monomers of the Piezo 2 cryo-EM structure with Pymol 2.5.0 (Schrödinger, LLC, 2015) and deleted overlapping lipids from the three envelops. We then placed the resulting lipid-enveloped Piezo 2 trimer into a coarse-grained, planar and asymmetric membrane of area 50 nm × 50 nm with the center of mass of the Piezo 2 ion channel about 4 nm below the membrane midplane (see Fig. 8(a)). In this placement, the Piezo 2 trimer remains oriented as in the cryo-EM structure, with the ion channel roughly oriented along the z-axis and, thus, perpendicular to the x-y plane of the membrane. To merge the curved lipid envelop of the Piezo 2 trimer and the planar membrane into a membrane nanodome, we first added water beads and sodium and chloride ions to neutralize the overall system charge at a physiological salt concentration of 150 mM with GROMACS. The system was then energy minimized in 10 000 steps of steepest descent and equilibrated in a 10 ns MD simulation in the NVT ensemble with an integration time step of 10 fs and a subsequent 500 ns MD simulation in the NPT ensemble with an integration time step of 20 fs. In these mimimization and equilibrium steps, the protein beads were constrained by a harmonic potential with force constant of 1000 kJ mol⁻¹ nm⁻² to maintain the protein shape. The temperature in the equilibration simulations was kept at 310 K by a velocity-rescale thermostat, and the pressure in the NPT simulations was kept at 1 bar with the Berendsen barostat. During these equilibration simulations, a continuous membrane nanodome is formed (see Fig. 8(b)). To embed the Piezo 1 trimer, we took advantage of the structural similarity to Piezo 2 and replaced the membrane-embedded Piezo 2 trimer by the modelled full-length Piezo 1 trimer after structural alignment in Pymol.

Figure 8.

Equilibration

After these membrane embedding and lipid equilibration steps, we slowly released the position restraints on the protein atoms by reducing the force constant to 1000, 100, 10 and 1 kJ mol⁻¹ nm⁻² in four MD simulation steps of 10 ns each in the NPT ensemble with an integration time step of 20 fs. In these simulations, protein, membrane, and solvent beads where coupled independently to an external bath to keep the temperature at 310 K using a velocity-rescale thermostat (Bussi et al., 2007). The lateral and normal pressure was maintained at 1 bar using semi-isotropic pressure coupling and the Berendsen barostat (Berendsen et al., 1984) with a collision frequency of 5 ps. We implemented bond constraints using the LINCS algorithm (Hess et al., 1997) and employed a cut-off of 1.2 nm for the non-bonded interactions. As last step, we removed the position restraints on protein beads and run 10 independent trajectories of length 200 ns to generate 10 different input structures for the production simulation runs.

Coarse-grained production simulations

In the production simulation runs, we switched to the Parrinello-Rahman barostat (Parrinello and Rahman, 1981) with a coupling constant of 20 ps to achieve more flexibility for box shape modulations with semi-isotropic pressure coupling (Braun et al., 2019). All other simulation parameters were identical to the parameters of the final equilibration simulations. To perform simulations at different membrane tensions, we kept the pressure in the normal z direction at P_z = 1 bar and applied a series of pressure values in the lateral x and y directions of the membrane plane. These pressure values were P_x = P_y = 1, 0.5, 0, −1, −3, −7, and −11 bar in our coarse-grained Piezo 2 simulations and P_x = P_y = 1, 0, −1, −3 and −11 bar in our coarse-grained Piezo 1 simulations. The pressure difference in the lateral and normal directions leads to the membrane tension (Yefimov et al., 2008)

where L_z is the average box height during the simulations. At each lateral pressure, we generated 10 independent production trajectories with a length of 4 μs each. For Piezo 2, we extended the 10 simulation trajectories at zero membrane tension to 8 μs.

Atomistic simulations

System setup

To set up our atomistic simulations of membrane-embedded Piezo 2, we used a coarse-grained Piezo 2 conformation obtained after 200 ns of final equilibration as starting point. The coarse-grained Piezo 2 trimer and surrounding membrane nanodome of this conformation were converted to atomistic resolution with the CHARMM-GUI server (Wassenaar et al., 2014; Jo et al., 2008; Brooks et al., 2009; Lee et al., 2016). We reduced the membrane area to about 33 × 33 nm² to decrease the total number of atoms in the system and removed atomic clashes of lipids and amino acids by energy minimization with 5000 steps of the steepest descent algorithm in GROMACS. We then solvated and neutralized the system at a total salt concentration of 150 mM of KCl by adding an appropriate number of potassium and chloride ion. In our atomistic simulations, we used the CHARMM36 force field (Huang and MacKerell Jr, 2013) with the TIP3P water model (Jorgensen et al., 1983).

System equilibration

We performed the equilibration and production simulations with the Amber20 software (Case et al., 2020) because of the efficient use of GPUs with this software (Salomon-Ferrer et al., 2013). To this end, we first generated Amber suitable input coordinates using the ParmEd software tool (Shirts et al., 2017) and energy-minimized the system with 10 000 steps of the steepest descent algorithm and subsequent 10 000 steps of the conjugate gradient algorithm. We then heated the system from temperature 0 K to 310 K in three simulation steps of 10 ns in the NVT ensemble using an integration time step of 1 fs and a Langevin thermostat with a collision frequency of 1 ps⁻¹. In these three heating simulations, the positions of the protein backbone atoms were restrained by harmonic potentials with a force constant of 10 kcal mol⁻¹ Å⁻². We subsequently released the harmonic restraints in 5 simulation steps of 10 ns in the NPT ensemble, decreasing the force constant by a factor of 2 in each of the steps. In these equilibration simulations with an integration time step of 2 fs, the temperature was maintained at 310 K by a Langevin thermostat with a collision frequency of 1 ps⁻¹, the pressure in normal and lateral directions was kept at 1 bar by a Berendsen barostat with a semi-isotropic pressure coupling and a pressure relaxation time of 2 ps, a cutoff length of 10 Å was used for non-bonded interactions, long-range electrostatic interactions were calculated with the Particle Mesh Ewald (PME) method (Darden et al., 1993; Essmann et al., 1995), and all bonds involving hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert et al., 1977).

Production simulations

In the production simulations starting from the equilibrated system conformation, we used hydrogen mass repartitioning (Hopkins et al., 2015) to increase the integration timestep to 4 fs for computational efficiency. Our production simulations consist of 5 independent trajectories at each of the membrane tension values 0, 1, 2, 4, 8, and 18 mN/M. In our simulations with the Amber20 software, the membrane tension was implemented by subjecting a constant surface tension at the two interfaces of the lipid bilayer in the x-y plane, using the Berendsen barostat with a reduced pressure relaxation time of 1 ps to ensure rapid pressure regulation under constant surface tension conditions. As in the equilibration simulations, the temperature was maintained at 310K by a Langevin thermostat with a collision frequency of 1 ps⁻¹, bonds containing hydrogen atoms were constrained using the SHAKE algorithm, a cutoff length of 10 Å was used for non-bonded interactions, and long-range electrostatic interactions were calculated with the PME method. At each of tension values 0, 2, 4, and 8 mN/M, we generated 5 independent trajectories with a trajectory length of 300 ns. At the tension values 1 and 18 mN/M, we generated 5 trajectories with a length of 180 ns each.

Analysis of trajectories

Membrane shape of simulation conformations

To determine the membrane midplane shape of a coarse-grained simulation conformation, we first translated the conformation so that the center of mass of the ion channel formed by the three helices 38 is located at the center of the simulation box. We next discretized the x-y plane of the simulation box into a square lattice with 25 × 25 equally sized squares (bins), separated the lipid head beads (head beads of POPC, POPE, POPS, POP2 and DPSM) into upper and lower membrane leaflet, and allocated these head beads to bins based on the x-y coordinates of the beads. In each bin with at least two lipid head beads in the upper and in the lower leaflet, we calculated the average z coordinate z_upper of the lipid head beads in the upper leaflet in this bin, the average z coordinate z_lower of the lipid head beads in the lower leaflet, and the midplane coordinate z_midplane of the bin as the average of z_upper and z_lower. We finally performed a bilinear interpolation of the calculated midplane coordinates of the square lattice with Mathematica 13 to determine z coordinates also for bins with less than two lipid head beads in the upper or lower leaflet, and to obtain continuous midplane shapes. We determined these membrane midplane shapes for simulations frames at intervals of 20 ns, so for e.g. 201 frames of a coarse-grained simulation trajectory with a length of 4 μs. The membrane midplane shape of atomistic simulation conformations was calculated similarly, for a square lattice of 17 × 17 equally sized bins because of the smaller membrane size in our atomistic simulations, using the coordinates of the common phosphorus atoms in the lipid heads for bin allocation and averaging. Along the atomistic simulation trajectories, we determined membrane midplane shapes for simulation frames at intervals of 1 ns.

Average membrane shapes and lipid densities

To calculate average membrane midplane shapes from the continuous midplane shapes of individual conformations, we first aligned the membrane shapes by rotations around an axis in z-direction through the center of mass of the ion channel in the center of the simulation box. In this rotational alignment, the rotation angles were chosen so that the x-y coordinates of the centres of mass of the three helices 38 forming the channel in the underlying simulation conformations were aligned. To generate average shapes that reflect the 3-fold symmetry of the Piezo proteins, we performed additional rotations by 120 and 240 degrees of all rotationally aligned membrane shapes and included the two additional resulting shapes per membrane shape in the averaging. In our figures and calculations, the average membrane shapes have a circular projected area with a diameter d_m that corresponds to the average x and y dimensions of the simulation boxes of the conformations, because only these ‘inscribed’ circular membrane segments are overlying in all rotated conformations. In an analogous way, we calculated average, 3-fold symmetric lipid densities from the interpolated lipid numbers in the bins of the square lattice of each conformation.

Excess area and bending energy

To numerically determine the excess area and bending energy of average membrane midplane shapes, we discretized the continuous average shapes and lipid densities obtained after the interpolation, rotation, and averaging steps described above by discretizing the reference x-y plane of the Monge parametrization into a square lattice with lattice constant a = 1 nm. The excess area at site (x, y) was then calculated as

where z(x, y) denotes the midplane z coordinate at site (x, y), and the overall excess area ΔA of the average shape was obtained by summing up the excess areas δA(x, y) at all sites of the shape, i.e. at all sites with x² + y² ≤ (d_m/2)². We determined the bending energy at site (x, y) based on the general expression for the mean curvature H in Monge parametrization as

with bending rigidity κ and the standard discretizations of the partial derivatives

The overall bending energy of the average membrane midplane shape of the membrane in the protein-membrane nanodome was then calculated as

where ρ(x, y) ≤ 1 is the relative average lipid density at site (x, y), calculated as the number of lipid head beads at site (x, y) divided by the mean number of lipid head beads for ‘pure’ membrane patches far away from the protein.

TM helices

In determining the average COMs of the 4-TM units in Figures 1(d) and 2(e), we discarded residues of TM helix ends that protrude out of the membrane. To identify these residues for the Piezo 1 and Piezo 2 TM helices, we calculated the average contact probability of all TM helix residues with lipid tails in the final two microseconds of the 10 coarse-grained trajectories with tensionless membranes. A TM helix residue was taken to be in contact with lipids tails in a simulation frame if the smallest distance between the beads of this residue to lipid tail beads was less than 0.6 nm. We discarded consecutive stretches of residues at the helix ends with an average contact probability of less than 10% as helix ends that protrude out of the membrane.

Acknowledgements

This research has been funded by the Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114/A4 and by the Max Planck Society.