Somatosensory neurons integrate the geometry of skin deformation and mechanotransduction channels to shape touch sensing

Touch sensation hinges on force transfer across the skin and activation of mechanosensitive ion channels along the somatosensory neurons that invade the skin. This skin-nerve sensory system demands a quantitative model that spans the application of mechanical loads to channel activation. Unlike prior models of the dynamic responses of touch receptor neurons in Caenorhabditis elegans (Eastwood et al., 2015), which substituted a single effective channel for the ensemble along the TRNs, this study integrates body mechanics and the spatial recruitment of the various channels. We demonstrate that this model captures mechanical properties of the worm’s body and accurately reproduces neural responses to simple stimuli. It also captures responses to complex stimuli featuring non-trivial spatial patterns, like extended or multiple contacts that could not be addressed otherwise. We illustrate the importance of these effects with new experiments revealing that skin-neuron composites respond to pre-indentation with increased currents rather than adapting to persistent stimulation.


Introduction
The sense of touch is a prime example of mechanotransduction in biology (Chalfie, 2009;Hoffman et al., 2011;Katta et al., 2015;Schneider et al., 2016) that culminates in the activation of mechanically-gated ion channels arrayed along sensory dendrites. These dendrites are not isolated from the tissues they innervate, nor are the relevant ion channels isolated from the extracellular matrix, plasma membrane, or the underlying cytoskeleton. In mammals, slowly-adapting mechanoreceptors depend on their association with keratinocyte-related Merkel cells for their response dynamics Maksimovic et al., 2014;Ikeda et al., 2014). Similarly, the response dynamics of rapidly-adapting mechanoreceptors are sensitive to their association with Pacininan corpuscles (Loewenstein and Mendelson, 1965). Extracellular links are present in the sensory neurons innervating hair follicles (Li and Ginty, 2014) and in dorsal root ganglion neurons in culture (Hu et al., 2010). Such protein tethers are also thought to be essential for mechanotransduction by vertebrate hair cells (Sakaguchi et al., 2009;Fettiplace, 2017). The NompC channels that mediate mechanosensation in Drosophila are directly linked to microtubules and, in campaniform receptors, this intracellular protein tether is essential for mechanosensitivity (Sun et al., 2019;Zhang et al., 2015;Liang et al., 2013). Thus, independent of their specific anatomy or mechanosensory function, the mechanically-gated ion channels that decorate sensory dendrites are intimately connected to surrounding tissues.
In C. elegans, the Touch Receptor Neurons (TRNs) are embedded in the skin and attached to a specialized extracellular matrix, and this structure is required for the proper distribution of MEC-4dependent Mechano-electrical Transduction (MeT) channels (Lumpkin et al., 2010;Emtage et al., 2004). The TRNs and the MEC-4 channels enable these roundworms to evade predatory fungi that trap nematodes in a noose-like structure (Maguire et al., 2011). This escape behavior can also be elicited manually by drawing an eyebrow hair across the animals body (Chalfie, 2014) or using mechanical stimulators (Petzold et al., 2013;Mazzochette et al., 2018). These observations suggest that laboratory stimuli are sufficiently good replicas of natural stimuli that they elicit the same behaviors. Forces in the nano-to micro-Newton range are sufficient to elicit this escape behavior (Petzold et al., 2013;Mazzochette et al., 2018) in wild-type animals. Sensitivity also depends on body stiffness such that larger forces are needed to trigger escape behaviors in stiffer animals (Petzold et al., 2013). Thus, touch sensitivity is a combined property of the skin-nerve tissue systems.
Delivering a touch by pushing a flexible probe against the worms body activates MEC-4-containing channels (O'Hagan et al., 2005), connecting touch stimulation directly to activation of a specific ion channel in living animals. This process of sensory mechanotransduction depends more on the depth of body indentation than it does on the force applied (Eastwood et al., 2015), reinforcing the importance of tissue mechanics in touch sensation. As found for the dendrites innervating Pacininan corpuscles, the TRNs depolarize in response to the application and removal of a simple touch (O'Hagan et al., 2005). In (Eastwood et al., 2015), we introduced a simplified, but quantitative description of sensory mechanotransduction that recapitulated this on/off response dynamic. This initial model introduced a hypothetical elastic tether connected to the channel that would be stretched in response to stimulus onset, relax during continued stimulation, and stretch in the opposite direction following stimulus offset. It was inspired by the tip-link model for auditory hair cells (Howard and Hudspeth, 1987;Hudspeth, 2014), but differs from this classical model in that it posits a tangential, rather than vertical, stretching of a tether. The picture emerging from this model replaces the hinged trapdoor of hair cell mechanotransduction with a sliding trapdoor. The tangential motion emerges from the mechanics of thin shells (Landau et al., 1986;Ventsel and Krauthammer, 2001;Audoly and Pomeau, 2011) and applies to the worm's body and its TRNs based on their anatomical position within the animal's skin (outer shell).
While appealing in its simplicity, the model in Eastwood et al. (2015) is incomplete: we replaced the ensemble of channels known to be distributed along TRN dendrites by a single effective channel. This simplification is similar in spirit to a mean-field approximation in physics and shares its utility for insight as well as its theoretical and predictive limitations. An important theoretical limitation was the neglect of the nonlinear mechanics of the worm's body. Additionally, the response to variations in contact areas or stimulus timing are not well described in this simplified model. The main object of the present study is to introduce a comprehensive and quantitative description linking focal mechanical stimuli to activation of single mechanically-gated ion channels, taking into account nonlinear mechanics and the spatial distribution of mechanically-gated ion channels. We evaluate the current model against prior experimental data, generating new insight into the contribution of internal hydrostatic pressure to touch sensitivity and as a major source of variation in experimental data. Additionally, we show that pre-indentation increases the response to subsequent indentation steps and use the model to reveal that this unexpected finding can be explained by the spatial distribution of tissue deformation and MeT channels. The approach underlying the present model and its evaluation by comparing simulated and experimental data could be adapted to other mechanosensory neurons that differ in their anatomy, their encapsulating tissues, and the distribution of MeT channels within the dendrites.

