<figure><img src="image_1.png" /></figure>

HAL
open science

# Dynamical study of Na v channel excitability under mechanical stress

Q. X. Ma, A. Arneodo, G. H. Ding, Françoise Argoul

► To cite this version:

Q. X. Ma, A. Arneodo, G. H. Ding, Françoise Argoul. Dynamical study of Na v channel excitability under mechanical stress. Biological Cybernetics (Modeling), 2017, 111 (2), pp.129-148. 10.1007/s00422-017-0712-3. hal-01502876

HAL Id: hal-01502876

https://hal.science/hal-01502876

Submitted on 6 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Dynamical study of Na$_v$ channel excitability under mechanical stress

Q. X. Ma<fnref n="1" /> · A. Arneodo<fnref n="2" />,<fnref n="3" /> · G. H. Ding<fnref n="1" /> · F. Argoul<fnref n="2" />,<fnref n="3" />

**Abstract** Alteration of Na$_v$ channel functions (channe-
lopathies) has been encountered in various hereditary muscle
diseases. Na$_v$ channel mutations lead to aberrant excitabil-
ity in skeletal muscle myotonia and paralysis. In general,
these mutations disable inactivation of the Na$_v$ channel, pro-
ducing either repetitive action potential firing (myotonia)
or electrical dormancy (flaccid paralysis) in skeletal mus-
cles. These “sick-excitable” cell conditions were shown to
correlate with a mechanical stretch-driven left shift of the
conductance factors of the two gating mechanisms of a frac-
tion of Na$_v$ channels, which make them firing at inappropriate
hyperpolarised (left-shifted) voltages. Here we elaborate on a
variant of the Hodgkin–Huxley model that includes a stretch
elasticity energy component in the activation and inactivation
gate kinetic rates. We show that this model reproduces fairly
well sick-excitable cell behaviour and can be used to predict
the parameter domains where aberrant excitability or paral-
ysis may occur. By allowing us to separate the incidences of
activation and inactivation gate impairments in Na$_v$ channel
excitability, this model could be a strong asset for diagnosing
the origin of excitable cell disorders.

**Keywords** Na$_v$ voltage-gated channels · Excitable cell disorders · Mechanical stress · Sick-excitable cell · Periodic firing and myotonia · Dynamical systems · Bifurcations

**1 Introduction**

Since the mid-1990s, voltage-gated sodium (Na<sub>v</sub>) channels have been recognised for their implication in the generation and conduction of action potentials in excitable cells (Hille 1992). Na<sub>v</sub> channels are macromolecular protein complexes containing pore-forming α subunits and smaller non-pore-forming β subunits which are embedded inside the plasma membrane. The vertebrate Na<sub>v</sub> channel α subunit is a single polypeptide chain (~260 kDa) that contains the ion-selective gating component. There are 10 genes encoding the mammalian α subunits Na<sub>v</sub>1.1–1.9, and an atypical channel Na<sub>x</sub>. The β subunits regulate channel expression and gating; these cell adhesion molecules (CAMs) belong to the immunoglobulin (Ig) superfamily (Namadurai et al. 2015). Many pieces of evidence suggest that, in addition to regulating electrical excitability, voltage-gated ion channels are also involved in numerous “non-conducting” signalling mechanisms (Kaczmarek 2006). Thus far, the non-conducting role of voltage-gated channels via the β subunits has been best characterised in cell adhesion (Brackenbury et al. 2008). The α subunits remain closed at rest but are activated during cell membrane depolarisation by a fast and transient structural change (activation particles) leading to a flow of sodium ions down their concentration gradient followed by another structural change (inactivation particles) that counterbalances this inward current (Ulbricht 2005). This inactivation mechanism is necessary for the channel to recover its initial excitability (Gribkoff and Kaczmarek 2009). Na<sub>v</sub> channels are not only expressed in excitable cells like nerve and muscle cells but

<figure><img src="image_3.png" /></figure>

F. Argoul
francoise.argoul@u-bordeaux.fr

G. H. Ding
ghding@fudan.edu.cn

<footnote n="1">1 Shanghai Key Laboratory for Acupuncture Mechanism and Acupoint Function, Department of Aeronautics and Astronautics, Fudan University, Shanghai 200433, China</footnote>

<footnote n="2">2 LOMA, CNRS, UMR 5798, Université de Bordeaux, 51 Cours de la Libération, 33405 Talence, France</footnote>

<footnote n="3">³ Laboratoire de Physique, ENS Lyon, CNRS UMR5672, Université de Lyon, 46 Allée d'Italie, 69364 Lyon, France</footnote>
also in cells ordinarily considered as non-excitable, includ-
ing glial cells, fibroblasts, vascular endothelial, immune and
metastatic cancer cells (Brackenbury et al. 2008). Neuronal
sodium channels Na<sub>v</sub> 1.1 – 1.3 and Na<sub>v</sub> 1.6 – 1.9 are found
in the central nervous system (CNS) and peripheral nervous
system (PNS) neurons and within glial cells, whereas Na<sub>v</sub> 1.4
and Na<sub>v</sub> 1.5 are expressed in skeletal and cardiac muscle cells,
respectively (Catterall 2012).

Channelopathies are caused by genetic or acquired defects
in ion channels, the former being the most common. Voltage-
gated channel defects are involved in a wide variety of dis-
eases like epilepsy, migraine, blindness, deafness (Spillane
et al. 2016), diabetes, hypertension, cardiac arrhythmia,
asthma, irritable bowel syndrome, and cancer (Bracken-
bury 2012). Na<sub>v</sub> channel mutations are associated to aber-
rant excitability in skeletal muscle myotonia and paralysis
(Jarecki et al. 2010). These channelopathies include hyper-
kalaemic periodic paralysis (hyperPP), potassium aggra-
vated myotonia (PAM), paramyotonia congenita (PMC)
and hypokalaemic periodic paralysis (hypoPP) (Davies and
Hanna 2001; Jarecki et al. 2010; Simkin and Bendah-
hou 2011). Hyperkalaemic periodic paralysis, paramyoto-
nia congenita and potassium-aggravated myotonia represent
myotonic disorders due to episodic membrane hyperex-
citability of skeletal muscles (Ptacek et al. 1993). These
diseases are all caused by mutations in the sodium chan-
nel genes SCN4A 3, 4, 5 and 6. They are all inherited in
an autosomal-dominant manner. In general, these mutations
slow down or impede inactivation of the Na<sub>v</sub> channel, pro-
ducing either repetitive action potential firing (myotonia) or
electrical dormancy (flaccid paralysis) in skeletal muscles.

Physical factors such as temperature, pressure, stretching
and bending can impact bilayer mechanics and the coupling
to their embedded voltage-gated channels (Markin and Sachs
2004). For instance, membrane stretching not only makes the
bilayer thinner, but also decreases the amplitude of its fluctu-
ations and can modify its viscoelastic properties depending
on its composition. In particular, if the membrane is highly
bent with a high density of inclusions (organic or proteic),
stretching may increase its lateral tension depending on the
tightness of the connection of the inclusions with the bilayer
amphiphilic molecules. In that regard, overstretching may
lead to the distension or even the failure of this connection
and to the inactivation of the embedded channel. In other
words, the coupling of an ion channel (whatever voltage-
gated or mechano-gated) to its surrounding membrane is of
fundamental importance for its good functioning.

The periodic contractility of the heart is one of the best
examples of the impact of a mechanical stimulation on the
organ function through stretch activation of ion channels
(Hu and Sachs 1997). It was early recognised that dis-
tending the right heart atria could increase the heart rate
(Bainbridge 1915). But it took almost a century (Hu and

Sachs 1997; Peyronnet et al. 2016) to distinguish between
mechanically modulated and mechano-gated ion channels
and their respective impacts on cardiac electrophysiology.
For instance, understanding heart rhythm disorder regula-
tion at the entire organ level from the knowledge of single
cell excitability remains a subject of intensive research (Bett
and Sachs 1997; Kohl et al. 1999; Despa and Vigmond 2016).
Stretch-activated channels (SACs) were first recognised from
experiments on chick skeletal muscles (Guharay and Sachs
1984) and further identified in many cell types of different
phylogenetic origins (animals, plants, fungi and even bacteria).
More interestingly, in eukaryotic cells, SACs were shown to
interact directly with the cortical cytoskeleton (Guharay and
Sachs 1984) and with the extracellular matrix (ECM) (Liu
et al. 1996). Stretch-inactivated channels (SICs) are less com-
mon, but have also been reported in several cell types (Sachs
1989; Morris 1990; Sackin 1995).

A single channel patch-clamp recording (Hamill et al.
1981; Sakmann and Neher 2009) implies a mild to strong
aspiration (suction) and bending of the membrane inside a
small capillary that likely increases its tension and activates
(inversely inactivates) the SACs (inversely the SICs). The
larger the stimulus, the higher the probability of the SAC to
open, as long as the channel to membrane connection does
not leak or break. Channels previously labelled as voltage-
gated (Laitko 2004, 2006; Lin et al. 2007; Morris and Juranka
2007) or ligand-gated (Casado and Ascher 1998; Maingret
et al. 2000) are also mechanically sensitive. The only require-
ment for mechanical sensitivity is that the channel changes
its shape between closed and open states and that stress on
the membrane can be transferred to the channel (Markin and
Sachs 2004).

Na<sub>v</sub>, K<sub>v</sub>, Ca<sub>v</sub> and hyperpolarisation-activated cyclic nucl-
eotide-gated (HCN) channels were all shown to exhibit
reversible gating changes with stretch (Shcherbatko et al.
1999; Tabarean et al. 1999; Gu et al. 2001; Calabrese et al.
2002; Lin et al. 2007). The impact of membrane deforma-
tions upon the Na<sub>v</sub> channel activity was suspected since the
late 1990s (Shcherbatko et al. 1999; Tabarean et al. 1999;
Pawson and Bolanowski 2002) and was recently elucidated
by voltage-clamp experiments (Hamill 2006; Morris and
Juranka 2007; Beyder et al. 2010). Na<sub>v</sub> 1.5 channel reversibly
responds to membrane stretch by an acceleration of the rate-
limiting voltage-dependent step leading to a Na<sup>+</sup> ion inward
flow. Inactivation may also be accelerated. Positive feedback
amplification is a fundamental property of excitable cells that
can further explain the very high sensitivity of Na<sub>v</sub> to even
small mechanical modulations. Given that minute deviations
in Na<sub>v</sub> channel kinetics may lead to cardiac arrhythmias,
epilepsy, neuropathic pain and general muscle disorders, it
is critical to determine to which extent Na<sub>v</sub> channel kinet-
ics can be altered by membrane deformations. In this work,
we perform a dynamical system analysis of a modified ver-
sion of the Hodgkin–Huxley model (Hodgkin and Huxley 1952a) for muscle fibres, aiming at quantifying how a modification of the bilayer membrane tension which surrounds a Na<sub>v</sub> channel may impair its normal function (excitability), either by driving it to a very excitable state (periodic spiking) or by blocking its excitability (paralysis).

