Biophysics and Modeling of Mechanotransduction in Neurons: A Review

: Mechanosensing is a key feature through which organisms can receive inputs from the environment and convert them into speciﬁc functional and behavioral outputs. Mechanosensation occurs in many cells and tissues, regulating a plethora of molecular processes based on the distribution of forces and stresses both at the cell membrane and at the intracellular organelles levels, through complex interactions between cells’ microstructures, cytoskeleton, and extracellular matrix. Although several primary and secondary mechanisms have been shown to contribute to mechanosensation, a fundamental pathway in simple organisms and mammals involves the presence of specialized sensory neurons and the presence of different types of mechanosensitive ion channels on the neuronal cell membrane. In this contribution, we present a review of the main ion channels which have been proven to be signiﬁcantly involved in mechanotransduction in neurons. Further, we discuss recent studies focused on the biological mechanisms and modeling of mechanosensitive ion channels’ gating, and on mechanotransduction modeling at different scales and levels of details.


Introduction
Mechanosensation is among the most important sensory functions in animals. Indeed, animals' navigation, communication, and survival are critically dependent on their ability to perceive the external and internal forces. To reach this purpose, animals have evolved specialized sensory cells that allow the conversion of mechanical stimuli into electrical signals to be processed by the nervous system. This conversion process is usually referred to as "mechanotransduction", and constitutes the essential basis of touch, hearing, proprioception, and of all the other biological functions that rely on force detection. Mechanotransduction at the single-cell level is allowed by the combination of the morphological and molecular features of the cell and by complex interactions between cell membrane structures, cytoskeleton, and extracellular matrix (ECM). Although several molecular mechanisms have been shown to be affected by the mechanical state of cells, only a small amount contributes to direct mechanical sensing of external stimuli. In particular, the main pathway through which organisms achieve mechanosensing accounts for specialized sensory neurons opportunely enriched with mechano-electrical transduction (MeT) complexes and interfaced with different types of biological structures. A typical example is represented by mammalian hair cells in the cochlea, equipped with hair bundles composed of several rows of variable-length stereocilia. The stereocilium tip is able to convert mechanical forces into electrical signals thanks to MeT complexes formed by stretch-gated ionic channels, allowing fast and effective transduction of the auditory stimulus. This general example shares common features with a large amount of mechanical-biological sensors based on neuronal transduction. In this review, we focus on a specific aspect of the mechanosensing investigation, that is the multiscale theoretical and computational description of systems and mechanisms. We first present an overview of the main mechanosensitive (MS) channels which have been shown to have an important role in mechanosensation both in simple organisms and in mammals, including the Degenerin/Epithelial Sodium Channels/Acid Sensing Channels (DEG/ENaC/ASIC), Two Pore domain Potassium (K2P), Piezo, Anoctamin, and Transient Receptor Potential (TRP) superfamilies, also discussing the most credited biological and physiological models of ion channel gating. Secondly, concerning the theoretical modeling, we discuss different scales and methods developed to investigate mechanotransduction. All presented methods have multiscale components. Starting from the most recent advances in mechanosensitive ion-channel modeling based on all-atom and/or coarse-grained molecular dynamics, we move to the description of continuum models for pressure-induced channel gating. Finally, we present methods coupling the biomechanics to the electrophysiological modeling of neuron response, with possible applications to biological tissues of increasing size and complexity. In summary, this contribution aims to dissect the molecular basis of mechanotransduction, to present the current opinions on its biological and physiological regulation, and to discuss the physico-mathematical frameworks introduced in the literature to built descriptive and predictive models to explain mechanosensation at the molecular, micro, and macro scales. Literature was searched primarily on PubMed and Google Scholar databases (and on the journals' webpages in some cases), including different specific keywords for each topic on mechanosensation afforded in this review, such as "mechanosensitive ion channels","mechanosensitive" + channel name, "mechano-gated ion channels", "mechano-gated" + channel name, "neuron mechanosensation/mechanosensing/mechanotransduction", "atomistic model" + channel name, "molecular dynamics" + channel name, "multiscale/multiphysics model mechanosensation/mechanosensing/mechanotransduction" (where the/symbol denotes all the possible combination of the words). We restricted our research to the papers published in the last twenty years.
The paper is organized as explained below. Section 2 deals with the mechanosensitive (MS) ion channels expressed in neurons, also discussing their functional role and gating mechanisms. In Section 3, we review different modeling approaches that have been applied to describe animals' mechanosensing, from the atomistic to the tissue scale. In Sections 3.1 and 3.2, we describe the recent advances in particle-based modeling of MS ion channels coupled to the lipid membrane and the continuum models of MS channels coupled to the membrane, respectively. Section 3.3 describes multiscale modeling studies merging biomechanics to neurons' electrical dynamics at different levels of detail. Finally, Section 4 is devoted to conclusions and outlooks.

Ion Channels in Neurons Mechanosensing
Specialized mechanosensory cells and neurons can sense a wide variety of stimuli, including shear stress, vibrations, and sound waves, that induce perturbations of the cellular environment. Mechanical stimuli are translated into electrical signals by a specific category of ion channels, the mechanosensitive ion channels, endowed with the ability to open or close in response to the applied force. The heterogeneity of the mechanical stimuli and of the channels involved in mechanotransduction have manifested the need for specific criteria to establish whether an ion channel is mechanosensitive or not. Four different criteria have been proposed and applied to classify MS channels [1].

1.
Temporal and spatial expression in mechanosensory cells. The channel must be expressed in specialized mechanosensory cells, and should not be necessary for cell maturation or integrity [1].

2.
Direct involvement in the mechanical response. The channel must be critical for mechanosensitivity, in the sense that its loss abolishes the ability of the cell to respond to mechanical stimulations. The fulfillment of this criterion is a necessary but not sufficient condition to establish the mechanosensitive nature of a candidate channel. Indeed, its participation in mechanosensation could be indirect, as in the case of channels involved in cell development or in the signaling downstream of the stimulus [1]. Furthermore, a cell can compensate for the loss of a subunit by forming heteromers with other members of the same family of channels [2].

