Ion channels in critical membranes: clustering, cooperativity, and memory effects

Much progress has been made in elucidating the inner workings of voltage-gated ion channels, but less understood is the influence of lipid rafts on gating kinetics. Here we propose that state-dependent channel affinity for different lipid species provides a unified explanation for the experimentally observed behaviors of clustering, cooperativity, and hysteresis. We develop models of diffusing lipids and channels engaged in Ising-like interactions to investigate the collective behaviors driven by raft formation in critical membranes close to the demixing transition. The model channels demonstrate lipid-mediated long-range interactions, activation curve steepening, and long-term memory in ionic currents. These behaviors likely play a role in channel-mediated cellular signaling and suggest a universal mechanism for self-organization of biomolecular assemblies.

Ion channels play a critical and ubiquitous role in cellular signaling by transmuting external forces into changes in ion permeability across membranes, a process known as gating.The discovery of on/off ionic currents catalyzed by a single pore [1,2] led to a reinterpretation of the seminal Hodgkin and Huxley equations [3] as stochastic gating by independent channels [4].As electrophysiological and biochemical techniques improved, gating schemes expanded to accommodate new data [5,6], but the assumption that channels gate independently has persisted.
The paradigm of the independent channel has been challenged by experiments that directly visualize the cellular membrane, such as electron microscopy, confocal fluorescence microscopy, and superresolution imaging [7][8][9][10].Channels are seen to aggregate in clusters ranging in size from a few channels to large groups containing tens of thousands of channels.Factors that contribute to clustering include direct channel-channel interactions [11] and connections to the cytoskeleton [12].
Here we consider another mechanism of clustering that is mediated by the lipid membrane.The fluid-like dynamics of a heterogeneous population of lipids can indirectly mediate channel interactions.The notion that conformational changes in the channel can affect and be affected by the thermodynamic properties of the membrane is supported by a range of evidence from cryoEM, electrophysiology, and computational studies [13][14][15][16][17][18][19][20][21][22][23].Lipid membranes under physiological conditions are close to the demixing transition, resulting in the formation of rafts and liquid disordered domains [24][25][26][27].Near the critical temperature, membranes are subject to fluctuations in local composition characterized by large correlation lengths [28].
In this work, we propose a unified explanation for the size distribution of clusters and gating anomalies based on the assumption that distinct conformational states of the channel (including open and closed) have affinities for different lipids species.To this end, we introduce two complementary phenomenological descriptions of ion channels embedded in a lipid mixture: i) a mesoscopic model of lipid membranes, which faithfully describes the membrane phase behavior and its coupling to the channel conformation; and ii) a discrete-state Markov model of interacting channel-lipid systems with a realistic description of channel gating kinetics.The first model is a minimalist representation of lipids and channels whose purpose is to demonstrate that these collective phenomena are independent of any specific molecular detail, but rather derive from the combination of state-dependent interactions and critical demixing.The second model builds on a more realistic representation of the channel structure and activation dynamics and allows us to ascertain that the effects investigated here occur on physiologically relevant length-and time-scales.In both models, we find that fluctuations in the underlying lipid medium give rise to long-range effective attractive forces leading to clustering, cooperative gating, hysteresis, and longmemory effects in the channel gating dynamics.