## 2 Methods

### 2.1 Modelling the dynamical response of excitable cells

#### 2.1.1 Characteristic features of excitable cells

Excitable cells (Hille 1992) are ubiquitous in animals and plants. They are distinguishable from non-excitable cells by their ability to sustain action potentials that manifest as short-lasting events in which the electrical cellular membrane potential rises and falls, separated by long-term recovery phases. These action potentials are produced by voltage-gated ion channels which can sustain the rapid upstroke and downstroke autocatalytic processes due to fast activation gates, followed by slow and rather silent phases due to slower inactivation and recovery mechanisms. The nonlinearity of these excitable systems underlies the fact that the action potential (spike) amplitude depends very little on the size of the perturbation, provided it is suprathreshold.

#### 2.1.2 General Hodgkin–Huxley (HH) formalism

The Hodgkin–Huxley (HH) models for nerve and muscle fibres (Hodgkin 1951; Hodgkin and Huxley 1952a, b; Hodgkin and Horowicz 1959; Huxley 1959) are constructed on the same basis as electrical model equations, combining a capacitance C and three ionic conducting units in parallel. The potential difference V between the inner and outer sides of the membrane depends on the different currents flowing from, towards or through the membrane (assisted or not assisted by ionic channels). The three ionic currents are respectively a sodium current I<sub>Na</sub>, a potassium current I<sub>K</sub> and a background or leak current I<sub>L</sub>. Currents are given as current densities in µA/cm<sup>2</sup>, ionic conductances as surface conductances in mS/cm<sup>2</sup> and capacitances as surface capacitances in µF/cm<sup>2</sup>. The applied (external) current I<sub>app</sub> is equal to the sum of the charging current and the ion transport currents through the membrane that depend on the membrane potential:

$$I_{app} = C \frac{dV}{dt} + I_{Na}(V, t) + I_{K}(V, t) + I_{L}(V, t). \quad (1)$$

Figure 1 illustrates the temporal variation of the different variables of a HH model (indeed its muscle variant defined

in Sect. 2.1.3) during a single spike. The leak current I<sub>L</sub> is purely ohmic with a fixed conductance g<sub>L</sub> and accounts for the passive permeability of the membrane to ions. The potassium current I<sub>K</sub> is produced by a voltage-gated channel (K<sub>v</sub>), a homotetramer with voltage-sensing domains that may move across the channel or change conformation upon gating (like lockers). Upon channel depolarisation and repolarisation (action potential), corresponding to a sharp increase and decrease of V, respectively (Fig. 1a), different steps can be distinguished. The raise of V produces an increase of m<sub>Na</sub> because the inward flow of Na<sup>+</sup> is facilitated (activation of the Na channel) (red curve in Fig. 1b, c). This sodium inward flux saturates within a few milliseconds. The maximum of the slope of V(t) correlates with a progressive increase in the potassium conductance (green curve in Fig. 1b, c), which reaches a maximum value slightly later than the maximum of V. The K<sub>v</sub> channel is modelled by a single gating mechanism, and its conductance factor m<sub>K</sub> is a monotonous and nonlinear function of the membrane potential and gives the K<sub>v</sub> conductance g<sub>K</sub>(V) = g<sub>K</sub>m<sub>K</sub><sup>4</sup>(V) (Fig. 2b); g<sub>K</sub> is the maximum conductance when m<sub>K</sub> = 1. The conductance of the Na<sub>v</sub> channel is modelled by the combination of two gating mechanisms: (i) a first gating mechanism (described by a variable conductance factor m<sub>Na</sub>) similar to the K<sub>v</sub> channel which activates the permittivity of the Na<sub>v</sub> channel to Na<sup>+</sup> ions (depolarisation step) and (ii) a second gating mechanism which is initially opened at low potential and which progressively closes during depolarisation and repolarisation of the membrane. The conductance factor for this second gating mechanism is written h<sub>Na</sub>. The whole conductance of the Na<sub>v</sub> channels g<sub>Na</sub>(V) = g<sub>Na</sub>m<sub>Na</sub><sup>3</sup>(V)h<sub>Na</sub>(V) was proposed by Hodgkin and Huxley (Hodgkin and Huxley 1952a) from experimental observations. Due to the counteracting effects of the two gates, g<sub>Na</sub>(V) reaches a maximum before the maximum of the voltage (Fig. 1c). The maximum of this conductance is decreased by the presence of the second gate h<sub>Na</sub>. This combination of two gates makes the negative spike of Na<sup>+</sup> ion current I<sub>Na</sub> (red curve in Fig. 1d) dramatically sharper than the voltage spike (black curve in Fig. 1a). Finally, after repolarisation, the voltage crosses more negative values that correspond to the hyperpolarisation of the membrane, before returning to its initial value. This slowest process is monitored by the slow recovery of the second sodium gating mechanism h<sub>Na</sub> (blue curve in Fig. 1b). The intersection of the two curves for m<sub>Na</sub> and h<sub>Na</sub> in Fig. 1b corresponds to V = -65 mV which is close to the resting potential of the membrane. From the current and voltage signals, 2D phase portraits can be reconstructed that shed light on the temporal transformation of the system trajectories plotted in phase space. In Fig. 1e are displayed the three ionic currents I<sub>Na</sub>(t), I<sub>K</sub>(t) and I<sub>L</sub>(t) and their sum I<sub>S</sub>(t) versus V(t). In Fig. 1f, I<sub>Na</sub>(t) is represented versus the capacitance charging current I<sub>C</sub>(t) = CdV(t)/dt. HH models [Eq. (1)] are constructed
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Dynamical signature of the different variables involved in a HH model derived for muscle cells (Hodgkin 1951; Adrian et al. 1970). **a** Voltage V(t) vs t. **b** Opening probabilities for the three gates involved in the HH model: m<sub>Na</sub>(t) (red), h<sub>Na</sub>(t) (blue) and m<sub>K</sub>(t) (green). **c** Conductance factors for the sodium m<sup>3</sup><sub>Na</sub>(t)h<sub>Na</sub>(t) (red) and the potassium m<sub>K</sub>(t) (green) channels. In **b** and **c** the rescaled voltage curve (black dashed-dotted curve) is also plotted for comparison. **d** Ionic currents:</figcaption></figure>

I<sub>Na</sub>(t) (red), I<sub>K</sub>(t) (green) and I<sub>L</sub> (grey). In black, the sum of these ionic currents I<sub>S</sub>(t) = I<sub>Na</sub>(t)+I<sub>K</sub>(t)+I<sub>L</sub>(t). e 2D phase plots of the three ionic current signals and their sum versus V. f 2D phase plot of I<sub>Na</sub>(t) versus the charging current I<sub>C</sub>(t) = CdV(t)/dt. The arrows indicate the orientation of the dynamics on the phase-space trajectories. For parameter values, see Sect. 2.1.3; I<sub>app</sub> = 7µA/cm<sup>2</sup> (colour figure online)

from four variables and these phase portraits are projections on a phase space with only two dimensions, which may not be sufficient if the dynamics is more complex than the one illustrated in Fig. 1. In the first phase portrait (Fig. 1e), the comparison of the three ionic current dynamics and their relative importance in the total ionic current is enlightened. For instance, it becomes clear that the leak current (grey dots) is not negligible for large V values and can even overtake the total ionic current for V values approaching zero. The interest of phase portrait representations is also to provide a direct separation of fast and slow dynamics from the observed density of sampling dots; the sparser the dots, the faster the dynamics.

When making explicit the different ionic conductances, Eq. (1) becomes:

$$C \frac{dV}{dt} = -g_{\text{Na}} m_{\text{Na}}^3 (V) h_{\text{Na}} (V - E_{\text{Na}}) \\ - g_{\text{K}} m_{\text{K}}^4 (V) (V - E_{\text{K}}) \\ - g_{\text{L}} (V - E_{\text{L}}) + I_{\text{app}}(t), \qquad (2)$$

where E<sub>Na</sub>, E<sub>K</sub> and E<sub>L</sub> are the reversal potentials of the membrane channels corresponding to vanishing currents I<sub>Na</sub>, I<sub>K</sub> and I<sub>L</sub>, respectively; these potentials are the Nernst equilibrium potentials for each ion. The conductance coefficients

g<sub>Na</sub>, g<sub>K</sub> and g<sub>L</sub> of the three ionic currents are constant parameters that do not depend on V.

The dependance of the activation and inactivation variables m<sub>Na</sub>, m<sub>K</sub> and h<sub>Na</sub> and of the corresponding ionic conductances on the voltage V is essential ingredients of the HH models. The activation and inactivation variables m<sub>Na</sub> and h<sub>Na</sub> have opposite voltage dependencies. Depolarisation causes m<sub>Na</sub> to increase and h<sub>Na</sub> to decrease, whereas repolarisation causes m<sub>Na</sub> to decrease and h<sub>Na</sub> to increase (Namadurai et al. 2015) (Fig. 2a). In HH models, the probability of Na<sup>+</sup> permissivity is considered as equal to the product m<sup>3</sup><sub>Na</sub>h<sub>Na</sub>, assuming that the two Na gates m<sub>Na∞</sub> and h<sub>Na∞</sub> act independently. These transition probabilities dm<sub>i</sub>/dt (i = 1, 2, 3) = (dm<sub>Na</sub>/dt, dh<sub>Na</sub>/dt and dm<sub>K</sub>/dt), from permissive to non-permissive states of the different ionic gates, are assumed to obey first-order rate equations with rate factors α<sub>i</sub>(V) and β<sub>i</sub>(V) that depend on the membrane potential V:

$$\frac{dm_i}{dt} = \alpha_i(V)(1 - m_i) - \beta_i(V)m_i. \qquad (3)$$

