Q X Ma

• A Arneodo

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

• F Argoul

Dynamical study of Na v channel excitability under mechanical stress

Keywords: Na v voltage-gated channels, Excitable cell disorders, Mechanical stress, Sick-excitable cell, Periodic firing and myotonia, Dynamical systems, Bifurcations

Alteration of Na v channel functions (channelopathies) has been encountered in various hereditary muscle diseases. Na v channel mutations lead to aberrant excitability in skeletal muscle myotonia and paralysis. In general, these mutations disable inactivation of the Na v channel, producing either repetitive action potential firing (myotonia) or electrical dormancy (flaccid paralysis) in skeletal muscles. 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 fraction 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 paralysis 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.

Introduction

Since the mid-1990s, voltage-gated sodium (Na v ) channels have been recognised for their implication in the generation and conduction of action potentials in excitable cells [START_REF] Hille | Ion channels of excitable membranes, 2nd edn. Sinauer Associates, sunderland Hodgkin AL (1951) The ionic basis of electrical activity in nerve and muscle[END_REF]. Na v channels are macromolecular protein complexes containing pore-forming α subunits and smaller non-poreforming β subunits which are embedded inside the plasma membrane. The vertebrate Na v 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 v 1.1-1.9, and an atypical channel Na x . The β subunits regulate channel expression and gating; these cell adhesion molecules (CAMs) belong to the immunoglobulin (Ig) superfamily [START_REF] Namadurai | A new look at sodium channel beta subunits[END_REF]. 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 [START_REF] Kaczmarek | Non-conducting functions of voltage-gated ion channels[END_REF]. Thus far, the non-conducting role of voltage-gated channels via the β subunits has been best characterised in cell adhesion [START_REF] Brackenbury | An emerging role for voltage-gated Na channels in cellular migration: regulation of central nervous system development and potentiation of invasive cancers[END_REF]. 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 [START_REF] Ulbricht | Sodium channel inactivation: molecular determinants and modulation[END_REF]. This inactivation mechanism is necessary for the channel to recover its initial excitability (Gribkoff and Kaczmarek 2009). Na v channels are not only expressed in excitable cells like nerve and muscle cells but also in cells ordinarily considered as non-excitable, including glial cells, fibroblasts, vascular endothelial, immune and metastatic cancer cells [START_REF] Brackenbury | An emerging role for voltage-gated Na channels in cellular migration: regulation of central nervous system development and potentiation of invasive cancers[END_REF]. Neuronal sodium channels Na v 1.1 -1.3 and Na v 1.6 -1.9 are found in the central nervous system (CNS) and peripheral nervous system (PNS) neurons and within glial cells, whereas Na v 1.4 and Na v 1.5 are expressed in skeletal and cardiac muscle cells, respectively [START_REF] Catterall | Voltage-gated sodium channels at 60: structure, function and pathophysiology[END_REF].

Channelopathies are caused by genetic or acquired defects in ion channels, the former being the most common. Voltagegated channel defects are involved in a wide variety of diseases like epilepsy, migraine, blindness, deafness [START_REF] Spillane | Genetic neurological channelopathies: molecular genetics and clinical phenotypes[END_REF], diabetes, hypertension, cardiac arrhythmia, asthma, irritable bowel syndrome, and cancer [START_REF] Brackenbury | Voltage-gated sodium channels and metastatic disease[END_REF]. Na v channel mutations are associated to aberrant excitability in skeletal muscle myotonia and paralysis [START_REF] Jarecki | Human voltage-gated sodium channel mutations that cause inherited neuronal and muscle channelopathies increase resurgent sodium currents[END_REF]. These channelopathies include hyperkalaemic periodic paralysis (hyperPP), potassium aggravated myotonia (PAM), paramyotonia congenita (PMC) and hypokalaemic periodic paralysis (hypoPP) [START_REF] Davies | The skeletal muscle channelopathies: basic science, clinical genetics and treatment[END_REF][START_REF] Jarecki | Human voltage-gated sodium channel mutations that cause inherited neuronal and muscle channelopathies increase resurgent sodium currents[END_REF][START_REF] Simkin | Skeletal muscle Na+ channel disorders[END_REF]. Hyperkalaemic periodic paralysis, paramyotonia congenita and potassium-aggravated myotonia represent myotonic disorders due to episodic membrane hyperexcitability of skeletal muscles [START_REF] Ptacek | Genetic and physiology of the myotonic muscle disorders[END_REF]). These diseases are all caused by mutations in the sodium channel 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 v channel, producing 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 [START_REF] Markin | Thermodynamics of mechanosensitivity[END_REF]. For instance, membrane stretching not only makes the bilayer thinner, but also decreases the amplitude of its fluctuations 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 voltagegated 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 [START_REF] Hu | Stretch-activated ion channels in the heart[END_REF]. It was early recognised that distending the right heart atria could increase the heart rate [START_REF] Bainbridge | The influence of venous filling upon the rate of the heart[END_REF]. But it took almost a century [START_REF] Hu | Stretch-activated ion channels in the heart[END_REF][START_REF] Peyronnet | Cardiac mechano-gated ion channels and arrhythmias[END_REF] to distinguish between mechanically modulated and mechano-gated ion channels and their respective impacts on cardiac electrophysiology. For instance, understanding heart rhythm disorder regulation at the entire organ level from the knowledge of single cell excitability remains a subject of intensive research [START_REF] Bett | Cardiac mechanosensitivity and stretchactivated ion channels[END_REF][START_REF] Kohl | Stretch-induced changes in heart rate and rhythm: clinical observations, experiments and mathematical models[END_REF][START_REF] Despa | From single myocyte to whole heart[END_REF]. 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 phylogenic 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) [START_REF] Liu | Interaction between a putative mechanosensory membrane channel and a collagen[END_REF]. Stretch-inactivated channels (SICs) are less common, but have also been reported in several cell types [START_REF] Sachs | Ion channels as mechanical transducers[END_REF][START_REF] Morris | Mechanosensitive channels[END_REF][START_REF] Sackin | Mechanosensitive channels[END_REF].

