Microscopic mechanism of PIEZO1 activation by pressure-induced membrane stretch

PIEZO1, one of only two vertebrate members of a small eukaryotic family of cation-selective mechanosensitive ion channels (Coste et al., 2012), plays essential roles across physiological systems and participate to the progression of numerous diseases including hereditary dehydrated stomatocytosis (xerocytosis), lymphoedema, hypertension, cancer, osteoporosis, and inflammation (Qin et al., 2021; Ranade et al., 2014; Song et al., 2022; Syeda, 2021). It is now well established that PIEZO1 senses mechanical forces directly transmitted by physical deformations of the cell membrane, such as increased membrane tension which results in elastic stretching of the lipid bilayer (Cox et al., 2016; Lewis and Grandl, 2015; Syeda et al., 2016). Mechanical activation of PIEZO1 can be quantitatively evaluated using the pressure-clamp technique, in which the tension of the membrane in the patch increases upon negative or positive pressurization of the path pipette according to the Young-Laplace law (Lewis and Grandl, 2015). When acutely evoked with a mechanical stimulus, including pipette pressure, the macroscopic relaxation currents produced by a population of PIEZO1 channels rise rapidly and decay within tens of milliseconds through one or more inactivation processes, suggesting the existence of multiple resting, open, and inactivated states (Bae et al., 2013; Lewis et al., 2017; Zheng et al., 2019a).

PIEZO structures solved in detergent micelles consist of homotrimers with a non-conducting pore located in the central region and three peripheral non-coplanar blades, giving the channel an upright bowl-like conformation. This non-planar structure is thought to correspond to a resting state populated in absence of mechanical stimulus (Ge et al., 2015; Saotome et al., 2018; Wang et al., 2019; Zhao et al., 2018). Interestingly, because of the non-zero bending rigidity of lipid bilayers, PIEZO1’s non-planar conformation is predicted to curve the lipid bilayer beyond the protein boundary, suggesting that PIEZO1’s conformation can be controlled by membrane tension through elastic membrane-protein interactions (Haselwandter and MacKinnon, 2018; Haselwandter et al., 2022a, 2022b). Cryo-EM images of PIEZO1 reconstituted into proteoliposomes and molecular dynamics simulations indeed show that the PIEZO1 structure can be made more curved or flatter by manipulating the curvature of the surrounding membrane (Jiang et al., 2021; Lin et al., 2019; Yang et al., 2022). Recently, a flattened PIEZO1 cryo-EM structure was captured by imposing a large curvature mismatch between PIEZO1 and the liposomal membrane (Yang et al., 2022). Although this flattened conformation could not be unambiguously assigned to a particular functional state due to a lower resolution in the pore region, molecular dynamics simulations show that flattening of PIEZO1 correlates with an open state (De Vecchis et al., 2021; Jiang et al., 2021).

In spite of these efforts, little is known about how membrane stretch energetically modulates PIEZO1’s functional states at the microscopic level. Here, we measured pressure-induced modulations of PIEZO1 single channel dwell time using 1-D and 2-D dwell time analysis and Markov-chain modeling.

Cell culture, reagents, and transfection methods described in our previous work are used in this project with slight modifications (Wijerathne et al., 2022). PIEZO1 knock out HEK293T^ΔPZ1 cells, gifted by Dr. Ardèm Patapoutian (Scripps Research, La Jolla, CA, USA), were cultured at 37°C and 5% CO2 in Dulbecco’s Modified Eagle Medium (Gibco, Thermo Fisher Scientific) supplemented with penicillin-streptomycin and 10% fetal calf serum (Life Technologies, Thermo Fisher Scientific). Cells were plated into 35-mm cell culture dishes 48 h prior to experiments, and transiently transfected with 1 µg pCDNA3.1-mPIEZO1-IRES-mCherry plasmid with Lipofectamine 3000 Transfection Reagent (Life Technologies) using the manufacturer’s recommended instructions 12–24 h before the experiments. Cells were seeded onto 9-mm round coverslips 3 h before experiments and used in patch clamp recordings 12–24 h after transfection. For MscL measurements, an MscL-IRES-GFP plasmid was created and generously gifted by David Clapham (Janelia Research Campus, Ashburn, VA, USA; Doerner et al., 2012) and a GFP plasmid was used as negative control.