V is measured in millivolts and is defined such that a positive V value corresponds to an increased inward voltage gradient across the membrane. The probability that a gate opens over a short interval of time is proportional to the probability (1 − m<sub>i</sub>) of finding the gate closed, multiplied by the opening rate
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Dependance on the voltage V of the steady-state conductances and currents in a HH model of muscle fibres ($\sigma = 0$) [Eqs. (5–8)]. **a** Probabilities of channel gate opening $m_{Na\infty}$ (red), $h_{Na\infty}$ (blue) and $m_{K\infty}$ (green). **b** Conductance factors in the steady state for the sodium channel $100 * m_{Na\infty}^3 h_{Na\infty}$ (red) and the potassium channel $m_{K\infty}$ (green). **c** Time constants for the three gates $\tau_{mNa}$ (red), $\tau_{hNa}$ (blue) and $\tau_{mK}$ (green). **d** Ionic currents in the steady state $I_{Na\infty}$ (red), $I_{K\infty}$ (green) and $I_{L\infty}$ (grey) and the sum of these three currents $I_{S\infty}$ (black). **e** Zoom in (c) on the range of V values corresponding to the local maxima of the time constants. **f** Zoom in (d) showing the sigmoidal shapes of $I_{Na\infty}$ and $I_{S\infty}$ (colour figure online)</figcaption></figure>

$\alpha_i(V)$. Likewise, the probability that a gate closes during this short interval of time is proportional to the probability $m_i$ of finding the gate open, multiplied by the closing rate $\beta_i(V)$. Dividing Eq. (3) by $\alpha_i(V) + \beta_i(V)$ brings out the time constants $\tau_i$ for each gate:

$$ \tau_i(V) \frac{dm_i}{dt} = m_{i\infty}(V) - m_i, \qquad (4) $$

where

$$ \tau_i(V) = \frac{1}{\alpha_i(V) + \beta_i(V)}, \quad m_{i,\infty}(V) = \frac{\alpha_i(V)}{\alpha_i(V) + \beta_i(V)}. \quad (5) $$

Equation (4) tells us that for a fixed voltage V, $m_i$ approaches the limiting value $m_{i,\infty}(V)$ exponentially with a time constant $\tau_i(V)$.

Experimental data have revealed that the $m_{i,\infty}(V)$ for Na and K activation gates ($i = 1, 3$) follow a sigmoidal increasing shape from 0, at large and negative V values, to 1 at large positive V values. Inversely, for the inactivation of Na channel ($h_{Na}$), the closing rate increases with V, whereas the opening rate decreases with V. These behaviours are nicely reproduced by the HH models (Fig. 2a).

### 2.1.3 Variants of the HH rate equations for muscle cells

The model equations for three gating subunits (Na activation, Na inactivation, K activation) was initially proposed in 1952 by Hodgkin and Huxley from squid giant axon voltage-clamp experiments (Hodgkin 1951; Hodgkin and Huxley 1952a, b). The interest of Hodgkin and Huxley and their coauthors for excitability in mammalian cells was not limited to
nerve cells, since they also devoted much of their efforts to characterise the excitability and conduction of muscle fibres (Hodgkin 1951; Hodgkin and Horowicz 1959; Huxley 1959). Similar to nerve action potentials, they confirmed older observations by Fenn and Cobb (1936) of the implication of sodium ion entry and potassium ion exportation in mammalian muscle excitation spikes. This similarity can be explained by the presence of a higher concentration of potassium inside these excitable cells, whereas sodium and chloride are relatively dilute. The mechanisms of excitation conduction elaborated for nerve cells were generalised to skeletal muscle cells (Butchal and Sten-Knudsen 1959). However, unlike for unmyelinated nerve fibres as the squid giant axons (Hodgkin and Huxley 1945), if the initial phase of membrane repolarisation to ~-70 mV was shown to occur within 3 ms at 20 °C, an additional 30 ms was necessary for the restitution of the resting membrane potential. The hyperpolarisation phase was also comparatively diminished, and repolarisation was observed to remain in progress after the peak of isometric twitch tension of the muscle fibre.

Thanks to the development of micro-electrode techniques for single cell membrane potential difference measurements, a variant of the original HH model was proposed by Adrian et al. (1970). The opening $\alpha_i$ and closing $\beta_i$ rates were adapted from the original HH equations (Hodgkin and Huxley 1952a) to match the excitability of muscle cells:

$$ 
\begin{aligned}
\alpha_{m_{Na}}(V) &= 0.288\gamma_1/(1 - e^{-\gamma_1/10}), \\
\beta_{m_{Na}}(V) &= 1.38 e^{-\gamma_1/18}, \\
\alpha_{h_{Na}}(V) &= 8.1 \ 10^{-3} e^{-\gamma_2/14.7}, \\
\beta_{h_{Na}}(V) &= 4.38/(1 + e^{-\gamma_2/9}), \\
\alpha_{m_{K}}(V) &= 0.0131\gamma_3/(1 - e^{-\gamma_3/7}), \\
\beta_{m_{K}}(V) &= 0.067e^{-\gamma_3/40},
\end{aligned}
\quad (6) $$

where $\gamma_1 = 44 + V$, $\gamma_2 = 45 + V$ and $\gamma_3 = 40 + V$. The $\alpha_i$ and $\beta_i$ are given in reciprocal ms ($T = 25$ °C) and the $\gamma_i$ in mV. Let us note that the rates of closing of the Na ($\beta_{m_{Na}}$) and K ($\beta_{m_K}$) activation gates and the rate of opening of the Na ($\alpha_{h_{Na}}$) inactivation gate are written as exponential decreasing functions of the channel voltage. Inversely, the rates of opening of the Na ($\alpha_{m_{Na}}$) and K activation ($\alpha_{m_K}$) gates and the rate of closing ($\beta_{h_{Na}}$) of the Na inactivation gate are written as the ratio of a linear function of V (shifted by a constant value) divided by an exponential decreasing function of V. This reproduces nicely the experimentally observed opposite variations of Na and K activation rates and of Na inactivation rate with the voltage V (Fig. 2a).

The reversal potentials corresponding to vanishing $Na^+$, $K^+$ and leak currents (Nernst potentials) for muscle cells are $E_{Na} = 47$ mV, $E_K = -93$ mV and $E_L = -85$ mV. The maximum transient conductances are $g_{Na} = 150$ mS/cm²

and $g_K = 21.6$ mS/cm², whereas the leak conductance is $g_L = 0.75$ mS/cm². The steady-state open probabilities for the activation and inactivation Na gates are combined in the total conductance factor for $Na^+$ ions:

$$ m_{Na\infty}^3 h_{Na\infty} = \left[ \frac{\alpha_{m_{Na}}(V)}{\alpha_{m_{Na}}(V) + \beta_{m_{Na}}(V)} \right]^3 \times \left[ \frac{\alpha_{h_{Na}}(V)}{\alpha_{h_{Na}}(V) + \beta_{h_{Na}}(V)} \right]. \quad (7) $$

The corresponding steady-state open probability for the K channel yields the total conductance factor for $K^+$ ions:

$$ m_{K\infty}^4 = \left[ \frac{\alpha_{m_K}(V)}{\alpha_{m_K}(V) + \beta_{m_K}(V)} \right]^4. \quad (8) $$

The time constants for each rate equations [Eq. (4)] are respectively given by Eq. (5).

The different activation and inactivation probabilities in the steady state ($m_{Na\infty}$, $h_{Na\infty}$ and $m_{K\infty}$) [Eq. (5)] and the conductance factors $m_{Na\infty}^3 h_{Na\infty}$ [Eq. (7)] and $m_{K\infty}^4$ [Eq. (8)] are calculated numerically from Eq. (6) and plotted in Fig. 2a, b. Note that the Na conductance factor $m_{Na\infty}^3 h_{Na\infty}$ is much smaller than the K conductance factor and has been multiplied by 100 for better visibility. The time constant of the $Na_v$ activation gate (fast upstroke and downstroke dynamics) is much smaller than the two other time constants for K activation gate and Na inactivation gate (Fig. 2c, e where $\tau_{m_{Na}}$ has been multiplied by 30). The range of potential values where these time constants are maximal shows a clear separation of the activation (Na and K) and inactivation (Na) processes. For membrane potentials close to the activation maximum ($V \sim -54$ mV), the time constant for the activation gate of the K channel (green line) is about 30 times larger than the one for the activation (red line) and inactivation (blue line) gates of Na channel. After the upstroke which takes a few milliseconds, the membrane potential rapidly returns to more negative potentials and then remains in a hyperpolarised state for several milliseconds at a non-null total current $I_S$ (Fig. 1a). This phenomenon is also confirmed by the variation of the stationary ion currents $I_i$ (Fig. 2d) in this range of potential values (Fig. 2f). In the steady state, $I_{Na\infty}(V)$ (red line) is smaller than $I_{K\infty}(V)$ (green line) and never crosses it (checked from numerical simulations) and both these currents vary much slower than the leak current $I_{L\infty}$ which dominates the total current. At membrane voltages below $V = -85$ mV where the sum of these three currents gets null, the capacitance charging current $I_C = CdV/dt$ dominates in Eq. (1).

Let us point out that both $I_{Na\infty}$ and $I_{S\infty}$ display a sigmoidal shape versus V (Fig. 2f), whereas in the nerve fibre HH model, $I_{S\infty}$ is a smoothly increasing function of V. Interestingly, the steady-state current $I_{S\infty}(V)$ shows a decreasing
regime between -69.3 and -58 mV sandwiched by two increasing regimes. This means that for a single value of the input current $I_{app}$, three values of V corresponding to three steady states are possible. The stability of these three steady solutions varies with the model parameters, and only the stable ones will be observable in experiments. The local stability of these steady states will be analysed in Sect. 3.1.

## 2.2 Modelling the mechanosensitivity of voltage-gated channels

Variants of HH model equations (Hodgkin and Huxley 1952a) were used for the interpretation and the prediction of impairments of sodium channel in myotonia and paralysis in humans (Cannon et al. 1993; Cannon and Corey 1993). The phenotype of myotonia muscle disorders corresponds to a delayed relaxation of tension after contraction and an enhanced electrical excitability that manifests in repetitive discharges occurring either spontaneously or after voluntary contraction or muscle percussion. HyperPP may be concomitant to myotonia. (Cannon et al. 1991) showed that functional defects in the sodium channels of hyperPP muscle fibres can impair the inactivation gate when the extracellular potassium is elevated. Partial disruption of the inactivation mechanisms of sodium channels by a mechanical disruption of muscle fibres and the formation of cell membrane blebs were also observed when adding a peptide toxin (ATX II) (Cannon and Corey 1993). All these defaults in sodium channel functions were broadly related to a persistent inactivation failure of a fraction of the sodium channels at strong depolarising potentials. More recently, the question of mechanical, ischaemic and inflammatory injuries of voltage-gated sodium channels and their harmful leakage was embodied in the term “sick-excitable cell” conditions (Morris et al. 2012), corresponding to a left shift of both the steady-state conductance factors $m_{Na\infty}(V)$ and $h_{Na\infty}(V)$ when plotted versus V (Fig. 2a). This left shift can be followed by the decrease of the half-voltage $V_{1/2}$ corresponding to $\beta_i/\alpha_i = 1$. Actually, rather than being inhibited, damaged and leaky Na$_v$ channels gate “too well” since they get activated at inappropriate hyperpolarised (left-shifted) voltages (Shcherbatko et al. 1999; Tabarean et al. 1999; Wang et al. 2009; Beyder et al. 2010). By virtue of molecular coupling between the fast-mode activation and inactivation (availability) processes in a Na$_v$ channel these two processes left-shift in synchrony; the kinetic factors $m_{Na\infty}(V)$ and $h_{Na\infty}(V)$ both shift of the same amount in voltage. This behaviour was shown for recombinant Na$_v$ 1.6 channels and termed “coupled left shift” (CLS) (Morris et al. 2012). CLS leads to maxima of the Na channel conductance factor $m_{Na\infty}^3 h_{Na\infty}$ in a range of membrane voltage which was subthreshold in normal situations. Even if only a small fraction of channel gates were altered, it would never-