Results
To improve understanding of mechanosensory transduction during touch at the systems level, we develop a comprehensive and quantitative description of how a focal mechanical stimulus or touch activates single MeT channels in C. elegans touch receptor neurons, taking their spatial distribution into account. The following sections include models of how touch is transformed into skin deformation and how deformation activates single MeT channels as well as comparisons between predicted and experimental mechanoreceptor currents.
Non-linear mechanics of the nematode body and its role in converting touch into mechanical strain within the skin A nematode's body is a tapering cylinder (Figure 1) that consists of an outer and an inner tube separated by a fluid-filled pseudocoelom. The outer shell includes the cuticle, skin (hypodermis), excretory system, neurons and body wall muscles, and the inner shell is formed by the pharynx, intestine and gonad (Altun and Hall, 2009). Adult C. elegans hermaphrodites are about 1 mm in length and 50m in diameter, at their widest point. This simple body plan has inspired models of its mechanics consisting of a cylindrical outer shell and internal pressure, which is conferred by the combined effects of internal organs and the pseudocoelom (Park et al., 2007). Small punctures in the cuticle and skin are thought to decrease, but not eliminate internal pressure. Indeed, this maneuver has been demonstrated to decrease stiffness inferred from force-indentation curves derived from experiments using self-sensing microcantilevers (Park et al., 2007;Eastwood et al., 2015). These and other experiments use glass probes or microbeads with a radius 5-10 mm to indent the body to a maximum depth of~10 mm, that is about half the radius of the shell and larger than its thickness 1 mm.
The simplest physical model consistent with the above observations is that the strain within the shell is small, so that Hookean elasticity applies, yet the displacement of material points can be of the same order or larger than the shell thickness (see Appendix 1 for details). The latter implies that the linear approximation of the strain is not appropriate and must be replaced by the nonlinear Green-Lagrange expression: z y x p L 2R+t Figure 1. Scheme of the geometry in our model for C. elegans mechanics. The figure shows the scheme of a worm in a natural posture (left), straightened (as in neurophysiology experiments), and the model (right) that we shall consider here: a cylinder of length L ' 1 mm and radius R þ t=2 ' 25m is indented by a spherical bead (with radius 10 mm unless stated otherwise), applied here at its center. R is the radius of the middle surface and t is the thickness of the shell. Only half of the cylinder is shown for clarity. DOI: https://doi.org/10.7554/eLife.43226.002 where q j is the derivative with respect to the spatial coordinate x j (j ¼ 1; 2; 3) and u i is the i-th component of the displacement. The equations for the stress tensor s ij (Piola-Kirchoff of the second type) read (Landau et al., 1986;Ventsel and Krauthammer, 2001;Audoly and Pomeau, 2011): where sum over repeated indices is implied. Quadratic terms in Equations (1) and (2) make them appropriate for large deformations. Nonlinearities are called 'geometric' due to their relation to the shape of material elements (Ventsel and Krauthammer, 2001;Audoly and Pomeau, 2011). Linear approximations that neglect geometric nonlinearities can lead to substantial discrepancies, as shown below.
The linear Hookean relation between stress and strain s ij ¼ E 1þn " ij þ n 1À2n " kk d ij À Á is consistent even for large deformations provided components of the strain tensor stay moderate, which will be the case here (see Appendix 1). Here, E is Young's modulus, n is the Poisson ratio, and s ij is energy conjugate to " ij , that is the elastic energy E el upon a deformation varies as dE el ¼ R s ij d" ij dV. In addition to E and n, parameters of the model are (see Figure 1): the length L and the thickness t of the shell, the radius R of its middle surface, and the internal pressure p. The pressure p is understood to be the difference between the internal and the external atmospheric pressure. To simplify the following dimensional analysis, we shall neglect the external atmospheric pressure, which has only minor effects for our thin cylindrical shells (see Appendices 2 and 3). Elastic parameters are effective quantities that subsume different contributions in the inner and outer tube. Previous estimates of those parameters are discussed in the Section on Results for experimental validations. Boundary conditions generated by the pressure p and the external forces are in the next Section. Finally, Appendix 4 discusses the reduction of the 3D Equation (2) to the 2D thin-shell limit.
The nonlinear structure of Equation (2) hampers analytical approaches, which pushed us to apply numerical finite-element methods (see Burnett, 1985 for an introduction). Numerical simulations of Equation (2) were performed by the open-source program code-aster (Électricité de France, 2001). An hexahedral element with eight standard nodes (HEXA8) was used in combination with a mesh sensitivity analysis to verify that results are minimally sensitive to the element size. The numerical procedure was benchmarked and tested by comparing its results to known elasticity problems. In particular, Appendix 5 reports on the comparison with the deformation field and the force-indentation relation produced in small indentations of cylindrical shells, where an analytical solution is available (Morley, 1960), as well as large indentations of pressurized spherical shells, where a simplified equation was derived in Vella et al. (2012a). In all cases, agreement was verified.
For the internal pressure p, active readjustments of the internal pressure are possible and could a priori be accommodated in our approach. Here, we shall make the simplest working hypothesis that p holds constant when the stimulus is exerted onto the body of the worm. Its justification is empirical, that is we argue that the simplest option is sufficient based on results reported below.
As for boundary conditions, neurophysiology experiments have worms glued onto a plate and limited in the vertical displacement of their body's lower half (see Materials and methods). Since a mathematical formulation is not obvious (and elasticity has long-range effects), we tested different boundary conditions and present two of them for comparison (see Appendix 6). For the first, the lower half of the body is vertically rigid, that is the upper half of the cylinder in Figure 1 is free to move, while the lower one is constrained to move only parallel to the plane onto which the worm is glued. For the second boundary condition, the lower half is fixed, that is allowed to move neither vertically nor laterally. Results presented in the main text are obtained with the latter condition. As for the two ends of the cylinder, we present results with vanishing lateral forces; Appendix 7 discusses the effects of plugs at the ends.
Single MeT channel gating and its relationship to total mechanoreceptor currents deformations described above. Because the touch-induced deformation and the channels are spatially distributed, we developed an approach that accounts for the activation of a single MeT channel by local deformation and subsequently computes their summation based on the spatial distribution of both the deformation and the channels. This model departs from previous work (Eastwood et al., 2015) by considering the contribution of individual channels to the total current and taking into explicit consideration the spatial features of the activation process.
The mechanism in Eastwood et al. (2015) posits that the dynamics of individual channels is the combination of an elastic and a relaxation (frictional) component. While various implementations may be contemplated, we shall refer for concreteness to a situation where each ion channel is connected to an elastic filament. We denote by r c;f the undeformed positions of the channel and the tip of its elastic filament; the corresponding displacements induced by the deformation of the embedding tissue are Dr c;f .
The elastic component reflects the Hookean response of the filament to its stretching. Elastic energy is VðxÞ ¼ k 2 x 2 , where x ¼ Dr f À Dr c is the elongation of the filament with respect to its undeformed configuration. The corresponding restoring force is F elastic ¼ Àkx : ( As for the frictional component, the TRN and its channels are embedded in the medium and expected to move with it, that is Dr c ¼ uðr c Þ where u is the displacement in Equation (2). Conversely, as the filament slides with respect to the medium, the friction force is where g is the friction coefficient and r 0 is the undeformed position of the material point that coincides with the location of the tip, that is r f þ Dr f ¼ r 0 þ uðr 0 Þ. Expanding uðr 0 Þ and using that gradients of u are small, we obtain r 0 À r f ' Dr f À u r f À Á , which is then inserted into Equation (4) to show that uðr 0 Þ can be replaced by uðr f Þ. Effects of inertia are negligible and the overdamped approximation holds at microscopic scales (Phillips et al., 2012), that is the sum of the forces F friction þ F elastic ¼ 0, which yields where t ¼ g=k is the relaxation time, x is the extension of the filament, and Dr c ¼ uðr c Þ was used. Equation (5) drives x to zero for G constant, which is the basis for adaptation. Equation (5) is supplemented by the constraints exerted by the neural membrane around the channel, which limit the motion in the vertical direction. The constraint can be written as x Áŵ 3 ! 0, where, for every channel, w 1;2 span the plane locally tangential to the neural membrane whileŵ 3 indicates the orthogonal direction. Specifically (see also Appendix 8), we define an orthonormal basisê 0 i as follows:ê 0 y is aligned with the local direction of the (deformed) axis of the cylinder running head-to-tail;ê 0 z is orthogonal to the neural membrane at the top of the TRN, and oriented outward;ê 0 x is tangential to the neural membrane, along the remaining direction of a right-handed system. The basesŵ i are constructed by rotating theê 0 i appropriately. For a channel placed at the top of the TRN, the local basiŝ w i coincides withê 0 i . If the channel is rotated by along the surface of the TRN, then (5), is obtained by the displacements u calculated using the mechanical model. The relation between stretching and " is (Landau et al., 1986;Audoly and Pomeau, 2011): G 2 ðtÞ À G 2 ð0Þ ¼ 2" ij ðtÞG i ð0ÞG j ð0Þ, where the G's are again assumed to be small.

The opening/closing dynamics of channels
Channels can be in multiple states: open, closed and several sub-conducting (Brown et al., 2008). Experiments and effects discussed here are captured by including a single sub-conducting state between the open and closed states (C* )S * )O). The respective probabilities P c;s;o obey the master equation dPc dt ¼ ÀR cs P c þ R sc P s ; where R ij are the respective transition rates, and R c;o ¼ R o;c ¼ 0 again to minimize free parameters. The channels are posited to work at equilibrium, so that where DG ij is the free energy difference between the states i and j, b ¼ 1=k B T, T is the temperature, and k B is the Boltzmann constant (see, e.g., Phillips et al., 2012). Channels are coupled to mechanics via their elastic filaments described by Equation (5). Namely, the extension of the filament modulates the free energy differences among the above states of the channels: where g 0 , g 1 are dimensional constants, and F is the amplitude of the tangential component of F elastic in Equation (3): where, for every channel,ŵ 1;2 span the plane locally tangential to the neural membrane whileŵ 3 indicates the orthogonal direction, as defined above.
Choices for the free energy other than Equation (8), for example a quadratic dependence on x, are discussed in Appendix 9. The free energy of the intermediate subconductance state has a priori its own parameters. However, to reduce free parameters, its free energy is assumed intermediate between the closed and the open state in Equation (8) DG os ¼ aDG oc ; DG sc ¼ ð1 À aÞDG oc ; with the only additional parameter 0 a 1. The ability of the model to quantitatively describe experimental data supports Equation (10) as a good empirical description of the free energy of the intermediate subconductance state. Ion channels are believed to be distributed in spots ('puncta') along the neural membrane (O'Hagan et al., 2005). Their distribution is consistent with uniformity in the angular and longitudinal directions, while spacings between successive puncta are distributed log-normally (Cueva et al., 2007). For simplicity, each punctum is assumed to contain a single channel.
The current along the TRN is the sum I ¼ P k i k of the currents of individual channels. Its mean is given by The non-linear elastic model estimates mechanical parameters that agree with experiments Next, we sought to determine the aspects of measured mechanoreceptor currents dynamics that are captured by this quantitative model incorporating body mechanics, single MeT channel gating, and the spatial distribution of MeT channels in touch receptor neurons. To achieve this goal, we compare simulations and experiments for both mechanics and neural responses. The first step for a proper comparison with experiments is an appropriate choice of the no-stress state: the corresponding length, thickness and radius should be such that the pressurization of the cylinder leads to the values relevant for experiments. In particular, if we want to keep the final (at pressure p) values fixed, the no-stress initial values should change as p varies. This point, as well as our below results, differ from Elmi et al. (2017), where an elastic model is discussed, yet the role of the no-stress state and pressure are not considered. Initial (no-stress) values are conveniently obtained by using perturbative analytical formulae in Appendix 2, which give the variation of various quantities with p.
The schematic of an indented shell in our numerical simulations is shown in Figure 1. Note that the size of the indenter is not negligible with respect to other dimensions, and the region of contact with the cylinder is expected to change with the indentation depth (Johnson, 1985).
The thickness of the shell t can be rescaled out, as discussed below in the analysis of bending and stretching contributions, and all geometric parameters appearing in Equation (2) are fixed. The variable factors are p and E, which can enter the deformation for a given indentation only via their nondimensional ratio p=E. We plot then in Figure 2B the dependence of the deformation profile along the longitudinal coordinate vsp=E. Vertical deformation is strongest at the center of the indenting bead, and its longitudinal extension decreases with p=E. The best least squares fit yields p=E ¼ 0:01. Figure 2C shows how the estimate p ' 40kPa for the internal pressure was obtained: we fix the ratio p=E ¼ 0:01, predict the relation force-indentation as p varies, and make a best fit to the experimental data. The value of 40 kPa is on the same order of magnitude as the range of pressures measured in the larger nematode Ascaris lumbricoides (Harris and Crofton, 1957). The two above estimates yield for the Young's modulus E~4MPa, which is the same range as the 1.3MPa obtained by measuring the bending stiffness of the nematode (Backholm et al., 2013) or values~10MPa obtained in Fang-Yen et al. (2010) and Petzold et al. (2011). Our estimate differs from the much higher values in Park et al. (2007), which were also obtained by indentation data, yet using formulae of linear elasticity that are only valid for indentations ( t.  Having fixed the parameters of our model, we can now independently test it against data on the mechanical response of C. elegans to changes in the external pressure (Gilpin et al., 2015). The variation DV of the initial volume V 0 was found to depend linearly on the variation of the external pressure Dp, and the resulting bulk modulus k ¼ Dp Performing the same operations in our simulations, we obtained estimated values of k ¼ 150 À 230kPa. This agreement is quite significant as we derived k, a global mechanical property, by using parameters inferred from local indentation measurements. Finally, Appendix 11 shows that mutations in the cuticle induced by disruptions of the lon-2 gene should modify the bulk modulus, contrary to suggestions in Gilpin et al. (2015).

Predictions and experiments for responses to pre-indented stimuli
With our model reflecting the mechanical properties of the worm, we can now estimate the forces transferred to each individual channel along the TRN. By summing these individual responses to calculate the total TRN response, we can test the model's ability to explain neural responses recorded in TRNs in vivo. Predictions developed using this model recapitulate experimental responses that could not be addressed by our prior model (Eastwood et al., 2015), highlighting the importance of the new elements introduced here.
A detailed discussion of micro-cantilever systems of stimulation, and the in vivo patch clamp system to record neural responses was presented in Eastwood et al. (2015). A summary, together with specific differences for the data first reported here, is in Materials and methods. Instances of stimuli and neural responses from Eastwood et al. (2015) are shown in Figure 3 and neural responses to new, pre-indented profiles are shown in Figure 4.
Let us then describe how model predictions are obtained. The pressure p and the parameters of the mechanical part are as in the previous Section. Positions of the channels along the TRNs were randomly generated according to the log-normal experimental distribution (Cueva et al., 2007), and results were averaged over those statistical realizations. As for the elastic filaments described by Equation (5), their length is rescaled to unity as discussed in Appendix 12, while their initial direction Gð0Þ is distributed randomly. Namely, directions were generated uniformly in the semisphere with a non-negative component along the local outward normalŵ 3 to the neural membrane. Based on the deformation field determined by Equation (2), we used Equation (5) to compute the force on the channels as a function of time, and obtained the dynamics of the channels via Equation (6). More details on the fits and the resulting values of the parameters are in Appendix 12.
Our results in Figure 3 manifestly capture the symmetric and rapidly adapting response of TRNs. Because of the onset-offset symmetry of the touch response, the response to sinusoidal stimuli oscillates at twice the input frequency. At high frequencies, inertia in the switch between open and closed states of the channels contributes to the reduced amplitude of the oscillations. The response to ramps is intuitive: the slower the indentation, the smaller the response because of adaptation.
Simulations also capture the empirical relationship between speed and current amplitude ( Figure 3c). Appendix 13 presents the histogram of the errors for individual realizations, which shows that neural responses are captured at that level as well (not just the mean, as in Figure 3). The histogram also shows that restricting filaments to be initially or permanently tangential (Gð0Þ Áŵ 3 ¼ 0 or xðtÞ Áŵ 3 ¼ 0), further improves results. The latter restrictions being speculative at this stage, we shall focus on the unrestricted model; we only note that tangential restrictions admit plausible molecular mechanisms, for example. by microtubules that run along TRNs, are attached to the neural membrane through filaments and are known to impact touch sensation (Bounoutas et al., 2009).  . Pre-indented steps yield stronger responses due to their more extended deformation profile. (A) Stimuli delivered for standard (blue) and pre-indented (purple) steps. Black arrows indicate the total displacements for the on-currents in the following panels (colors match). Green arrows indicate relative displacements. Experimental stimuli and neural responses from ALM neurons were obtained as detailed in Materials and methods. We recorded from 11 separate worms with 3 À 11 presentations of each stimulus per recording. Recordings were only included if they met the criteria outlined in the Data Analysis section of the experimental methods, which led to a final number of biological replicates per displacement point that varied from 5 to 11. Representative traces shown here are from one biological replicate. (B) The on-current vs the total displacement (the preindentation for the purple points is 5 mm). Dotted curves and diamonds (in this panel and the following ones) report the prediction of our previous Mean Field (MF) model in Eastwood et al. (2015). The goal is to stress the importance of spatial integration effects, which constitute the main contribution of this paper and were neglected in Eastwood et al. (2015). (C) Off-currents are statistically indistinguishable, as expected since offsteps are identical and adaptation erased the memory of the pre-step. (D) The on-current vs the relative displacements. Note the stronger response for pre-indented stimuli. (E) Changes in the profile of deformation: Dz is the difference between the deformations after and before the (relative) stimuli. Note the greater extension for the pre-indented case, which is the reason underlying results in panel D. Open circles in panels B-D were normalized to the maximal currents detected and show the mean ± s.e.m. DOI: https://doi.org/10.7554/eLife.43226.006 The following source data is available for figure 4: Source data 1. Experimental and simulated neural responses to pre-indented steps. DOI: https://doi.org/10.7554/eLife.43226.007 Additional insight is gained by delivering stimuli alternative to the classical profiles in Figure 3, namely pre-indented stimuli in Figure 4. Panel A contrasts standard and pre-indented steps, that is where an initial step (5 mm in our data) is delivered. The neural response to two steps of equal amplitude, one pre-indented and the other not, is substantially stronger for the former. That is surprising at first, since the amplitude of the steps is identical, and enough time between the successive half steps was left for adaptation. The explanation was obtained by using our model: it is indeed the case that channels adapt and return to their rest state; however, the tissue is deformed by preindentation, which leads to a more extended region of stimulation and more channels activated, as shown in Figure 4E. The previous mean-field model in Eastwood et al. (2015) was unable to account for this increase in the number of channels reached by stimulation with pre-indentation ( Figure 4B,D). The resulting predictions reproduce experimental trends, highlighting the importance of the coupling between mechanics and channel activation that constitutes the main focus of our paper.

Variance in residual internal pressure following dissection accounts for variation in responses among individual worms
The full model not only allows us to predict neural responses to complex stimuli, but also to delve further into how body mechanics can explain the variation in experimental responses. The dissection procedure required for recording from TRNs in vivo necessarily alters body mechanics: a small incision allows some portion of the tubular internal organs (intestine and gonads) to be released outside the animal. The cuticle is largely re-sealed by the remaining internal pressure pushing large organs over the hole, and a second incision is then made to release the TRN cell body without other organs. This standard procedure results in 'soft' worms with varying fractions of their internal organs released. A modified dissection also used in Eastwood et al. (2015) omits the first incision and release of organs, resulting in 'stiff' worms. Names stem from the force-indentation curves in Figure 5A, which evidences that the latter procedure better preserves the body's integrity. Most experimental recordings reported here are for 'soft' worms while data for 'stiff' worms in Figure 5A were used to predict the ratio p=E ( Figure 2C). Eastwood et al. (2015) empirically showed that neural responses for soft and stiff worms are similar for displacement-clamped stimulations, while they strongly differ for force-clamped protocols. This Section analyzes the mechanical consequences of the above procedures, their effects on neural responses, and explains the empirical observation in Eastwood et al. (2015).
Since internal organs of soft worms are removed away from the stimulation point, it is plausible that the dissection affects the internal pressure and has weaker effects on the indented external shell. This suggests to conservatively keep (in our model) the Young's modulus E and the thickness t fixed, and modify p. Results of the corresponding simulations are shown in Figure 5A. The slope of the force-indentation relation decreases with p; the best fit for soft worms is p = 1.6 kPa, which is 4% of the value for stiff worms. Finally, the scatter around the mean shows that the more invasive dissection procedure results in stronger variability, with the corresponding p ranging from 0.04 to 16 kPa.
To further analyze the effects of the dissection procedure, Figure 5B shows the longitudinal profiles of vertical deformation for soft and stiff worms. The point is that the curves differ much less when displacement is clamped, rather than force. That is translated into predictions for currents as follows. We assume that the dissection procedure does not affect the channels and calculate their respective currents as described previously. Specifically, we fix the distribution of the channels, change p to the value corresponding to soft or stiff worms, and compute the neural response to force or displacement-clamped stimuli. Results are shown in Figure 5C: responses for force-clamped stimuli widely differ for soft and stiff animals, yet they are similar for displacement-clamped stimuli. In sum, empirical observations reported in Eastwood et al. (2015) are explained by the mechanics of the nematode and its coupling to neural activation.
Finally, we address the variability in Figure 5C among soft worms, which we tentatively related to p varying over three orders of magnitude. For further support, we tested whether the observed variability could indeed be reproduced by keeping all parameters fixed but p. Results in Figure 5D show that the peak current increases systematically with p for displacement-clamp and has the opposite behavior for force-clamped stimuli. The predicted change is larger in the former case, which is consistent with differences between soft and stiff worms. The model predictions are compared with an experimental dataset of 21 worms obtained in Eastwood et al. (2015). For each worm we inferred p from the force-indentation relation, pooling together animals with similar trends. Figure 5D supports the initial hypothesis that differences in p among dissected animals are a major component in the observed variability.

Testable model predictions for future experiments
In the previous sections, we were able to validate many of our predictions using existing or new data, increasing our confidence in the model. The following subsections illustrate how the model can be used to make predictions to be tested in future experiments.
Shell bending is weak compared to stretching; stiffness is dominated by internal pressure The mechanics of pressurized shells relies on the balance between the internal pressure p, bending and stretching of the shell. Contradictory results have left undecided the previous balance for C. elegans (Park et al., 2007;Gilpin et al., 2015). Here, we exploit our model to clarify this issue.
We computed the vertical deformation for different values of p=E and t, as previously done for the validation of the mechanical model. The longitudinal extension is quantified by the distance y h for the deformation to reduce to half of its maximum value (at the center of the bead).  shows that y h decreases, that is the deformation is more localized, when p=E increases. Conversely, as t reduces, the deformation is wider if p=E < 10 À4 and narrower if p=E > 10 À4 . To gain insight regarding the consequences of Figure 6, we can use the thin-shell limit of Equation (2) in Appendix 4. Reducing the limit equation to non-dimensional form as previously done for spheres (Vella et al., 2012a), we find that the bending term is multiplied by the factor 1=t 2 E 2 t 4 =p 2 R 4 . If t ) 1, the bending term is small, and (with the possible exception of boundary layer regions) the only remaining dependence on t is via the stiffness S Et. These arguments suggest to plot y h vs p=Et as in Figure 6: curves with different t indeed collapse for the values of p=Et that are relevant for experiments. We conclude that internal pressure and stretching of the shell provide the dominant balance.
We next compared the elastic energy of the shell with the work by the external forces. Results in Figure 6C show that their ratio reduces as p=E increases, and the elastic energy tends to become marginal, which illustrates the dominance of the internal pressure in the body stiffness.

Mechanical and neural responses depend on the radius of the indenting bead
Previous research has noted that the amplitude of neural currents depends on the radius R b of the indenting bead (O'Hagan et al., 2005), but no systematic study has been made of this relationship. We fill this knowledge gap with simulations of how bead size and internal pressure interact to affect the deformation of the worm and thus the neural currents produced.
Results are shown in Figure 7: R b influences both the deformation profile ( Figure 7A) and the force-indentation curve ( Figure 7B) for p=E in the experimentally relevant range. The curves are intuited as follows. The curvature of the deformation field at the indentation point increases with p=E until it matches the radius R b . As p=E increases further, the shell cannot become any steeper (the bead is rigid), so that it adapts to the bead in the contact region (see Figure 2B), which widens with R b . The radius R b also controls the deformation outside of the contact region, namely the mid- maximum extension of the deformation (/ R b , data not shown). As p=E reduces, the deformation becomes shallower at the indentation point and the role of R b vanishes (see Figure 2B). Similarly, Figure 7A shows that the volume of the body to be deformed increases with R b , and more work is needed for a given maximum indentation w 0 . In formulae: the work Fdw 0 by the indenter roughly balances the contribution of the internal pressure ÀpdV (the elastic energy is small, as discussed previously) and the force-indentation relation is then given by F / Àp dV=dw 0 .
Qualitatively, larger dV's associated with larger beads yield the nonlinear dependence in Figure 7B. Quantitatively, we write F ¼ p, where ¼ ÀdV=dw 0 has the dimension of a length squared and depends on w 0 , the shell and bead radii R, R b , and the ratio p=E. Keeping the latter fixed, we investigated numerically the behavior of and observed that the ratio Fðw 0 Þ=F max does not depend on R b (see Figure 7C). It follows that the dependence on R b should factorize out: ¼ G 1 ðR b =RÞG 2 ðw 0 ; RÞ where G 1 depends on R b =R for dimensional reasons. The function G 2 brings the length squared dimensionality and, in the limit of small R b and large p, behaves as R w 0 (Vella et al., 2012b). It follows that G 2 ¼ Rw 0 G 3 ðw 0 =RÞ.
The above functions G 3 and G 1 are determined as follows. We computed numerically the force indentation relation of cylinders of different radii (R = 25, 40, and 50 mm) to stimulations produced by beads of different size (R b from 3 to 10 mm); results of the simulations are then used to fit coefficients of the Taylor expansions of G 1 ðxÞ and G 3 ðxÞ. Using this approach we find that the functional form with a 1 = 0.76, a 2 = 2.1, and a 3 = 0.66, captures quantitatively the behavior of the force indentation relation ( Figure 7B, R 2 = 0.995). Variations between R b = 3 mm and R b = 10 mm are on the order of few mN, hence they should be accessible experimentally. Equation (12) generalizes the linear relationship, valid in the limit of very small R b , obtained in Vella et al. (2012b).
Consequences for neural responses are in Figure 7D. deformation decays rapidly as one moves circumferentially from the north pole (where the bead is indenting the body) toward the equator. That implies an appreciable dependence on the angular position of the TRN with respect to the bead, which will be stronger for smaller beads as the extension of their deformation is reduced.
Similar tangential forces at stimulus onset and offset drive symmetric on/off responses Thus far, we have treated the fact that TRNs respond to both the application and release of a step stimulus as a given. Yet channels in many other mechanosensitive systems respond preferentially to stimuli in a particular direction -either on or off -rather than responding symmetrically to both (see Katta et al., 2015 for review). Here, we further analyze the origin of this symmetry, by calculating the stimuli upon the channels and analyzing differences among microscopic gating mechanisms that are consistent with the symmetry. A first key remark, which generally applies to thin shells (Landau et al., 1986;Audoly and Pomeau, 2011), is that the off-diagonal components " xz and " yz are small compared to the rest of the components of the tensor, namely the tangential ones. Indeed, those two off-diagonal terms are proportional to the corresponding components of s, which vanish due to the thinness of the shell (see Landau et al., 1986;Audoly and Pomeau, 2011). It follows from the definition of " (see Appendix 8) that vectors initially tangential and perpendicular to the surface of the cylinder, remain orthogonal even after deformation.
In addition to the general above property, the component xy is also negligible when the indenting bead is applied on top of the cylinder. The strain tensor is then diagonal, as confirmed by Figure 8B.
The force acting on a single channel, as defined by Equation (9) and calculated using Equation (5), is shown in Figure 8C. The force is maximal if the elastic filament is initially in the tangential plane while orthogonal filaments generate negligible forces. Notably, forces for tangential filaments have opposite signs yet very similar amplitudes at the onset and offset. The relation between vertical and tangential directions is key to the onset-offset symmetry and stems from the above discussion on thin shells. That constitutes the physical reason for our positing that tangential stimuli gate the channels: the orthogonal dynamics is indeed affected by the neural membrane, which a priori prevents any symmetry between inward and outward extensions.

Models with symmetric or directional channel populations could be distinguished experimentally
Though the forces reaching the channel at stimulus onset and offset are similar in amplitude, they are opposite in direction ( Figure 8B). Two alternative mechanical models could then explain the observed symmetry: a 'symmetric' model in which each individual channel responds to force in both directions, and a 'directional' model in which individual channels respond preferentially to force in one direction, but the population as a whole responds to both. Namely, alternatively to the isotropic choice in Equation (8), we could consider the 'directional' model with the preferential direction v: Contrary to Equation (8), Equation (13) breaks the symmetry for individual channels (see Figure 8D), which can be restored though for the total current if channels along the TRN have their directions v in Equation (13) independently and isotropically distributed ( Figure 9A).
A more quantitative analysis leans toward the symmetric model. Indeed, for the experimental value (O'Hagan et al., 2005;Brown et al., 2008) of the single-channel current i 0 ¼ À1:6 AE 0:2pA, Equation (13) underestimates the mean current ( Figure 9A). More generally, we can optimize parameters and calculate the errors in the fits of the experimental datasets: the probability distribution for Equation (13) is broader and shifted to higher errors with respect to the isotropic model (see Figure 9B). Similar conclusions hold if i 0 is allowed to vary ( Figure 9C).
In sum, the analysis of available experimental data favors symmetric channels, but is not fully conclusive. New data will be needed, which is our motivation to describe hereafter two possible experiments.
A first approach relies on the noise level of currents. The intuition is that anisotropy reduces (for a given density of channels) the number of active channels along the TRN, and thereby leads to more noise. Specifically, the number of active channels could be inferred (Appendix 10 includes a generalization of the noise analysis in O'Hagan et al., 2005 to non-equally-stimulated channels) and compared to the number of channels measured by fluorescent tags (Cueva et al., 2007). Figure 9D presents the Coefficient of Variation (CV) of the TRN current vs the stimulus strength, calculated over repetitions of a given stimulus. Differences in CVs are poised to permit discrimination and the approach described in Appendix 10 estimates that~100 trials suffice for their reliable measurement.
A second alternative exploits the architecture of C. elegans neurons: TRNs extend longitudinally for about half of the nematode's length, leaving a region around its center that is relatively insensitive to touch (see Mazzochette et al., 2018). Figure 9E indicates the range over which effects of indentation are felt by individual channels; panel F shows the differences between microscopic models as the indenting bead slides along the longitudinal direction. An additional relevant statistic is the asymmetry between on-and off-currents. The logic is that, as the number of stimulated channels The sketch on the right illustrates that directional channels respond only to stimuli properly aligned with respect to their preferential direction while symmetric channels respond isotropically. DOI: https://doi.org/10.7554/eLife.43226.012 decreases, asymmetries should become more substantial if the channels are not isotropic, see Figure 9G. An appealing possibility is the stimulation of the worm in its center (negative coordinates in panels F,G). There, the number of activated channels is small, which could indeed bring microscopic insight.

Discussion
We present a quantitative description of the response of C. elegans touch receptor neurons to simple and complex mechanical cues. This work combines modeling and simulations of how touch deforms the skin and its embedded TRNs with a detailed model of the activation of single MeT channels, linking skin indentation to neuronal strain and MeT activation. This model explains several facets of the coupling between tissue mechanics and neural responses that were not previously understood.
Our model explains several aspects of the coupling between mechanics and neural responses. First, the model replicates experimental observed currents evoked by a wide range of indentation profiles. The prior model (Eastwood et al., 2015) relied on a mean-field approximation of channel activation, which could not properly account for the fact that pre-indentation increases the number of channels that contribute to the total current. Second, the model explains how the mechanics of skin deformation contribute to the empirical observation previously made in Eastwood et al. (2015) that neural responses are variable for force-clamped stimulation but less so for displacement clamp. It shows that variation in internal pressure resulting from the dissection procedure is a major source of variability in the neural responses of individual nematodes. Third, the model predicts that the neural response should increase with the size of the indenter in a manner dependent on internal pressure. Finally, the variation in deformation around the circumference of the worm suggests that the angular position of the TRN with respect to the indentation bead affects the response. These findings can help design future experiments to include measurements of these important parameters.
Our analysis also provides insight on two points related to the mechanics of the worm and the underlying biology. These findings reconcile conflicting results (Park et al., 2007;Gilpin et al., 2015) on mechanical properties in wild type and mutants. Specifically, we showed that the mechanical response of the nematode is captured by an elastic cylindrical thin shell in a pressure-dominated regime. The theory makes testable predictions on the dependence of the force-indentation relation on the indenter size, as well as the effects of mutations in the cuticle on the bulk modulus. An experimental verification of those predictions would further support the major role of internal pressure in C. elegans body mechanics. The fact that our model yields best results for boundary conditions with no force at the two ends of the cylinder, suggests that the body of the nematode might relax longitudinal components of the stress. One possible mechanism is through the annular structure of the cuticle (Altun and Hall, 2009), which may be effectively described as a shell with anisotropic Young's moduli. This relaxation of longitudinal stresses could facilitate the bending motions required for nematode motility.
We have shown how neural responses have similar amplitude when the stimulus is applied or released, which addresses the long-standing puzzle of the onset-offset symmetry. The picture emerging from our work is that the most likely gating model involves displacement tangential to the neuronal membrane. This insight can guide the search for the biophysical basis of in vivo activation of MeT channels and differs fundamentally from other models that involve orthogonal components (Howard and Hudspeth, 1987;Chalfie, 2009;Hudspeth, 2014). At the microscopic level, we inquired whether the symmetry holds for individual channels or at the population level only (individual channels are asymmetric yet their preferential directions are randomly distributed and their cumulative effect is again symmetric, as suggested for Drosophila sound receiver [Albert et al., 2007]). We showed, however, that a model with symmetric channels gives a better description of existing data, hinting at symmetry for single channels. Definitive evidence could be obtained by the experiments that we suggested in Figure 9C and G, with the noise levels better controlled and the stimulation point moved along the longitudinal axis so as to assay a variable number of channels. The ideal experiment would be to precisely assay the neural response to stimuli in the central deadzone of the body, where few channels are likely to be directly stimulated (see Mazzochette et al., 2018).
In our current description, we assumed that the material composing the shell is purely elastic and the dependencies on frequency result from the gating of the channels. While this procedure successfully captures many experimental observations, it is known that tissues do feature viscous effects (Backholm et al., 2013). Future developments will address viscoelastic effects, which should be relevant to the understanding of touch sensation at high frequencies.
Finally, it is worth noting that our modeling ultimately relies on the fact that touch receptor neurons are close to the surface of the skin's thin layers. This leads to physical effects peculiar to thin shells, namely the importance of tangential forces, which are at the basis of the gating mechanism discussed here. Since the above features are common in touch sensation, we expect results and methods that we developed to be widely relevant.

Materials and methods
We incorporated most of our modeling methods in Results and Discussion. Numerical simulations were performed as discussed in Results by the open-source program code-aster (Électricité de France, 2001). For additional details, please see Appendices. Experimental methods are found hereafter.

Imaging Cuticle Deformation
TP12 worms were immobilized with 0.1 mm polystyrene beads on a 6% NGM agarose pad. 10 mm glass beads (Duke Scientific) for indenting the worms were spread onto a coverslip, which was inverted to cover the agarose pad holding the worms. To image the worms, we used a high-magnification camera (Orca-R2, Hamamatsu) on an inverted microscope (Leica) with an EGFP filter set and a high-numerical aperture 63x oil immersion lens, to yield a shallow depth of field » 0.1 mm for optical sectioning. When glass beads were trapped between the cuticle of the animal and the coverslip, we were able to capture fluorescence images of COL-19::GFP in the cuticle at >10 different focal planes. At each focal plane, we measured the radius of the bead and the radius of the cuticle deformation (by identifying where the cuticle was in focus). We then calculated the depth of the plane based on the radius of the bead at the focal plane. Experimental data shown are a combination of all focal planes for two adult animals.
Membrane current and voltage were amplified and acquired with an EPC-10 USB amplifier and controlled through Patchmaster software (HEKA/Harvard Biosciences). The liquid junction potential between the extracellular and intracellular solutions was À14 mV and was accounted for by the Patchmaster software. Data were sampled at 10 kHz and filtered at 2.9 kHz.
Electrophysiology source data from Eastwood et al. (2015) are available upon request.

Mechanical stimulation
For mechanical stimulation during patch-clamp electrophysiology, previous studies used either open-loop systems with a piezoelectric bimorph (O'Hagan et al., 2005) or stack (Bounoutas et al., 2009;Arnadó ttir et al., 2011;Geffeney et al., 2011;Chen et al., 2015) with no measurement of actual displacement or a closed-loop system with a stimulus bead at the end of a piezoresistive cantilever for force detection, driven by a piezoelectric stack (Eastwood et al., 2015). Here, we use an open-loop system adapted from the piezoelectric stack system with a photodiode motion detector described in Peng et al. (2013). This enables faster stimulation than the force-clamp system (Eastwood et al., 2015;Petzold et al., 2013) at the expense of control over and measurement of exact force and indentation. The photodiode detector allows for a readout of the time course of the displacement of the stimulator. An open-loop piezoelectric stack actuator with 20 mm travel distance (PAS-005, ThorLabs) was attached with marine epoxy (Loctite) to a 0.5'' diameter, 8'' length tungsten rod, and mounted on a micromanipulator (MP-225, Sutter) at a 17˚angle to allow the stimulator to fit beneath the microscope objective.
For detecting probe motion at the 0.5-10 mm scale, we adapted the system from Peng et al. (2013) to use the SPOT-2D segmented photodiode (OSI Optoelectronics), and mounted it in an XY translator on top of a rotation stage (ST1XY-D, LCP02R, ThorLabs) to enable alignment of the photodiode gap perpendicular to the direction of probe motion. This was affixed above a secondary camera port on the microscope (Eclipse E600FN, Nikon) with no additional magnification.
To create a defined and reproducible contact surface for the stimulation probe, we adapted the bead gluing technique used previously for the force-clamp system (Petzold et al., 2013;Eastwood et al., 2015), but with an opaque bead that allowed for a clear signal from the photodiode motion detector. Borosilicate glass pipettes (Sutter, BF150-86-10) were pulled and polished to a tip diameter of 10-15 mm, and 20-23 mm diameter black polyethylene beads (BKPMS-1.2, Cospheric) were attached with UV-curable glue (Loctite 352, Henkel). Pipettes with attached beads were trimmed to a length of 1-2 cm, placed in the pipette holder, and waxed in place with sealing wax (Bank of England wickless, Nostalgic Impressions). A high-resolution 3D-printed acrylic pipette holder (custom design) was attached with marine epoxy to a steel tip (PAA001, ThorLabs) mounted on the piezo stack.
After cell dissection, but before making a gigaseal for patch clamp, the front edge of the stimulator bead was moved into place and visually aligned under the 60X objective with the highest visible edge of the worm's cuticle at a distance of 108 ± 36 mm anterior to the ALM cell body.

Stimulus control and data acquisition
All systems described here were controlled through HEKA Patchmaster software with a 10 kHz sampling frequency. The voltage output from the EPC-10 amplifier (HEKA) was adjusted based on the total range of the stack for a relationship of 0.418 V/mm. This command signal was filtered at 2.5 kHz on an 8-pole Bessel filter (LPF-8, Warner Instruments) and then amplified with a high-voltage, highcurrent Crawford amplifier (Peng and Ricci, 2016) to achieve a signal between 0-75V which was sent to the stack. The stack was biased with a starting offset of 3-4 mm, and the largest displacement used was 3-4 mm less than the upper limit of the stack's travel distance, ensuring that stack motion was linear. The analog signal from the photodiode circuit was digitized at a rate of 10 kHz by the EPC-10 amplifier and Patchmaster software, for temporal alignment of the probe motion signal with the evoked current response.
Only recordings with holding current <-10pA at -60mV and series resistance <210MW were included in the analysis. Since the voltage was not changed during the course of these experiments, we did not correct for voltage errors due to uncompensated series resistance. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. In the main text, we exploited the property that gradients of the deformation field u (which are non-dimensional quantities) are small compared to unity. This was used first to justify a Hookean relation between stress and strain, and then in the derivation of the friction force acting on the elastic filaments. The goal of this Appendix is to provide an empirical a posteriori validation of this assumption.

Author contributions
Appendix 1- figure 1 shows the various gradients components for an indentation of 10 mm, which is the greatest value in experiments. To provide an upper bound, we focus on the region of maximum deflection, that is along the longitudinal direction at the top of the cylinder. Over a wide range of positions, both in soft and stiff worms, the gradients are indeed smaller than unity. The only exception is the component du z =dy, which approaches unity in a small region around y = 5 m, which is where the bead detaches. However, it is sensible to neglect even this contribution for predictions of the neural current. Indeed, values in Appendix 1- figure 1 are an upper bound, and the region is just a few m wide, hence only a small number of channels are possibly concerned (we remind that the average inter-channel distance is~2 m). Note also that duy dz has a negative value that compensates for duz dy in the yz component of the strain tensor, which could otherwise potentially be large. In this Appendix we discuss modifications in the geometry of a cylindrical shell upon application of internal pressure. We first compute analytically the deformation field in the linear approximation, and then extend the derivation in the nonlinear regime with the assistance of numerical analysis. The 3D elasticity equations (Equation (2) in the main text) in the linear approximation and cylindrical coordinates read where x ¼ r cosðÞ, z ¼ r sinðÞ; non-diagonal terms of s vanish due to the symmetry of the problem. Boundary conditions are where R in ¼ R À t=2 and R out ¼ R þ t=2. The condition of zero longitudinal stress at the two ends of the cylinder is motivated by the results in Appendix 7.
The third line of Equation (14) and the boundary conditions imply s yy ¼ 0 along the shell. By using the constitutive Hookean relations between stress and strain tensor in the main text, we obtain For small deformations, the diagonal components of the strain tensor (see Equation (1) in the main text) are given by Due to the geometry of the problem qu =q ¼ 0; by using Equations (14), (16) and (17) we find the following equation for the radial deformation d 2 u r dr 2 þ 1 r du r dr À u r r 2 ¼ 0 ; whose general solution is u r ¼ a=r þ br. Using the boundary conditions we obtain It follows from the solution Equation (19) that the relative change in radius, thickness and length, in the limit t=R ( 1, are given by Appendix 2-figure 1 demonstrates good agreement of these predictions with numerical simulations for small values of p=E. The corresponding expressions for the change in the volume V ¼ pL R þ t=2 ð Þ 2 of the cylinder and its bulk modulus read As p=E increases, nonlinear behaviors beyond the linear description of Equations (14) and (17) become important. An analytical description of the corresponding deformations would require to take into account nonlinear terms in the original Equations (1) and (2) of the main text.
For our purpose here, the following empirical approach will suffice. The linear solutions Equation (21) depend on two dimensionless small parameters: pR=Et and t=R. In fact, Equation (21) has t=R appearing only multiplied by pR=Et; that makes its contribution small, and implies that the functions depends on pR=Et only, at the dominant order. Appendix 2- figure 1 shows that this property extends in the nonlinear regime: indeed, the curves for relevant values of t=R collapse when plotted vs pR=Et. Using this empirical observation, we looked for a power series in pR=Et and found numerically that the following functional forms describe the deformations in the regime relevant for our problem (pR=Et 0:4): where a R ¼ À0:6182, a t ¼ À0:0626, a L ¼ 0:0479. Finally, by using Equation (23), we obtain the volume V and the bulk modulus k of the cylinder as a function of p=E, as shown in Appendix 2- figure 1G,H.
In the main text, we studied the mechanical properties of C. elegans body as a function of the pressure parameter p, that is the difference between the internal and the external atmospheric pressure P atm , neglecting the latter. This Appendix shows that our results hold also when P atm is taken into account.
To gain insight, we first adapted the calculation in Appendix 2 to the case where atmospheric pressure P atm is considered. The equations are not modified yet the boundary conditions on the internal and external surface of the shell become P in ¼ P atm þ p ; P out ¼ P atm : By following the same procedure detailed in the previous Section B, we obtain that the radial deformation of a pressurized shells is given by Equation (25) suggests that effects of atmospheric pressure on the mechanical response of the shell are negligible as long as P atm t=pR ( 1. To test this hypothesis, we simulated indentation experiments with and without atmospheric pressure for different values of p; results are shown in Appendix 3-figure 1. As expected, for p=E » 10 À2 , both the force indentation relation and the deformation profile are not modified by the atmospheric pressure. It follows that all mechanical properties and neural responses derived in the main text for stiff worms (where the internal pressure is high and p=E » 10 À2 ) are not modified by the inclusion of atmospheric pressure.
Soft worms were shown in the main text to have smaller values of p=E because of the dissection procedure. Appendix 3- figure 1 shows that the force indentation relation does not change significantly yet the deformation profile becomes wider as p=E reduces. This modification is not relevant to describe mechanical properties of intact animals but it might a priori influence the neural response in soft worms. Therefore, we computed numerically the neural response to indentation experiments in soft worms with and without atmospheric pressure. As shown in Appendix 3-figure 1, the performance of the model in describing neural data is not modified substantially, which shows that results of the main text are generally valid even when atmospheric pressure is taken into account.  Figure 3 of the main text. Green curves are experimental data (as in Figure 3 of the main text); purple and black curves are model predictions with and without atmospheric pressure, respectively. The ratio p=E is 7‰ of the value for stiff worms.
By using the numerical approach described in the main text, we determined the mechanical response of a thin cylindrical shell, and compared it to the above solution. Note that we are solving directly the equations of three-dimensional elasticity, at variance with the simplified set of equations in Morley (1960). Results for different mesh sizes are shown in Appendix 5- figure 1: the deformation profile is indeed captured by our code, even with a relatively coarse mesh; a finer mesh is needed to capture quantitatively the force-indentation relation. The agreement improves as the numerical mesh becomes finer (the mesh length is t (blue), t=2 (green), t=3 (red), where t is the thickness of the shell). Parameters of the simulations are:

Large indentation of spherical shells
Large indentations of pressurized spherical shallow shells were investigated in Vella et al. (2012a). In response to a point indentation (and in the absence of buckling [Vaziri and Mahadevan, 2008] Here, r indicates the distance in the plane orthogonal to the indentation point, p is the internal pressure, wðrÞ is the deformation field, and is related to the components of the stress tensor as s ¼ 0 and s rr ¼ =r. The nonlinear term in Equation (32) is due to geometrical nonlinearities generated by large deformations. Boundary conditions are: where the prime denotes differentiation with respect to r. A test of our numerical scheme is provided by the comparison of its results for a thin spherical shell to the solution of the simplified Equation (32). Examples of force-indentation relations and deformation profiles given by Equation (32)  Appendix 5-figure 2. Mechanical response of a pressurized thin spherical shell to a point indentation at the north pole. (A) The force required to produce a deformation w 0 =R ¼ À10 À2 for different values of p. (B) The deformation profile for p=E ¼ 10 À6 (blue), 10 À5 (green), 10 À4 (red). Black lines and colored dots correspond to the solutions of Equation (32) and the results by our code, respectively. Note that, as for a pressurized cylinder (see main text), the deformation profile narrows as p increases. Parameters of the simulations are : E ¼ 1 MPa, n ¼ 0:3, t=R ¼ 10 À3 , t ¼ 1m. Standard touch sensation experiments have the nematode glued onto a plate. As mentioned in the main text, the displacement of its body is strongly limited at locations where the glue is applied. This Appendix will analyze the effect of the gluing onto mechanical responses.
We consider the limiting case where only the line of contact with the plate is glued, that is the limit opposite to the one considered in the main text. There, the entire lower half of the body was glued, which was motivated by the experiments reported in the paper. In our model, the south pole of the cylinder corresponds to the line of contact with the plate in the absence of indentation. Appendix 6- figure 1 shows the corresponding response of the shell to indentation.
It is of interest that the stiffness in Appendix 6- figure 1A is smaller for the south-pole gluing, even though its longitudinal deformation is more extended than for the lower-half gluing, as visible in Appendix 6-figure 1B. That is accounted by the deformation along the orthogonal direction, which expands in the entire lower half of the body (see Appendix 6- figure 1C). That deformation is forbidden for the lower-half gluing, which is the reason for the increased stiffness. Generally, the stiffness is expected to decrease as we reduce the region where the nematode is glued to the plate, which could be verified experimentally. The main text discusses a pressurized cylindrical shell with free conditions at the two ends of the cylinder. Here, we analyze how mechanical properties are modified when the ends are closed by plugs, which leads to an additional component of stress.
We computed numerically the response to the indentation of a closed pressurized cylindrical shell. The numerical procedure is quite analogous to the main text, with the only difference of the boundary conditions. The action of the pressure on the plugs produces a longitudinal force on the shell, whose magnitude does not depend on their structure (but for a boundary layer close to the ends). Without loss of generality, we used semi-spherical plugs.
Results of the simulations are shown in Appendix 7-figure 1. The deformation profile for a given value of p=E is more extended than for free conditions at the sides (Appendix 7- figure 1A). Since the associated change in volume is bigger, and the stiffness is dominated by internal pressure, the stiffness of the shell is greater (Appendix 7- figure 1B). The additional stiffness stems from the longitudinal stress introduced by the lateral plugs, as confirmed by the decrease of the difference between the closed and the open conditions with the internal pressure (Appendix 7- figure 1A,B).
Let us now show that closed ends cannot reproduce experimental data, which is the reason why the main text focuses on open cylinders. We first fix the range p=E 2 ½0; 0:02 considered so far. The extension of the deformation decreases with p=E but, even at the largest value p=E ¼ 0:02, it remains too wide to account for experimental data. Further increase in p=E does reduce the extension, yet it runs in conflict with experimental data on the bulk modulus (Gilpin et al., 2015). Indeed, the estimate for p obtained in the main text depends only on the deformation profile and the force-indentation relation, which are both given by the experiments. We can therefore fix p ¼ 40kPa, and predict the bulk modulus for the corresponding various values of E. Results in Appendix 7-figure 1C are systematically smaller than experiments (and even further increase of p=E would not help as the bulk modulus decreases with p=E).
In summary, we showed that the description of the nematode body as an elastic shell with closed ends cannot reproduce experimental data on the mechanics of C. elegans. Conversely, the main text showed that the same elastic model with free sides does capture main features of the mechanical response. Our results suggest that the longitudinal stress generated by the plugs is somehow relaxed in the worm, which may relate to the annular structure of the cuticle. In our model, the free energy of a channel is modulated by the deformation of its elastic filament. A general rotationally-symmetric form of the free energy reads where x ¼ jxj. In the main text, we used the linear form of Equation (35); here, we repeat the analysis of the experimental data for g 1 ¼ 0, that is a quadratic form. Results are shown in Appendix 9-figure 1, with results for the linear model included for the sake of comparison. The quadratic dependence limits the sensitivity range, and leads to a worse description of the data (see Appendix 9-figure 1); the same also holds for individual realizations of the responses (data not shown). That is witnessed by the step response in Appendix 9-figure 1, where the current predicted by the model saturates at smaller values than the peaks observed in the data, even though the best-fitted baseline activity is stronger. We also repeated the analysis of the directional model replacing Equation (13) of the main text by where Q is the Heaviside function. Results were comparable to those presented in the main text.
Here, we compute mean and variance of the current as a function of the statistics of the ion channels. We use these results in the main text to compare model predictions to experimental data. The point is to generalize standard noise analysis to the case where the channels are not equivalent. That can be the case either because they are not identical or, as it the case here, because their stimulation differs due to their location with respect to the stimulation. We also calculate the scaling of the level of fluctuations expected in a finite sample of N measurements. We start by deriving the expression of the mean and variance of the neural current from the statistics of the single channels. The total ion current I through the neuron is the sum where i k is the current flowing through the k-th channel and the sum runs over the channels along the neuron. Channels can take three distinct states in our model : open, with maximal current i 0 and probability P o ðkÞ; sub-conducting, with intermediate current i s and probability P s ðkÞ; closed, with no current and probability P c ðkÞ ¼ 1 À P o ðkÞ À P s ðkÞ. Each channel follows a generalized Bernoulli distribution : the associated mean and variance are hi k i ¼ i o P o ðkÞ þ i s P s ðkÞ ; s 2 k ¼ i 2 o P o ðkÞ 1 À P o ðkÞ ð Þ þ i 2 s P s ðkÞ 1 À P s ðkÞ ð Þ À 2i o i s P o ðkÞP s ðkÞ: Assuming that the gatings of the channels are independent random variables, we obtain In the main text we suggested that the variance of the current could be used to study experimentally the microscopic properties of ionic channels. Since the variance should be inferred from a finite sample, we want to quantify the scaling of fluctuations in the sample variance with the number of measurements. Given N measurements I n of the TRN current, the sample mean m 1 and sample variance m 2 are defined as The expected sample variance and its variance are (Kenney and Keeping, 1951;Rose and Smith, 2002): To determine the relation between the moments 2 ðIÞ, 4 ðIÞ, . . . of the total current I and the statistics of the single-channel currents, it is convenient to use cumulants and the cumulant-generating function of I, defined as Q I ðtÞ ¼ loghexp tI ð Þi (Papoulis, 1991). Its advantage is that the function is additive : where we have exploited independence among the i k 's. Since the cumulant q n ðIÞ of order n (see Papoulis, 1991) is q n ðIÞ d n dt n Q I ðtÞj t ¼ 0 , Equation (42) implies the additivity q n ðIÞ ¼ P k q n ði k Þ of the cumulants.
The cumulants q n ði k Þ are calculated using the fact that the channels obey a generalized Bernoulli distribution : Q ik ðtÞ ¼ log 1 À P s ðkÞ À P o ðkÞ ð Þ þ P s ðkÞ e tis þ P o ðkÞ e ti0 Â Ã ; whence we obtain q 1 ði k Þ ¼ i s P s ðkÞ þ i o P o ðkÞ ; q 2 ði k Þ ¼ i 2 s P s ðkÞ 1 À P s ðkÞ ð Þ þ i 2 o P o ðkÞ 1 À P o ðkÞ ð Þ À 2i s i o P s ðkÞP o ðkÞ ; q 4 ði k Þ ¼ i 3 s i s À 4q 1 ði k Þ ½ P s ðkÞ þ i 3 o i o À 4q 1 ði k Þ ½ P o ðkÞ þ 6 q 1 ði k Þ ð Þ 2 q 2 ði k Þ þ 3 q 1 ði k Þ ð Þ 4 À3 q 2 ði k Þ ð Þ 2 : The derivation is completed by relating the central moments ðIÞ to its cumulants via standard formulae (Papoulis, 1991), for example 2 ðIÞ ¼ q 2 ðIÞ ; 4 ðIÞ ¼ q 4 ðIÞ þ 3 q 2 ðIÞ ð Þ 2 ; . . . ; by expressing q n ðIÞ as P k q n ði k Þ, and finally using Equation (44 Mutations are modeled via changes of the stretching stiffness S with respect to the wild type S wt . (A) As S increases, the geometry of the pressurized cylinder modifies: its length L increases and its radius R decreases. (B) As S increases, the bulk modulus k (blue) and the stiffness f (red) of the force-indentation relation increase. (C-D) Our theoretical predictions and experimental measurements (Park et al., 2007;Gilpin et al., 2015) for f and k vs R. Since the radius of the mutants was not reported in Gilpin et al. (2015), we used values in Park et al. (2007). The upshot is that lon-2 mutants should have a bulk modulus significantly different than the wild type.

Dependence of neural responses on the properties of elastic filaments
This Appendix will investigate the dependence of the TRN neural response on the interaction between the channels and the surrounding medium, which we have embodied here into an elastic filament. In our model, those interactions are described by two parameters: the elastic constant and its friction coefficient with the medium. The neural response depends on the ratio t ¼ g=k, which appears in Equation (5) of the main text. To understand the effects of t, we computed the TRN current predicted by our model in response to a step of fixed amplitude as a function of t, keeping fixed all other parameters detailed in the main text.
Appendix 13- figure 1 shows the peak current and the decay time, i.e. the time for the current to reach half of its peak value, both averaged over the statistical realizations of the distributions for the channels. Both the peak current and the decay time increase with t: filaments relax more slowly, which provides a stimulus on the associated channel that lasts longer and thereby allows for a higher current. detailed in the main text. Appendix 13- figure 2 shows the histogram of least-squares errors for individual responses to the stimuli in Figure 3 of the main text for different geometries of the filaments attached to the channels. The black curve shows that individual realizations as well (and not just the average as in the main text) reproduce neural responses for the various profiles, strengths and frequencies. The figure also presents results for a model with filaments initially or permanently restricted to be tangential. Graphs indicate that our predictions are further improved by introducing those additional assumptions.