A. A mesoscopic model of ion channels in phase-separating membranes
To study the dynamics of channels embedded in a lipid bilayer, we introduce a simplified mesoscale representation of the membrane and of the embedded ion channels.
The membrane is described by an assembly of arXiv:2401.07660v2[cond-mat.soft]22 Jan 2024 molecules of two types, each representing saturated or unsaturated lipids.The assembly is meant to mimic the phase behavior of a three component membrane in physiological conditions [44,45].In cholesterol containing ternary mixtures with saturated and unsaturated lipids, saturated lipids tend to segregate into liquid ordered (L o ) domains, while unsaturated lipids promote the formation of a liquid disordered phase (L d ).L o domains are characterized by hexatic order (see [17,27,46,47] and SM for a definition) and a larger membrane thickness compared to the L d phase, see Fig. 1b.At fixed temperature, the transition from a well-mixed liquid disordered membrane to a phase-separated coexistence of L o and L d domains is driven by a change in lipid composition (i.e. a change in the cholesterol and lipids molar fractions).Importantly, at fixed composition, this mixing/demixing transition can be driven by a change in temperature: higher temperatures promote mixing of the two lipids into a liquid disordered (L d ) phase.
We model each ion channel as a multi-domain inclusion in the membrane in which each domain is allowed to fluctuate across a single-barrier conformational landscape.The stable basins are meant to recapitulate the two-state character of ion channels that cycle through open/closed (in the case of the central pore unit) and resting/activated (in the case of ancillary sensory domain) conformations .We further assume that the open and closed or activated and resting conformations have affinity for saturated and unsaturated lipids, respectively.This feature of the model is meant to mimic the dependence of ion channel activation on the properties of the embedding membrane: this has been observed for KvAP channels, which have been shown to respond to the membrane potential differently if immersed in a bilayer of saturated or unsaturated lipids [36] and for the Eag Kv channel whose activation entails a disorder/order transition in the surrounding lipid membrane [22].These observations raise the intriguing possibility that local lipid composition (i.e. the ratio between saturated and unsaturated lipids) and degree of orientational order (i.e.L o vs L d ) can potentially affect the activation properties of channels.
We built a model for the lipid bilayer by extending the particle-based Ising model described in Ref. [48], Fig. 1a and SM.In the particle-based Ising, basic building blocks are connected pairs of beads (i.e. a dumbbell) orthogonally oriented with respect to the membrane plane.The bottom beads are constrained to lay on the plane, and interact only through a volume exclusion short-range repulsion, while the top beads interact via a Lennard-Jones potential (i.e. it also includes an attractive tail).The z position of the top bead can assume two values (the two minima of a quartic potential); we showed in Ref. [48] that these two top-bead heights can be mapped to the spin values -1 and +1, and the dumbbells, fixed in a triangular lattice, have a ferromagnetic/paramagnetic transition varying the temperature.This work extends the previous model by allowing dumbbells to freely diffuse (xy-directions), and the dumbbell vertical orientation (described by the tilting angle θ) can vary as well through a harmonic potential with force constant k (SM).All physical quantities referring to this model are expressed in Lennard-Jones units of mass m, energy ϵ and length σ.
Considering the fact that unsaturated lipids form thinner membranes and occupy on average a larger area per lipid than saturated lipids, we modelled the former by "unsaturated" dumbbells (d u ) which occupy preferentially the s=-1 spin state and have a small force constant k u = 70, and the latter by "saturated" dumbbells (d s ) which occupy preferentially the s=+1 state and have a larger value of k s = 10 4 (Fig. 1b).Each dumbbell type has thus a typical height, although both the s=-1 and s=+1 states are accessible and occasionally populated by all dumbbells.At any fixed temperature, constant pressure molecular dynamics simulations of single-type dumbbells show that due to the different values of k, d u occupy on average a larger area a compared to d s (Fig 1c are simulated together in a 50-50 mixture, we have the additional tendency of dumbbells of same height to interact more strongly, and thus to segregate.This effect depends on the dumbbells height fluctuation, and thus is expected to be temperaturedependent.In general, within the d s -d u mixture, we have a non trivial combination of the demixing tendency and the L o /L d phase transition occurring for each dumbbell type (Fig 1d).Our choice of parameters led to a behavior which is consistent with the one of ternary mixtures.We find in fact that for T = 0.41 ∼ T * u < T * s , the two dumbbell types demix, with d s forming a highdensity, hexatically-ordered phase, identified as L o domains, while d u form a disordered (L d ) phase.For T = 0.436 ∼ T * s , we still observe a nearly phase-separated system where mixing starts to occur but clusters still fluctuate in size, and both phases are disordered (L d ).These fluctuations are large as expected for binary mixture near the demixing transition [49].For T = 0.49 > T * s we have a completely mixed and disordered phase (L d ).
We now introduce a simplified representation of ion channels, to model the observation that each distinct configuration (open or closed) interacts preferentially with a specific lipid species.The goal is to understand how the coupling between the internal degrees of freedom of the channel and those of the lipids can give rise to the collective phenomena mentioned in the introduction.Our goal is to introduce a simple non-trivial model representing a voltage-gated ion channel.A single dumbbell has only two accessible states, up and down, while a collection of dumbbells increases the combinatorial landscape.The presence of multiple states makes it possible to study cooperativity as an enhanced two-state character of the system.
The channel is a collection of dumbbells arranged in a hexagonal shape, with co-planar neighboring beads connected via springs (Fig. 2a).This shape does not perturb the hexatic order of dumbbells and gives a sufficiently large interaction surface with lipids.The z coordinate of top beads is subjected to an external bias force, F p , mimicking the membrane potential.Top beads move on z in a two-well potential, with equally probable states without external bias (F p = 0) and with a smaller barrier than lipid dumbbells, to facilitate transitions between the two states.When biased, one state becomes more populated than the other.By assigning the role of "pore domain" to the central dumbbell, we define a proxy for the two conductance states, closed and open, based on the dumbbell spin state.Accordingly, the remaining six dumbbells can be interpreted as voltage sensor domains.Importantly, due to the attractive Lennard-Jones potential, each channel state interacts preferentially with lipid species having the same height (or spin state).Channels interact with one another or with lipids via the same potential (SM).
We first look at the channel-channel clustering tendency due to lipid-mediated interactions and the resulting segregation effects.We then turn our attention to cooperativity and dynamical and out-of-equilibrium effects.
Lipid-mediated interactions.Figure 2b demonstrates the effect of channels within a membrane at the three temperatures considered (T=0.41,0.436, 0.49) with a fixed channels-to-lipids ratio of f p = 0.05 (we will refer from now on to d s and d u dumbbells as lipids).The top row shows the effect without an applied external potential on the channels, F p = 0, while the bottom row shows the effect with F p = 1.
Comparing the first row of Fig. 2b to Fig. 1d, we see that channels promote lipid mixing at T = 0.41 and F p = 0, in the sense that the L d and L o interface margins are less well-defined.This can be explained by the fact that unbiased channels can act as lineactants between the two phases [50].On the other hand, the lipid-channel attraction strength is the same as for channel-channel, and thus channels remain in a diluted state.Across all temperatures, unbiased channels can diffuse into both lipid phases and tend to assume their height or spin state.
When an external field is applied to channels, they are biased towards one or the other spin state based on the direction of the force (Fig. 2a), assuming either an open or closed conformation.In the case of Fig. 2b, we apply a force of F p = 1 which biases the channels to adopt the spin-up (open) conformation.At equilibrium, channels migrate towards the phase with the same height, which in this case is the phase with spin-up lipids.As a result, channels effectively behave like lipids of the same height, causing the system to become more similar to one made up of only lipids (compare second row of Fig. 2b to 1d).
Of particular interest is the determination of whether channels show an attraction or repulsion between each other, due exclusively to the fact that channels are im-mersed in lipids with the same or opposite spin.To investigate this effect, we calculated the potential of mean force between channels when they are both open or both closed and when one is open and the other one is closed (Fig. 2c).The potential of mean force is computed as G ss (r) = −k B T ln P ss (r), with P ss (r) being the probability of having two channel pores with spin values ss at distance r.We find that channels in the same conformation display long-range lipid-mediated attraction, which is stronger near the demixing transition and decreases as the lipids begin to mix.Moreover, channels in the open conformation will feel a stronger attraction than channels in the closed conformation.In contrast, channels in different conformations show long-range lipid-mediated repulsion.These effects can be understood by considering that channels in the same phase are constrained to have the same conformation, regardless of their distance.Open channels are embedded in a denser environment and, as a result, lipid-mediated interactions are stronger.The range of the potential of mean force reflects the typical size of the underlying lipid patches.Note that the existence of these lipid-mediated long-range interactions occurring under critical conditions (known as Casimir forces) was already predicted theoretically [15,51] and shown in atomistic simulations [17,51].
Conversely, the effect of channels on the lipid patch can be understood by considering the lipid-lipid spincorrelation function (Fig. 2d), defined as C ss (r) = ⟨s i s j ⟩, with s i the spin assigned to each lipid, and r the distance between two lipids.At F p = 0 and temperatures near the demixing, the function decays much faster than the corresponding decay calculated for systems without channels (dashed line), while at F p = 1 the correlation becomes longer.The longer-range correlation for F p = 1 is a consequence of the fact that channels are all in the same state and thus attract same-height lipids around them, locally increasing the lipid spin-spin correlation.Conversely, at F p = 0 channels are free to fluctuate between the two states and adapt to the local environment.This lineactant-like behavior causes a roughening of the domain boundaries and thus a decrease in lipid spin-spin correlation.
Percolation transition.We then turned our attention to the effects of channel density on the system, Fig. 3, fixing T = 0.436.At F p = 0 and for small values of f p , we observe that the size distribution of channel clusters has the typical shape of a power law with exponential cut-off correction, P (n) ∼ n −τ e −n/n * , with n * a typical cluster size (Fig. 3a, first row).By increasing the channel density, the distribution first tends to a pure power-law for f p ∼ 0.2, P (n) ∼ n −τ , and then it develops a peak at large n (not shown).This behavior is in general associated with a percolation transition [52], which in two dimensions is characterized by the critical Fisher exponent τ = 187/91.Notably, the same exponent is observed in the experiments shown in Fig. S1 from Ref. [8].
An interesting role in this percolation transition is played by the external field F p (Fig. 3a, second row).
When applied, the field shifts the transition towards lower densities, f p ∼ 0.07.The reason is that, when activated, channels are encouraged to migrate to the most favorable lipid phase.In this way, the surface area occupied by channels decreases (by a factor of 2) and, accordingly, the effective density increases (Fig. 3b).Overall, the percolation transition is observed for a smaller channel-lipid ratio.
Channel-Channel coupling and cooperativity We now investigate whether the presence of lipidmediated interactions has any effect on the ion channel response to the trans-membrane potential.Experimentally, this may be investigated in a voltage-clamp setting by recording the steady-state ion current, or pore conductance, G, as a function of the applied voltage V .Voltage-gated ion channels contain specialized domains (voltage sensors) that increase the voltage sensitivity of the conduction pore through allosteric coupling; this increased responsiveness to the transmembrane potential is manifested as a steepening of the "G-V" curve.We expect lipid-mediated interactions between voltage sensors in the same channel and between neighboring channels to enhance voltage sensitivity of conductance in a predictable temperature-dependent manner.
To reproduce the "G-V" curve, we computed the equilibrium probability, P up , that the central "pore" particle of each 7-dumbbell channel is in a spin-up state (open conformation).The conductance is related to P up through G = N gP up , where N is the number of channels and g is the single pore conductance (assumed constant).We simulated a normalized "G-V" curve by calculating P up as a function of the bias force, F p , as shown in Fig. 4a, where we compare these curves in the case of an isolated (bare) channel, and of a channel embedded in the lipid bilayer, for three different temperatures.As anticipated, we notice a temperature dependence, which is stronger in the case of embedded channels.
However, one must distinguish between the effect due to a trivial temperature dependence of Boltzmann weights and the temperature dependence that arises from subtle modulations of intra-and inter-channel interactions.To discriminate between these two effects, we considered the "conductance" Hill plot, which eliminates the trivial T -dependence by reporting on the free energy of pore activation W = k B T log(P up /(1 − P up )) as a function of the external field [53].In particular, we are interested in the vertical separation ∆W between the two linear asymptotic trends observed at extreme field strengths.This amounts to removing from the free energy the linear field-coupling term, thereby isolating the effective interaction energy between the pore and the voltage-sensitive particles [53].The entire procedure eliminates the weak inverse-temperature dependence of the "G-V" slope intrinsic to the Boltzmann distribution, which exists even for the bare channel (Fig. 4a), and reveals a more interesting lipid-mediated temperature sensitivity.
An exact expression of W is found for the isolated (bare) 7-dumbbell channel, which exists as a system of connected spins with ferromagnetic coupling constant J 0 between the pore and each voltage sensor (SM).In this system, the coupling energy of the bare channel is found to be proportional to J 0 , ∆W 0 = 24J 0 .Embedding the channels in a lipid environment affects the coupling energy ∆W , see Fig. 4b.This modulation of lipid-channel and channel-channel interactions can be captured by an effective, temperature-dependent coupling constant, J ≡ ∆W/24, which takes into account the effect of the environment and provides a quantitative measure of the system's cooperativity.Increasing J , and thus ∆W , steepens the closed-to-open transition emphasizing the two-state character of the channel.Fig. 4c shows the fitted ∆W for different temperatures and channel densities, compared to the temperatureindependent value ∆W 0 for the bare channel (SM for the fitting procedure).Two major observations can be made from the figure.
First, there is a crossover between J and J 0 around the critical temperature (T ∼ 0.436) with J > J 0 or J < J 0 depending on whether the channel is embedded in a demixed or mixed lipid membrane.This dependence on temperature highlights the crucial role of membrane phase behavior.In the lower temperature demixed state, channels can migrate toward the lipid patch that stabilizes the instantaneous conformational state (i.e.closed or open).This is not possible in the higher temperature mixed state, in which the fluctuating local lipid composition gives rise to a more heterogeneous ensemble that decreases the overall two-state character of the system.
Second, higher channel densities (f p ) result in larger J at a given temperature, suggesting that channels are coupled to one another and tend to activate in a cooperative fashion; the biggest gap occurs around the demixing temperature, T ∼ 0.436.We infer that channelchannel interactions (both direct and lipid-mediated) increase cooperative gating beyond what is expected from independently-gated channels.Notably, similar channelchannel gating cooperativity was found in mechanosensitive channels [54] where their attraction promotes channels closure.
Out-of-equilibrium effects and hysteresis.The demonstration that channel-lipid interactions enhance cooperative gating raises an interesting question: since the response of any given channel depends crucially on the activation state of neighbor channels (and on the thermodynamic phase of the surrounding membrane), is the closed-to-open transition quasi-static or dominated by non-equilibrium effects?
To answer this question we characterized the response of channels to a time-dependent external field with a strength that increases linearly with time.Fig. 4d illustrates P up as a function of the instantaneous value of the applied F p .The system starts at equilibrium with F p = −1 and the field is gradually increased to F p = 1 for a total ramp duration τ (red curve, ON ramp protocol); the black curve shows the mirror-image OFF ramp proto-col.A microscopically reversible process should give rise to overlapping curves.Instead, there is definite hysteresis even for relatively long ramp times τ .Close inspection of a sample trajectory highlights the reason underlying this behavior.Inverting the polarity of the applied field results in two combined effects: it causes channels to migrate from one phase to the other, and it promotes a phase transformation in the region surrounding channels.If τ is small compared to both the typical diffusion time L 2 /D (where L is the typical patch size and D is the channel diffusion constant) and the phase-separation relaxation time, then the process is out of equilibrium and hysteresis appears.To quantitatively characterize this behavior, we report the field value resulting in the activation of half of the channels, namely F 1/2 , as a function of τ for ON and OFF ramps (Fig. 4e).We note that the ON and OFF protocols tend asymptotically to the same equilibrium value and that the rate of convergence is markedly temperature-dependent, with hysteresis being stronger for lower temperatures.
Regarding finite size effects, one expects that the system size imposes a cutoff on the maximum patch size encountered, and in turn the latter changes the characteristic time needed for channels to cross different patches.We thus expect that dynamical properties are affected by finite size effects more than equilibrium ones.Hysteresis, which depends on the time required for channels to migrate from one patch to the other is indeed enhanced by considering larger systems (see SM).

B. Bridging membrane physics and physiology: a lattice model
We shift gears from the molecular dynamics model that is chiefly concerned with the phase structure of the membrane to a system with more physiological relevance.Physiologists have traditionally modeled ion channel gating kinetics using continuous-time discrete-state Markov models [5,6].Gating schemes with a finite number of states adequately fit macroscopic and single-channel data in response to a variety of external stimuli [55], but may be inadequate in describing types of anomalous gating explored in this paper-clustering-induced cooperativity, hysteresis, modal gating, and long-term memory.The crux of our thesis is that a complex membrane absorbs and distributes information like a field, altering the kinetics of bare channels.We thus set out to include additional degrees of freedom to the traditionally used discrete-state models to account for lipid-channel coupling.We demonstrate that this paradigm can provide a unified explanation for anomalous gating.
We consider an established kinetic model for potassium channel gating and extend it to include the effect of channel-lipid interactions.A channel system suitable for our purposes is the large-conductance calcium-activated potassium (BK) ion channel, which, under conditions of zero-or saturating calcium, is regulated by four iden-tical voltage-sensing domains in a manner closely described by the Monod-Wyman-Changeux (MWC) model of allosterism [7,56,57].The kinetic model corresponding to an MWC scheme contains ten states (Fig 5a).Each forward step corresponds to the activation of a single domain, either the pore or one of the voltage sensors.This discrete-state model results in a more realistic description of gating kinetics and a more faithful structural representation than could be achieved with the hexagonal model studied with molecular dynamics.
In order to study the effect of a dynamic lipid environment, we constructed a square-lattice system that best adapts to the four-fold symmetry of the BK channel and possesses Ising-like interactions between neighboring binary cells, see Fig 5b for a schematic representation.The lattice contains channels and lipids, of which there is a 50:50 mix of two types, "saturated" and "unsaturated".The complete lattice is large enough (128 2 cells) to contain 100 channels.An individual channel is composed of a central pore in contact with four symmetrically arranged voltage sensors, occupying a total of 20 cells.The activation of each voltage sensor corresponds to the transport across the membrane of an elementary charge q.Macroscopically, this gives rise to a current whose time integral Q is a typical measurement done in patch-clamp experiments.The response of Q to the applied voltage V is encoded by the "Q-V" curve.Each lipid, which occupy a single cell, is permanently fixed in the "up" or "down" state depending on whether it represents a saturated or unsaturated species.Channels can translate and rotate, while lipids only exchange their state with neighbouring lipids (right side of Fig 5b), so that the saturated/unsaturated ratio is preserved.Lipids are mobile and preferentially associate with their own type via an energy penalty incurred for oppositely aligned pairs.Importantly, lipids interact with the pore and the voltage sensor in a state-dependent fashion with the same energy penalty, where unsaturated/saturated lipids preferentially bind the resting/activated states.
The model parameters were chosen to be consistent with broad stroke features of ion channel dynamics, namely i) the steepness and relative positions of equilibrium gating charge (Q-V) and conductance (G-V) curves [53,58], ii) out-of-equilibrium decay rates of these same quantities [7], and iii) translational diffusion coefficients of lipids and channels.The dynamics of the system were simulated using the kinetic Monte Carlo method of Gillespie [59] based on the chemical master equation whose rate constants from two states are assumed to have an Arrhenius form with a different rate coefficient for each specific Monte Carlo move.A detailed description of parameter choices and methods for computation can be found in the SM.
Fig. 5c shows a snapshot of the lattice system near the mid-point of channel activation (V = −40 mV) at subcritical (10 • C) and supra-critical (80 • C) temperatures.At 10 • C the system is phase-separated (demixed).Similar to the mesoscopic model, the activation state of a channel generally aligns with the lipid phase in which it resides.At 80 • C, phases are mixed.Fig. 5d shows a long (20 sec) single channel tracing characterized by periods of substantially different open probabilities linked to slow transitions between different lipid phases, increasing in frequency with rising temperature.Because these phase-like transitions are slow compared to the millisecond Markovian kinetics of the isolated channel (bottom trace, at 30 • C), they can be considered modal events [38,39].
Cooperativity and hysteresis To test for hysteresis, we implemented a double ramp protocol with rising (ON) and falling (OFF) phases acquired at a ramp speed of ±0.1 mV/ms.With this slow-ramp protocol, entire (quasi-static) ON and OFF activation curves were acquired in a single 4-second simulation.We measured the conductance Hill plots and Q-V curves of averaged activation curves from 100 embedded channels (channel density of f p = 7.8 • 10 −3 ), Fig. 6a-b.The vertical separation of skew asymptotes in the Hill plot was used to measure the pore-voltage sensor coupling energy ∆W at T = 10 • C, for isolated (bare) and embedded channels.
The bare channel yielded overlapping ON and OFF activation curves while Q-V hysteresis was seen with embedded channels, accompanied by temperature-sensitive cooperativity (Fig. 6c-d).These behaviors are consistent with the idea developed in the previous section that membrane demixing promotes channel-channel cooperativity and hysteresis arises from the relatively long time it takes a channel to migrate between phases.Increasing the temperature, we again observed lessening hysteresis from the mixing of lipid phases.
Long-term memory and Hurst analysis.We questioned whether lipid-embedded channels can exhibit long-term memory of the type observed in single channel records from BK channels [30][31][32][33][34].This phenomenon differs from short-term memory in Markov models, in which dwell time distributions demonstrate a finite number of exponential decays [60].Long-term memory follows a power law characterized by self-similarity in a process B(t) satisfying B(at) ∼ |a| H B(t), where H is the Hurst exponent [61].If H > 1/2 (H < 1/2), the process has long-term positive (negative) autocorrelation, while H = 1/2 for uncorrelated processes.
We applied Hurst analysis to consecutive open dwell times from 100 embedded channels (Fig. 6e).Simulations were run under equilibrium conditions at the halfactivation potential V = −40 mV and at temperatures T = 30 • C and 80 • C, until every channel experienced at least 2 14 opening events.Bare channels governed by the ten-state Markov model demonstrated longitudinal homogeneity in their dwell-time sequence, whereas dwelltime sequences of embedded channels were more erratic.It was possible to fit the broadened dwell-time distribution of the embedded channel to the original ten-state Markov model by adjusting kinetic parameters.The fitted Markov model of the embedded channel was used as a control in the Hurst analysis.
We performed the standard determination of the Hurst exponent from the slope of log 2 R(n)/S(n) versus log 2 (n), averaged over subsets of data of size n, where R(n) is the range of the summed deviation from the mean dwell time and S(n) is the standard deviation [61].The slope was estimated through linear regression for n ranging from 2 to 14.This method overestimates H due to sampling bias in small bins [62] but is consistent with previous analysis of experimental data [35].In practice, determining the absence of long-term correlations was achieved by comparing to internal controls.
Typical Hurst plots and fitted exponents are shown in Fig. 6f-g.The average Hurst exponent for the bare channel (Markov case) is equal to 0.542 ± 0.018 (mean ± std.dev.), a result that is independent of temperature.Although H is slightly greater than 0.5, it is statistically no different from an uncorrelated series of exponential deviates (H = 0.539 ± 0.019, p = 0.24, Student's t-test), demonstrating the absence of long-term memory.
On the other hand, in embedded channels we find H = 0.701 ± 0.018, consistent with outside experimental values of H in the BK channel that vary from 0.61 to 0.71 [30][31][32]34].There is lipid-induced broadening of dwelltime distributions that accompanies increased H values, but this does not explain long-term memory, as the fitted Markov model to the embedded channel distribution is near the uncorrelated value (H = 0.561 ± 0.016).Increasing the temperature to 80 • C reduces but does not abolish long-term memory (H = 0.642 ± 0.016), consistent with the earlier hysteresis results.