theless imply a subthreshold persistent current which could depolarise the adjacent (intact) Na$_v$ membrane channels.

These modifications not only correspond to a drastic change of the $I_S(t)$ curves of stretch-traumatised nerve or muscle cells, but also to a transformation of the type of temporal dynamics that such cells can sustain on longer terms. A way to understand how this modification of the upstroke and downstroke dynamics of muscle fibres can induce ectopic action potential spiking is to use a strategy inspired from dynamical system theory (Izhikevich 2007) and which consists in reconstructing the so-called bifurcation diagrams in both normal and altered sodium-gated channels. A first attempt of bifurcation diagram reconstruction for single ion channel cluster was recently proposed by Assmann and Lenz (2014), based on a modification of the squid axon model of Hodgkin and Huxley (1952a). However, the classification of the different dynamics encountered in the normal and modified HH models was limited to the observation that when the membrane surface tension was fixed at values larger than its equilibrium values, stable limit cycle oscillations could replace single, isolated excitability spikes for biologically relevant $I_{app}$ values. They predicted that the membrane potential would undergo a transition from non-spiking to stable periodic spiking regions upon increasing the surface tension. In this work, we follow the same modelling rationale, and we investigate how a modification of the kinetic rates of the Na activation ($m_{Na\infty}$) and inactivation ($h_{Na\infty}$) gates can change the whole channel dynamics. We start from a different HH model that describes muscle cell excitability and presents multiple steady solutions at normal tensions. From this model, we perform a full analysis of the different dynamics with and without tension and we reconstruct the corresponding bifurcation diagrams (Sect. 3.1).

### 2.2.1 Membrane tension controls channel activation and inactivation rates

The easiest way to account for the elasticity of the bilayer membrane and its modifications upon mechanical stretch is to write its tension as the product of an elasticity constant $\sigma$ and the fractional increase of its area $\Delta A/A$ (Helfrich 1973; Guharay and Sachs 1984; Assmann and Lenz 2014):

$$ \mathcal{T} = \sigma \Delta A/A. \qquad (9) $$

$\sigma$ represents the amount of tension required for doubling the area ($\Delta A = A$). The energy of membrane deformation is:

$$ U = \frac{A}{2\sigma} \mathcal{T}^2 = \frac{\sigma}{2} \left( \frac{\Delta A^2}{A} \right). \qquad (10) $$

It was demonstrated experimentally (Guharay and Sachs 1984) that stretching modifies the rate constants of the
channel gating mechanism according to a simple Eyring formulation of the rate constant between an initial state i and an intermediate activated complex i*:

$$k_{i \to i^*}^{\sigma=0} = k_0 e^{-\Delta G^{\ddagger}/kT}, \quad (11)$$

where k is the Boltzmann's constant, T the absolute temperature and $\Delta G^{\ddagger}$ the free energy of transition from the initial state to the activated complex. Equation (11) for the rate constant $k_{i \to i^*}^{\sigma=0}$ follows an Eyring formulation (transition state theory), which does not mean that the system is strictly at equilibrium (the reactants and the products may not be in equilibrium) but in local quasi-equilibrium. Given that the probability distribution functions of the channel conductance factors in the initial and transition states were unknown a priori, we assumed them as standard distributions, putting aside more complex distributions such as power laws. Within this assumption we related the free energy of transition between the initial state and the intermediate activated complex to an elastic energy. Replacing the membrane energy term defined in Eq. (10) into Eq. (11) yields:

$$k_{i \to i^*} = k_{i \to i^*}^{\sigma=0} e^{\frac{\sigma}{2kT} \left(\frac{\Delta A^2}{A}\right)} = k_{i \to i^*}^{\sigma=0} e^{\sigma \mathcal{B}}, \quad (12)$$

where $\mathcal{B} = \left(\frac{\Delta A^2}{2kT A}\right)$ describes the amount of stretching (assumed fixed) implied by the experimental measurement. $\Delta G^{\ddagger} = \Delta G_{\sigma=0}^{\ddagger} - U$ because the intermediated activated complex is stabilised upon stretching. The open probability of the channel gate and the corresponding current flow will also follow this same relationship.

Here we have restricted our modelling to the linear elasticity regime. In real situations, such as mechanically damaged membrane bilayers, the membrane elasticity may become a nonlinear function of the deformation, such as $\sigma = \sigma_0(\Delta A/A)^\zeta$ with $\zeta > 0$ for strain-stiffening membranes and $\zeta < 0$ for strain-softening membranes.

### 2.2.2 SLS HHm model

Along the same line as previous works on ionic channel mechanosensitivity (Guharay and Sachs 1984; Sachs 1988; Assmann and Lenz 2014; Morris et al. 2015), we introduce this modification on the opening/closing rate of the Na activation gate only (single left shift—SLS), to study more specifically how it changes the excitability and periodic firing of the HH muscle (HHm) model. The $\sigma$ coefficient can be considered as a tension characterising the mechanical sensitivity of the gate. An increase of $\sigma$ (in positive direction) speeds up the opening of the Na activation gate. Conversely, a decrease of $\sigma$ (in negative direction) slows down the opening of the Na activation gate. Intuitively, increasing the membrane tension surrounding the channel favours its opening.

The closing of the Na activation gate works in the opposite, an increase of $\sigma$ (in positive direction) slows down the closing of the Na activation gate and a decrease of $\sigma$ (in negative direction) speeds up its closing. The Na activation gate closing rate is therefore favoured when reducing the membrane tension:

$$\alpha_{m_{Na}} = \alpha_{m_{Na,0}} e^{\sigma \mathcal{B}}, \quad \beta_{m_{Na}} = \beta_{m_{Na,0}} e^{-\sigma \mathcal{B}}. \quad (13)$$

Introducing the parameter $\sigma$ in the model for Na channel activation, we get the following expression for $m_{Na\infty}$:

$$m_{Na\infty} = \frac{1}{1 + \frac{\beta_{Na\infty}(V)}{\alpha_{Na\infty}(V)} * e^{-2\sigma\mathcal{B}}} = \frac{1}{1 + \frac{\beta_{\sigma}}{\alpha_{\sigma}}}. \quad (14)$$

Equation (14) shows that when introducing a mechanical energy term inside the Na channel activation kinetics, this not only implies a nonlinear change of the half-voltage $V_{1/2}$ (corresponding to $\frac{\beta_{\sigma}}{\alpha_{\sigma}} = 1$) but also induces a nonlinear distortion of the $m_{Na\infty}(V)$ curves for $\sigma$ values away from zero.

We take $\mathcal{B}$ in the range of values expected from real situations. For $\mathcal{B} = 129.65 \text{ m}^2/\text{J}$ (Assmann and Lenz 2014), we compute the different probability $m_i$, the conductance factors, the time constants and the steady currents versus the membrane voltage with and without tension $\sigma$ (Fig. 3). Note that the conclusions of our dynamical analysis of the HHm model with mechanical tension do not depend on the specific value of $\mathcal{B}$. By playing on the tension parameter $\sigma$, one can compensate a change in $\sigma$ and recover similar bifurcation diagram changes. For $\sigma > 0$, the Na activation gate opening probability displays a left shift towards negative V values (Fig. 3a) (SLS HHm model). This stretch-activated left shift of the conductance factors (Fig. 3b) reproduces nicely previous experimental measurements on adult rate ventricular myocytes (Banderali et al. 2010). This simple mechanical modification of the HHm model (limited to the Na activation gate) embodies not only an important feature observed in experiments, namely the shift of the currents to negative V values, but also their global amplification (Fig. 3d) (Wang et al. 2009; Banderali et al. 2010; Beyder et al. 2010; Yu et al. 2012; Morris et al. 2015).

Given that the opening probabilities of the other gates are unchanged, we can examine precisely how the shift in voltage of the Na activation gate probability impacts the conductance factor for sodium, the Na activation gate constant and the steady-state currents. The introduction of a positive tension $\sigma$ in the conductance factor $m_{Na\infty}^3(V)h_{Na\infty}(V)$ not only shifts its maximum to lower voltages, but also increases its maximum (Fig. 3b). Hence, one suspects an enlargement of the voltage interval where spiking and/or periodic firing occurs for larger positive $\sigma$ values, coming along with a sharp increase of the amplitude of the ion currents (Fig. 3d). The
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Dependance on the voltage V of the conductances and currents of a SLS HHm model for muscle fibres. **a** Probabilities of Na channel gate opening $m_{Na\infty}$ (red), $h_{Na\infty}$ (blue) and $m_{K\infty}$ (green). **b** Conductance factors in the steady state for the sodium $100 * m_{Na\infty}^3 h_{Na\infty}$ (red) and the potassium $m_{K\infty}$ (green) channels. **c** Time constants for the three gates $\tau_{mNa}$ (red), $\tau_{hNa}$ (blue) and $\tau_{mK}$ (green). **d** Ionic currents in the steady state $I_{Na\infty}$ (red), $I_{K\infty}$ (green), $I_{L\infty}$ (grey) and their sum $I_{S\infty}$ (black). The symbols have the following meanings: $\sigma = 0$ mN/m (solid line), 2 mN/m (empty circles), 4 mN/m (filled circles), -2 mN/m (empty triangles), -4 mN/m (filled triangles) (colour figure online)</figcaption></figure>

