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

  1. Alessandro Sanzeni
  2. Samata Katta
  3. Bryan Petzold
  4. Beth L Pruitt
  5. Miriam B Goodman
  6. Massimo Vergassola  Is a corresponding author
  1. University of California, San Diego, United States
  2. National Institute of Mental Health Intramural Program, National Institutes of Health, United States
  3. Stanford University School of Medicine, United States
  4. Stanford University, United States


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.


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 (Woo et al., 2014; 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-4-dependent 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.


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 50μm 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 μm to indent the body to a maximum depth of ~10 μm, that is about half the radius of the shell and larger than its thickness ~1 μm.

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 L1mm and radius R+t/225μm is indented by a spherical bead (with radius 10 μm 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.

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:

(1) 2εij=iuj+jui+(iuk)(juk),

where j is the derivative with respect to the spatial coordinate xj (j=1,2,3) and ui is the i-th component of the displacement. The equations for the stress tensor σij (Piola-Kirchoff of the second type) read (Landau et al., 1986; Ventsel and Krauthammer, 2001; Audoly and Pomeau, 2011):

(2) j[σij+σkjuixk]=0,

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 σij=E1+ν(εij+ν1-2νεkkδ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, ν is the Poisson ratio, and σij is energy conjugate to εij, that is the elastic energy el upon a deformation varies as δel=σijδεij𝑑V.

In addition to E and ν, 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

The TRN dendrites are decorated by MeT channels along their entire length (Chelur et al., 2002; Emtage et al., 2004; Cueva et al., 2007). We propose that these channels are activated by the 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 𝐫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 𝚫𝐫c,f.

The elastic component reflects the Hookean response of the filament to its stretching. Elastic energy is V(𝐱)=k2x2, where 𝐱=𝚫𝐫𝐟-𝚫𝐫𝐜 is the elongation of the filament with respect to its undeformed configuration. The corresponding restoring force is

(3) 𝐅elastic=-k𝐱.

As for the frictional component, the TRN and its channels are embedded in the medium and expected to move with it, that is 𝚫𝐫c=𝐮(𝐫c) where 𝒖 is the displacement in Equation (2). Conversely, as the filament slides with respect to the medium, the friction force is

(4) 𝐅friction=-γd(𝚫𝐫f-𝐮(𝐫))dt,

where γ is the friction coefficient and 𝐫 is the undeformed position of the material point that coincides with the location of the tip, that is 𝐫f+𝚫𝐫f=𝐫+𝐮(𝐫). Expanding 𝐮(𝐫) and using that gradients of 𝐮 are small, we obtain 𝐫-𝐫f𝚫𝐫f-𝒖(𝐫f), which is then inserted into Equation (4) to show that 𝐮(𝐫) can be replaced by 𝐮(𝐫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 𝐅friction+𝐅elastic=0, which yields

(5) d𝐱dt+1τ𝐱=d𝚪dtd(𝐮(𝐫f)-𝐮(𝐫c))dt,

where τ=γ/k is the relaxation time, 𝐱 is the extension of the filament, and 𝚫𝐫c=𝐮(𝐫c) was used. Equation (5) drives 𝐱 to zero for 𝚪 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 𝐱𝐰^30, where, for every channel, 𝐰^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 𝐞^i as follows: 𝐞^y is aligned with the local direction of the (deformed) axis of the cylinder running head-to-tail; 𝐞^z is orthogonal to the neural membrane at the top of the TRN, and oriented outward; 𝐞^x is tangential to the neural membrane, along the remaining direction of a right-handed system. The bases 𝐰^i are constructed by rotating the 𝐞^i appropriately. For a channel placed at the top of the TRN, the local basis 𝐰^i coincides with 𝐞^i. If the channel is rotated by θ along the surface of the TRN, then 𝐰^1=cos(θ)𝐞^x-sin(θ)𝐞^z, 𝐰^2=𝐞^y, and 𝐰^3=sin(θ)𝐞^x+cos(θ)𝐞^z.

The dynamics of 𝚪 in Equation (5), is obtained by the displacements 𝒖 calculated using the mechanical model. The relation between stretching and ε is (Landau et al., 1986; Audoly and Pomeau, 2011): Γ2(t)-Γ2(0)=2εij(t)Γi(0)Γj(0), where the Γ’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 (CSO). The respective probabilities Pc,s,o obey the master equation

(6) {dPcdt=-RcsPc+RscPs,dPsdt=-(Rsc+Rso)Ps+RosPo+RcsPc,dPodt=-RosPo+RsoPs,

where Rij are the respective transition rates, and Rc,o=Ro,c=0 again to minimize free parameters. The channels are posited to work at equilibrium, so that

(7) Rcs/Rsc=e-βΔGsc;Rso/Ros=e-βΔGos,

where ΔGij is the free energy difference between the states i and j, β=1/kBT, T is the temperature, and kB 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:

(8) βΔGoc=g0-g1,

where g0, g1 are dimensional constants, and is the amplitude of the tangential component of 𝐅elastic in Equation (3):

(9) 1=𝐅elastic𝐰^1;2=𝐅elastic𝐰^2,

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 𝒙, 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)

(10) ΔGos=aΔGoc;ΔGsc=(1-a)ΔGoc,

with the only additional parameter 0a1. 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=kik of the currents of individual channels. Its mean is given by

(11) I=iokPo(k)+iskPs(k),

and its variance is calculated similarly (see Appendix 10). Here, io and is are the channel current in its open/subconducting state. P(k) obey Equation (6), and their rates depends on the position along the TRN via Equation (3), (5), (8). The single-channel current i0=-1.6±0.2pA, as measured in O'Hagan et al. (2005) and Brown et al. (2008); other parameters will be inferred from experimental data.

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 formulæ 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 non-dimensional 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.

Deformation profiles and force indentation relations.

(A) Representative photomicrograph of a transgenic animal with GFP-tagged cuticular annuli being pressed into a glass bead. Experimental deformation profiles in (B) were derived from a stack of images at different focal planes. (B) Experimental and numerical deformation profiles along the longitudinal axis (the generatrix of the cylinder). Data were obtained as described in Materials and methods by using 2 biological replicates (adult animals). (C) Experimental (from Eastwood et al., 2015; Eastwood et al., 2019) and numerical force-indentation relationships. Length is L=1mm and the Poisson coefficient ν=0.3.
Figure 2—source data 1

Measurements of cuticle deformation by beads.

Figure 2C shows how the estimate p40kPa 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 E4MPa, 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 formulæ 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 ΔV of the initial volume V0 was found to depend linearly on the variation of the external pressure Δp, and the resulting bulk modulus κ=ΔpΔVV0=140±20kPa. Performing the same operations in our simulations, we obtained estimated values of κ=150-230kPa. This agreement is quite significant as we derived κ, 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.

Our model captures experimental neural responses to various stimuli.

(A) The applied experimental indentation (top); TRN’s response (bottom, green) and average predictions (solid black). Dot-dashed black lines correspond to one standard deviation above/below the mean. Experimental stimuli and neural responses are from Eastwood et al. (2015) and Eastwood et al. (2019). (B) A typical ramp-like profile of indentation (top) and the corresponding current (TRN’s response in green; black lines as in panel A). (C) The predicted peak current vs the slope of the ramp for a total fixed indentation of 8 μm. Red circles indicate experimental data from Eastwood et al. (2015) and Eastwood et al. (2019).
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 pre-indentation for the purple points is 5 μm). 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 off-steps 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: Δz 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.
Figure 4—source data 1

Experimental and simulated neural responses to pre-indented steps.

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 𝚪^(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 (𝚪(0)𝐰^3=0 or 𝐱(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).

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 μm 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 pre-indentation, 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).

Residual internal pressure accounts for current amplitude in soft and stiff worms.

(A) Experimental data (dots) and the average theoretical prediction (lines) for force-indentation relations. Best fit of pooled data for soft worms gives p = 1.6 kPa; individual values are variable, with estimated p in the range 0.04–16kPa. (B) The vertical deformation profiles uz vs the position y along the longitudinal axis for stiff (red) and soft (blue) animals. Note the widely differing profiles for the force-clamped curves. (C) Experimental (dots) and theoretical (mean value as continuous lines) peak current for force (top) and displacement-clamped (bottom) stimuli. The current is normalized by the mean peak in soft and stiff worms, respectively. (D) Peak current vs the pressure p, which shows that the model (continuous lines are the mean; dot-dashed lines are above/below one standard deviation) captures experimental trends (dots). Experimental data reproduced from Eastwood et al. (2015) and Eastwood et al. (2019) and derived from 4 and 21 recordings in the stiff (red) and soft (blue) conditions, respectively.
Figure 5—source data 1

Experimental and simulated neural responses to force-clamped or displacement-clamped stimuli.

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 yh for the deformation to reduce to half of its maximum value (at the center of the bead). Figure 6 shows that yh decreases, that is the deformation is more localized, when p/E increases. Conversely, as t reduces, the deformation is wider if p/E10-4 and narrower if p/E10-4.

The mechanical balance for our model of pressurized shell.

(A) The longitudinal extension yh vs p/E for various thicknesses t. Different trends at small and large values of p/E reflect the contributions of bending to the elastic energy. (B) yh vs p/Et. The collapse of the curves at the right end reflects the small value of the bending term coefficient (see the text). The value of yh found experimentally (black line) is well inside that asymptotic region. (C) The various contributions to energy, and the work done by the indenter for increasing values of the ratio p/E.

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/τ2E2t4/p2R4. If τ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 SEt. These arguments suggest to plot yh 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 Rb 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: Rb 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 Rb. 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 Rb. The radius Rb also controls the deformation outside of the contact region, namely the mid-maximum extension of the deformation (Rb, data not shown). As p/E reduces, the deformation becomes shallower at the indentation point and the role of Rb vanishes (see Figure 2B).

Effects of the radius Rb of the indenting bead.

(A) In pressurized shells, the deformation profile depends on Rb (colored lines) while the dependence disappears in the absence of internal pressure (black line). (B) The force-indentation relation for p=40kPa and various Rb. In strongly pressurized shells, the relation (colored dots) follows Equation (12) (solid lines). For shells with p=0, the curves collapse onto a unique curve (black dots). (C) The ratio F(w0)/Fmax of the force normalized by its maximum value is essentially independent of Rb. (D) The peak current increases with Rb in a p-dependent manner (the same holds for the sensitivity to the stimulus). The current is normalized by its value for Rb=5μm.

Similarly, Figure 7A shows that the volume of the body to be deformed increases with Rb, and more work is needed for a given maximum indentation w0. In formulæ: the work Fdw0 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-pdV/dw0.

Qualitatively, larger dV’s associated with larger beads yield the nonlinear dependence in Figure 7B. Quantitatively, we write F=pχ, where χ=-dV/dw0 has the dimension of a length squared and depends on w0, the shell and bead radii R, Rb, and the ratio p/E. Keeping the latter fixed, we investigated numerically the behavior of χ and observed that the ratio F(w0)/Fmax does not depend on Rb (see Figure 7C). It follows that the dependence on Rb should factorize out: χ=𝒢1(Rb/R)𝒢2(w0,R) where 𝒢1 depends on Rb/R for dimensional reasons. The function 𝒢2 brings the length squared dimensionality and, in the limit of small Rb and large p, behaves as Rw0 (Vella et al., 2012b). It follows that 𝒢2=Rw0𝒢3(w0/R).

The above functions 𝒢3 and 𝒢1 are determined as follows. We computed numerically the force indentation relation of cylinders of different radii (R = 25, 40, and 50 μm) to stimulations produced by beads of different size (Rb from 3 to 10 μm); results of the simulations are then used to fit coefficients of the Taylor expansions of 𝒢1(x) and 𝒢3(x). Using this approach we find that the functional form

(12) F=α1pRw0(1+α2RbR)(1+α3w0R),

with α1 = 0.76, α2 = 2.1, and α3 = 0.66, captures quantitatively the behavior of the force indentation relation (Figure 7B, R2 = 0.995). Variations between Rb = 3 μm and Rb = 10 μm are on the order of few μN, hence they should be accessible experimentally. Equation (12) generalizes the linear relationship, valid in the limit of very small Rb, obtained in Vella et al. (2012b).

Consequences for neural responses are in Figure 7D. In agreement with O'Hagan et al. (2005), the peak current increases by ~20% as Rb goes from 3 to 10 μm, hence our prediction could be tested experimentally. A quantitative comparison with data in O'Hagan et al. (2005) was hampered by the lack of force-indentation measurements in O'Hagan et al. (2005), preventing us from inferring p.

Finally, it is worth remarking that the bead size also affects the dependence of the response on the circumferential position of the TRN. Indeed, Appendix 6—figure 1 evidences that the profile of 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 σ, 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.

Stimuli on the channels due to a step.

(A) Indentation profile for a step stimulus. (B) The diagonal components of the strain tensor εijvs the longitudinal position along the cylinder. The overlap of those components (color) and the eigenvalues of εij (black) show that the tensor is essentially diagonal, which leads to the conservation of angles under deformation discussed in the main text. (C) The two components (green and blue curves) tangential to the neural membrane of the force acting upon on a channel (computed using Equations (5), (9)) for the stimulus in panel A. The panels refer to the different directions of the elastic filament (the first two are tangential and the third orthogonal to the neural membrane). (D) Gating probability for an individual symmetric or directional channel, as produced by the two tangential extensions in panel C. Parameters are: y=1μm, θ=0, 𝑣=cos(π/3)𝐰^1-sin(π/3)𝐰^2. 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.

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 𝐯:

(13) βΔGoc=g0-g2𝐯.

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 𝐯 in Equation (13) independently and isotropically distributed (Figure 9A).

Symmetry of the single channel response.

(A) Mean neural current response to a step, for symmetric (black) vs directional (red) channels. (B) Histogram of the errors in reproducing the data of Figure 3 obtained with the two above models for different realizations of the channels’ distribution. Symmetric channels (black) give a better description. (C) The mean error as the maximum current i0 per channel is varied. (D) The Coefficient of Variation (CV) of the TRN current vs the stimulus strength, calculated over many repetitions of a given stimulus. (E) The average current for a channel (normalized by its maximum value i0) as a function of its distance to the center of the indenting bead. (F) The current flowing along the TRN vs the position of the indenting bead. The origin indicates the end of the TRN; negative coordinates correspond to the relatively insensitive zone in the middle of the body of the worm. (G) The colored curves show the predicted current for symmetric and directional channels, for a given distribution of the channels. The black curves show the expected level of asymmetry between onset and offset, as quantified by the standard deviation (Ion-Ioff)21/2 between the peak responses Ion,off at the onset/offset of the stimulus, averaged over the distributions of the channels. Dashed-dotted curves show the range of expected asymmetries in individual realizations.

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 i0=-1.6±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 i0 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 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.


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 dead-zone 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.

Experimental methods

Nematode strains

Request a detailed protocol

The following transgenic C. elegans nematodes were used: TU2769 uIs31[mec-17p::gfp] III (O'Hagan et al., 2005) and TP12 kaIs12[col-19::gfp] (Thein et al., 2003). The corresponding identifiers are RRID:WB-STRAIN:TU2769 and RRID:WB-STRAIN:TP12, respectively. The uIs31 transgene expresses GFP exclusively in the TRNs, enabling in vivo recordings from these neurons and the kaIs12 transgene encodes a fusion between the COL-19 collagen protein and GFP, labeling cuticular annuli. Animals were grown on OP50 at either 15°C (TU2769) or 20°C (TP12) and used as well-fed L4 larvae or young adults.

Imaging Cuticle Deformation

Request a detailed protocol

TP12 worms were immobilized with 0.1 µm polystyrene beads on a 6% NGM agarose pad. 10 µm 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 µm 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.


Request a detailed protocol

Worms were immobilized on 2% agarose pads with WormGlu (GluStitch), dissected, and patch-clamped as described in Eastwood et al. (2015). Recordings were performed on the ALMR neuron due to geometric constraints of the stimulator system; ALMR is bilaterally symmetric to the previously used ALML neuron. The extracellular solution contained (in mM): NaCl (145), KCl (5), MgCl2 (5), CaCl2 (1), and Na-HEPES (10), adjusted to pH 7.2 with NaOH. Before use, 20 mM D-glucose was added, bringing the osmolarity to ~325mOsm. The intracellular solution contained (in mM): K-Gluconate (125), KCl (18), NaCl (4), MgCl2 (1), CaCl2 (0.6), K-HEPES (10), and K2EGTA (10), adjusted to pH to 7.2 with KOH. Before use, 1 mM sulforhodamine 101 (Invitrogen) was added to help visualize successful recording of the neuron.

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

Request a detailed protocol

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 μm 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 μm 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 μm, and 20–23 μm 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 μm anterior to the ALM cell body.

Stimulus control and data acquisition

Request a detailed protocol

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/μm. 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, high-current 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 μm, and the largest displacement used was 3–4 μm 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.

Data analysis

Request a detailed protocol

Whole-cell capacitance and series resistance were measured as previously described (Goodman et al., 1998). Data analysis was performed with MATLAB from Mathworks (data import and analysis functions are available online at: and Igor Pro (Wavemetrics). The identifier of the MATLAB-Patchmaster analysis code is (copy archived at

Only recordings with holding current <-10pA at -60mV and series resistance <210MΩ 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.

Appendix 1

Numerical evaluation of the amplitude of the deformation gradients

In the main text, we exploited the property that gradients of the deformation field 𝐮 (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.

Appendix 1—figure 1 shows the various gradients components for an indentation of 10 μm, 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 duz/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 duydz has a negative value that compensates for duzdy in the yz component of the strain tensor, which could otherwise potentially be large.

Appendix 1—figure 1
The gradients of the deformation along the longitudinal direction of an indented cylinder.

Gradients of 𝐮 computed by using numerical simulations for soft (p/E=4×10-4, blue) and stiff (p/E=0.01, red) shells. Note their moderate amplitude even for the indentation of 10 μm shown here, which is the strongest that we consider.

Appendix 2

How does the radius of the cylinder change with the internal pressure?

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

(14) σrrr+σrr-σθθr=0,1rσθθθ=0,σyyy=0,

where x=rcos(θ), z=rsin(θ); non-diagonal terms of σ vanish due to the symmetry of the problem. Boundary conditions are

(15) σrr(r=Rin)=-p,σrr(r=Rout)=0,σyy(y=±L/2)=0,

where Rin=R-t/2 and Rout=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 σyy=0 along the shell. By using the constitutive Hookean relations between stress and strain tensor in the main text, we obtain

(16) σrr=E1-ν2(ϵrr+νϵθθ),σθθ=E1-ν2(ϵθθ+νϵrr).

For small deformations, the diagonal components of the strain tensor (see Equation (1) in the main text) are given by

(17) ϵrr=urr,ϵθθ=1ruθθ+urr,ϵyy=uyy.

Due to the geometry of the problem uθ/θ=0; by using Equations (14), (16) and (17) we find the following equation for the radial deformation

(18) d2urdr2+1rdurdr-urr2=0,

whose general solution is ur=a/r+br. Using the boundary conditions we obtain

(19) ur=p/ERout2/Rin2-1((1+ν)Rout2r+(1-ν)r),uy=-νp/ERout2/Rin2-1y,uθ=0,
(20) σrr(r)=pRout2/Rin2-1(1-Rout2r2),σθθ(r)=pRout2/Rin2-1(1+Rout2r2).

It follows from the solution Equation (19) that the relative change in radius, thickness and length, in the limit t/R1, are given by

(21) ΔRR=pREt(1+(ν-1)t2R),Δtt=-pREt(ν+(1-ν)t2R),ΔLL=-ν2pREt.

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=πL(R+t/2)2 of the cylinder and its bulk modulus read

(22) ΔVV(2-ν/2)pREt,klinearEtR(2-ν/2).

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/Et0.4):

(23) ΔRR=pREt(1+αRpREt),Δtt=pREt(-ν+αtpREt),ΔLL=pREt(-ν2+αLpREt),

where αR=-0.6182, αt=-0.0626, αL=0.0479.

Appendix 2—figure 1
Deformations of a cylindrical shell due to internal pressure.

Relative change in the radius R (A), thickness t (B) and length L (C) of the shell as a function of p/E and various t/R, as obtained from our numerical simulations. As p increases, R increases whilst t and L decrease. (D,E,F) While the curves in the previous panels change with t/R, the curves are collapsed by plotting them against pR/Et, which suggests that the contribution of terms in t/R is negligible. The collapsed behavior agrees with the linear prediction Equation (21) for pR/Et<0.1 (dash-dotted lines), and is well captured by the empirical Equation (23) in the moderately nonlinear regime. Using the nonlinear Equation (23), we computed the change in volume (G) and bulk modulus (H) of the shell as a function of p/E (dashed lines are the linear predictions Equation (22) valid for small values of p/E).

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.

Appendix 3

Effects of external atmospheric pressure on mechanical and neural response

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 Patm, neglecting the latter. This Appendix shows that our results hold also when Patm is taken into account.

To gain insight, we first adapted the calculation in Appendix 2 to the case where atmospheric pressure Patm is considered. The equations are not modified yet the boundary conditions on the internal and external surface of the shell become

(24) Pin=Patm+p,Pout=Patm.

By following the same procedure detailed in the previous Section B, we obtain that the radial deformation of a pressurized shells is given by

(25) ΔRRpREt[1+(ν-1)t2R(1+2Patmp)].

Equation (25) suggests that effects of atmospheric pressure on the mechanical response of the shell are negligible as long as Patmt/pR1. 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/E10-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/E10-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.

Appendix 3—figure 1
Effect of atmospheric pressure on mechanical response.

Deformation profile (first row) and force indentation relation (second row) for a cylindrical shell (with the same properties as in the main text) with (purple) and without (black) atmospheric pressure.
Appendix 3—figure 2
Effect of atmospheric pressure on neural response.

Neural responses of the model presented in the main text to the experimental stimuli in 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.

Appendix 4

The thin shell limit

In the limit when the shell is thin and the surface is shallow, the 3D elasticity Equation (2) reduce to (Ventsel and Krauthammer, 2001; Audoly and Pomeau, 2011):

(26) {B4uz+k2ϕ-[ϕ,uz]=p-F(x,y),1S4ϕ-k2uz=-12[uz,uz],

which is the limit equation that was used in the main text for analyzing the energetic balance between bending, stretching and internal pressure. The brackets in Equation (26) are defined as [f,g]2fx22gy2-22fxy2gxy+2gx22fy2, while the derivative k2=1R2y2. Equation (26) is two-dimensional, with the variables z=z(x,y) and uz=uz(x,y) representing the middle surface and the deformation field of the cylindrical shell, respectively. The deformed surface is z+uz and we chose the axes so that the plane z=0 is tangent to the top of the cylinder. The Airy stress function ϕ is the scalar function that parametrizes the in-plane components of the stress tensor as σxx=yy2ϕ, σyy=xx2ϕ and σxy=-xy2ϕ.

The above simplifications are due to the thinness of the shell and the resulting small vertical components of σ. The parameters B=Et3/12(1-ν2) and S=Et are the bending and stretching stiffness, respectively. Finally, p and F are the internal pressure and the external force applied by the indenter. In the limit R, the k2 term is negligible, and Equation (26) reduce to the Föppl-von Karman equations for a thin plate (Ventsel and Krauthammer, 2001; Audoly and Pomeau, 2011).

Appendix 5

Validation of the numerical scheme used to determine the mechanical response

In the main text, we described the mechanics of the nematode as an elastic cylindrical shell under pressure. Because of the geometrical nonlinearities involved, numerical simulations are the main tool available for determining the resulting mechanical response. The goal of this Appendix is to give more details on the tests that we employed to validate the numerical procedure discussed in the main text. Tests rely on elasticity problems that were previously investigated, namely small indentations of cylindrical shells (where an analytical solution is available [Morley, 1960]), and large indentations of pressurized spherical shells (where a simplified framework was derived in Vella et al., 2012a). In all cases, agreement was obtained.

1 Small indentations of cylindrical shells

A thin cylindrical shell subject to equal and opposite concentrated radial loads was investigated in Morley (1960). For indentations w0t, where t is the shell thickness, the equations of 3D elasticity (Equation (2) in the main text) reduce to (Morley, 1959):

(27) 4(2+1R2)2w+4K4R4w,xxxx=1D4q,

where q is the applied load, x and θ are the cylindrical longitudinal and angular coordinates, w is the radial displacement and

(28) K4=3(1ν2)(Rt)2,D=Et312(1ν2), w,xxxθ=1R4wx3θ,w,xxxx=4wx4,2=2x2+1R22θ2.

By ‘equal and opposite concentrated radial loads’, it is meant that two equal, spatially localized forces are applied at (x,θ)=(0,0) and (x,θ)=(0,π).

The solution to the above problem was obtained (Morley, 1960) by writing q(x,θ) as

(29) q(x,θ)=Fπ2R2limδ00[sin(λδ/R)λδ/R+2n=2,4,sin(nδ/R)sin(λδ/R)(nδ/R)(λδ/R)cos(nθ)]cos(λx/R)𝑑λ,

and solving Equation (27) order by order in n. The resulting deformation profile reads (Morley, 1960) :

(30) w(x,θ)=R2F2πD{exK4K4[RK+cos(xK+)+RKsin(xK+)]2n=2,4,Im[λ1n(λ1n2+n2)eiλ1nx/R(λ1n2+n21)(λ1n4n4+n2)+λ2n(λ2n2+n2)eiλ2nx/R(λ2n2+n21)(λ2n4n4+n2)]cos(nθ)},


(31) K±=1RK2±12,λ1n,2n2=-n2+12iK2[11i(2n2-1)K2].

The force-indentation relation is obtained from Equation (30) by determining the deformation w at either one of the loading points, that is w(0,0) in the above formulæ, as a function of F. An example of the force-indentation relation and the deformation profile predicted by Equation (30) are shown in Appendix 5—figure 1.

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.

Appendix 5—figure 1
Mechanical response of a thin cylindrical shell to equal and opposite concentrated radial loads.

(A) Force-indentation relation. (B–C) Radial deflection along the longitudinal and angular directions. The analytical solution Equation (30) (black line) is well approximated by numerical solutions (colored lines) obtained by using our numerical code. 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: t/R=10-2, t=1μm, E=1MPa, ν=0.3.

2 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]), the 3D equations of nonlinear elasticity reduce to

(32) {F2π=pr22+ψ(dwdr-rR),rddr[1rddr(rψ)]=Et[rRdwdr-12(dwdr)2].

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 σθθ=ψ and σrr=ψ/r. The nonlinear term in Equation (32) is due to geometrical nonlinearities generated by large deformations. Boundary conditions are:

(33) w(0)=-w0,limr0(rψ-νψ)=0,w()=0,ψ()=pR2,

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) for different values of p are shown in Appendix 5—figure 2. Results of our code are also compared in Appendix 5—figure 2: both the force-indentation relation and the deformation profile are well reproduced.

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 w0/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=1MPa, ν=0.3, t/R=10-3, t=1μm.

Appendix 6

How the gluing of the nematode onto the plate influences its mechanical response

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.

Appendix 6—figure 1
Influence of the gluing of the worm on its mechanical response.

Blue and red curves correspond respectively to gluing of the entire lower half of the cylinder or the line of contact with the plate (south pole) only. (A) Force-indentation relations; note that the stiffness is greater for the blue curve. The deformation profiles along the longitudinal (B) and orthogonal (C) directions are wider for the south-pole gluing, which has the lower half of the shell deformed in the orthogonal direction as well. The undeformed geometry is the black dash-dotted line.

Appendix 7

Influence of the boundary conditions at the ends of the cylinder on its mechanical response

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[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.

Appendix 7—figure 1
Boundary conditions at the ends of the shell influence mechanical properties.

(A) Comparison of the mechanical response for a closed (dash-dot line) and open (continuous line) cylinder. For a given value of p/E, the deformation profile is more extended (left) and the shell stiffer (right) if the two ends are closed; the difference increases with p/E. (B) Experimental and numerical deformation profiles along the longitudinal axis. None of the values p/E[0,0.02] for the closed cylinder (colored line) captures experimental data. (C) Experimental and predicted values for the bulk modulus. The value predicted for closed conditions at the ends decreases with p/E and is too small to account for the data. Results for free lateral conditions are shown for comparison. Parameters of the simulations are as in Figure 2 of the main text.

Appendix 8

The local geometry of the channels along the TRNs

We treat the TRN as a (small) cylinder running (at rest) in the upper part of the (big) cylinder in Figure 1, at x=0 along the y axis. The TRN is placed right above the internal part of the shell. Upon indentation, the orthogonal basis 𝐱^i (i=1,2,3) in Figure 1 is deformed into a triad of vectors 𝐱^i in a way that depends on its original location. The separation between pairs of neighboring material points (𝐫,𝐫+𝐝𝐫) is calculated using the gradients of the displacement field 𝐮(𝐫). In particular, the variation of the squared distance dr2-dr2=2εijdridrj, and the angle between two vectors 𝐝𝐫𝟏 and 𝐝𝐫𝟐 (see Landau et al., 1986; Audoly and Pomeau, 2011)

(34) 𝐝𝐫1𝐝𝐫2-𝐝𝐫1𝐝𝐫2=2εijdr1,idr2,j.

A convenient orthonormal basis 𝐞^i to analyze the dynamics of the channels is defined as follows : 𝐞^y is aligned with the local direction of the (deformed) axis of the cylinder running head-to-tail; 𝐞^z is orthogonal to the neural membrane at the top of the TRN, and oriented outward; 𝐞^x is tangential to the neural membrane, along the remaining direction of a right-handed system.

For every channel, we define a local base 𝐰^i such that 𝐰^1,2 span the plane locally tangential to the neural membrane while 𝐰^3 indicates the orthogonal direction. The bases 𝐰^i are constructed by rotating the 𝐞^i appropriately. For a channel placed at the top of the TRN, the local basis 𝐰^i coincides with 𝐞^i. If the channel is rotated by θ along the surface of the TRN, then 𝐰^1=cos(θ)𝐞^x-sin(θ)𝐞^z, 𝐰^2=𝐞^y, and 𝐰^3=sin(θ)𝐞^x+cos(θ)𝐞^z.

Appendix 9

Comparison of different functional forms for the activation of the channels

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

(35) βΔGoc(𝐱)=g0+g1x+g2x2

where x=|𝐱|. In the main text, we used the linear form of Equation (35); here, we repeat the analysis of the experimental data for g1=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.

Appendix 9—figure 1
The linear form of Equation (35) outperforms its quadratic counterpart in describing experimental data.

(A) The average predictions for the experimental data (green) given by the linear (black) and quadratic (purple) forms of Equation (35).

We also repeated the analysis of the directional model replacing Equation (13) of the main text by

(36) βΔGoc=g0g2vΘ(v),

where Θ is the Heaviside function. Results were comparable to those presented in the main text.

Appendix 10

Noise analysis for non-equally stimulated channels

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 I=k=1Kik, where ik 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 i0 and probability Po(k); sub-conducting, with intermediate current is and probability Ps(k); closed, with no current and probability Pc(k)=1-Po(k)-Ps(k). Each channel follows a generalized Bernoulli distribution : the associated mean and variance are

(37) ik=ioPo(k)+isPs(k),σk2=io2Po(k)(1Po(k))+is2Ps(k)(1Ps(k))2ioisPo(k)Ps(k).

Assuming that the gatings of the channels are independent random variables, we obtain

(38) I=k=1Kik,σI2=(II)2=k=1Kσk2.

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 In of the TRN current, the sample mean m1 and sample variance m2 are defined as

(39) m1=1Nn=1NIn,m2=1Nn=1N(In-m1)2.

The expected sample variance and its variance are (Kenney and Keeping, 1951; Rose and Smith, 2002):

(40) m2=N1Nμ2(I),var(m2)=(N1)2N3μ4(I)(N1)(N3)N3μ2(I)2,
(41) μ2(I)=(II2=σI2,μ4(I)=(II4.

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 QI(t)=logexp(tI) (Papoulis, 1991). Its advantage is that the function is additive :

(42) QI(t)=logetI=logetkik=klogetik=kQik(t),

where we have exploited independence among the ik’s. Since the cumulant qn(I) of order n (see Papoulis, 1991) is qn(I)dndtnQI(t)|t=0, Equation (42) implies the additivity qn(I)=kqn(ik) of the cumulants.

The cumulants qn(ik) are calculated using the fact that the channels obey a generalized Bernoulli distribution :

(43) Qik(t)=log[(1Ps(k)Po(k))+Ps(k)etis+Po(k)eti0],

whence we obtain

(44) q1(ik)=isPs(k)+ioPo(k),q2(ik)=is2Ps(k)(1Ps(k))+io2Po(k)(1Po(k))2isioPs(k)Po(k),q4(ik)=is3[is4q1(ik)]Ps(k)+io3[io4q1(ik)]Po(k)+6(q1(ik))2q2(ik)+3(q1(ik))43(q2(ik))2.

The derivation is completed by relating the central moments μ(I) to its cumulants via standard formulæ (Papoulis, 1991), for example

(45) μ2(I)=q2(I),μ4(I)=q4(I)+3(q2(I))2,,

by expressing qn(I) as kqn(ik), and finally using Equation (44).

Appendix 11

Mechanical effects of mutations in proteins of the cuticle

Mutations of proteins composing the cuticle have been used as a probe to investigate mechanical properties of C. elegans (Park et al., 2007; Gilpin et al., 2015). This Appendix will explore the consequences of the simplest possible assumptions on the effects of those mutations, obtain predictions for our mechanical model, and compare them to experiments.

We shall describe the effects of mutations in the cuticle through variations in the stretching stiffness S, which is the only parameter related to the thin external layers (see Appendix 4). We assume that mutations do not affect the internal pressure of the nematode. The argument in Park et al. (2007) is that the genes involved, for example dpy-5 and lon-2, do not affect transport proteins likely to regulate osmotic pressure. It is possible, though, that mutations in the cuticle affect the development of the nematode, hence its body structure and the internal pressure. To constructively advance the issue, we assume that that does not happen and explore consequences.

Appendix 11—figure 1A shows the dependence of geometrical properties on S : the radius R of the pressurized cylinder decreases, whilst its length L increases with S. This reflects changes with respect to the unpressurized condition brought by the internal pressure p. Appendix 11—figure 1B shows that the bulk modulus increases with S. Finally, the stiffness S enters the response to indentation experiments. Appendix 11—figure 1B shows the ratio fF/w0, that is the force F needed to reach an indentation w0=5μm, vs S (we verified that results do not depend on the choice of w0). The increase with S is due to two contributions that affect the change in volume via the deformation field. First, if the external radius of the shell is kept constant, the deformation field is wider (see Appendix 11—figure 1B). Second, the radius of the shell becomes smaller as S increases (see Appendix 11—figure 1A) which leads to a larger deformation field (data not shown).

The parameter S is not measured experimentally, which forces us to use R as a proxy. In this formulation, the model predicts that mutations in the cuticle which increase R will decrease L, stiffness and bulk modulus. These features are in qualitative agreement with experimental observations (Park et al., 2007; Gilpin et al., 2015).

Quantitatively, we observe that : (i) the length L of the mutants in Park et al. (2007) is systematically larger than our predictions; (ii) the radius of lon-2 mutants in Park et al. (2007), which is about 25% smaller than the wild type, cannot be obtained in our simulations, as seen in Appendix 11—figure 1C. (i) may be due to the likely non-isotropy of the Young’s modulus generated by the annular structure of the external layers. That is expected to make the stiffness smaller in the longitudinal than the orthogonal direction. (ii) is due to the fact that R cannot be smaller than the value for the unpressurized shell, which is only 20% smaller than the wild type. Since relative changes in R grow with p/E, it is likely that uncertainties on R are due to fluctuations in p/E among worms (see main text). As already mentioned above, both R and L could also be affected by our neglecting developmental effects.

We finally compare in Appendix 11—figure 1C-D our predictions to experimental data from Park et al. (2007) and Gilpin et al. (2015). The agreement is notable and leads to the prediction that the bulk modulus in lon-2 mutants should significantly deviate from the wild type. This differs from the conclusion that the bulk modulus is not affected by the cuticle (Gilpin et al., 2015), which was based on mutants other than lon-2.

Appendix 11—figure 1
Changes in the mechanics caused by mutations of proteins in the cuticle.

Mutations are modeled via changes of the stretching stiffness S with respect to the wild type Swt. (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.

Appendix 12

Inference of the model parameters

The amplitude of 𝚪(0), which controls the scale of the elongation in Equation (5), is set to unity by redefining the parameters gh and gs in Equation (8). As for the rates R, a parsimonious form that respects Equations (7) and (10) is

(46) Rsc=rcse(1+b)(1a)βΔGoc;Rcs=rcseb(1a)βΔGoc;
(47) Ros=rsoe(1+d)βΔGoc;Rso=rsoedaβΔGoc;

Here, rcs (rso) controls the rate of the transitions between the closed and the subconductance states (the subconductance and the open states) and the parameters b,d[-1,0] control their global shift with respect to variations of the free-energy difference. Parameters are inferred from experimental curves by using the MATLAB optimization function ‘lsqcurvefit’, based on the least-squares distance between the predicted and the observed current profiles. For every realization of the quenched disorders (distribution of the channels and the initial direction of their filaments), we obtain the best parameters, which are then averaged. Their means are τ=1.4msgh=1.4103; gs/gh=0.09; rcs=1/69.5msb=0.75 ; rso=1/18ms; d=0.56; a=0.50; is/i0=0.71.

Appendix 13

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 τ=γ/k, which appears in Equation (5) of the main text. To understand the effects of τ, we computed the TRN current predicted by our model in response to a step of fixed amplitude as a function of τ, 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 τ: 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.

Appendix 13—figure 1
Dependence of the neural response on the filament-medium interaction.

Peak current (A) and decay time (B) as a function of the relaxation time τ of the elastic filament connected to the channel. The model predicts that both the peak current and the decay time increase with τ.
Appendix 13—figure 2
Dependence of the neural response on the geometry of the filaments.

The error for the profiles of stimulation in Figure 3 of the main text. The black histogram is built from individual realizations for the unconstrained model discussed in the main text. The purple and the blue curves refer to the corresponding histograms for filaments initially or permanently restricted to be tangential to the neural membrane.

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files. The electrophysiology source data from Eastwood et al. is available on Dryad ( See also transparent reporting form.

The following previously published data sets were used
    1. Eastwood AL
    2. Sanzeni A
    3. Petzold BC
    4. Park S
    5. Vergassola M
    6. Pruitt BL
    7. Goodman MB
    (2019) Dryad Digital Repository
    Data from: Tissue mechanics govern the rapidly adapting and symmetrical response to touch.


    1. Altun ZF
    2. Hall DH
    Introduction, WormAtlas, CSHL Press, 10.3908/wormatlas.1.1.
  1. Book
    1. Audoly B
    2. Pomeau Y
    Elasticity and Geometry: From Hair Curls to the Non-Linear Response of Shells
    Oxford: Oxford University Press.
  2. Book
    1. Burnett DS
    Finite Element Analysis: From Concepts to Applications
    Reading MA: Addison-Wesley.
    1. Électricité de France
    (2001) Finite element code_aster, Analysis of Structures and Thermomechanics for Studies and Research
    Finite element code_aster, Analysis of Structures and Thermomechanics for Studies and Research,
    1. Harris JE
    2. Crofton HD
    Structure and function in the Nematodes: internal pressure and cuticular structure in Ascaris
    The Journal of Experimental Biology 34:116–130.
  3. Book
    1. Kenney JF
    2. Keeping ES
    Mathematics of Statistics
    Princeton NJ: Van Nostrand.
  4. Book
    1. Landau LD
    2. Lifshitz EM
    3. Kosevich AM P
    4. Elasticity T
    Course of Theoretical Physics
    1. Papoulis A
    McGraw-Hill Series in Electrical Engineering
    Probability, random variables, and stochastic processes, McGraw-Hill Series in Electrical Engineering, New York, McGraw-Hill.
    1. Peng AW
    2. Ricci AJ
    Auditory and Vestibular Research. Methods in Molecular Biology
    Glass Probe Stimulation of Hair Cell Stereocilia, Auditory and Vestibular Research. Methods in Molecular Biology, 1427, New York NY, Humana Press, 10.1007/978-1-4939-3615-1_27.
  5. Book
    1. Phillips R
    2. Kondev J
    3. Theriot J
    4. Garcia H
    Physical Biology of the Cell
    Taylor & Francis Group.
  6. Book
    1. Rose C
    2. Smith MD
    Mathematical Statistics with Mathematica
    New York: Springer-Verlag.

Article and author information

Author details

  1. Alessandro Sanzeni

    1. Department of Physics, University of California, San Diego, La Jolla, United States
    2. National Institute of Mental Health Intramural Program, National Institutes of Health, Bethesda, United States
    Conceptualization, Software, Formal analysis, Validation, Investigation, Methodology, Writing—original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8758-1810
  2. Samata Katta

    Neuroscience Program, Stanford University School of Medicine, Stanford, United States
    Conceptualization, Formal analysis, Validation, Investigation, Methodology, Writing—original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9748-1653
  3. Bryan Petzold

    Department of Mechanical Engineering, Stanford University, Stanford, United States
    Validation, Investigation, Visualization, Writing—original draft
    Competing interests
    No competing interests declared
  4. Beth L Pruitt

    1. Department of Mechanical Engineering, Stanford University, Stanford, United States
    2. Department of Bioengineering, Stanford University, Stanford, United States
    Supervision, Validation, Investigation, Writing—original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4861-2124
  5. Miriam B Goodman

    Department of Molecular and Cellular Physiology, Stanford University School of Medicine, Stanford, United States
    Conceptualization, Formal analysis, Supervision, Validation, Investigation, Methodology, Writing—original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5810-1272
  6. Massimo Vergassola

    Department of Physics, University of California, San Diego, La Jolla, United States
    Conceptualization, Formal analysis, Supervision, Validation, Investigation, Methodology, Writing—original draft
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7212-8244


National Institutes of Health (R01-NS-047715)

  • Miriam B Goodman

National Institutes of Health (R35-NS-105092)

  • Miriam B Goodman

National Institutes of Health (R01-EB-006745)

  • Beth L Pruitt

National Institutes of Health (F31-NS-093825)

  • Samata Katta

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


We thank Z Liao for technical assistance with C. elegans animals. Experimental research was supported by NIH grants (R01-NS-047715, R35-NS-105092 to MBG; R01-EB-006745 to BLP and MBG) and fellowships (F31-NS-093825 to SK).

Version history

  1. Received: October 30, 2018
  2. Accepted: July 17, 2019
  3. Version of Record published: August 13, 2019 (version 1)


© 2019, Sanzeni et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 1,564
  • 206
  • 13

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Alessandro Sanzeni
  2. Samata Katta
  3. Bryan Petzold
  4. Beth L Pruitt
  5. Miriam B Goodman
  6. Massimo Vergassola
Somatosensory neurons integrate the geometry of skin deformation and mechanotransduction channels to shape touch sensing
eLife 8:e43226.

Share this article

Further reading

    1. Cell Biology
    2. Physics of Living Systems
    Ivan Castello-Serrano, Frederick A Heberle ... Ilya Levental
    Research Article

    The organelles of eukaryotic cells maintain distinct protein and lipid compositions required for their specific functions. The mechanisms by which many of these components are sorted to their specific locations remain unknown. While some motifs mediating subcellular protein localization have been identified, many membrane proteins and most membrane lipids lack known sorting determinants. A putative mechanism for sorting of membrane components is based on membrane domains known as lipid rafts, which are laterally segregated nanoscopic assemblies of specific lipids and proteins. To assess the role of such domains in the secretory pathway, we applied a robust tool for synchronized secretory protein traffic (RUSH, Retention Using Selective Hooks) to protein constructs with defined affinity for raft phases. These constructs consist solely of single-pass transmembrane domains (TMDs) and, lacking other sorting determinants, constitute probes for membrane domain-mediated trafficking. We find that while raft affinity can be sufficient for steady-state PM localization, it is not sufficient for rapid exit from the endoplasmic reticulum (ER), which is instead mediated by a short cytosolic peptide motif. In contrast, we find that Golgi exit kinetics are highly dependent on raft affinity, with raft preferring probes exiting the Golgi ~2.5-fold faster than probes with minimal raft affinity. We rationalize these observations with a kinetic model of secretory trafficking, wherein Golgi export can be facilitated by protein association with raft domains. These observations support a role for raft-like membrane domains in the secretory pathway and establish an experimental paradigm for dissecting its underlying machinery.

    1. Microbiology and Infectious Disease
    2. Physics of Living Systems
    Chi Zhang, Rongjing Zhang, Junhua Yuan
    Research Article

    Bacteria in biofilms secrete potassium ions to attract free swimming cells. However, the basis of chemotaxis to potassium remains poorly understood. Here, using a microfluidic device, we found that Escherichia coli can rapidly accumulate in regions of high potassium concentration on the order of millimoles. Using a bead assay, we measured the dynamic response of individual flagellar motors to stepwise changes in potassium concentration, finding that the response resulted from the chemotaxis signaling pathway. To characterize the chemotactic response to potassium, we measured the dose–response curve and adaptation kinetics via an Förster resonance energy transfer (FRET) assay, finding that the chemotaxis pathway exhibited a sensitive response and fast adaptation to potassium. We further found that the two major chemoreceptors Tar and Tsr respond differently to potassium. Tar receptors exhibit a biphasic response, whereas Tsr receptors respond to potassium as an attractant. These different responses were consistent with the responses of the two receptors to intracellular pH changes. The sensitive response and fast adaptation allow bacteria to sense and localize small changes in potassium concentration. The differential responses of Tar and Tsr receptors to potassium suggest that cells at different growth stages respond differently to potassium and may have different requirements for potassium.