II. CONCLUSIONS
The aim of our study was to investigate the consequences of state-dependent channel affinity for different lipid species in demixing membranes.We found that when lipids are phase-separated, ion channels exhibit clustering, cooperativity, hysteresis, and long-term memory effects.
In particular, we have shown that fluctuations in local lipid composition give rise to effective attractive interactions and to correlations between ion channels even when they are far apart.These effects are responsible for the observed clustering tendency and cooperativity of activation.Moreover, since the typical relaxation times of these fluctuations is longer than any time-scale involved in channel activation, channel dynamics shows memory effects, including hysteresis of the activation curve.
The membrane influence provides explanations for previously puzzling experimental observations and generates testable predictions.In particular, fluorescence microscopy techniques could be used to measure the size distribution of clusters of ion channels for different lipid composition to correlate cluster properties to the phase behavior of the underlying membrane.Our model predicts that the largest clustering tendency occurs for near-critical membranes with a characteristic size distribution showing a power-law tail with a precise exponent.From the point of view of electrophysiology, we expect the slope of the G-V and Q-V curves (and their hysteresis effects) to show distinct dependencies on temperature and lipid composition, reflecting the fact that critical membranes enhance cooperativity in a predictable fashion.Similarly, the hypothesis of a long-range coupling mediated by lipids could be tested by measuring the influence of channel density on the slope of the activation curve.Finally, single channel recording experiments could be designed to characterize long-term memory in ionic currents for different temperatures and lipid compositions.
Experimental confirmation of the paradigm proposed here would pave the way for a new class of effective models bridging membrane physics and physiology.Coupling ion channel conformational states and local lipid composition introduces a new level of complexity to the quantitative description of channel activity by taking into account the cellular context.For instance, the activity of several members of voltage-gated ion channels family is regulated by PIP2, a signaling lipid molecule that can facilitate or inhibit channel activation [63] and that has been suggested to localize preferentially in lipid rafts [64].The coupling between channel activation and rafts formation/disruption hypothesized here could have consequences on the PIP2 local concentration dynamics and thus potentially underlie some of the complex mechanisms employed by cells to modulate membrane excitability.Although the generality of these mechanisms is not yet fully understood, it is likely that state-dependent lipid coupling plays a significant physiological role, particularly in neurons where channels are abundant and interconnected [8].Specifically, clustering and cooperativity [37] appear to have important functional consequences [9].Importantly, some of these mechanisms are likely at work in other types of membrane proteins with similar state-dependent affinity [23], suggesting a universal class of self-organizing biomolecular systems.