A single channel patch-clamp recording [START_REF] Hamill | Improved patch-clamp techniques for high-resolution current recording from cells and cell-free membrane patches[END_REF][START_REF] Sakmann | Single channel recording[END_REF] 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 voltagegated [START_REF] Laitko | Membrane tension accelerates rate-limiting voltagedependent activation and slow inactivation steps in a shaker channel[END_REF][START_REF] Laitko | Membrane stretch slows the concerted step prior to opening in a Kv channel[END_REF][START_REF] Lin | Dual stretch responses of mHCN2 pacemaker channels: accelerated activation, accelerated deactivation[END_REF][START_REF] Morris | Nav channel mechanosensitivity: activation and inactivation accelerate reversibly with stretch[END_REF] or ligand-gated [START_REF] Casado | Opposite modulation of NMDA receptors by lysophospholipids and arachidonic acid: common features with mechanosensitivity[END_REF][START_REF] Maingret | Lysophospholipids open the two-pore domain mechano-gated K+ channels TREK-1 and TRAAK[END_REF] are also mechanically sensitive. The only requirement 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 [START_REF] Markin | Thermodynamics of mechanosensitivity[END_REF].

Na v , K v , Ca v and hyperpolarisation-activated cyclic nucleotide-gated (HCN) channels were all shown to exhibit reversible gating changes with stretch [START_REF] Shcherbatko | Voltage-dependent sodium channel function is regulated through membrane mechanics[END_REF][START_REF] Tabarean | Membrane stretch affects gating modes of a skeletal muscle sodium channel[END_REF]Gu et al. 2001;[START_REF] Calabrese | Mechanosensitivity of N-type calcium channel currents[END_REF][START_REF] Lin | Dual stretch responses of mHCN2 pacemaker channels: accelerated activation, accelerated deactivation[END_REF]). The impact of membrane deformations upon the Na v channel activity was suspected since the late 1990s [START_REF] Shcherbatko | Voltage-dependent sodium channel function is regulated through membrane mechanics[END_REF][START_REF] Tabarean | Membrane stretch affects gating modes of a skeletal muscle sodium channel[END_REF][START_REF] Pawson | Voltage-gated sodium channels are present on both the neural and capsular structures of Pacinian corpuscles[END_REF] and was recently elucidated by voltage-clamp experiments [START_REF] Hamill | Twenty odd years of stretch-sensitive channels[END_REF][START_REF] Morris | Nav channel mechanosensitivity: activation and inactivation accelerate reversibly with stretch[END_REF][START_REF] Beyder | Mechanosensitivity of Nav1.5, a voltage-sensitive sodium channel[END_REF]. Na v 1.5 channel reversibly responds to membrane stretch by an acceleration of the ratelimiting voltage-dependent step leading to a Na + 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 v to even small mechanical modulations. Given that minute deviations in Na v channel kinetics may lead to cardiac arrhythmias, epilepsy, neuropathic pain and general muscle disorders, it is critical to determine to which extent Na v channel kinetics 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 v channel may impair its normal function (excitability), either by driving it to a very excitable state (periodic spiking) or by blocking its excitability (paralysis).

Methods

Modelling the dynamical response of excitable cells

Characteristic features of excitable cells