time constant $\tau_{mNa}(V)$ for the Na activation gate is also left-shifted in V for $\sigma > 0$ (Fig. 3c), as a consequence of the left shift of $m_{Na\infty}$ (Fig. 3a). Less intuitive is the fact that the time constant maximum of this shifted curve is not changing much with $\sigma$, meaning that the slower affordable rate for this gate is quasi-invariant. For spiking events to occur (which is the fundamental property of excitable cells), one major ingredient is to have well-separated time-scales for the different gating mechanisms. This is still fulfilled in this SLS HHm model. The time constant of the K activation rate is so large as compared to the one of Na activation rate that a left shift of the later for $\sigma > 0$ does not change much the spiking ability of the model. Indeed, it simply changes the range of $I_{app}$ values where these voltage and current large amplitude incursions occur. The modifications of the sodium and total steady-state currents with $\sigma$ (Fig. 3d) are also indicative of the changes observed in the temporal ionic current signals. Actually, we note in Fig. 3d that changes in the sign and amplitude of the steady-state currents $I_{Na\infty}$ and $I_{S\infty}$ are very different for positive and negative values of $\sigma$. For $\sigma < 0$, the sigmoidal shape of $I_{Na\infty}$ disappears for $\sigma \lesssim -2$ mN/m, likely corresponding to a loss of excitability (this aspect will be discussed later), since spiking and periodic firing are conditioned by the presence of several steady states for a given value of $I_{app}$. If nonetheless the system keeps on oscillating, the small amplitude oscillations remain subthresholded without large amplitude spiking. The steady current curves no longer change for arbitrarily large negative $\sigma$ values, which is explained by the dramatic decrease of the conductance factor $m_{Na\infty}^3(V)h_{Na\infty}(V)$ (Fig. 3b). Positive $\sigma$ values lead to a more important change of the steady-state

current curves, since their sigmoidal shape is significantly reinforced (Fig. 3d).

The description of the temporal evolution of the different opening probabilities and conductance factors when varying the parameter $\sigma$ is a preliminary step towards the bifurcation diagram analysis reported in Sect. 3.1. In Fig. 4a, we compare two large amplitude voltage oscillations obtained with the SLS HHm model for $I_{app} = 11~\mu\text{A/cm}^2$, $\sigma = 0$ (thin dashed-dotted line) and $\sigma = 2$ mN/m (thicker solid line). Whereas the relaxational (non-harmonic) character of these oscillations is not affected, the range of V is enlarged, both towards negative and positive V values. To better distinguish the shapes of the two series of curves, we rescale the time t to the period of the large amplitude oscillations. Codimension two diagrams ($I_{app}$, $\sigma$) visualising how the amplitude and period of large amplitude periodic dynamics change with the parameters $\sigma$ and $I_{app}$ is presented in Sect. 3.1.

More interesting are the changes of the opening probability signals $m_{Na}(t)$ (red), $h_{Na}(t)$ (blue) and $m_K(t)$ (green) when varying $\sigma$ in the SLS HHm model (Fig. 4b). We note that a positive $\sigma$ value does not change the sharp rise of $m_{Na}(t)$ since the two red curves are quite superimposed for $t/T_c \in [-0.1; 0]$ and $t/T_c \in [0.4; 1]$, but the length of the plateau at high $m_{Na}(t)$ values for $t/T_c \in [0; 0.1]$ is clearly increased. Since the inactivation gate dynamics is almost unchanged (Fig. 4b), the retardation of the decay of $m_{Na}(t)$ for $\sigma > 0$ means that the inward flux of Na$^+$ ions should increase, leading to a sharper current excursion to negative $I_{Na}(t)$ and $I_S(t)$ values, as confirmed in Fig. 4d. The impact of the tension on the Na inactivation gate opening probability $h_{Na}(t)$, although less important, leads also to a global post-
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Comparison of the dynamical signatures of the different variables involved in the SLS HHm model. **a** Membrane voltage $V(t/T_c)$ vs $t/T_c$, where $T_c$ is the period of the large amplitude oscillations. **b** Opening probabilities for the three gates $m_{Na}(t)$ (red), $h_{Na}(t)$ (blue) and $m_K(t)$ (green). **c** Conductance factors for the sodium $m_{Na}^3(t)h_{Na}(t)$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(red) and the potassium $m_K(t)$ (green) channels. **d** Ionic currents $I_{Na}(t)$ (red), $I_K(t)$ (green), $I_L(t)$ (grey) and their sum $I_S(t)$ (black). $\sigma = 0$ (thin dashed-dotted line), 2 mN/m (thicker solid line). $I = 11 \mu A/cm^2$ (colour figure online)</figcaption></figure>

poning of the dynamics of this gate with time, favouring again larger negative current spike amplitudes for $Na^+$ and also for the total ionic current $I_S(t)$. Surprisingly, the $K_v$ activation gate $m_K(t)$ is shifted in time although its equation does not contain a mechanical term (Fig. 4b). Indeed, the potassium current $I_K$ is a “slaved” variable to the other currents through the variations of $V$ [Eq. (1)]. As shown in Fig. 3b, c, its steady opening probability $m_K$ and its time constant $\tau_K$ change when $V$ exceeds –50 mV. The $K_v$ gate is also much slower than the other ones, and its temporal evolution does not show the drastic change from fast to slow dynamics of the $Na_v$ gates. The strong dynamical amplification (amplitude and relative duration during one oscillation) of the $Na^+$ conductance (~×3) (Fig. 4c) is quite impressive although expected from the SLS HHm model modification. The dynamical amplification of the $K$ conductance factor is more amazing since it means that the outward flux of potassium should also dramatically change, although its gating properties are assumed unchanged under tension. This means that from the measurement of sodium and potassium ion conductances, we could misleadingly conclude of some impact of the stretching of the membrane upon both these ionic gates, whereas the modification of the membrane tension has impacted only one of these gates. This prediction of the SLS HHm model therefore questions the conclusions commonly drawn from simple experimental measures of ion fluxes, even when performed with highly sophisticated and sensitive methods. The different current signals (Fig. 4d) confirm this effect since both

negative and positive incursions of the currents are amplified and extended in time for $\sigma > 0$. The shape of the increasing part of $I_{Na}(t)$ for $t/T_c \in [0.05; 0.25]$ is modified since the exponential decay of the current amplitude (negative value) is replaced around $t/T_c \sim 0.05$ by a linear decay which ends on a small plateau for $t/T_c \sim [0.14–0.2]$ that vanishes through another exponential decay. The small plateau of $I_{Na}(t)$ likely represents the postponement of the Na channel return to its initial state; it exists with and without tension and is simply postponed for $\sigma > 0$. Note that increasing $\sigma$ from 2 mN/m to 4 mN/m enhances further the effects outlined above.

### 2.2.3 CLS HHm model

The CLS HHm model for muscle fibres ensues by introducing a mechanical factor $e^{\sigma \mathcal{B}}$ on both $\alpha_{m_{Na}}$, $\beta_{m_{Na}}$ (Eq. 13) and $\alpha_{h_{Na}}$, $\beta_{h_{Na}}$, to reproduce the experimentally observed coupled left shifts of both the steady-state conductance factors $m_{Na\infty}(V)$ and $h_{Na\infty}(V)$. To obtain a similar left shift of the opening and closing rates for the inactivation gate, $\alpha_{h_{Na}}$ and $\beta_{h_{Na}}$ must be changed in an opposite manner than for the Na activation gate [Eq. (13)]:

$$ \alpha_{h_{Na}} = \alpha_{h_{Na,0}} e^{-\sigma \mathcal{B}_h}, \quad \beta_{h_{Na}} = \beta_{h_{Na,0}} e^{\sigma \mathcal{B}_h}, \qquad (15) $$

where the stretching factor $\mathcal{B}_h = 200 m^2/J$. Note that this choice is completely different from that proposed by Ass-
<figure><img src="image_9.png" /><figcaption>**Fig. 5** Dependance on the voltage *V* of the conductances and currents of a CLS HHm model for muscle fibres. **a** Probabilities of Na channel gate opening *m*<sub>Na∞</sub> (*red*), *h*<sub>Na∞</sub> (*blue*) and *m*<sub>K∞</sub> (*green*). **b** Conductance factors in the steady state for the sodium 100 *m*<sub>Na∞</sub><sup>3</sup> *h*<sub>Na∞</sub> (*red*) and the potassium *m*<sub>K∞</sub> (*green*) channels. **c** Time constants for the three gates</figcaption></figure>

*τ*<sub>mNa</sub> (*red*), *τ*<sub>hNa</sub> (*blue*) and *τ*<sub>mK</sub> (*green*). **d** Ionic currents in the steady state *I*<sub>Na∞</sub> (*red*), *I*<sub>K∞</sub> (*green*), *I*<sub>L∞</sub> (*grey*) and their sum *I*<sub>S∞</sub> (*black*). The symbols have the following meanings: *σ* = 0 mN/m (*solid line*), 2 mN/m (*empty circles*), 4 mN/m (*filled circles*), −2 mN/m (*empty triangles*), −4 mN/m (*filled triangles*) (colour figure online)

mann and Lenz (2014), where all the channel rates were
modulated identically.

We test the impact of such a modification on the con-
ductance factor for sodium, the Na activation gate constant
and the steady-state currents (Fig. 5) for the same values of
the tension σ as previously considered with the SLS HHm
model (Fig. 3). In the CLS HHm model, the introduction
of a positive tension σ leads to a concerted left shift of the
Na activation and inactivation rates (Fig. 5a), as observed
in experiments (Morris et al. 2012). Interestingly, the con-
ductance factor for the Na channel not only shifts to the left
(lower voltage values) but also significantly decreases for
σ > 0 (Fig. 5b), whereas the time constant for the inactiva-
tion gate increases (Fig. 5c). All these evidences should lead
to a further extension of the depolarised state, favouring the
activation gating mechanism at the expense of inactivation.
The *I*<sub>S</sub>(*V*) curve (Fig. 5d) shows also a global shift of the
range of *V* values where multi-steady states exist when vary-
ing σ. Contrarily to the SLS HHm model (Fig. 2d), the *I*<sub>app</sub>
values where this occurs are restricted to a narrower band
close to zero.

The computation of the dynamical signatures of the large
amplitude oscillations with the CLS HHm model for I<sub>app</sub> =
3.95 µA/cm<sup>2</sup> (Fig. 6) shows that the range of V values is
now limited to negative values. The overall V amplitude
is decreased when both m<sub>Na</sub> and h<sub>Na</sub> are left-shifted by the
mechanical stretching energy. It is also very important to
remark that the length of the plateau of maximal m<sub>Na</sub>(t) val-
ues is dramatically reduced (Fig. 6a), meaning that the inward