Cells were patched using unpolished 5–6 MΩ pipettes to increase reproducibility of pipette tip geometry. Negative pressure was delivered to the backside of patch pipettes using a Clampex-controlled high-speed pressure clamp (ALA Scientific Instruments). Pipettes were filled with 140 mM KCl, 10 HEPES, 10 mM TEA, and 2 mM EGTA (pH 7.4 with NaOH). HBSS with calcium, magnesium, and sodium bicarbonate (#14025; Gibco) was used as a bath solution. Currents were acquired using an Axopatch 200B amplifier and an analog-digital converter Digidata 1550B (both from Molecular Devices) as individual 1-min continuous gap-free recordings using a sampling rate of 500 kHz and a filtering frequency of 5 kHz. Data were subsequently filtered at 500 Hz for display.

All-point current histograms were generated using Clampfit software (Molecular Devices) and fitted with OriginPro2018’s Multiple Peak Fit tool. Single channel time series were idealized using Clampfit using a dead time of 10 µs. A MATLAB (MathWorks) script posted to Github (https://github.com/LacroixLaboratory) was used to ignore events with duration shorter than 0.05 ms and with change of current amplitude <25% of the preceding event. Adjacent open-open and shut-shut events produced by this correction were then combined by running the corrected time series one more time with the script.

The sequence of alternating open/shut events was scrambled using a home-made protocol in OriginPro (2018). Briefly, the sequence of open/shut intervals (taken from idealized traces) was split in two sets: one having all open intervals, the other having all shut intervals. In both sets, each interval was associated with a uniform random number. Each set of open and shut intervals was then scrambled by re-ordering them by ascending random numbers. Next, these random numbers were erased and replaced with continuous odd positive integers for open intervals starting from 1 (1, 3, 5, 7, …) and with continuous even positive integers for shut intervals starting from 2 (2, 4, 6, 8, …). All open and shut intervals were then recombined and sorted by ascending integer. This operation effectively scrambles the sequence of open/shut events while keeping the strict alternance between open and shut intervals.

2-D dwell time histograms, used to quantify the likelihood of a shut event of length a following an open event of length b (and vice versa), were estimated using the krn software, available to download from the online supplementary material of an earlier paper (Rosales et al., 2004). Two open states and three shut states were defined for kernel analysis, corresponding to the two and three distinct exponential distributions observed in both open and shut 1-D dwell time histograms. A sampling period of 0.01 ms was used for all calculations. The number of total evaluations was set to 100. Prior to data entry, traces were idealized into general open and shut event states via MATLAB as described above. Dependency plots were obtained by first calculating the dwell time distribution assuming random transitions (i.e., no preference for any shut state to selectively transition into an open state of a specific duration, and vice versa). This value was then subtracted from the observed distribution to produce ∆-dependency plots. Output data sets were imported to OriginPro 2018 for plotting. The x and y axes, generated in natural log scale by krn software, were log-converted prior to plotting.

The maximum interval likelihood function of the online software suite Qub was used for the determination of rate constants between O1, O2, S1, S2, and S3 states. Prior to data entry, input traces were idealized by MATLAB as described above. These idealized current traces were converted to .DWT files and uploaded to Qub. For each round of rate calculations, 10⁵ iterations were run with a step size of 10⁻⁶ s⁻¹ and 10⁵ steps per iteration. To achieve model convergence and avoid local minima, initial rate constants were first determined by simultaneously fitting all current traces from −40 mmHg experiments using an unconstrained model that allow all interstate transitions with a nominal rate constant of 1 s⁻¹ for all transitions. If convergence was not achieved, both step size and iterations were changed. Upon successful convergence the data were saved and model was reset to initial rates of 1 s⁻¹ for all transitions and rerun until convergence. Rates calculated by the unconstrained model with best Bayesian Information Criterion (BIC) score was used as the initial rates for the BKU model (Bauer et al., 1987; Kienker, 1989). The BKU model connects every open state to every shut state and vice versa. Upon successful convergence, the data was saved, and rates were reset to initial rates from the unconstrained model and the model was reoptimized. The rates from this BKU model were used as the initial rates in all other models. When a new model required a silent transition which are absent in BKU model, the transition was added with initial rates set to 1 s⁻¹. Detailed balance was enabled during Markov modeling.

Data simulation was done by online software suite Qub. 240 s long single-channel traces were simulated at 100 kHz sampling rate. Output traces were downloaded as idealized dwelt time files (.DWT).

Log Likelihood, BIC, and AIC scores were computed by the in-built functions of Qub and Clampfit softwares. We used “Delta BIC+1” to compare the models fitted from the same data set. Delta BIC was obtained by subtracting the BIC from the current model from the lowest BIC obtained for the same data set. All the delta BIC values were added a +1 to bring the minimal BIC values to one instead of zero, facilitating BIC comparison using logarithmic scale.

Cell-attached pressure-clamp single channel recordings represent an ideal approach to test how membrane stretch modulates discrete open and shut states of mechanosensitive channels in the physiological context of a living cell (McBride et al., 1992). Here, we recorded single channel currents from PIEZO1-deficent HEK293T^ΔPZ1 cells (Lukacs et al., 2015) transiently transfected with a plasmid encoding mouse PIEZO1, as described earlier by us and others (Nosyreva et al., 2020; Wijerathne et al., 2022). To increase the chance of single vs. multiple channel capture, we used high-resistance patch pipettes (5–6 MΩ, ∼0.8 µm tip diameter) to record from a small membrane area and patched cells <24 h after transfection.

We acquired 60-s-long single channel current traces at a steady-state negative pipette pressure ranging from 0 to 90 mmHg (relative to ambient air pressure; Fig. 1 A). We limited our most negative pressure to −90 mmHg because more negative pressures tend to cause premature patch rupture. We occasionally applied a short −90 mmHg pressure pulse at the end of low-pressure recordings to promote simultaneous opening of multiple channels in case more than one channels were trapped in the patch. Most recordings displayed either no current or discrete currents consistent with non-overlapping single channel openings (recordings with overlapping openings were eliminated from further analysis). Because PIEZO1 is cation-selective (Coste et al., 2010), its reversal potential is near 0 mV in physiological solutions. A large patch voltage was thus necessary to resolve single channel current above noise. To mimic the negative membrane environment, we used a constant voltage of −90 mV in all our patches.

The durations of conducting (open) events and non-conducting (shut) events were extracted from raw current trajectories displaying non-overlapping single channel opening events by performing a standard time series idealization using Clampfit, followed by a homemade event correction script to eliminate transitions with amplitude change <25% and duration <0.05 ms (see Materials and methods; Wijerathne et al., 2022). We calculated PIEZO1’s apparent open probability (NPo) for each pressure using (1) the area under the curve from all-point histograms, (2) the sums of individual open and shut durations, and (3) the mean duration of open and shut dwell times (see Materials and methods). Scatter plots of Npo values against pressures were well fitted by a two-state Boltzmann function, with half-maximal pressure values (P₅₀) ranging from −51.6 ± 3.9 to −54.4 ± 5.7 mmHg, giving an error-propagated average of −52.9 ± 2.8 mmHg (Fig. 1, B–D). Plotting the mean duration of open and shut events as a function of pressure reveals that pressure-induced increase of PIEZO1 open probability correlates with a large (approximately >100-fold) decrease of the mean shut duration and with a comparatively smaller (approximately <10-fold) increase of the mean open duration (Fig. 1 E).

Because the stimulus that PIEZO1 senses are the membrane stretch produced by the transmembrane pressure gradient, we sought to estimate the half-maximal tension (T₅₀) that corresponds to the observed P₅₀ values. Although a direct determination of membrane tension is possible by optically measuring the patch dome radius (Cox et al., 2016; Lewis and Grandl, 2015), this approach was not possible due to the narrowness of our pipettes. Instead, we measured the macroscopic pressure-activation curve of MscL, a bacterial mechanosensitive ion channel with a known T₅₀ value of ∼10.4 to ∼11.8 mN m⁻¹ (Chiang et al., 2004; Sukharev et al., 1999). When measured with our pipettes, MscL has a P₅₀ value of −172 ± 6 mmHg (Fig. 2, A and B). Assuming tension linearly scales with pressure for identical cells and pipette geometry (Sukharev et al., 1999), the average PIEZO1 single channel P₅₀ value would yield a T₅₀ value of 3.4 ± 0.3 mN m⁻¹, which falls within the range, but near the upper limit, of T₅₀ values observed from non-equilibrium macroscopic relaxation recordings (∼1.4 to ∼4.5 mN m⁻¹; Cox et al., 2016; Lewis and Grandl, 2015).

The relative distribution of discrete conformational states is exponentially related to their relative energies (Boltzmann, 1868). Hence, open and shut dwell times extracted from single channel trajectories tend to be exponentially distributed, with one or more exponential components depending on the number of discrete functional states and on the complexity of conformational routes connecting them.

We plotted open and shut dwell time histograms using the log-transformed Sigworth et al. representation (Sigworth and Sine, 1987) and determined the optimal number of exponential components from 1 to 6 by comparing the values of the BIC obtained from multi-exponential fitting (see Materials and methods). Fig. 3 A shows that, regardless of pressure, 2–4 and 3–4 exponential components yield approximately similar BIC values for open and shut dwell time histograms, respectively. However, a careful examination of the fitted histograms reveals that adding a third and a fourth exponential terms for open and shut histograms, respectively, leads to slow components with negligible areas for most pressures and with time constant (Taus) that have large errors or values beyond the length of observation time (Fig. 3, B–D). This shows that adding these additional terms did not improve fitting. Hence, open and shut dwell time histograms are best described with sums of two and three exponentials components, which we, respectively, call open components 1 (fast, OC1) and 2 (slow, OC2), and shut components 1 (fast, SC1), 2 (intermediate, SC2), and 3 (slow, SC3).

Increasing the amplitude of the pressure stimulus correlates with an ∼5- and ∼8-fold increase of OC1 and OC2 time constants, respectively (two terms plot in Fig. 3 C), and with an ∼3-, ∼4-, and ∼15-fold decrease of SC1, SC2, and SC3 time constants, respectively (three terms plot in Fig. 3 D). These reciprocal changes of time constants are within the same orders of magnitude between open and shut components. These modulations thus seem unable to fully explain the strongly asymmetric effect imparted by pressure on shut and open mean duration (Fig. 1 E). Interestingly, while increasing pipette suction does not significantly affect the relative amplitude of OC1 (>80%) and OC2 (<20%) components to open dwell times (Fig. 3 E), it fully reversed the relative amplitudes of SC1 (from ∼10 to ∼60%) and SC3 (from ∼65 to ∼10%) components to shut dwell times (Fig. 3 F). This amplitude reversal explains why the effect of pressure is more pronounced on the mean shut dwell times: at pressures near 0 mmHg, the majority of shut events belong to the slowest component SC3 and thus have a long duration, whereas at more negative pressures, the majority of shut events belong to the fastest component SC1 and thus have a short duration. When fitted using a two-state Boltzmann function, the pressure-dependent open probability plot derived from the amplitude-weighted time constants of exponential components yields a P₅₀ value of −62.6 ± 6.1 mmHg (Fig. 3 G), slightly more negative than the P₅₀ calculated using fit-independent methods (Fig. 1, B–D).

A crucial information tractable from single channel data is the correlation between adjacent open and shut intervals (Magleby and Weiss, 1990; Rosales et al., 2004; Rothberg and Magleby, 1998). This correlation can be visualized by plotting the relative frequency of paired adjacent open/shut events as a function of their individual dwell times. The top panels in Fig. 4 A show the corresponding 2-D dwell time histograms for each pressure. These plots show that, near 0 mmHg, most paired events have a long shut duration and a short open duration. As pressure increases, the dominant correlation shifts toward pairs having shorter shut and longer open durations. Several clear patterns seem to emerge at all pressures: (1) Long open intervals are strongly associated with short shut intervals; (2) long shut intervals are strongly associated with short open intervals; and (3) long open intervals are not associated with long shut intervals.

Because these correlations are necessarily influenced by the relative frequency of individual open and shut intervals, 2-D dwell time plots should also be normalized by subtracting the observed 2-D distribution by the distribution that would be observed if open and shut events were randomly paired based on their relative abundance. When differences exist between these two distributions, positive correlations show up as positive (dark gray) and negative (dark red) peaks in 2-D dependency plots (bottom panels in Fig. 4 A; Magleby and Song, 1992; Rosales et al., 2004; Rothberg and Magleby, 1998). Across most pressures, these dependency plots confirm the three patterns seen in non-normalized 2-D histograms and further suggest that short shut intervals do not statistically pair with short open intervals.

To test the meaningfulness of these correlations, we scrambled the sequence of alternating open/shut events from our −40 mmHg data (refer to Materials and methods). This modification had no effect on 1-D dwell time histograms, as expected since 1-D histograms are only sensitive to the distribution of intervals, not to the temporal sequence in which they occur during the course of the recordings (Fig. 4 B). By contrast, scrambling this sequence eliminates the characteristic L-shape from the 2-D dwell time plot at −40 mmHg, and completely eliminates the four correlations peaks seen in almost all dependency plots (Fig. 4 C). This shows that these correlations are highly specific to the observed temporal sequence of open/shut events. Overall, our 2-D dwell time analysis suggests that, at low pressures, open/shut transitions mostly occur between long shut and short open events, whereas at more negative pressures, they mostly occur between shorter shut and longer open events. This information will be used in a subsequent paragraph to help narrow down the number of kinetic topologies that can realistically recapitulate PIEZO1’s single channel behavior.

To generate 1-D dwell time distributions with two open and three shut exponential components with significant amplitude, PIEZO1 must populate a minimum of two open (O1 and O2) and three shut (S1, S2, and S3) microscopic states (Colquhoun and Hawkes, 1981). To determine how the channel transits between that many minimal number of states, we used the Markov modeling software Qub (Milescu et al., 2005; Wijerathne et al., 2022) to fit our −40 mmHg single channel data using five canonical models (M1–M5) with distinct topologies: one authorizing all possible transitions (M1), one having every possible O-S links but no silent O-O or S-S link (M2, also called Bauer-Kienker Uncoupled [BKU] model; Bauer et al., 1987; Kienker, 1989), a box model (M3), and two minimal models (i.e., models in which any two states are connected by no more than one path) with a single gateway state (i.e., with only one independent O-S link; M4 and M5; Fig. 5 A).

Comparison of optimal BIC values obtained by fitting our data to each canonical model (see Materials and methods) shows that our −40 mmHg data is best fitted with either M1, M2, or M3, whereas M4 and M5 yield relatively poor fits (Fig. 5 B and Table 1). All models predict visually similar 1-D dwell time distributions (Fig. 5 C). However, whereas single channel data simulated with M1, M2, or M3 recapitulates the L-shaped 2-D dwell time plot and the dependency pattern seen in unmodified experimental data (Fig. 4 A and Fig. 5 D), data stimulated with M4 or M5 fail to reproduce the L-shaped 2-D dwell time plot and yield “flat” dependency plots (Fig. 5 D). This is due to the fact that both M4 and M5 harbor a single gateway state: in these models, adjacent open and shut durations necessarily pair independently of their durations (Rothberg et al., 1997).

A graphical enumeration theorem demonstrates that there are precisely 728 ways to uniquely connect five states with bidirectional links (Harary and Palmer, 1973). Our initial Markov modeling results suggests that a minimal model with only four links could be sufficient to fit our data, as models with five (M3) or six (M2) links converged to models with four dominant links; the remaining ones having negligible rate constants (red links in Fig. 5 A). Cayley’s formula shows that the number of minimal models connecting n states is equal to nⁿ⁻² (Harary and Palmer, 1973). This means that there are 125 minimal five-states models, which is considerably less than 728 but still fairly large to test one by one without an automatic approach (Sukharev et al., 1999).

Our 2-D dwell time analysis and dependency plots suggest that long open events pair with short shut events and that short open events pair with long shut events (Fig. 4). Furthermore, the 2-D dwell time analysis suggests the existence of more than one gateway state, implying that a realistic kinetic model should have more than one independent open-shut link (Rothberg et al., 1997). Based on these observations, we made three simple assumptions for our kinetic models: (1) O2 is connected to S1; (2) O1 is connected to S2 and/or S3; and (3) O1 and O2 must be connected to different shut states. Although the dependency plots suggest that short open events do not pair with short shut events, we decided to include models with the O1-S1 link, as short open/shut intervals are among the most frequently observed intervals across all tested pressures. These constraints led us to 3 families of minimal models, encompassing a total of 18 new possible models, named M6 to M23 (Fig. 6, A and B). Some models have silent open-open and shut-shut links, thus allowing long open or shut events to arise from sojourns into more than one state.

In the Qub software, input states are only defined by their open/shut class, not by their duration. Hence, from the perspective of Qub, these 18 models correspond to only 8 distinct input topologies (Fig. 6 C). We fitted our −40 mmHg data with each of these input topologies multiple times, using initial rates obtained from the converged BKU model (see Materials and methods). For each of eight possible input topologies, the fits consistently converged to a unique model (Fig. 6 D). All these unique models recaptured the exponential distributions of 1-D dwell time histograms (Fig. 6 E), but M21 had a much better (lower) BIC score compared to other models (Fig. 6 F).

We fitted data obtained at different pressures with M21 (Fig. 7 A) using the fitted rates obtained at −40 mmHg as initial rates. The converged models recapitulate the exponential components of 1-D dwell time distributions (Fig. 7 B). Because M21 authorizes silent transitions between shut states (S3-S2 and S2-S1), the mean dwell times of shut states predicted by M21 may be shorter than the time constant of exponential components identified from 1-D dwell time fitting. This is because silent microscopic transitions between shut states increase the apparent duration of observed shut events. This seems to be the case for the highly connected S2 state, whose dwell time is consistently shorter than the time constant of the intermediate component SC2 at pressures between −30 and −90 mmHg (Fig. 7 C). In contrast, the dwell times of O1 and O2 states are very close to the time constants of the OC1 and OC2 components at all pressures, as expected since no silent transition connecting O1 to O2 exits in M21 (Fig. 7 D).

Increasing pipette suction imparts profound effects on most rate constants in the M21 scheme. One of such effect is to increase the two opening rates S1→O2 and S2→O1, the largest effect being an ∼40-fold increase of the S2→O1 opening rate (Fig. 7 E). These changes overall decrease shut dwell times and increase the frequency of transitions into open states. The two closure rates O2→S1 and O1→S2 also decrease, thus prolonging the duration of open intervals (Fig. 7 F). Many of the rates between silent shut-shut transitions are also affected by pressure. As pressure decreases, the rates of S3→S2 and S2→S3 transitions increase ∼10-fold, while the rate of S2→S1 transitions increases by >20-fold (Fig. 7 G). These changes result in a shift of the proportion of shut events toward the S1 population, thus explaining the reversal of relative amplitudes between the SC1 and SC3 components (Fig. 3 F).

We next simulated single channel data rates obtained from M21 fitting. The Po values calculated from the sum of occupancy of O1 and O2 states display a Boltzmann curve when plotted as a function of pressure, yielding a P₅₀ nearly identical to the P₅₀ calculated from the amplitude-weighted time constants of exponential components of dwell time distributions (−62.98 ± 4.09 vs. −62.6 ± 6.1 mmHg; Fig. 3 G and Fig. 7 H). Regardless of pressure, the channel spends most of the time into the long shut S3 state, whose relative occupancy only decreases from ∼0.73 at 0 mmHg to ∼0.5 at −90 mmHg. This seems at odds with the fact that only a minority of shut events are long at negative pressures (Fig. 3 F). This conflict is resolved by realizing that the dwell times of the S3 state is at least one order of magnitude longer than any other state (Fig. 7, C and D). In addition, single channel data simulated with M21 yields 2-D dwell times distributions and dependency plots (Fig. 8) comparable to those obtained from model-independent analysis (Fig. 4). Overall, modeling our data with M21 suggests that pipette suction increases open probability by (1) increasing opening rates, (2) decreasing closure rates, and (3) shifting the distribution of shut intervals such that the relative frequency of short shut intervals increases while the relative frequency of long shut intervals decreases.

At constant physiological temperatures and under true equilibrium conditions (i.e., with no entropy production), thermal agitation is the only source of energy enabling ion channel molecules to continuously transit between open and shut states, generating a seemingly endless sequence of open/shut transitions. Because this energy cannot be converted into organized molecular motions, microscopic transitions are purely stochastic. Agreement with these conditions implies that, for sufficiently large samples, the number of forward and reverse microscopic transitions between any pair of connected states should be the same, producing time reversibility in the single channel record.

Microscopic reversibility can, however, be violated in the single channel record, for instance if an inconspicuous source of free energy influences channel gating (Richard and Miller, 1990), or when data is not obtained in true equilibrium conditions. We used two known approaches to test whether our PIEZO1 single channel data obeys microscopic reversibility. The first approach uses a χ² test to evaluate statistical differences between 2-D dwell time distributions obtained using single channel traces in either chronological or time-reversed order (Song and Magleby, 1994). Table 2 shows that, for all tested pressures, these differences are not significant at the 5% level (Z < 1.95; Snedecor and Cochran, 1989). The second method uses the detailed balance feature of the Qub software which was used by default during Markov modeling. To check if the data are truly time-reversible, we tested whether the Log Likelihood value obtained by fitting our data with M21 would improve without detailed balance. Table 2 shows that, at all tested pressures, these values remain unchanged. We thus concluded that, in steady-state conditions, PIEZO1 obeys microscopic reversibility, suggesting that thermal agitation acts as sole source of energy to stochastically move the channel between open and shut conformations.

In this study, we showed that pipette suction increases PIEZO1 open probability through an asymmetric modulation of open and shut dwell times. Although the exponential time constants describing the distribution of both open and shut interval vary in opposite direction and with similar magnitude as a function of pressure, shut dwell times shift from a distribution dominated by long shut events near zero pressure to a distribution dominated by short shut events at more negative pressures. Overall, these modulations impart a large reduction of the mean shut duration and a relatively small increase of the mean open duration.

Other mechanosensitive channels, including TWIK-related K⁺ channels type 2 (TREK-2) and MscL, also respond to stretch predominantly by reducing mean shut dwell times (Clausen et al., 2017; Sackin, 1989; Sukharev et al., 1999). This similarity seems to contrast with the bewildering diversity of stretch-activated channels, which have been identified in evolutionary-distant protein families with no significant homology in primary amino acid sequences or tri-dimensional structures (Kefauver et al., 2020). In addition, these channels seem to use widely different molecular mechanisms to sense stretch. For instance, non-planar PIEZOs are thought to sense tension through local changes in membrane curvature (Haselwandter et al., 2022a; Jiang et al., 2021; Yang et al., 2022), whereas MscL has been proposed to sense lateral tension due to hydrophobic mismatch between transmembrane helices and the lipid bilayer, opening the pore by tilting these helices in an iris-like fashion (Perozo et al., 2002; Wang et al., 2014). MscS seemingly employs allosteric lipid-protein interactions to couple lipid stretch to pore opening (Reddy et al., 2019), whereas TRAAK (TWIK-related arachidonic acid activated K⁺) and TREK channels are thought to sense tension through minute rearrangements of an extended transmembrane helix (Brohawn et al., 2014; Lolicato et al., 2014).

Our minimal model M21 recapitulates all important results obtained from model-free analysis of single channel data. This model was, however, selected based on model complexity-corrected maximum likelihood (BIC) comparison from data obtained at a single pressure of −40 mmHg. We thus cannot exclude the possibility that another minimal model could better fit our data at other pressures. In addition, both M21 and the more complex canonical BKU model (M2) yield identical BIC values, showing that both models are equally likely to describe our data. The presence of silent transitions between shut states (S3-S2 and S2-S1) in the M21 model may, however, confer M21 an advantage over the BKU model: these transitions could provide a microscopic kinetic basis for the recovery from inactivation observed from macroscopic relaxation recordings. Indeed, this kinetic process is predicted to be driven by silent transitions from inactivated (unable to open) to resting (ready to open) macroscopic states (Lewis and Grandl, 2020; Wijerathne et al., 2022). However, before this microscopic basis can be established, one would first need to determine whether some of the shut states identified from the single channel data can be specifically assigned to an inactivated state. This could be done in future studies, in principle, by conducting single channel analysis in non-activating mutants, which are available for PIEZO1 (Bae et al., 2013; Zheng et al., 2019a).

An important similarity between M21 and the BKU models is the presence of two gateway states, consistent with our independent observation of correlations between adjacent open and shut intervals. We recently used a simpler, linear four-states model to study the effects of small molecules on PIEZO1 (Wijerathne et al., 2022). Consistent with the present study, this linear model also includes two gateway states. Both linear and nonlinear (branched) models encompassing multiple gateway states have also been proposed for other mechanosensitive channels, including TRAAK (Sorum et al., 2021) and stretch-activated channels endogenously expressed in Xenopus laevis oocytes (Gil et al., 2001). In contrast, linear models with a single gateway state have been proposed for MscL (Sukharev et al., 1999) and for stretch-activated channels present in embryonic chick skeletal muscles (Guharay and Sachs, 1984), in Necturus proximal tubules (Sackin, 1989), and in molluscan heart cells (Sigurdson et al., 1987). Because Markov models are not always identifiable (Fredkin et al., 1985) or distinguishable based on single channel data (Kienker, 1989), comparing kinetic models proposed for these different channels remains difficult. In this context, 2-D dwell time analysis enables rapid identification of dominant microscopic transitions and informs about the number of gateway states, which should help validate or invalidate kinetic schemes independently of Markov modeling (Gil et al., 2001; Magleby and Weiss, 1990).

An important fundamental question is to know how the microscopic transitions between multiple open and shut states relate to structural changes in the pore and in the putative mechanosensory blade domains. A pore-only PIEZO1 construct lacking the entire blade domain undergoes steady-state open/shut transitions and yields the same number of exponential dwell time components as full-length channels (Nosyreva et al., 2020). These observations suggest that the pore domain itself populates functional states that are normally visited by wild-type PIEZO1. They also suggest that pressure-induced conformational changes, predicted to be initiated in the blade domain, could control the pore by energetically influencing the distribution of open and shut intervals rather than by switching the pore between open and closed conformations in a strictly binary fashion. Experiments aiming at correlating the functional dynamic of the pore with the structural dynamic of the blade would be interesting to pursue in this regard. To this aim, we recently used site-specific macroscopic fluorimetry measurements to establish temporal correlations between blade motions and pore opening (Ozkan et al., 2022 Preprint). Single molecule fluorescence techniques should help establish this correlation at the microscopic level in future studies.

PIEZO1 senses an unusually wide range of mechanical stimuli, such as shear stress, cellular indentations, and membrane stretch, and detects forces transmitted from both lipids (Cox et al., 2016; Lewis and Grandl, 2015; Syeda et al., 2016) and tethered extracellular and/or intracellular filaments (Ellefsen et al., 2019; Gottlieb et al., 2012; Mylvaganam et al., 2022; Wang et al., 2022). In addition, PIEZO1 seems to possess distinct mechanisms to sense distinct mechanical forces (Ozkan et al., 2022,Preprint; Verkest et al., 2022; Wang et al., 2018). It is thus possible that the activatory effects of distinct mechanical stimuli could be additive, at least partially. If that is the case, co-stimulating PIEZO1 with multiple mechanical stimuli might increase open probability above the saturation limit observed here using pipette suction (∼0.25 at a V = −90 mV). This hypothetical additivity means that the channel could potentially use a larger range of open probability in certain physiological contexts in which multiple forms of mechanical stress are present, potentially making PIEZO1 an effective multimodal integrator of mechanical forces. Future studies will be needed to test this hypothesis.

Beside mechanical forces, many other parameters modulate PIEZO1 activity. These include temperature (Zheng et al., 2019b), the membrane potential (Moroni et al., 2018; Nosyreva et al., 2020), small molecules (Syeda et al., 2015; Wang et al., 2018; Wijerathne et al., 2022), extracellular pH (Bae et al., 2015), and dietary lipids (Romero et al., 2019). Why PIEZO1 responds to such non-mechanical stimuli is not fully understood. These modulations could exist as mere side effects of the complex channel structure. They could also arise as bona fide mechano-modulatory functions selected during evolution to provide a selective physiological advantage. Future studies might shed light on this interesting aspect of PIEZO1 channel regulation.

Crina M. Nimigean served as editor.

This work was supported by National Institutes of Health grant GM130834 to J.J. Lacroix.

The authors declare no competing financial interests.

Author contribution; T.D. Wijerathne: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing - Review & Editing, Visualization; J.J. Lacroix: Conceptualization, Methodology, Writing - Original Draft, Review & Editing, Supervision, Project administration, Funding acquisition; A.D. Ozkan: Methodology, Software, Validation, Formal analysis.