Excitable cells [START_REF] Hille | Ion channels of excitable membranes, 2nd edn. Sinauer Associates, sunderland Hodgkin AL (1951) The ionic basis of electrical activity in nerve and muscle[END_REF]) 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 voltagegated 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.

General Hodgkin-Huxley (HH) formalism

The Hodgkin-Huxley (HH) models for nerve and muscle fibres (Hodgkin 1951;Hodgkin and Huxley 1952a, b;[START_REF] Hodgkin | Movements of Na and K in single muscle fibres[END_REF][START_REF] Huxley | Local activation of muscle[END_REF]) 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 Na , a potassium current I K and a background or leak current I L . Currents are given as current densities in µA/cm 2 , ionic conductances as surface conductances in mS/cm 2 and capacitances as surface capacitances in µF/cm 2 . The applied (external) current I app 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 dV dt + I Na (V, t) + I K (V, t) + I L (V, t). ( 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 L is purely ohmic with a fixed conductance g L and accounts for the passive permeability of the membrane to ions. The potassium current I K is produced by a voltage-gated channel (K v ), 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 Na because the inward flow of Na + 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 v channel is modelled by a single gating mechanism, and its conductance factor m K is a monotonous and nonlinear function of the membrane potential and gives the 2b); g K is the maximum conductance when m K = 1. The conductance of the Na v channel is modelled by the combination of two gating mechanisms: (i) a first gating mechanism (described by a variable conductance factor m Na ) similar to the K v channel which activates the permissivity of the Na v channel to Na + 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 Na . The whole conductance of the Na v channels g Na (V ) = g Na m 3