flux of Na⁺ ions is diminished (Fig. 6d) in the CLS HHm
model. This phenomenon is also visible on the sharpness
of the conductance factor for Na⁺ ions which increases and
decreases much more rapidly when the excitable cell fires
(Fig. 6c). We also note that the voltage *V*(*t*) (Fig. 6a) and
the leak current *I*<sub>L</sub>(*t*) (Fig. 6d) signals have a longer slowly
decreasing tail in the CLS HHm model, as compared to the
SLS HHm model. The corresponding 2D phase plots (*I*<sub>ions</sub>,
*V*) (Fig. 6e) and (*I*<sub>Na</sub>, *I*<sub>C</sub>) (Fig. 6f) have to be compared to
the ones obtained in the original HHm model in Fig. 1e, f,
respectively.

3 Results

**3.1 Dynamical analysis of the HHm models**

Following the terminology used by Izhikevich (Izhikevich
2000, 2007), we consider an excitable cell as *quiescent* if
its membrane potential is at rest or if it exhibits only small
amplitude (subthreshold) oscillations. In terms of dynamical
system theory, these two cases correspond to a stable steady
state or to a small amplitude limit cycle attractor, respectively.
Excitability occurs when a small perturbation can drive the
system from its quiescent state to a large excursion (much
larger than the small amplitude perturbation), also called
a spike, before returning to its initial quiescent state. Such
excitable behaviour does occur when the quiescent state is
close to a bifurcation that allows the system to visit a large
<figure><img src="image_10.png" /><figcaption>**Fig. 6** Comparison of the dynamical signatures of the different variables involved in the CLS HH model. **a** Voltage $V(t/T_c)$ vs $t/T_c$, where $T_c$ is the period of the large amplitude oscillations. **b** Opening probabilities for the three gates included in the HH model $m_{Na}(t)$ (red), $h_{Na}(t)$ (blue) and $m_K(t)$ (green). **c** Conductance factors for the sodium $m_{Na}^3(t)h_{Na}(t)$ (red) and the potassium $m_K(t)$ (green) channels. In **(b)**</figcaption></figure>

and **(c)** the scaled voltage curve is also plotted for comparison. **d** Ionic currents $I_{Na}(t)$ (red), $I_K(t)$ (green), $I_L$ (grey) and their sum $I_S(t)$ (black). **e** 2D phase plots of the three ionic current signals and their sum versus V. **f** 2D phase plot of $I_{Na}(t)$ versus the capacitance charging current $I_C(t) = CdV(t)/dt$. $I_{app} = 3.95 \ \mu A/cm^2$, $\sigma = 2.3 \ mN/m$ (colour figure online)

amplitude periodic pseudo-orbit passing near the quiescent state. The type of bifurcations that occur nearby the quiescent state actually defines the nature of the cell excitability (Izhikevich 2000). When the quiescent state is near a saddle-node on invariant circle bifurcation, the duration of the spike may be arbitrary long and it is rather easy to define a threshold voltage above which excitability occurs for an initial resting state. Below the threshold, the channel features damp the perturbations and relax the cell to quiescence. Above the threshold, the sharper the perturbation in time, the quicker the spike occurrence. Izhikevich (2000, 2007) qualified this type of behaviour as an “integrator”. When the quiescent state is near a Hopf bifurcation, the cell does not necessarily respond with a single spike, nor it has a well-defined threshold for spiking and it can even fire for a subthreshold perturbation. Such a situation is qualified as “resonator” because the cell responds preferentially to precise perturbation frequencies and increasing the perturbation frequency ultimately leads to a loss of cell excitability. In this section we use both analytical and numerical tools to locate these two types of bifurcations in the parameter space ($I_{app}$, $\sigma$) and to evaluate the impact of the mechanical stress upon cell excitability in the SLS and CLS HHm models described

above. The existence of stable periodic firing (large amplitude spikes that repeat indefinitely) are also characterised and circumscribed and the evolution of their period with $I_{app}$ and $\sigma$ are compared in the SLS and CLS HHm models.

### 3.1.1 Computation of the steady states and their stability

A dynamical system analysis (Wiggins 1988; Guckenheimer and Labouriau 1993; Govaertz 2000a, b) of the HHm model with or without mechanical stress requires that we first compute analytically the steady states ($V_S$, $m_{Na_S}$, $h_{Na_S}$, $m_{K_S}$) from the four ODEs of the HHm model [Eqs. (2) and (3)], by solving [$dV/dt = 0$, $dm_{Na}/dt = 0$, $dh_{Na}/dt = 0$, $dm_K/dt = 0$].

The local stability of these steady solutions is obtained by linearising the four ODES of the HHm model around these steady solutions

$$ \begin{pmatrix} dV/dt \\ dm_{Na}/dt \\ dh_{Na}/dt \\ dm_K/dt \end{pmatrix} = \begin{pmatrix} J_{11} & J_{12} & J_{13} & J_{14} \\ J_{21} & J_{22} & J_{23} & J_{24} \\ J_{31} & J_{32} & J_{33} & J_{34} \\ J_{41} & J_{42} & J_{43} & J_{44} \end{pmatrix} \begin{pmatrix} V \\ m_{Na} \\ h_{Na} \\ m_K \end{pmatrix}. \quad (16) $$
The Jacobian matrix in Eq. (16) (matrix of first derivatives) is then diagonalised analytically, and its eigenvalues are computed numerically by solving a characteristic fourth degree equation. The corresponding eigenvectors are also computed numerically. The saddle-node (fold) and Hopf bifurcations are codimension one bifurcations corresponding, respectively, to a simple real eigenvalue ($\lambda$) or a pair of complex conjugate eigenvalues ($\lambda \pm i\omega(\omega > 0)$) that cross the imaginary axis ($\lambda = 0$). The saddle-node bifurcation corresponds to a collision between a stable and an unstable steady states leading to the disappearance of the former when changing $I_{app}$. The Hopf bifurcation corresponds to the emergence of a periodic orbit from the destabilisation of a stable steady state.

### 3.1.2 Steady states and their stability in the HHm model without tension

Hopf ($H_1$ and $H_2$) and saddle-node (limit points $LP_1$ and $LP_2$) bifurcations are reported in the bifurcation diagram of the HHm model without tension shown in Fig. 7. We notice that with the original HHm model, three steady states do exist over a $I_{app}$ interval bounded by the two saddle-node bifurcation points $LP_2$ and $LP_1$ (Fig. 7a). When increasing $I_{app}$ from 0, the lower-branch steady state exhibits a subcritical Hopf bifurcation $H_1$ giving rise to an unstable periodic orbit (brown vertical dashed-dotted line) for $I_{app} < I_{H_1}$ which disappears via a homoclinic bifurcation (homoclinic to hyperbolic saddle $HHS_3^*$) corresponding to a collision with

the middle-branch hyperbolic saddle steady state (Fig. 7b). The upper-branch steady state also exhibits a subcritical Hopf bifurcation $H_2$ (Fig. 7a) from which emerges an unstable periodic cycle ($I_{app} > I_{H_2} = 8.745$ µA/cm$^2$) (pink vertical dashed-dotted line). Global bifurcation like homoclinic bifurcations or bifurcations of periodic orbits cannot be computed from the linearised ODEs but indeed requires to solve numerically the whole HHm ODEs. Alternatively, we used Matlab ODE solver libraries (ODE15s) and MATCONT (Dhooge et al. 2003) which provides a graphical interface for interactive numerical study of dynamical systems. The overall bifurcation analysis was based on continuation methods, tracing out solution manifolds of various dynamical system states while varying our model parameters $I_{app}$ and $\sigma$. As shown in Fig. 7a, the numerical bifurcation analysis reveals the existence of large amplitude periodic oscillations (pink vertical solid line) over a rather wide range of $I_{app}$ values: $I_{LPC_1} < I_{app} < I_{LPC_2}$. When increasing $I_{app}$, this periodic firing (Izhikevich 2000, 2007) originates from a saddle-node bifurcation of cycle (limit point bifurcation of cycle $LPC_1$). Whereas the new born unstable periodic cycle is going to disappear via a homoclinic collision with the middle-branch hyperbolic saddle steady state at $I_{HHS_3} > I_{LPC_1}$ (Fig. 7b), the stable large amplitude periodic cycle develops up to rather high $I_{app}$ values and finally disappears via a saddle-node bifurcation of cycle at $I_{app} = I_{LPC_2}$ when colliding the unstable periodic cycle originating from the subcritical Hopf bifurcation $H_2$ of the upper branch ($I_{H_2} < I_{LPC_2}$) (Fig. 7a).

<figure><img src="image_11.png" /><figcaption>**Fig. 7** Bifurcation diagram of the original HHm model without tension. Codimension one Hopf ($H$) and saddle-node (fold—LP) bifurcations of steady states are reported together with saddle-node bifurcation of cycle (LPC). The *thick brown* (resp. *green*) lines correspond to stable focus (resp. stable steady state) solutions. The *dashed light brown* (resp. *green*) lines correspond to unstable focus (resp. hyperbolic saddle steady state) solutions. The transitions points $T_0$ (resp. $T_2$) correspond to neutral saddle points where two complex eigenvalues with positive (resp. negative) real parts become real. The *pink vertical solid lines* illustrate the large amplitude stable limit cycle. **b** Zoom in (**a**) showing the two</figcaption></figure>

unstable limit cycles arising from either the saddle-node bifurcation of cycle $LPC_1$ (*pink vertical dashed-dotted line*) or the subcritical Hopf bifurcation $H_2$ of the lower-branch state (*brown vertical dashed-dotted line*). Both these unstable cycles disappear via homoclinic bifurcations $HHS_3$ and $HHS_3^*$ by colliding with the middle-branch hyperbolic steady state (*green dashed line*). $I_{H_1} = 7.050$ µA/cm$^2$, $I_{H_2} = 8.745$ µA/cm$^2$, $I_{LPC_1} = 7.181$ µA/cm$^2$, $I_{LPC_2} = 4.335$ µA/cm$^2$, $I_{T_0} = 4.728$ µA/cm$^2$, $I_{T_1} = 7.176$ µA/cm$^2$, $I_{T_2} = 4.564$ µA/cm$^2$, $I_{HHS_3} = 6.914$ µA/cm$^2$, $I_{HHS_3^*} = 6.98$ µA/cm$^2$ (colour figure online)
<figure><img src="image_12.png" /><figcaption>**Fig. 8** Steady-state manifolds computed from the HHm models with tension. **a** SLS HHm model. **b** CLS HHm model. The *dark* (resp. *light*) *green areas* correspond to stable (resp. hyperbolic saddle) steady states. The *dark* (resp. *light*) *brown areas* correspond to stable (resp. unstable) steady-state solutions.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