3.
Channel alterations alter the mechanical response. Alterations of the channel biophysical properties, such as ion selectivity, activation, or inactivation, result in alterations of the physiological response to mechanical stimuli. In addition, in this case, the satisfaction of the criterion does not guarantee that the considered channel is mechanosensitive; indeed, multiple auxiliary subunits could modify the biophysical properties of the channel [2].

4.
Heterologous expression induces mechanical responses in the host cell. This condition is one of the most difficult to meet. Various pore-forming subunits do not show mechanosensitive properties when heterologously expressed, as their ability to sense mechanical forces could be determined by the specific lipid composition of the bilayer [2] or by auxiliary subunits that anchor the channel to the intracellular matrix [2].
The majority of the identified MS channels belongs to the following families: TRP, DEG/ENaC/ASIC, K2P, Piezo, and Anoctamin.
Two models have been proposed to explain how the channel responds to mechanical stimuli: the "force-from-tether" and the "force-from-lipid" model. The main difference between the two models is that in the first, the channel has to be bound to intra/extracellular structures such as the cytoskeleton and/or the extracellular matrix [3] ( Figures 1B and 2B), while in the second, the force is transmitted to the channels directly through the lipid bilayer ( Figure 3B,D).  [4,5]). The conversion of mechanical stimuli into electrical signal happens at the tip of the stereocilia where MeT complexes are localized. MeT complexes have specular stoichiometry and are composed of two Transmembrane Channels (TMCs) that interact with the tip link (PCDH15-CD2, in pink) through the TMIE and the LHFPL5 proteins [5][6][7][8][9][10][11]. In addition, the channel directly interacts with the Calcium Binding Protein-2 (CIB-2) [9,12]. A recent study on the Caenorhabditis elegans (C. elegans) model organism showed that the connection of the complex with the cytoskeleton is realized through the binding of the CIB protein with ankyrin protein UNC-44 [9]. However, in mammals, this aspect is still poorly explored, and the target of CIB-2 proteins is unknown.
In the following, we briefly describe the MS ion channels that have been identified in the literature, with a particular focus on their functions in mechanosensory neurons. The complete list of the mechanosensitive channels here reviewed is reported in Tables 1 and 2.