K v conductance g K (V ) = g K m 4 K (V ) (Fig.
Na (V )h Na (V ) was proposed by Hodgkin and Huxley (Hodgkin and Huxley 1952a) from experimental observations. Due to the counteracting effects of the two gates, g Na (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 Na . This combination of two gates makes the negative spike of Na + ion current I Na (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 Na (blue curve in Fig. 1b). The intersection of the two curves for m Na and h Na 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 Na (t), I K (t) and I L (t) and their sum I S (t) versus V (t). In Fig. 1f 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 dV dt = -g Na m 3 Na (V )h Na (V )(V -E Na ) -g K m 4 K (V )(V -E K ) -g L (V -E L ) + I app (t), ( 2
)
where E Na , E K and E L are the reversal potentials of the membrane channels corresponding to vanishing currents I Na , I K and I L , respectively; these potentials are the Nernst equilibrium potentials for each ion. The conductance coefficients g Na , g K and g L of the three ionic currents are constant parameters that do not depend on V . The dependance of the activation and inactivation variables m Na , m K and h Na and of the corresponding ionic conductances on the voltage V is essential ingredients of the HH models. The activation and inactivation variables m Na and h Na have opposite voltage dependencies. Depolarisation causes m Na to increase and h Na to decrease, whereas repolarisation causes m Na to decrease and h Na to increase [START_REF] Namadurai | A new look at sodium channel beta subunits[END_REF]) (Fig. 2a). In HH models, the probability of Na + permissivity is considered as equal to the product m 3 Na h Na , assuming that the two Na gates m Na ∞ and h Na ∞ act independently. These transition probabilities dm i /dt (i = 1, 2, 3) = (dm Na /dt, dh Na /dt and dm K /dt), from permissive to nonpermissive states of the different ionic gates, are assumed to obey first-order rate equations with rate factors α i (V ) and β i (V ) that depend on the membrane potential V :

dm i dt = α i (V )(1 -m i ) -β i (V )m i .
(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 (1m i ) of finding the gate closed, multiplied by the opening rate (f)

α 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 β i (V ). Dividing Eq. ( 3) by α i (V ) + β i (V ) brings out the time constants τ i for each gate:

τ i (V ) dm i dt = m i ∞ (V ) -m i , (4)
where

τ i (V ) = 1 α i (V ) + β i (V ) , m i,∞ (V ) = α i (V ) α i (V )+β i (V ) . (5)
Equation ( 4) tells us that for a fixed voltage V , m i approaches the limiting value m i,∞ (V ) exponentially with a time constant τ i (V ).

Experimental data have revealed that the m i,∞ (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).

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 voltageclamp 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;[START_REF] Hodgkin | Movements of Na and K in single muscle fibres[END_REF][START_REF] Huxley | Local activation of muscle[END_REF]. Similar to nerve action potentials, they confirmed older observations by [START_REF] Fenn | Electrolyte change in muscle during activity[END_REF] 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 [START_REF] Butchal | Impulse propagation in striate muscle fibers and the role of the internal current in activation[END_REF]. However, unlike for unmyelinated nerve fibres as the squid giant axons [START_REF] Hodgkin | Resting and action potentials in single nerve fibres[END_REF], 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 [START_REF] Adrian | Slow changes in potassium permeability in skeletal muscle[END_REF]. The opening α i and closing β i rates were adapted from the original HH equations (Hodgkin and Huxley 1952a) to match the excitability of muscle cells:

α m Na (V ) = 0.288V 1 /(1 -e -V 1 /10 ), β m Na (V ) = 1.38 e -V 1 /18 , α h Na (V ) = 8.1 10 -3 e -V 2 /14.7 , β h Na (V ) = 4.38/(1 + e -V 2 /9 ), α m K (V ) = 0.0131V 3 /(1 -e -V 3 /7 ), β m K (V ) = 0.067e -V 3 /40 , ( 6
)
where 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).

V 1 = 44 + V , V 2 = 45 + V and V 3 = 40 + V .
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 2 and g K = 21.6 mS/cm 2 , whereas the leak conductance is g L = 0.75 mS/cm 2 . The steady-state open probabilities for the activation and inactivation Na gates are combined in the total conductance factor for Na + ions:

m 3 Na ∞ h Na ∞ = α m Na (V ) α m Na (V ) + β m Na (V ) 3 × α h Na (V ) α h Na (V ) + β h Na (V ) . ( 7
)
The corresponding steady-state open probability for the K channel yields the total conductance factor for K + ions:

m 4 K ∞ = α m K (V ) α m K (V ) + β m K (V ) 4 . ( 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 ∞ , h Na ∞ and m K ∞ ) [Eq. ( 5)] and the conductance factors m 3 Na ∞ h Na ∞ [Eq. ( 7)] and m 4 K ∞ [Eq. ( 8)] are calculated numerically from Eq. ( 6) and plotted in Fig. 2a,b. Note that the Na conductance factor m 3 Na ∞ h Na ∞ 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 τ 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 ∼ -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 ∞ (V ) (red line) is smaller than I K ∞ (V ) (green line) and never crosses it (checked from numerical simulations) and both these currents vary much slower than the leak current I L ∞ 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 ∞ and I S ∞ display a sigmoidal shape versus V (Fig. 2f), whereas in the nerve fibre HH model, I S ∞ is a smoothly increasing function of V . Interestingly, the steady-state current I S ∞ (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.

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. [START_REF] Cannon | A sodium channel in hyperkalemic periodic paralysis:potassium-induced failure of inactivation[END_REF] 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 "sickexcitable cell" conditions [START_REF] Morris | Left-shifted Nav channels in injured bilayer: primary targets for neuroprotective Nav antagonists?[END_REF], corresponding to a left shift of both the steady-state conductance factors m Na ∞ (V ) and h Na ∞ (V ) when plotted versus V (Fig. 2a). This left shift can be followed by the decrease of the halfvoltage V 1/2 corresponding to β i /α 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 [START_REF] Shcherbatko | Voltage-dependent sodium channel function is regulated through membrane mechanics[END_REF][START_REF] Tabarean | Membrane stretch affects gating modes of a skeletal muscle sodium channel[END_REF][START_REF] Wang | Membrane trauma and Na+ leak from Nav1.6 channels[END_REF][START_REF] Beyder | Mechanosensitivity of Nav1.5, a voltage-sensitive sodium channel[END_REF]. 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 ∞ (V ) and h Na ∞ (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) [START_REF] Morris | Left-shifted Nav channels in injured bilayer: primary targets for neuroprotective Nav antagonists?[END_REF]. CLS leads to maxima of the Na channel conductance factor m 3 Na ∞ h Na ∞ 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 [START_REF] Izhikevich | Dynamical systems in neuroscience[END_REF]) 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 [START_REF] Assmann | Membrane tension influences the spike propagation between voltage-gated ion channel clusters of excitable membranes[END_REF], 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 nonspiking 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 ∞ ) and inactivation (h Na ∞ ) 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).

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 σ and the fractional increase of its area A/A [START_REF] Helfrich | Elastic properties of lipid bilayers: theory and possible experiments[END_REF]Guharay and Sachs 1984; Assmann and Lenz 2014):

T = σ A/A. (9)
σ represents the amount of tension required for doubling the area ( A = A). The energy of membrane deformation is:

U = A 2σ T 2 = σ 2 A 2 A . ( 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 σ =0 i→i * = k 0 e -G ‡ /kT , ( 11
)
where k is the Boltzmann's constant, T the absolute temperature and G ‡ the free energy of transition from the initial state to the activated complex. Equation ( 11) for the rate constant k σ =0 i→i * 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→i * = k σ =0 i→i * e σ 2kT A 2 A = k σ =0 i→i * e σ B , (12)
where

B = A 2
2kT A describes the amount of stretching (assumed fixed) implied by the experimental measurement.

G ‡ = G ‡ σ =0 -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 σ = σ 0 ( A/A) ζ with ζ > 0 for strain-stiffening membranes and ζ < 0 for strain-softening membranes.

SLS HHm model

Along the same line as previous works on ionic channel mechanosensitivity (Guharay and Sachs 1984;[START_REF] Sachs | Mechanical transduction in biological membranes[END_REF][START_REF] Assmann | Membrane tension influences the spike propagation between voltage-gated ion channel clusters of excitable membranes[END_REF][START_REF] Morris | Mechanosensitive gating of Kv channels[END_REF], 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 σ coefficient can be considered as a tension characterising the mechanical sensitivity of the gate. An increase of σ (in positive direction) speeds up the opening of the Na activation gate. Conversely, a decrease of σ (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 σ (in positive direction) slows down the closing of the Na activation gate and a decrease of σ (in negative direction) speeds up its closing. The Na activation gate closing rate is therefore favoured when reducing the membrane tension:

α m Na = α m Na,0 e σ B , β m Na = β m Na,0 e -σ B . ( 13
)
Introducing the parameter σ in the model for Na channel activation, we get the following expression for m Na ∞ :

m Na ∞ = 1 1 + β Na∞ (V ) α Na∞ (V ) * e -2σ B = 1 1 + β σ α σ . ( 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 β σ α σ = 1) but also induces a nonlinear distortion of the m Na ∞ (V ) curves for σ values away from zero.

We take B in the range of values expected from real situations. For B = 129.65 m 2 /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 σ (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 B. By playing on the tension parameter σ , one can compensate a change in σ and recover similar bifurcation diagram changes. For σ > 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 [START_REF] Banderali | Impaired stretch modulation in potentially lethal cardiac sodium channel mutants[END_REF]. 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) [START_REF] Wang | Membrane trauma and Na+ leak from Nav1.6 channels[END_REF][START_REF] Banderali | Impaired stretch modulation in potentially lethal cardiac sodium channel mutants[END_REF][START_REF] Beyder | Mechanosensitivity of Nav1.5, a voltage-sensitive sodium channel[END_REF][START_REF] Yu | Spontaneous excitation patterns computed for axons with injury-like impairments of sodium channels and Na/K pumps[END_REF][START_REF] Morris | Mechanosensitive gating of Kv channels[END_REF].

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 σ in the conductance factor m 3 Na ∞ (V )h Na ∞ (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 σ values, coming along with a sharp increase of the amplitude of the ion currents (Fig. 3d). The time constant τ m Na (V ) for the Na activation gate is also leftshifted in V for σ > 0 (Fig. 3c), as a consequence of the left shift of m Na ∞ (Fig. 3a). Less intuitive is the fact that the time constant maximum of this shifted curve is not changing much with σ , 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 σ > 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 σ (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 ∞ and I S ∞ are very different for positive and negative values of σ . For σ < 0, the sigmoidal shape of I Na ∞ disappears for σ -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 σ values, which is explained by the dramatic decrease of the conductance factor m 3 Na ∞ (V )h Na ∞ (V ) (Fig. 3b). Positive σ 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 σ 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 µA/cm 2 , σ = 0 (thin dashed-dotted line) and σ = 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 , σ ) visualising how the amplitude and period of large amplitude periodic dynamics change with the parameters σ 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 σ in the SLS HHm model (Fig. 4b). We note that a positive σ value does not change the sharp rise of m Na (t) since the two red curves are quite superimposed for t/T c ∈ [-0.1; 0] and t/T c ∈ [0.4; 1], but the length of the plateau at high m Na (t) values for t/T c ∈ [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 σ > 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- 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 and its time constant τ 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 σ > 0. The shape of the increasing part of I Na (t) for t/T c ∈ [0.05; 0.25] is modified since the exponential decay of the current amplitude (negative value) is replaced around t/T c ∼ 0.05 by a linear decay which ends on a small plateau for t/T c ∼ [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 σ > 0. Note that increasing σ from 2 mN/m to 4 mN/m enhances further the effects outlined above.

CLS HHm model

The CLS HHm model for muscle fibres ensues by introducing a mechanical factor e σ B on both α m Na , β m Na (Eq. 13) and α h Na , β h Na , to reproduce the experimentally observed coupled left shifts of both the steady-state conductance factors m Na ∞ (V ) and h Na ∞ (V ). To obtain a similar left shift of the opening and closing rates for the inactivation gate, α h Na and β h Na must be changed in an opposite manner than for the Na activation gate [Eq. ( 13)]:

α h Na = α h Na,0 e -σ B h , β h Na = β h Na,0 e σ B h , ( 15
)
where the stretching factor B h = 200 m 2 /J. Note that this choice is completely different from that proposed by Ass- mann and [START_REF] Assmann | Membrane tension influences the spike propagation between voltage-gated ion channel clusters of excitable membranes[END_REF], where all the channel rates were modulated identically.

We test the impact of such a modification on the conductance 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 [START_REF] Morris | Left-shifted Nav channels in injured bilayer: primary targets for neuroprotective Nav antagonists?[END_REF]. Interestingly, the conductance 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 inactivation 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 S (V ) curve (Fig. 5d) shows also a global shift of the range of V values where multi-steady states exist when varying σ . Contrarily to the SLS HHm model (Fig. 2d), the I app 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 app = 3.95 µA/cm 2 (Fig. 6) shows that the range of V values is now limited to negative values. The overall V amplitude is decreased when both m Na and h Na are left-shifted by the mechanical stretching energy. It is also very important to remark that the length of the plateau of maximal m Na (t) values 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 L (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 ions , V ) (Fig. 6e) and (I Na , I C ) (Fig. 6f) have to be compared to the ones obtained in the original HHm model in Fig. 1e,f, respectively.

Results

Dynamical analysis of the HHm models

Following the terminology used by Izhikevich [START_REF] Izhikevich | Neural excitability, spiking and bursting[END_REF][START_REF] Izhikevich | Dynamical systems in neuroscience[END_REF], 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 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 [START_REF] Izhikevich | Neural excitability, spiking and bursting[END_REF]. When the quiescent state is near a saddlenode 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. [START_REF] Izhikevich | Neural excitability, spiking and bursting[END_REF][START_REF] Izhikevich | Dynamical systems in neuroscience[END_REF] 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 , σ ) 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 σ are compared in the SLS and CLS HHm models.

Computation of the steady states and their stability

A dynamical system analysis [START_REF] Wiggins | Global bifurcations and chaos. Analytical methods[END_REF][START_REF] Guckenheimer | Bifurcations in Hodgkin and Huxley equations: a new twist[END_REF]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

⎛ ⎜ ⎜ ⎝ dV /dt dm Na /dt dh Na /dt dm K /dt ⎞ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎝ 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 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ V m Na h Na m K ⎞ ⎟ ⎟ ⎠ . (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 (λ) or a pair of complex conjugate eigenvalues (λ ± iω(ω > 0)) that cross the imaginary axis (λ = 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.

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 [START_REF] Dhooge | MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs[END_REF]) 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 σ . 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 [START_REF] Izhikevich | Neural excitability, spiking and bursting[END_REF][START_REF] Izhikevich | Dynamical systems in neuroscience[END_REF] 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 saddlenode 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). 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 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 H1 = 7.050 µA/cm 2 , I H2 = 8.745 µA/cm 2 , I LPC1 = 7.181 µA/cm 2 , I LPC2 = 4.335 µA/cm 2 , I T0 = 4.728 µA/cm 2 , I T1 = 7.176 µA/cm 2 , I T2 = 4.564 µA/cm 2 , I HHS3 = 6.914 µA/cm 2 , I HHS * 3 = 6.98 µA/cm 2 (colour figure online)

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 , σ, V ) in Fig. 8a, b, respectively. Increasing the parameter σ > 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 LP 1 ≤ I app ≤ I LP 2 ) widens (resp. shrinks) in the SLS HHm model for positive (resp. negative) values of σ (Fig. 8a). This interval of I app values does not increase for positive σ 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 σ 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 CLC HHm models by varying σ .

Excitability of the HHm models with tension

To investigate the excitability in the two-dimensional parameter space (I app , σ ) of the SLS and CLS mechanical variants of the HHm model, we proceed through the following numerical strategy. For each parameter pair (I app , σ ), the steady-state solutions of the four ODEs are computed first and the one with the lower V S value is chosen (lowerbranch 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 σ > 0 and that the amplitude of the single spikes also increases with σ . 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 multisteady-state window (Fig. 9a). The rightmost bottom part of the excitability domain (larger I app and smaller σ values) corresponds to single spikes with smaller voltage amplitude, 14 which ultimately resemble subthreshold oscillations. Beyond the light coloured boundary domain of these oscillations, the SLS HHm model has lost its excitability and remains in a resting state whatever the perturbation amplitude. This phenomenon is also simultaneous to the disappearance of large (issued from saddle-node LPC bifurcation) or small (issued from supercritical Hopf bifurcation) amplitude limit cycles.

The excitability diagram of the CLS HHm model (Fig. 10) bears some resemblance to the one of the SLS HHm model (Fig. 9). The excitable domain still covers a whole area on the left part of the Hopf (H 1 ) and saddle-node (LP 1 ) bifurca-tion lines on the lower steady-state fold. The fact that the two saddle-node LP 1 and LP 2 lines (green) remain very close to each other in I app (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 existence of the resting state excitability has not shrunk. On the contrary, the Hopf bifurcation H 1 on the lower branch (H 1 ) 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 voltage 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 v 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.

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 v channels gating kinetics, we analyse the range of parameters where periodic firing (large amplitude oscillations) occurs. The large amplitude limit cycles of the two SLS and CLS mechanical variants of the HHm model are computed numerically in the two-dimensional parameter space (I app , σ ). For each parameter pair (I app , σ ), the steady-state solutions of the four ODEs are computed first and the one with the largest V S value is chosen (in multi-steady-state situations). When increasing I app , the large amplitude limit cycle emerges from a saddle-node bifurcation of cycle, as previously observed in the original HHm model without tension (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 perturbation of the channel voltage, we choose the second alternative to build these diagrams. We perturb the steady state S S [V S , m Na S , h Na S , m K S ] by a negative drop (-50 mV) in V S 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).

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 app starts widening when increasing σ > 0 and shifts to lower I app 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 + and K + 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) resting 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 v channel is speeded up by tension, not only the cells become more excitable but they can also fire periodically at I app 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 σ > 0, the larger ones being observed at the left border (LPC 1 ) of the periodic firing domain. When decreasing σ 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 σ ≤ -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 σ 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 1 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 σ values, but mainly expands towards negative σ values (Fig. 12a). Moreover, the large amplitude oscillations are rather observed at negative σ 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-16 tive tensions, the periodic firing regime becomes robust and insensitive to extra perturbations. Indeed, both the right LPC 2 boundary line and the domain of existence of a stable focus on the upper branch are shifted to larger I app values, getting rid off bistability situation between the large amplitude periodic cycle and the upper branch saddle-focus steady state. This situation resembles more the myotonia muscle disorders (non-dystrophic). Myotonia is described by delayed muscle relaxation after contraction; its manifestations are stiffness, cramp or locking of muscles together with an inability to immediately relax muscles after forceful contraction or contraction induced by direct percussion (percussion myotonia). Percussion is a local and sharp mechanical stress that could be mimicked by our HHm models with tension. Actually in the CLS HHm model, the fact that the periodic firing domain spreads very far at large I app values is an indication that if the gated channels reach such a periodic firing state (which would correspond to a persistent muscle contraction), they may stay there without being able to recover to normal (uncontracted muscle cell) situation. The extension of both excitability and periodic firing domains in the CLS HHm model is indeed an indication of the inability of the system to stay in the rest state for long terms. Stochastic simulations (adding external noise with or without specific frequential signature) of the four HHm ODEs confirm that the CLS HHm model is much more sensitive to external perturbations than the SLS HHm model [START_REF] Cannon | Stochastic ion channel gating in dendritic neurons: morphology dependence and probabilistic synaptic activation of dendritic spikes[END_REF].

In Fig. 13 are compared the distributions of amplitude and period for the periodic firing regimes encountered in both SLS and CLS HHm models. The voltage amplitude distributions (Fig. 13a) are not very different in the range from 60 to 120 mV; they differ only at small and large amplitudes. In particular, the SLS HHm model predicts higher amplification of the periodic oscillations for σ > 0 than the CLS HHm model; these large amplitude oscillations occur in the LPC 1 leftmost part of the periodic firing domains (Fig. 11a). The firing period distributions (Fig. 13b) are not drastically different for the large periods; they differ noticeably at low periods (< 10 ms) that correspond to the rightmost part of the periodic firing domains in Figs. 11b and12b, respectively. The SLS HHm model can afford much faster dynamics close to the rightmost LPC 2 saddle-node of cycle bifurcation curve (Fig. 7a) at positive tension 3 mN/m σ 6 mN/m that has no counter part in the SLS HHm model.

To summarise, we show that the response of the Hodgkin-Huxley model for muscle cells to a mechanical perturbation is completely asymmetric when the Na v activation gate is changed solely (SLS HHm model). The enlargement of both the excitability domain and the window of periodic firing in I app for σ > 0 indicates that single spike and ultimately periodic firing are likely to occur more frequently when the cell is under tension. This hyperexcitability predicted by the SLS HHm model mimics HPP in damaged muscle cells and their stalling (paralysis) far from the rest state. For σ -2 mN/m, the SLS HHm model no longer fires and progressively looses its excitability. The passivity of the SLS HHm model for slightly negative σ values is an important feature that could be used to check to which extent the Na v activation gate has been impacted compared to the Na v inactivation gate. For instance, if in a hyperkalaemic periodic paralysis, lowering the tension by relaxing the muscle drives it back or close to normal situation, this could be a strong indication that the inactivation gate has not been impacted. On the contrary when both activation and inactivation Na v gates are mechanically perturbed to resemble a sick-excitable cell (CLS HHm model), the excitability and periodic firing processes persist for positive σ values and extend to negative σ values. In such a situation a tension relaxation treatment would be inefficient to restore the muscle resting state as observed in non-dystrophic myotonia disorders where the muscle remains contracted (keeps on firing) upon percussion without any ability to relax.

The dynamical system analysis proposed here opens new perspectives to physiologists as regards the possibility of distinguishing those of the Na v 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.

Fig. 2

 2 Fig. 2 Dependance on the voltage V of the steady-state conductances and currents in a HH model of muscle fibres (σ = 0) [Eqs. (5-8)]. a Probabilities of channel gate opening m Na∞ (red), h Na∞ (blue) and m K∞ (green). b Conductance factors in the steady state for the sodium channel 100 * m 3 Na∞ h Na∞ (red) and the potassium channel m K∞ (green). c Time constants for the three gates τ mNa (red), τ hNa (blue) and τ mK (green). d Ionic currents in the steady state I Na∞ (red), I K∞ (green) and I L∞ (grey) and the sum of these three currents I S∞ (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∞ and I S∞ (colour figure online)

  The α i and β i are given in reciprocal ms (T = 25 • C) and the V i in mV. Let us note that the rates of closing of the Na (β m Na ) and K (β m K ) activation gates and the rate of opening of the Na (α h Na ) inactivation gate are written as exponential decreasing functions of the channel voltage. Inversely, the rates of opening of the Na (α m Na ) and K activation (α m K ) gates and the rate of closing (β 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 .

Fig. 3

 3 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∞ (red), h Na∞ (blue) and m K∞ (green). b Conductance factors in the steady state for the sodium 100 * m 3 Na∞ h Na∞ (red) and the potassium m K∞ (green) channels. c Time constants for the three gates τ mNa (red), τ hNa (blue) and τ mK (green). d Ionic currents in the steady state I Na∞ (red), I K∞ (green), I L∞ (grey) and their sum I S∞ (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)

Fig. 5

 5 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 Na∞ (red), h Na∞ (blue) and m K∞ (green). b Conductance factors in the steady state for the sodium 100 * m 3 Na∞ h Na∞ (red) and the potassium m K∞ (green) channels. c Time constants for the three gates

Fig. 6

 6 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 3 Na (t)h Na (t) (red) and the potassium m K (t) (green) channels. In (b)

Fig. 7

 7 Fig. 7Bifurcation 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

Fig. 8

 8 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)

Fig. 9 Fig. 10

 910 Fig. 9 Excitability domain in the two parameter space (I app , σ ) of the SLS HHm model. a Excitable states are shown as coloured pixels according to the amplitude of the single spike (in mV) obtained by a numerical integration of the four ODEs. The lines for saddle-node

Fig. 11

 11 Fig. 11 Large amplitude oscillations (periodic firing) domain in the two parameter space (I app , σ ) of the SLS HHm model. a Colour coding of the voltage amplitude of the periodic oscillations (mV). The lines for saddle-node bifurcation of cycle (LPC-pink) are also represented

Fig. 12 Fig. 13

 1213 Fig. 12 Large amplitude oscillations (periodic firing) domain in the two parameter space (I app , σ ) of the CLS HHm model. a Colour coding of the voltage amplitude of the periodic oscillations (mV). The bifurcation lines for saddle-node bifurcation of cycle (LPC-pink) are also represented because they delimit the domain of existence of the

  Na (t) (red), I K (t) (green) and I L (grey). In black, the sum of these ionic currents I S (t) = I Na (t)+ I K (t)+ I L (t). 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 charging current I C (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 app = 7µA/cm 2 (colour figure online)

	(a)	(b)	(c)
	(d)	(e)	(f)
	Fig. 1 Dynamical signature of the different variables involved in a HH	I
	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 Na (t) (red), h Na (t) (blue) and m K (t) (green). c Con-
	ductance factors for the sodium m 3 Na (t)h Na (t) (red) and the potassium
	m K (t) (green) channels. In (b) and (c) the rescaled voltage curve (black
	dashed-dotted curve) is also plotted for comparison. d Ionic currents:
			,
			I Na (t) is represented versus the capacitance charging current
			I C (t) = CdV (t)/dt. HH models [Eq. (1)] are constructed

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.