foci. The boundaries of the *dark* and *light brown areas* correspond to Hopf bifurcation points. The boundaries of the *dark* (resp. *light*) *red* and *green* domains correspond to neutral saddle steady states (colour figure online)

### 3.1.3 Global change of the steady-state manifolds with tension

The system sensitivity and displacement of the bifurcation points (H, LP, LPC) when introducing tension in the SLS and CLS HHm models are shown in a 3D representation ($I_{app}$, $\sigma$, $V$) in Fig. 8a, b, respectively. Increasing the parameter $\sigma > 0$ leads to significantly different changes in the original HHm bifurcation diagram (Fig. 7), whether only $m_{Na}$ is sensitive to tension (SLS) or whether both $m_{Na}$ and $h_{Na}$ are sensitive to tension (CLS). For instance, the range of $I_{app}$ values where the three steady states coexist ($I_{LP1} \le I_{app} \le I_{LP2}$) widens (resp. shrinks) in the SLS HHm model for positive (resp. negative) values of $\sigma$ (Fig. 8a). This interval of $I_{app}$ values does not increase for positive $\sigma$ values but rather stabilizes to a fixed width that drifts towards negative values for the CLS HHm model (Fig. 8b). Another feature which changes from the original HHm model is the range of $I_{app}$ values visited by the boundary between brown (stable focus) and light brown (unstable focus) areas corresponding to the Hopf bifurcation ($H_1$) of the lower-branch steady state. Consistent with its proximal saddle-node bifurcation $LP_1$, $H_1$ extends to higher $I_{app}$ values when decreasing $\sigma$ to negative values. Since, as discussed above, Hopf and saddle-node bifurcations of steady states are favourable situations for excitability, we expect that it will be observed over a much wider range of $I_{app}$ values in the SLS and CLS HHm models by varying $\sigma$.

## 3.2 Excitability of the HHm models with tension

To investigate the excitability in the two-dimensional parameter space ($I_{app}$, $\sigma$) of the SLS and CLS mechanical variants of the HHm model, we proceed through the following numerical strategy. For each parameter pair ($I_{app}$, $\sigma$), the

steady-state solutions of the four ODEs are computed first and the one with the lower $V_S$ value is chosen (lower-branch when multiple steady states). This steady state $S_S[V_S, m_{Na_S}, h_{Na_S}, m_{K_S}]$ is perturbed by +40 mV on $V_S$ and used as initial condition for numerical integration of the ODEs equations with the Matlab ODE solver ODE15s. When the dynamics displays a single large amplitude oscillation (i.e. larger than the perturbation) before returning to the steady state, then the steady state $S_S$ is considered as excitable. Of course the excitability diagrams so obtained are not strictly exact but qualitative in the sense that their boundaries are dependent on the amplitude of the perturbation; the larger the perturbation amplitude the wider the expected excitability domain. We choose the perturbation amplitude larger than the voltage drop between the upper and lower steady states but smaller than the amplitude of the large limit cycle (pink vertical lines in Fig. 7).

The excitability diagram of the SLS HHm model (Fig. 9) shows that the system keeps its excitability when increasing $\sigma > 0$ and that the amplitude of the single spikes also increases with $\sigma$. The $I_{app}$ interval where excitability occurs shifts to lower (more negative) values of $I_{app}$, following nicely the Hopf bifurcation $H_1$ curve (blue dots) of the lower branch (right side of the excitable domain). This diagram illustrates the excitability of the resting state (corresponding to low $V$ values) specifically and not of the upper-branch steady state which does not correspond to realistic experimental situations. Interestingly, the voltage amplitude of the excitability spike does not depend on the amplitude of the perturbation and is larger at larger distance from the $H_1$ Hopf bifurcation line. This comes along with a widening to lower negative $I_{app}$ values of the multi-steady-state window (Fig. 9a). The rightmost bottom part of the excitability domain (larger $I_{app}$ and smaller $\sigma$ values) corresponds to single spikes with smaller voltage amplitude,
tion lines on the lower steady-state fold. The fact that the two
saddle-node LP₁ and LP₂ lines (green) remain very close to
each other in *I*<sub>app</sub> (they look quite parallel for *σ* > 0) is the
signature of a smaller voltage drop between the upper and
lower steady-state folds, which explains that the amplitude
of the excitability spike in the CLS HHm model is smaller
than in the SLS HHm model. However, the domain of exis-
tence of the resting state excitability has not shrunk. On
the contrary, the Hopf bifurcation H₁ on the lower branch
(H₁) persists for *σ* < 0 and the CLS HHm model remains
excitable for much more negative values of *σ* than the SLS
HHm model. Importantly, in the CLS HHm model, the volt-
age spike (Fig. 10b-g) is triggered very fast, independently of
the special location in the excitability domain. In that sense,
the excitability is stronger in the CLS HHm model than in
the SLS HHm model. If a Na<sub>v</sub> channel behaves as predicted
by the CLS HHm model, we can thus expect that it will be
more sensitive to its neighbouring channels and that the cell
will be globally more excitable.

3.3 Periodic firing in the HHm models with tension

To further understand and predict the dynamics of injured
muscle cells produced by some mechanical alteration of
the Na<sub>v</sub> channels gating kinetics, we analyse the range of
parameters where periodic firing (large amplitude oscilla-
tions) occurs. The large amplitude limit cycles of the two SLS
and CLS mechanical variants of the HHm model are com-
puted numerically in the two-dimensional parameter space
(I<sub>app</sub>, σ). For each parameter pair (I<sub>app</sub>, σ), the steady-state
solutions of the four ODEs are computed first and the one
with the largest V<sub>S</sub> value is chosen (in multi-steady-state
situations). When increasing I<sub>app</sub>, the large amplitude limit
cycle emerges from a saddle-node bifurcation of cycle, as
previously observed in the original HHm model without ten-
sion (Fig. 7). We use continuation methods to follow this
limit cycle that could also be reached from stable steady
state via a large amplitude perturbation. To illustrate what
could be performed in experimental situations by a perturba-
tion of the channel voltage, we choose the second alternative
to build these diagrams. We perturb the steady state S<sub>S</sub>[V<sub>S</sub>,
m<sub>Nas</sub>, h<sub>Nas</sub>, m<sub>Ks</sub>] by a negative drop (−50 mV) in V<sub>S</sub> and use
it as initial condition for numerical integration of the ODEs
equations with the Matlab ODE solver ODE15s. Choosing
a negative voltage drop (inversely to the excitability study
just above) is important because it is directed towards the
resting state voltage level. Again it is important to choose a
perturbation amplitude in V large enough to get out from the
attraction basin of the stable focus that exists on the upper
branch. Note that if we perturb the stable large amplitude
oscillations, depending on the direction and amplitude of the
perturbation the system may return back to the upper-branch
steady state (at high voltage).

4 Discussion

Periodic firing could provide some understanding of what
is observed in hyperPP which is characterised by muscle
hyperexcitability or weakness which can lead to uncontrolled
shaking followed by paralysis. In the SLS HHm model,
the window of periodic firing in I<sub>app</sub> starts widening when
increasing σ > 0 and shifts to lower I<sub>app</sub> values (Fig. 11a).
Periodic firing is likely to occur more often when the cell is
under tension. If the muscle cells fire periodically, both Na<sup>+</sup>
and K<sup>+</sup> ion inward fluxes increase. Moreover, if the channels
which have stepped to large amplitude oscillations are again
perturbed, they may stall on a stable upper-branch steady
state which is different from their normal (lower-branch) rest-
ing state, without being able to fire or to recover their initial
state, mimicking a sort of paralysis. The SLS HHm model
demonstrates that if the activation gate of the Na<sub>v</sub> channel
is speeded up by tension, not only the cells become more
excitable but they can also fire periodically at I<sub>app</sub> values
where they would have relaxed to the resting state in normal
situations.

The periodic firing diagram of the SLS HHm model (Fig. 11a) shows that the amplitude of periodic oscillations increases when increasing $\sigma > 0$, the larger ones being observed at the left border (LPC₁) of the periodic firing domain. When decreasing $\sigma$ to negative values, the oscillation amplitude diminishes progressively to become subthreshold (Fig. 11f) due to the change from subcritical to supercritical of the original Hopf (H) bifurcation and finally the oscillations disappear. For $\sigma \le -2$ mN/m, the system no longer fires and remains in its resting state. This characteristic passivity of the SLS HHm model for slightly negative $\sigma$ values is a hallmark of this system. For instance, if a SLS HHm-type cell has encountered a hyperkalaemic periodic paralysis phenomenon at higher tension and is stalled in a non-resting state, lowering the tension by relaxing the muscle cell will likely help it to recover a normal resting state, driving the whole system back or close to normal situation. In the diagram of Fig. 11b where the period of the large amplitude oscillations is colour coded, we note that this period does not change much and remains in the interval 5 to 20 ms, except in the neighbourhood of the left LPC₁ saddle-node of cycle bifurcation line where it increases sharply.

The periodic firing diagram of the CLS HHm model (Fig. 12) is very different from the one of the SLS HHm model (Fig. 11). If the periodic firing domain is still delimited by the LPC saddle-node of cycle bifurcation lines, it no longer develops very far to positive $\sigma$ values, but mainly expands towards negative $\sigma$ values (Fig. 12a). Moreover, the large amplitude oscillations are rather observed at negative $\sigma$ values and the $I_{app}$ interval where they occur spreads dramatically to larger $I_{app}$ values. Contrarily to the SLS HHm model, when switching the system to slightly nega-
The dynamical system analysis proposed here opens new perspectives to physiologists as regards the possibility of distinguishing those of the Na<sub>v</sub> gates which are mostly impacted by a mechanical damage. It should therefore be interesting to reconstruct excitability and/or periodic firing domains (such as those reported in this work) from experiments (e.g. modulating the cell tension) to estimate the level of mechanical damage (percentage of impacted activation and inactivation gates) and to propose an adapted treatment. If we believe that the CLS HHm model is relevant for myotonia disorders, then better not to try to relax the muscle but rather to increase its tension to take it away from the periodic firing regime.

**Acknowledgements** The authors gratefully acknowledge financial support from China Scholarship Council. This project was supported by the National Science Foundation of China (81590953, 81574053, 8157150277), the National Program on Key Basic Research Project (973 Program) (2012CB518502), the Shanghai Key Laboratory of Acupuncture Mechanism and Acupoint Function (14DZ2260500) and the National Science Foundation of Shanghai (15441903800). This work was also supported by the Centre National de la Recherche Scientifique, l'Ecole Normale Superieure of Lyon and the University of Bordeaux.