FIG. 1 .
FIG. 1.(a) Mesoscopic model of lipids, represented by a dumbbell which can diffuse on the xy-plane (bottom bead constrained in plane), has a varying z-height of top bead, and a deviation (θ angle) from the vertical orientation restrained by a harmonic potential of force constant k (SM).(b) Bilayer-dumbbell mapping.Unsaturated (saturated) lipids are mapped onto du (ds) dumbbells, with a global minimum of the height potential (see V (z)) corresponding to short (long) dumbbell or spin s=-1 (+1), and a small (large) tilting constant, ku = (ks = 10 4 ), which entails a large (small) surface area a. (c) Plots of temperature versus a for du and ds, obtained via heating/cooling the system (two curves for each type).Both dumbbell types exhibit a liquid ordered (Lo, associated to hexatic ordering) to liquid disordered (L d ) phase transition, see snapshots colored with local hexatic modulus |ψ6| (SM).As curves do not coincide due to hysteresis [17](ordering from the disordered phase is much slower than disordering from the ordered phase, we considered as transition temperature the average value in between the heating/cooling curves plateau, T * u ∼ 0.41 and T * s ∼ 0.436.The melting temperature is shifted towards higher T as k is increased.(d) 50-50 mixture of du and ds dumbbells at three temperatures, T = 0.41, 0.436, 0.49.Each column shows an equilibrated conformation for the indicated temperature, colored by the hexatic modulus and spin value (top and bottom row, respectively).

2 FIG. 3 49 FIG. 4 .FIG. 5 .FIG. 6 .
FIG. 2. (a) Mesoscopic model of ion channel, with seven dumbbells assembled in a hexagonal shape and co-planar neighbour beads connected via springs (SM).Top beads are subject to an external field Fp in the z direction, which models the membrane potential.Average top bead height fluctuates around two equally probable states when Fp = 0, while one of the two states becomes favored when Fp ̸ = 0 (see potential V (z) of z, the average channel top beads height).The central dumbbell has the role of the pore, with its spin mapped to a closed-open state.(b) Snapshots of embedded channels (darker colors) in the binary dumbbell mixture of Fig. 1d, at three temperatures (T = 0.41, 0.436, 0.49) and Fp = 0, 1 (top and bottom row).Channel/lipid ratio is fp = 0.05, and dumbbells are colored according to their spin value (color scheme in Fig. 1d).(c) Potential of mean force as a function of distance between two channel pores in the open-open, closed-closed, and open-closed configuration, same T as in (b) and Fp = 0. (d) Lipids spin-spin correlation with embedded channels and compared to lipid only mixture, for the cases shown in (b).
). Importantly, both d u and d s show a transition between the L o and L d phase at a critical temperature T *