DEG/ENaC/ASIC
The DEG/ENaC/ASIC is a superfamily of voltage-insensitive, Na + permeable, and amiloride-sensitive channels. These channels are found in both vertebrates and invertebrates and include the nematode degenerin (DEG), the Drosophila Pickpocket (ppk), Ripped Pocket (rpk) and Balboa, and the vertebrate Epithelial Sodium Channels (ENaC) and Acid-Sensitive Channels (ASIC). Channels of the DEG/ENaC family share a similar structure with two transmembrane domains and a large extracellular domain [13][14][15][16][17] (Figure 2A) and assemble in homo-or hetero-trimers with a chalice-like shape [13,15,18,19]. DEG/ENaC channels are involved in vertebrate and invertebrate mechanosensation.
Drosophila expresses three DEG/ENaC channels encoded by the pickpocket, ripped pocket, and balboa genes [13,26]. These channels are expressed in nociceptive class IV multidendritic (md) neurons of the fly peripheral nervous system. PPK and Balboa can form heterotrimeric channels that play a critical role in mechanonociception [13]. In addition, PPK mutant flies display impaired locomotion [27].
There are several pieces of evidence that DEG/ENaC/ASIC channels give rise to mechanically activated currents. However, the gating mechanisms are still not completely understood. The nematode DEG channels are proposed to rely on a "force-from-tether" model to transduce mechanical stimuli [42]. In this scheme, the pore-forming subunits (MEC-4 and MEC-10) may interact with the extracellular proteins MEC-1, MEC-5, and MEC-9 [43] and with the special tubulins MEC-7 and MEC-12 [42,44]. However, evidence of in vivo co-localization of MEC-4 with MEC-5 proteins and of MEC-7 and MEC-12 with the complex is still missing [42,45]. Moreover, MEC-4 and MEC-10 subunits associate with the stomatin-related protein MEC-2 and the paraoxonase-like protein MEC-6, which are critical regulators of the channel activity [20,[46][47][48]. Concerning the mammalian ASIC channel, both the "force-from-lipid" and the "force-from-tether" models have been proposed to explain the mechanical gating [38]. Similarly to their nematode homologs, ASIC channels associate with stomatin-domain proteins, STOLM1 and STOLM3 [49][50][51][52][53][54], suggesting that they could rely on a force-from-tether mechanism. However, this aspect still deserves further investigation, in particular, to highlight eventual interactions with extracellular and intracellular matrix proteins. The wrist connects the transmembrane domain to the seven stranded β-sheet palm and the α-helical thumb, while the β-ball is surrounded by the fingers and knuckle [14][15][16] (adapted from [15,55]). (B) Hypothetical gating mechanism of C. elegans DEG channels. The molecular bases of mammalian ASIC channels mechanical gating are still poorly explored. In the case of C. elegans, it has been established that the pore-forming subunits MEC-4 and MEC-10 interact with the paraoxonase-like and the stomatin-like proteins MEC-2 and MEC-6 [46][47][48] and with the special tubulins MEC-7 and MEC-12 [44]. In addition, the extracellular matrix components, MEC-1, MEC-5, and MEC-9, might be critical for force transduction [42,43] by directly associating with the pore-forming subunits (adapted from [38,55]).

K2P
The two-pore domain K2P channels are a family of tetrameric channels with two-pore forming loops [56] (Figure 3A,B). Some members of the family, such as TREK-1, TREK-2, and TRAAK, form thermosensitive and mechanosensitive channels expressed in different organs and tissues, including the central and peripheral nervous system [57][58][59][60]. In particular, TREK and TRAAK channels are expressed in dorsal root ganglia neurons where they regulate the threshold for mechanical responses [56,61]. Mice lacking K2P channels do not show significant sensory function impairment; rather they show hypersensitivity to mechanical and thermal stimuli (mechanical and thermal allodynia) [58,59]. When heterologously expressed, TREK and TRAAK show stretch-activated currents [57,60,62,63]. These currents are equally induced by the application of positive and negative pressure and are not altered by the disruption of intracellular components (e.g., the cytoskeleton) [60,62,63]. Moreover, these currents do not require the integrity of intracellular components, suggesting that K2P channels mechanical gating relies on force-from-lipid mechanism [60,62,63] ( Figure 3B). Drosophila orthologue of TREK channels, ORK1, is critical for the regulation of sleep and cardiac rhythm, similarly to its mammalian counterpart [64][65][66]. However, in contrast to mammalian K2P, there is still no report of a direct activation by mechanical stimuli. and two pore-forming regions (in black). They assemble in tetramers to form potassium channels such as the mechanosensitive TREEK-1/2 and TRAAK (adapted from [3]). (B) Model of mechanosensitive K2P gating. Mechanosensitive K2P channels rely on the force-from-lipid model to sense mechanical stimuli. The channel opens in response to membrane stretch thanks to its direct interaction with membrane lipids (fuchsia) (adapted from [67]). (C) Cartoon representation of mouse Piezo2 channel monomers. Piezo2 channels are predicted to have up to 38 transmembrane domains [3]. In addition to the transmembrane domains, Piezo channels possess also a C-Terminal-Extracellular domain (CED) that constitutes the central dome, as well as a beam and an anchor that connect the central pore of the channel with the mechanotransduction module (adapted from [68]). (D) Model of Piezo channel gating. Piezo channels assemble in homotrimeric structures that sense the stretch of the lipid bilayer (adapted from [69]).

PIEZO
Piezo channels are large homotrimeric channels with a peculiar structure that resembles a three-bladed propeller [70] ( Figure 3C). Vertebrates have two Piezo channels, Piezo-1 and Piezo-2, that, despite the similar structure, share only 42% of sequence homology [68]. Piezo-1 is mainly found in non-neuronal cells [55]. In contrast, Piezo-2 channels are primarily expressed in mechanosensory cells (Merkel cells, hair follicles, and hair cells of the auditory system) and neurons (dorsal root ganglion neurons) [71][72][73]. The two channels have only one homologue in Drosophila and C. elegans [74]. The Drosophila Dmpiezo gene shows an expression pattern similar to that of mice (it is found in non-neuronal and neuronal tissues), and the corresponding channel plays a critical role in noxious mechanosensation [75]. In C. elegans, Piezo channels are encoded by the gene pezo-1, which is widely expressed throughout the reproductive tissues [76].
Piezo channels generate mechanically-activated currents in response to different mechanical stimuli, including pipet pressure [74,77], Myosin-II mediated internal forces [78], optical tweezers [79], nanoparticles [80], and laminar shear stress [81,82]. There are several recent works devoted to the study of the gating mechanisms of Piezo channels. All these works report a direct activation of Piezo channels with membrane indentation, explained by the force-from-lipid model [70,71,74,75,[83][84][85] (Figure 3D). Moreover, an important contribution to Piezo gating comes from the actin cytoskeleton [86,87]. Works by Cox et al. [86] and Retailleau et al. [87] suggest that it plays a mechanoprotective role, regulating the activation of the channels and the membrane tension [86][87][88]. In addition, Zheng et al. [85] have recently shown that Piezo-1 and Piezo-2 differ in the activation mechanisms. Indeed, while Piezo-1 channels are extremely sensitive to cold-induced membrane stiffness alterations, Piezo-2 channels seem to be unaffected by such changes, suggesting that the force-from-lipid is not the principal activation pathway for Piezo-2 [85].
TMC channels are expressed both in vertebrates and invertebrates, and they play a critical role in mechanotransduction [55,90]. The TMC family is composed of two channels, TMC-1 and TMC-2, that are essential for hearing in different organisms. In mammals, TMC-1 and TMC-2 are expressed in the stereocilia of the inner ear hair cells, where they are part of an MeT complex ( Figure 1) [91,92]. TMC channels are important also for Drosophila proprioception and locomotion. In Drosophila larvae, DmTMC-1 and DmTMC-2 are expressed in the sensory neurons that provide sensory feedback for locomotion, such as the class I and class II dendritic arborization neurons and bipolar dendrites neurons [93,94]. TMC channels are critical also for Danio rerio water motion detection and hearing, with a prominent role of TMC-2 instead of TMC-1 [8,95]. C. elegans possesses two TMC genes: tmc-1 and tmc-2 that are expressed in mechanosensory neurons [9,96]. Recent studies, based on cryo-EM and size-exclusion chromatography, confirmed that both TMCs are pore-forming subunits that assemble in dimers [10,90,97]. However, TMC channels alone are not sufficient to transduce mechanical stimuli into electrical currents and couple with other proteins forming MeT complexes [6,[8][9][10] (Figure 1).

TRP Superfamily
Transient Receptor Potential channels are a superfamily of channels firstly discovered in Drosophila. TRP channels are found in yeast, fungi, and animals, where they are involved in processing a wide variety of sensory stimuli, such as temperature, sound, light, chemicals, osmotic changes, and mechanical forces [98].
Based on structural homology, TRP channels are divided into nine subfamilies: TRPC (canonical), TRPV (vanilloid), TRPM (melastatin), TRPP (polycystin), TRPML (mucolipin), TRPA (ankyrin), TRPN (no mechanoreceptor potential C, NOMPC), TRPS (soromelastin), and TRPY (yeast) [99,100]. Except for TRPN and TRPS, which are exclusively expressed in invertebrates, at least one member per subfamily is expressed in both vertebrates and invertebrates [99,101]. Most TRP channels are non-selective cation channels, structurally similar to voltage-gated channels (VGCs) [100,101]. The channels have 6 transmembrane domains (TM), a pore loop domain between TM5 and TM6, a cytosolic (N), and a carboxy (C) termini. In addition, some classes of TRP channels (i.e., TRPN, TRPC, TRPA) possess ankyrin repeats that may be important to define their mechanosensitive properties [102] ( Figure 4). TRP channels assemble in homo-or hetero-tetramers [3,101] in which the subunits belong either to the same TRP subfamily or to a different subfamily, giving rise to a wide variety of channels [101].  [99]). (B) TRPV channels. Representative diagram of human TRPV1 channels. In addition to the six transmembrane domains, TRPV channels have six ankyrin repeats (adapted from [99]). (C) TRPC channels. Representative diagram of TRPC channels (e.g., hTRPC5). In addition to the six transmembrane domain, TRPC5 channels have four ankyrin repeats, a PRE-S1 domain, as well as a coiled-coil domain (adapted from [99,103]). (D) TRPM channels. Representative structure of the TRPM channels (e.g., hTRPM8). The channels have six transmembrane domain, a SLOG domain, a TRP domain, a coiled-coil domain, and a Variable C-terminal domain. The Variable C-terminal domain could encode an alpha-kinase domain or a nudix-hydrolase domain. (adapted from [99,104]). (E) TRPA channels. Structure of human TRPA-1 channels (e.g., hTRPA1, adapted from [99,105]). TRPA channels have a high number of ankyrin. The number of ankyrin domains is variable between species and between TRPA1 channels of the same species. (F) TRPM-L channels. (adapted from [99]). (G) Structure of the vacuolar TRPY channels (adapted from [99,106]). (H) TRPS channels. The structure of TRPS channels needs further investigation. Their structure is similar to that of TRPML channels, and they are expressed among mollusks, nematodes, tardigrades, myriapods, and chelicerates [99]. Their structure includes a SLOG domain, a TRP domain, and a nudix-Hydrolase domain. (adapted from [99,107]). (I) TRPP channels. Structure of the TRPP channels (e.g., hPDK-1 and hPDK-2). PDK-1 and PDK-2 are expressed together and form a protein complex (adapted from [108]). Channels from almost every family have been shown to be involved in mechanosensation [1,101,109]. The only two families whose members are not involved in mechanosensation are the TRPML and the TRPS [99,110]. The complete list of mechanosensitive TRP channels expressed in neuronal cells is reported in Table 2. In the following, we discuss the current knowledge of the role of TRP channels in mechanosensation and their gating mechanisms, with some specific examples.
The mechanosensory ability of TRP channels and their gating mechanism by force are still debated. Based on how they react to mechanical stimuli, mechanosensitive TRP channels can be divided into two main categories: the mechanically-gated (MG) channels and the mechanically sensitive channels [129]. The first class includes the force sensors themselves that are directly activated by the stimulus [130,131], while the others are activated by second messengers whose release is triggered by the activation of the primary force sensor [128][129][130][131][132][133][134].
A particular example of indirect mechanical activation of TRP channels is given by the photomechanically activated TRPs in Drosophila. In this scenario, the phospholipase C (PLC) pathway activated by the G-protein signaling results in the conversion of phosphatidyli-nositol4, 5-biphosphate (PIP 2 ) into diacylglycerol (DAG). DAG occupies significantly lower space in the membrane than PIP 2 ; thus, PIP 2 hydrolysis might induce modifications of the membrane curvature in the proximity of the TRP channel leading to its mechanical activation [128]. In other cases, the indirect activation of TRP channels might be realized through the binding of specific ligands, such as DAG, PIP 2 [131], reactive oxygen species (ROS) [134], or serotonin [155].

Anoctamin superfamily
Mammalian TMC-1/2 Human and mice hearing Stereocilia of inner ear hair cells Pore forming subunits in MeT complexes in inner ear hair cells.

Modeling
We here deal with the modeling of mechanosensing at multiple scales and levels of detail. In particular, we focus on three main approaches that have been pursued in the literature to explain mechanosensing functioning and mechano-electrical coupling in neurons: • Atomistic modeling. Suitable for simulating single channels embedded in the membrane lipid bilayer and short timescales processes (tens of microseconds). These models have the potential to describe at the finest detail the molecular mechanisms at the basis of mechanosensing and to analyze channel kinetics in controlled environments from a mechanical and chemical point of view. The principal drawback is the high computational cost inherent in many-body dynamics. • Mechanosensing continuum modeling. Appropriate to describe single channel as well as channel pools kinetics. Two approaches can be pursued in this context. The first is the continuum-version of channel-membrane systems, designed to overcome computational limitations of atomistic modeling and span different length scales. The second is the detailed modeling of channels kinetics describing states and transitions and their dependency on mechanical strains and stresses, disregarding the spatiality of the system. • Multiscale mechano-electrical modeling. Fundamental to simulate the macro and microscale mechanical problem in soft biological tissues and couple the mechanical state to the electrical neuronal response. These models can take advantage of the information derived from atomistic and continuum modeling of channels' dynamic, defining suitable phenomenological and biophysical coupling laws.

Particle Dynamics Modeling of Ion Channel Mechanosensing-Lagrangian Description
A detailed atomistic molecular description, performed via a proper statistical sampling, is in principle able to capture the molecular functioning and the coupling mechanism at the channel-membrane interface. Limitations of this kind of modeling are related to size limits on the simulations, both in terms of space and time. The atomistic molecular modeling is able to describe at most isolated channels inserted in the lipid membrane and a timescale of a few tens of microseconds in the largest simulations. This level of description is able, on the other side, to disentangle the very molecular mechanism of mechanotransduction, including effects of mutations and drugs functioning, and it can provide parameters for higher-level or multiscale methods, as described in the following sections.
In this section, the very recent outcomes of molecular dynamics (MD) simulations in the investigation of lipid membrane regulation of membrane proteins are described [192], focusing on eukaryotes channels. We refer the reader to [193] for more details on molecular dynamics simulation results for the mechanosensing proteins MscS and MscL in prokaryotes.
In general, lipids of the cellular membrane may interact with membrane proteins as a whole, i.e., the membrane exerting lateral pressure on the protein transmembrane domain, or instead acting as ligands, binding in protein pockets and exercising allosteric effects. Some atomistic simulations have been performed specifically concerning this last point, which can be viewed as related to mechanosensing, in the sense that upon concerted motion of channel protein and membrane lipids, isolated lipid chains can occupy specific binding pockets and, therefore, favor or disfavor channel opening. In particular, Kir2 and PC2 have been investigated [194][195][196], and specific lipid-binding sites have been identified for these channels (see Figure 5). However, focusing on the modeling of membrane tension and its effects on channel proteins, molecular dynamics simulations of mechanosensing channels are, overall, a quite recent success for two main reasons. First of all, the X-rays or cryo-EM structure of eukaryotes mechanosensing channels only recent became available [63,67,[197][198][199]. Furthermore, the simulation of mechanosensing mechanism requires more than a standard MD simulation to model the effect of applied membrane tension, with simulation sizes that only recently become doable. For this reason, while several standard molecular dynamics simulations are available, characterizing various ion channels which are, among other properties, also mechanosensitive, very few studies are focused on the simulation of the effects of membrane tension on the gating mechanism or on the role of isolated lipids in controlling the gating.
In particular, the TREK-2 channel, from the K2P family, has been recently investigated in depth [200][201][202][203] via a combination of experimental and modeling tools, and in particular via extended molecular dynamics simulations. The used MD protocol has been borrowed from previous studies on prokaryotic mechanosensitive channels to simulate the lateral tension [204]. The increasing symmetric membrane tension is simulated by increasing the lateral dimensions of the lipid bilayer. For various lateral pressure values, down to -50 bar (that is, at increasing tensions), all-atom MD simulations were performed on the TREK-2 channel in down (closed) conformation, embedded in a 1,2-palmitoyl-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer. In addition, to simulate asymmetric stresses, applied only to a single leaflet of the membrane, a slightly different protocol has been used [202], keeping fixed the lateral size of the membrane and selectively removing lipids from either inner or outer leaflets. This protocol simulates the effects of membrane curvature or of the insertion of noncylindrical lipids in the membrane. Many interesting outcomes result from this investigation [200][201][202]. The TREK-2 channel, subjected to increasing membrane tension, in particular to a symmetric stretch of the bilayer, undergoes a large conformational change, from the so-called 'down' conformation to the 'up' conformation, with an overall expansion of the protein section. It is worth noting that the change in the cross-sectional area of the channel is limited to the inner leaflet, while the selectivity filter region is almost unchanged upon tension. Notably, the selectivity filter is the putative location of gating in TREK-2 and TRAAK channels. Indeed, when the asymmetric stretch is applied, the increasing tension within the inner leaflet is able to induce the down-up conformational transition, similarly to the symmetric protocol, while tension applied to the outer face does not induce conformational rearrangements [202]. The relatively low values of pressure inducing the down-up conformational change are consistent with the fact that K2P channels are much more sensitive than MscL to pressure variations. Another interesting result provided by all-atom simulations is the lateral, non-homogeneous profile of the pressure, which can be characterized as a function of the depth into the membrane and of the applied tension. A comparison with a lipid bilayer producing a different pressure profile demonstrated its key role, with respect to other putative ingredients as bilayer thickness and hydrophobic mismatch, at least for the considered channel. The same MD protocol used for TREK-2, applied to the homologous but non-mechanosensitive K2P channel TWIK-1, could not obtain conformational changes, assessing the specificity of the TREK-2 channel and the quality of the atomistic modeling [200] (Figure 6). The connection between conformational states (up or down) due to applied tension and the conductive state of the channel deserves more attention [203] because the transmembrane domain is intrinsically mechanosensitive, but it is also allosterically coupled to the gating mechanisms. A striking difference between MscL and TREK-2, concerning the channel pore behavior upon membrane tension, is that the protein expansion drives the direct pore opening in MscL, while the mechanism related to the channel conductivity is more complex in TREK-2 and it retains the K + selectivity [201]. Indeed, the TREK-2 selectivity filter is partially affected by membrane stretch [200] because ion occupancy changes depending on applied tension, while the ionic selectivity, interestingly, is not altered by membrane stretch [201]. Overall, the MD results for TREK-2 confirm that the model force-fromlipid can successfully explain the TREK-2 behavior, that the direct lipid occlusion has a non-prevalent role in the gating mechanism, that the ionic selectivity is retained at the selectivity filter, and that the asymmetric sensitivity of TREK-2 helps in maintaining the filter structure [200][201][202]. Further investigations have been devoted to the relationship between conformational and conductive states of the channel by simulating TREK-2 up and down states with physiological membrane potential, together with applied membrane tension [203]. Apart from gaining even more details of the gating mechanisms, those last MD simulations could identify the down conformation as a state most probably non-conductive, the up conformation as mostly conductive, and verify that membrane tension increases conductivity by favoring the up conformation. Figure 6. Membrane tension applied to the TREK-2 channel via molecular dynamics simulations. The transmembrane domain modifications induced by lateral membrane tension are shown. When the membrane is not strained, the channel is mostly in the down conformation. When tension is applied to the membrane, the channel moves towards the up conformation. The selectivity filter region, outside of the membrane, is not structurally altered by the applied tension. Adapted from [200] under the CC BY license.
The TRAAK channel, pertaining to the same K2P family, has also been modeled and studied [205] via molecular dynamics and Brownian dynamics [206]. The steric occlusion of the channel produced by lipids has been in this channel related to the loss of conduction in the closed conformation. Therefore, the direct lipid occlusion can also play a role in the complex regulation of channel proteins by the lipid membrane.
As said, the membrane-channel protein system is a typical case of multiscale systems, for which proper multiscale approaches should be used. The multiscale description may be serial, in the sense that results from a certain level of description can be used in a sequential protocol as input data for a different level of description. An example regarding the discrete description is given by the combination of coarse-grained molecular dynamics (CG-MD) and all-atom molecular dynamics. CG-MD simulations may efficiently explore membrane protein/lipid interactions, producing starting configurations to be converted to an atomistic resolution to refine and characterize the coupling [207]. The mechanosensitive ASIC channel has been used as a test case to verify the quality of the protein-lipid interface description at the coarse and atomistic level [207].
The combined CG-MD and all-atom MD approach has also been used in a recent case study, similar to TREK-2, provided by Piezo1 [208]. Its activation mechanism, related to the transition from closed to open conformation, is due to structural rearrangements related to cell membrane tension. Piezo1 has been investigated within a two-step, multiscale protocol. An initial coarse-grained molecular dynamics (CG-MD), using the Martini force field, has been performed to obtain the initial conformation, with a complex asymmetric membrane conformation. Piezo1 has an initial stable conformation characterized by a curved shape, with indentated membrane at rest. Once stabilized, the system has been converted to an all-atom description. The membrane stretching protocol is similar to the one used for TREK-2 [200]. Upon membrane stretching, both membrane and Piezo1 flatten, and Piezo1 expands, remaining embedded in the lipid bilayer. The extracellular domain becomes more exposed upon Piezo1 flattening, and the water-filled region inside the pore increases. Therefore, the applied membrane tension drives the channel opening. Due to the structural similarities between Piezo1 and Piezo2, it is hypothesized that the same mechanism holds for Piezo2 [208] (Figure 7). The molecular dynamics is, therefore, a powerful tool to describe and understand at the molecular level the functioning of mechanosensing channels, and its desirable application to other families of MS proteins would provide much more information that could also be included in multiscale continuum descriptions.

Continuum Modeling of Ion Channel Mechanosensing-Eulerian Description
The multiscale description necessary to study mechanotransduction can be afforded via continuum models that can mimic various mechanisms contributing to channel activation. The continuum models have been so far only applied, to the author's best knowledge, to scL channels, systems in which the pressure-driven gating mechanism is directly related to the channel opening, therefore less complex than in the case of TREK-2 or Piezo1 described above in terms of all-atom models.
A possible continuum mechanics-based simulation method, called Molecular Dynamicsdecorated Finite Element Method (MDeFEM) [209] is based on the description of the membrane as an elastic/viscoelastic medium, whose mechanical properties are obtained from molecular dynamics, as well as the coupling parameters with the embedded channel protein. The channel proteins are modeled as combinations of elastic continuum structures, whose motion (collective behavior) is described by the lowest modes obtained from the principal component analysis. The method can, in principle, couple different length scales and manage corresponding loadings, and it has been successfully applied to MscL [209]. A similar method, with some differences in the details of the system description, has been used more recently to study MscL [210], with a special focus on the combined effect of membrane bending and hydrophobic mismatch. The representative volume element (RVE) model of MscL is able to describe the channel pore corresponding to various functional states (in detergent, resting, and open [210]). The hydrophobic pairwise interaction between protein and bilayer is also included. The study identified the most favorable condition for MscL gating in the local bending together with zero hydrophobic mismatch and found that inward bending (negative bending) is more effective for MscL opening compared to outward (positive) bending.
An important modeling step is the proper evaluation of microscopic stresses [211,212], that connect the membrane continuum [211] to the protein described, via statistical mechanics, at the atomistic level. In some mathematical definitions of microscopic stress, some ambiguity may emerge, in particular when applied to complex materials, including protein-lipid assemblies. A detailed mathematical definition of microscopic stress tension has been recently developed, based on an accurate decomposition of forces even for many-body potentials [211,212]. The most detailed description of continuum microscopic stresses obtained from discrete systems mechanics is based on the Irving-Kirkwood (IK) stress definition for two-body interactions, extended to include the central force decomposition (CFD) for treating the many-body potentials and satisfying the mechanical equilibrium laws. The method, applied to a lipid bilayer, successfully incorporates the details of membrane layered and torquing stress. Therefore, it is important to use a reliable definition to evaluate microscopic stress to derive microscopic properties from MD simulations to be implemented in continuum models.
A generalized model for the activation of MeT channels in C. elegans [213] is based on the main idea of coupling, in a suited manner, MeT channels to the viscoelastic medium subject to external stimuli, and including in the description the dynamics of MeT gating, modeled in a mean-field theory framework. The coupling to the viscoelastic medium is modeled with each MeT channel linked to a mechanical anchor via an elastic filament. Both the filament and the anchor are embedded in the viscoelastic extracellular matrix. External strain induces a displacement between the channels and their anchors. Two other forces act on the filament/anchor systems, the elastic force due to the elongation of the elastic filament and the friction forces between the extracellular matrix and the filament/anchor system. Both forces are suitably represented with springs. Figure 8 represents a sketch of the proposed MeT channels gating mechanism.  [213]. In clockwise order from top to bottom: a deformation applied to the extracellular matrix causes a lateral shift between the anchor and the matrix and induces a stretch of the filament that activates the channel. Viscoelastic forces allow the filament to return into its relaxed state, inducing the channel adaption and closure. When the external deformation is released, a second lateral shift of the anchor with respect to the extracellular matrix occurs, resulting in the channel opening. Finally, the initial configuration is restored by viscoelastic forces. Figure adapted from [213] under the exclusive License to Publish.
Concerning the gating dynamics, the opening/closing of a single MeT channel is driven by two factors: the filament elongation, inducing the channel opening, and the statistical probability of the channel being in a conductive state. The gating dynamic is described by three states (open, closed, and subconducting) and by their transition rates. A realistic model should describe the different effects of membrane stretch at different positions in the viscoelastic medium. To overcome the complexity of a point-like description to describe the response of many MeT channels in different spatial locations, the membrane tension and its effect on elastic filaments are described via a mean-field method. Overall, the population of MeT channels is collapsed in a unique, effective MeT channel, whose gating dynamics overall takes into account the statistical response of a region of space with embedded MeT channels. The obtained effective channel, based on available experimental data, describes the dynamics of 15-30 actual channels, in good agreement with activated channels as estimated from experimental data. The inclusion of higher-order mechanical effects in the model shall require more detailed experimental data. Notably, biophysical modeling of mechano-activated channels can take inspiration from other detailed models built for different transduction mechanisms, e.g., [214,215].

Multiscale Mechano-Electrical Modeling
In this section, we discuss modeling studies dealing with the multiscale coupling of cells or tissue mechanics to the electrophysiological response of neurons. Those applications range from descriptions of mechano-electrical coupling at the single neuron level in simplified loading conditions to more realistic modeling investigations, including complex mechanical problems at the tissue level.
In the context of simplified loading mechanical conditions, Jérusalem et al. [216] adopted a multiscale mathematical model to investigate the effect of spinal cord injury and mechanical damage on myelinated axons' electrical response. In particular, they formulated a mechanical model of the myelinated axon-extracellular matrix system, comprising a Maxwell material for the axon in parallel to a viscoelastic material for the surroundings. They also included irreversible plastic damage in the model and assumed a small strain regime and quasi-static conditions. The microscopic axial strain is opportunely split into two contributions from Nodes of Ranvier (NRs) and Internodal Regions (IRs). Axon incompressibility is further used to couple the axial strain to cell membrane strain and diameter deformations of both surface membrane and myelin layers. Finally, axon electrical response is modeled via Hodgkin-Huxley and Cable theories, where electrical parameters of myelin layers and cell membranes depend on the corresponding mechanical state by means of the strain. Among the electrical parameters affected by the strain, we can find membrane capacitance, axial and membrane resistance, ion conductances, Nernst potentials, and steady-state activation/inactivation functions of ion channels' dynamics. The model is calibrated on tensile tests on guinea pig spinal cord nerves, performed at variable axial strain rates and different maximum applied strains. Model results show action potentials (AP) impairment after the removal of damaging loads consistent with experiments.
Another interesting contribution in the same direction comes from Tian et al. [217], who analyze alterations in action potential (AP) propagation in neurons subjected to stretching and plastic deformations. Compared to Jérusalem et al., their model is based on a slightly modified mechanical description, relying on power-law responses of neuronal deformations and a similar electrical compartment and mechano-electrical coupling. Simulations results are consistent with experimental observations and show AP suppression, increased spiking frequency, and altered conduction velocity at increasing axial strains.
Additional efforts in coupling mechanics to electrophysiology in nerves focused instead on the analysis of pressure waves traveling along the axon as induced by action potential propagation (Engelbrecht et al. [218]). In this study, it is considered both the electromechanical and the mechano-electrical coupling. In particular, pressure, longitudinal and transversal mechanical waves along the cell membrane, modeled as a two-dimensional elastic cylinder, are assumed to be triggered by forces induced by local excitations of the nerve, and the longitudinal displacement due to mechanical waves affects back electrical activation through proper forcing terms. Longitudinal and transversal waves are described based on the Heimburg and Jackson model for biomembranes and theory of rods, respectively. Electrical activation is reproduced with a FitzHugh-Nagumo model, generalized to include mechanical activation.
An interesting application of multiscale modeling, coupling mechanics and neuronal dynamics, concerns the mechanotransduction in skin. In this direction, Sripati et al. [219] formulated a model to reproduce the stimulus-response coupling in cutaneous slowly adapting 1 (SA1) and rapidly adapting afferent fibers, innervating Merkel and Meissner corpuscles, respectively. The model includes a continuum mechanical description of skin, allowing for predictions of strain profiles in response to applied loads. The skin is modeled as a linear elastic, homogeneous, and isotropic medium, and its mechanical state is described based on the Timoshenko and Goodier stress model for a point load, on a decomposition of the stress on a basis of such unit point-load responses, and on the Hooke's constitutive law. The skin's mechanical state is phenomenologically coupled to the afferent fibers' discharge rate by fitting models' parameters to reproduce experimental data recorded on monkeys. Interestingly, a combination of mechanical variables satisfying specific rotational invariance properties is included in the coupling term, among which the maximum compressive and tensile strains (and stresses), the vertical strain (and stress), the maximum horizontal strain (and stress), the maximum deformative strain (and stress), the sum and product of the strain eigenvalues, and the strain energy density. In regards to mechanotransduction in SA1 fibers, a further contribution by Lesniak et al. [220] introduced an interesting model coupling the skin mechanics to an integrate-and-fire electrical compartment to reproduce neuronal discharges. In this study, the authors built an FEM linear elastic model of fingertip skin reproducing skin's microstructure and multilayer composition. A strain energy density for nearly incompressible material, defined by means of the octahedral shear stress, is used to couple the mechanical state of the fingertip to the electrical response of the afferent fibers, whose neuronal dynamics are reproduced with a classic leaky integrate-and-fire (LIF) description. The coupling term is represented by a specific receptor current whose activation is shaped via sigmoidal function, dependent on the strain energy density. In this case, also, model parameters are fitted to reproduce experimental data recorded on monkeys. A similar and more novel contribution dealing with multiscale modeling of Merkel cell-SA1 fiber complex comes from Gerling et al. [221]. In this study, the authors formulated a model based on three components: a skin FEM model, a current-generator function, and a leaky integrate-and-fire electrical model.
The skin model is based on a multilayer description and on a quasi-linear viscoelastic material, with parameters fitted on mouse skin data, and permits to compute compressive stresses within the skin in response to ramp-and-hold displacement stimuli. The generator function basically couples the stress' dynamic to a stimulating current composed of fastinactivating, slowly-inactivating, and ultra-slowly-inactivating components. The coupling is performed through a convolution of the time derivative of the stress profile with different exponential functions describing the output for each current component in response to a step-stress stimulation. The total induced current is taken as input in an LIF model of neuronal response. Figure 9 shows the general modeling concept.
Detailed multiscale models were also developed for mechanotransduction in Pacinian corpuscles (PCs). In a specific application, Quindlen et al. [222] developed a comprehensive description of the PC based on finite-element models of the PC's outer and inner regions, and a detailed model of a neurite included in the PC. In this study, the mechanical description of the PC's outer shell accounts for both lamellar spherical shells and an interlamellar fluid component. In addition, the inner core is modeled as a linear elastic isotropic sphere, including a cylindrical neurite with a spherical ending, enriched with smaller cylinders to reproduce the neurite's filopodia. The outer shell model is stimulated by periodic pressure profiles applied at the outer surface at different frequencies, predicting the displacement of each spherical shell. The displacement of the inner shell is then used as a boundary condition for the inner core mechanical model, which permits the computation of the mechanical state of the neurite. Finally, the mechanical state of the neurite is used to modulate a stretch-gated ion current in a biophysical Hodgkin-Huxley (HH) neurite implemented in the NEURON environment. This sub-model includes classic HH potassium and sodium currents and point current sources to simulate stretch-activated channels, whose maximum current value is modulated by the neurite's strain through a sigmoidal activation function.
Another interesting multiscale model of mechanotransduction concerns the touch sensing in C. elegans. In this particular scenario, Sanzeni et al. [223] presented a detailed model coupling the macroscopic nonlinear mechanics of the nematode body to the activation of MeT channels in touch receptor neurons. The authors develop a model accounting for the mechanical deformation of the worm in response to the application of glass probes or microbeads and a generalized model of the activation of MeT channels, based on the work by Eastwood et al. [213]. In particular, the worm is modeled as a thin elastic cylindrical shell in a regime of large deformations (Figure 10). At the microscale level, MeT channels' activation is coupled to the displacement field computed by the mechanical model by assuming that each channel is connected to an elastic filament giving rise to elastic and friction forces. The computation of the filaments extension and the resulting elastic forces along the membrane tangential component is further used to modulate the change in free energy between MeT ion channels' states and compute the activation of ion channels' populations and the induced mechanosensitive current.

Discussion and Conclusions
In this contribution, we reviewed the main ion channels expressed in neurons directly involved in mechanosensation and in transducing mechanical signals from the environment in specific electrical outputs and provided a comprehensive description of available theoretical frameworks successfully applied to the investigation of mechanosensing at different scales.
The identified MS ion channels superfamilies include DEG/ENaC/ASIC, K2P, PIEZO, Anoctamin, and TRP channels. These channels are expressed both in specialized cells and in sensory neurons and share similar properties in the underlying mechanisms of channels' gating. Those specific properties permit mechanical sensing by means of stretch and stress thanks to interactions between ion channels, cell membrane, cytoskeleton, and ECM. The gating of all those channels is regulated by complex mechanisms, including the two main contributions given by lipid-and tether-transmitted forces, and it is reasonable to assume that those features are opportunely combined to achieve optimal mechanical sensitivity.
In the last years, several efforts have been devoted to the investigation of MS channels' gating and to the development of accurate models able to describe channels' activation in response to mechanical stimulation. Those models were developed on different levels of detail, including the molecular and the cell/tissue scales. Our analysis on atomistic and molecular dynamics approaches to the investigation of MS channels showed that modeling could enlighten key features in MS channel gating, and, up to now, such approaches were able to highlight important mechanisms involved in the lipid-force-driven gating of specific MS channels types. MD derived properties can also be used to develop the continuum description of mechanotransduction, based on viscoelastic modeling of the membrane coupled to mechanic models of ion channels. It is worth noting that the increasing complexity of gating mechanisms, from MscL to mammalian Piezo1, that can be described at the atomistic level, will contribute to developing more refined continuum models, with the advantage of simulating larger systems and a wider range of mechanical stimulations. Such descriptions can be eventually used to obtain macroscopic parameters to be included in coarse electrophysiological models of whole-cell mechano-gated ion currents. Concerning modeling of mechanosensation in neurons at the whole-cell scale and the tissue scale, we further showed that the microscale activation of MS channels in neurons could be coupled to the embedding tissue's biomechanics to achieve a comprehensive description of the whole mechanotransduction process. However, most of these multiscale descriptions are based on simplistic models of neuronal dynamics, mechano-electrical coupling, and tissue biomechanics and still need to be refined and generalized to achieve a higher degree of realism and descriptive/predictive power. In particular, other statistical mechanics tools, multiphysics couplings, and more detailed biomechanical and electrophysiological models can be introduced, as performed in modeling studies on different biological systems [224][225][226][227][228][229][230][231][232][233][234]. In addition, to achieve a better description of mechanotransduction, further studies focused on the electrophysiological dissection of mechano-gated current should be performed. In principle, this information can be merged with atomistic single-channels simulations and is crucial to built reliable biophysical models of neuronal dynamics at the whole-cell scale. Indeed, our research identified only a small amount of papers presenting those data but limited only to MS macro-currents or not reporting a full characterization of the ion currents with respect to the mechanical state. On this basis, future efforts should be devoted both to the generalization and development of atomistic, electrophysiological, and multiscale models, not only to grasp key features in neuron mechanotransduction but also to translate this knowledge to the design of bioinspired tools and devices.