## References

Adrian R, Chandler W, Hodgkin A (1970) Slow changes in potassium permeability in skeletal muscle. J Physiol 208:645-668

Assmann MA, Lenz P (2014) Membrane tension influences the spike propagation between voltage-gated ion channel clusters of excitable membranes. Phys Biol 11(4):046,006

Bainbridge FA (1915) The influence of venous filling upon the rate of the heart. J Physiol 50(2):65-84

Banderali U, Juranka PF, Clark RB, Giles WR, Morris CE (2010) Impaired stretch modulation in potentially lethal cardiac sodium channel mutants. Channels 4(1):12-21

Bett GCL, Sachs F (1997) Cardiac mechanosensitivity and stretch-activated ion channels. Trends Cardiovasc Med 7(1):4-8

Beyder A, Rae JL, Bernard C, Strege PR, Sachs F, Farrugia G (2010) Mechanosensitivity of Nav1.5, a voltage-sensitive sodium channel. J Physiol 588(Pt 24):4969-4985

Brackenbury W (2012) Voltage-gated sodium channels and metastatic disease. Channels 6(5):352-361

Brackenbury WJ, Djamgoz MBA, Isom LL (2008) An emerging role for voltage-gated Na channels in cellular migration: regulation of central nervous system development and potentiation of invasive cancers. Neuroscientist 14(6):571-583 NIHMS150003

Butchal F, Sten-Knudsen O (1959) Impulse propagation in striate muscle fibers and the role of the internal current in activation. Ann N Y Acad Sci 81:422-445

Calabrese B, Tabarean IV, Juranka P, Morris CE (2002) Mechanosensitivity of N-type calcium channel currents. Biophys J 83(November):2560-2574

Cannon RC, O'Donnell C, Nolan M (2010) Stochastic ion channel gating in dendritic neurons: morphology dependence and probabilistic synaptic activation of dendritic spikes. PLoS Comput Biol 6(8):e1000886

Cannon SC, Corey DP (1993) Loss of Na+ channel inactivation by anemone toxin (ATX II) mimics the myotonic state in hyperkalaemic periodic paralysis. J Physiol 466:501-520

Cannon SC, Brown RH, Corey DP (1991) A sodium channel in hyperkalaemic periodic paralysis:potassium-induced failure of inactivation. Neuron 6(4):619-626

Cannon SC, Brown RH, Corey DP (1993) Theoretical reconstruction of myotonia and paralysis caused by incomplete inactivation of sodium channels. Biophys J 65(1):270-288

Casado M, Ascher P (1998) Opposite modulation of NMDA receptors by lysophospholipids and arachidonic acid: common features with mechanosensitivity. J Physiol 513(2):317-330

Catterall WA (2012) Voltage-gated sodium channels at 60: structure, function and pathophysiology. J Physiol 590(11):2577-2589

Davies NP, Hanna MG (2001) The skeletal muscle channelopathies: basic science, clinical genetics and treatment. Curr Opin Neurol 14(5):539-551

Despa S, Vigmond E (2016) From single myocyte to whole heart. Circ Res 118(2):184-186

Dhooge A, Govaerts W, Kuznetsov YA (2003) MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans Math Softw 29:141-164

Fenn W, Cobb D (1936) Electrolyte change in muscle during activity. Am J Physiol 115:345-356

Govaerts W (2000a) Numerical bifurcation analysis for OEDs. J Comput Appl Math 125:57-68

Govaerts W (2000b) Numerical methods for bifurcation of dynamical systems. SIAM, Philadelphia

Gribkoff VK, Kaczmarek LK (2009) Structure, function and modulation of neuronal voltage-gated ion channels. Wiley, Hoboken

Gu CX, Juranka PF, Morris CE (2001) Stretch-activation and stretch-inactivation of shaker-IR, a voltage-gated K+ channel. Biophys J 80(6):2678-2693

Guckenheimer J, Labouriau IS (1993) Bifurcations in Hodgkin and Huxley equations: a new twist. Bull Math Biol 55(5):937-952

Guharay BYF, Sachs F (1984) Stretch-activated single ion channel currents in tissue-cultured embryonic chick skeletal muscle. J Physiol 352:685-701

Hamill OP (2006) Twenty odd years of stretch-sensitive channels. Pflugers Arch Eur J Physiol 453(3):333-351

Hamill OP, Marty A, Neher E, Sakmann B, Sigworth FJ (1981) Improved patch-clamp techniques for high-resolution current recording from cells and cell-free membrane patches. Pflugers Arch Eur J Physiol 391:85-100

Helfrich W (1973) Elastic properties of lipid bilayers: theory and possible experiments. Z Naturforsch 28c:693-703

Hille B (1992) Ion channels of excitable membranes, 2nd edn. Sinauer Associates, sunderland

Hodgkin AL (1951) The ionic basis of electrical activity in nerve and muscle. Biol Rev 26(4):339-409

Hodgkin AL, Horowicz P (1959) Movements of Na and K in single muscle fibres. J Physiol 148(2):405-432

Hodgkin AL, Huxley AF (1945) Resting and action potentials in single nerve fibres. J Physiol 104(2):176-195

Hodgkin AL, Huxley AF (1952a) A quantitative description of membrane current and its application to conduction and excitation of nerve. J Physiol 117:500-544

Hodgkin AL, Huxley AF (1952b) Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo. J Physiol 116:449-472

Hu H, Sachs F (1997) Stretch-activated ion channels in the heart. J Mol Cell Cardiol 29(6):1511-1523

Huxley AF (1959) Local activation of muscle. Ann N Y Acad Sci 81:446-452

Izhikevich EM (2000) Neural excitability, spiking and bursting. Int J Bifurc Chaos 10(6):1171-1266

Izhikevich EM (2007) Dynamical systems in neuroscience. MIT Press, Cambridge
Jarecki BW, Piekarz AD, Jackson JO, Cummins TR (2010) Human voltage-gated sodium channel mutations that cause inherited neuronal and muscle channelopathies increase resurgent sodium currents. *J Clin Investig* 120(1):369-378

Kaczmarek LK (2006) Non-conducting functions of voltage-gated ion channels. *Nat Rev Neurosci* 7(10):761-771

Kohl P, Hunter P, Noble D (1999) Stretch-induced changes in heart rate and rhythm: clinical observations, experiments and mathematical models. *Prog Biophys Mol Biol* 71(1):91-138

Laitko U (2004) Membrane tension accelerates rate-limiting voltage-dependent activation and slow inactivation steps in a shaker channel. *J Gen Physiol* 123(2):135-154

Laitko U (2006) Membrane stretch slows the concerted step prior to opening in a Kv channel. *J Gen Physiol* 127(6):687-701

Lin W, Laitko U, Juranka PF, Morris CE (2007) Dual stretch responses of mHCN2 pacemaker channels: accelerated activation, accelerated deactivation. *Biophys J* 92(5):1559-1572

Liu J, Schrank B, Waterston RH (1996) Interaction between a putative mechanosensory membrane channel and a collagen. *Science* 273(5273):361-364

Maingret F, Patel AJ, Lesage F, Lazdunski M, Honoré E (2000) Lysophospholipids open the two-pore domain mechano-gated K+ channels TREK-1 and TRAAK. *J Biol Chem* 275(14):10,128-10,133

Markin VS, Sachs F (2004) Thermodynamics of mechanosensitivity. *Phys Biol* 1:110-124

Morris CE (1990) Mechanosensitive channels. *J Membr Biol* 107:93-107

Morris CE, Juranka PF (2007) Nav channel mechanosensitivity: activation and inactivation accelerate reversibly with stretch. *Biophys J* 93(3):822-833

Morris CE, Boucher PA, Joos B (2012) Left-shifted Nav channels in injured bilayer: primary targets for neuroprotective Nav antagonists? *Front Pharm* 3(19):1-13

Morris CE, Prikryl EA, Joos B (2015) Mechanosensitive gating of Kv channels. *PLoS ONE* 10(2):1-24

Namadurai S, Yereddi NR, Cusdin FS, Huang CLH, Chirgadze DY, Jackson AP (2015) A new look at sodium channel beta subunits. *Open Biol* 5(1):140,192

Pawson L, Bolanowski SJ (2002) Voltage-gated sodium channels are present on both the neural and capsular structures of Pacinian corpuscles. *Somatosens Mot Res* 19(3):231-237

Peyronnet R, Nerbonne JM, Kohl P (2016) Cardiac mechano-gated ion channels and arrhythmias. *Circ Res* 118(2):311-329

Ptacek LJ, Johnson KJ, Griggs R (1993) Genetic and physiology of the myotonic muscle disorders. *N Engl J Med* 328(7):482-489

Sachs F (1988) Mechanical transduction in biological membranes. *CRC Crit Rev Biomed Eng* 16(2):141-169

Sachs F (1989) Ion channels as mechanical transducers. In: Stein WD, Bronner F (eds) Cell shape: determinants, regulation and regulatory role. *Academic Press*, San Diego, pp 63-92

Sackin H (1995) Mechanosensitive channels. *Annu Rev Physiol* 57:333-353

Sakmann B, Neher E (2009) Single channel recording, 2nd edn. Springer, New York

Shcherbatko A, Ono F, Mandel G, Brehm P (1999) Voltage-dependent sodium channel function is regulated through membrane mechanics. *Biophys J* 77(4):1945-1959

Simkin D, Bendahhou S (2011) Skeletal muscle Na+ channel disorders. *Front Pharm* 2(63):1-10

Spillane J, Kullmann DM, Hanna MG (2016) Genetic neurological channelopathies: molecular genetics and clinical phenotypes. *J Neurol Neurosurg Psychiatry* 87:37-48

Tabarean IV, Juranka P, Morris CE (1999) Membrane stretch affects gating modes of a skeletal muscle sodium channel. *Biophys J* 77(2):758-774

Ulbricht W (2005) Sodium channel inactivation: molecular determinants and modulation. *Physiol Rev* 85:1271-1301

Wang JA, Lin W, Morris T, Banderali U, Juranka PF, Morris CE (2009) Membrane trauma and Na+ leak from Nav1.6 channels. *Am J Physiol Cell Physiol* 297(4):C823-C834

Wiggins S (1988) Global bifurcations and chaos. *Analytical methods. Springer, New York*

Yu N, Morris CE, Joós B, Longtin A (2012) Spontaneous excitation patterns computed for axons with injury-like impairments of sodium channels and Na/K pumps. *PLoS Comput Biol* 8(9):e1002,664