Abstract
A mechanistic description of the generation of whisker movements is essential for understanding the control of whisking and vibrissal active touch. We explore how facial-motoneuron spikes are translated, via an intrinsic muscle, to whisker movements. This is achieved by constructing, simulating, and analyzing a computational, biomechanical model of the motor plant, and by measuring spiking to movement transformations at small and large angles using high-precision whisker tracking in vivo. Our measurements revealed a supralinear summation of whisker protraction angles in response to consecutive motoneuron spikes with moderate interspike intervals (5 ms < Δt < 30 ms). This behavior is explained by a nonlinear transformation from intracellular changes in Ca2+ concentration to muscle force. Our model predicts the following spatial constraints: (1) Contraction of a single intrinsic muscle results in movement of its two attached whiskers with different amplitudes; the relative amplitudes depend on the resting angles and on the attachment location of the intrinsic muscle on the anterior whisker. Counterintuitively, for a certain range of resting angles, activation of a single intrinsic muscle can lead to a retraction of one of its two attached whiskers. (2) When a whisker is pulled by its two adjacent muscles with similar forces, the protraction amplitude depends only weakly on the resting angle. (3) Contractions of two adjacent muscles sums up linearly for small amplitudes and supralinearly for larger amplitudes. The model provides a direct translation from motoneuron spikes to whisker movements and can serve as a building block in closed-loop motor–sensory models of active touch.
Introduction
Most sensory organs in mammals, including whiskers of rodents, are embedded within rich muscular systems. For example, eyes are controlled by muscles that allow a wide range of motor behavior (Haslwanter, 2002), and hands contain both the highest density of tactile receptors and a complex motor apparatus (An et al., 1989). Active control of even a minute sensory movement is believed to be essential for accurate sensation (Ahissar and Arieli, 2001; Kleinfeld et al., 2006), perhaps at the level of movements induced by single motor spikes (Simony et al., 2008). Understanding of sensation-targeted motor control at a single spike level has not yet been possible so far in any of the active senses, including vision; as far as we know, oculomotor control has been modeled only at a level of population firing rate, and not at a level of spiking neurons. The vibrissal system of the rat provides an excellent model for the investigation of sensation-targeted motor control due to its robustness in reduced dimensionality; rodents exert large whisker movements primarily along the posterior–anterior axis (Bermejo et al., 2002; Kleinfeld et al., 2006; Hill et al., 2008; Knutsen et al., 2008). Moreover, the detailed anatomy of the vibrissa musculature is largely known (Dörfl, 1982; Haidarliu et al., 2010). Intrinsic muscles, which connect two adjacent whiskers in a row, are sling muscles that wrap around the more anterior follicle above the plate, and are connected to the more posterior follicle near the skin (Dörfl, 1982). Extrinsic muscles are attached to bones and move substantial areas of the pad (Haidarliu et al., 2010). A whisking cycle during exploratory movement is composed of three phases. Contraction of the protracting extrinsic muscle m. nasalis, followed by the contraction of intrinsic muscles cause protraction. Subsequent contraction of the retracting extrinsic muscles, together with passive properties of the pad, cause retraction (Berg and Kleinfeld, 2003; Hill et al., 2008; Haidarliu et al., 2010).
While the general muscular scheme of whisking seems to be understood, the principle of temporal and spatial integration by the muscles and their effect on whisker movement is still mostly unknown. This issue was addressed experimentally by Herfst and Brecht (2008) using a juxtacellular method (Pinault, 1996) to evoke single spikes in individual motoneurons (MNs) in the lateral facial nucleus of the anesthetized rat, and a high-resolution video camera to measure resulting whisker movements. In some of their experiments, a supralinear temporal summation of whisker protraction was observed in response to two consecutive spikes of an MN. Spikes of MNs that stimulate putative intrinsic muscles caused measurable protraction of no more than two adjacent whiskers in a row, with variable ratios of protraction angles between the whiskers and various patterns of temporal summation. The source of this variability is still not understood.
In this work, we analyze, simulate, and characterize transformations between spiking of facial MNs and the movement of whiskers generated by the contraction of intrinsic muscles these neurons innervate. A combination of experimental, modeling, and theoretical techniques were used to investigate the following: (1) the force developed by an intrinsic muscle, and the resulting vibrissae movements generated, in response to consecutive spikes fired by an MN; (2) how two neighboring whiskers protract in response to the contraction of the intrinsic muscle connecting them; and (3) how movement of a single whisker is affected by the contraction of two muscles attached to it.
Materials and Methods
Animal preparations and artificial whisking.
Experiments were performed on six male albino Wistar rats weighing 200–300 g. Rats were anesthetized with urethane (1.5 g/kg, i.p.), secured in a stereotaxic device (SR-6; Narishige), and placed on a servo-controlled heating blanket. Supplemental doses (10% of initial) of anesthesia were administered when required. During surgery and experimental manipulations, body temperature was maintained at 37°C and monitored rectally.
Artificial whisking was induced by stimulating the upper buccal motor branch of the facial nerve (Fig. 1A,B) (Semba and Egger, 1986; Szwed et al., 2003) as follows. The facial nerve was exposed and cut, its distal end was mounted on bipolar silver electrodes, and the nerve was kept moist by frequent brief washes with warm saline between periods of stimulation. Single-, two-, and three-pulse biphasic rectangle electrical stimuli (amplitude 0.5–1.5 V, duration of 40 μs) were applied through an isolated pulse stimulator (STG2000; Multi Channel Systems). The interpulse intervals, Δts, varied from 5 to 40 ms with 2 ms resolution, and from 40 to 80 ms with 10 ms resolution. Additional interpulse intervals from 1 to 4 ms with 1 ms resolution were applied to three of nine whiskers. Artificial whisking was induced by 5 Hz (50% duty cycle) trains for 2 s for 10 cycles. We recorded whisker movements at 1000 frames/s with a fast digital video camera (MotionScope PCI 8000; Redlake;), upon stimulation of the facial nerve. The resulting whisker movements were quantified using the protraction angle θ(t) (Fig. 1C). Negative values of θ(t) indicate retraction.
Juxtacellular stimulation of vibrissa motoneurons.
Herfst and Brecht (2008) used a juxtacellular method (Pinault, 1996) to evoke single spikes from individual vibrissa motoneurons (vMNs) located in the lateral facial nucleus (Fig. 1A), and high-resolution video to measure the resulting whisker movements. We reanalyzed whisker trajectories (32 trials) from four experiments of their study to characterize the properties of temporal summation of spikes by the intrinsic muscle, and the biomechanical plant. In each of the four experiments, a single MN was stimulated. The temporal resolution of the recording camera was different for each of the four experiments used (1500, 670, 165, and 1000 fps). The response of adjacent whiskers to motoneuron spikes was evaluated by analyzing data from three experiments in which the following was true: (1) several neighboring whiskers in a row, and the angles between those whiskers and the skin, were tracked; (2) movements were confined for two adjacent whiskers, not for the entire row; and (3) peak whisker amplitude could be extracted. All three experiments were recorded at 50 Hz deinterlaced.
Model of force generation in muscle cells.
Control of whisker movements by MN spikes is modeled as a two-stage process (Fig. 1B). First, the spikes are translated to force in the muscle cell. Second, the resulting muscle force causes whisker movement. We will start by describing our model for generation of force in muscle cells. Suppose that the MN fires spikes at times tk, k = 1, … , j. MN spikes induce spikes in muscle cells via neuromuscular junctions, and muscle spikes trigger the excitation–contraction cycle in the muscle cell. A simple model for transformation from motoneuron spiking to muscle force production, consisting of five steps, is presented here (Fig. 1D).
An action potential in the muscle activates a voltage-sensor protein, the dihydropyridine receptor (Rhoades and Pflanzer, 2003). This protein opens ryanodine receptors (RyRs) on the membrane of the sarcoplasmic reticulum (SR), and Ca2+ is released from the SR to the cytoplasm. If the muscle cell has not fired for a long time (say, >50 ms), the fraction of the open RyRs, r, is 0. We assume that r increases to 1 immediately following the firing of an action potential, which is assumed to be very brief (<1 ms), and decays exponentially with a time constant τr. This assumption is based on several studies with skeletal muscles (Miledi et al., 1982; Stein et al., 1988), in which intracellular Ca2+ transients, in response to repetitive action potentials, exhibited sublinear processes. The value of r (fraction of open RyRs) between spikes k and k + 1 is therefore Because of the saturation effect of r (to 1) by each spike, the process r(t) exhibits sublinear temporal summation. The sublinearity is substantial for brief interspike intervals (ISIs), with time difference Δt between spikes on the order of τr or less.
Depending on the fraction of open RyR receptors r(t), Ca2+ ions are released from the SR, and are sequestered to longitudinal elements of the SR by an active transport mechanism (Rhoades and Pflanzer, 2003). This mechanism is assumed to depend linearly on the normalized intracellular Ca2+ concentration, [Ca2+]i, and to have a time scale τc. Therefore, [Ca2+]i evolves according to the following equation: where r0 is a constant. The solution of Equation 2 with the initial condition [Ca2+]i(t = tk) = Ck is where Ck = 0 for k = 1 and for k ≥ 1, where Δtk = tk+1 − tk.
The muscle force (F) depends on the level of [Ca2+]i (Rhoades and Pflanzer, 2003), and this dependency is denoted by FC([Ca2+]i). The function FC([Ca2+]i) is nonlinear because several (n) Ca2+ ions are needed to bind to troponin to enable the excitation–contraction cycle (Bobet and Stein, 1998). Following Bobet and Stein (1998), we assume that n = 4 and hence where A is a force scaling factor and A0 = 1 mg · mm/ms2. Equation 4 accounts for the supralinear temporal summation for moderate ISIs, on the order of τc (see Discussion).
The muscle force depends also on muscle length. This dependency is known as the “force–length curve” (Rhoades and Pflanzer, 2003), FL(z), where lm(t) is the instantaneous muscle length, and l0m is the muscle length at rest. Following Hill et al. (2008), FL(z) is modeled as where 0 < zh < zl < 1. Note that FL(z) = 1 for small protraction angles (Fig. 1D). Following Hill et al. (2008), we ignore the velocity dependence of the muscle force.
The total muscle force Fm(t) generated by consecutive MN spikes is a product of FC([Ca2+]i) and FL(z): We use the following parameter set, referred to as the “reference parameter set,” throughout the paper unless otherwise stated: r0 = 1.9, τr = 5 ms, τc = 6 ms, zh = 0.1, zl = 0.45, A = 1.33. This parameter set was chosen to fit the model to data from stimulation of a particular single MN (see Fig. 4 below). These parameters were chosen to minimize the sum over responses to one spike and to two spikes, of the integrals over time of the mean-square difference between the experimentally measured θ(t) and the function θ(t) obtained from the model. The integrals were computed from t = 0 to 17 ms after the time of the last stimulating spike.
Biomechanical model.
Our model of a row of whiskers in the motor plant consists of equations of motion for the whiskers (follicles and hairs) as a functions of the forces generated by the muscles and by the viscoelastic properties of the mystacial pad. Since the goal here is to quantify whisker dynamics in response to MN spikes that cause contractions of intrinsic muscles, only intrinsic muscles are included in the model. The viscoelastic tissue is modeled as a system of linear springs and dampers, an approximation good for small protraction angles. Deviations from the linear representation are expected for large protraction angles (see Discussion). A scheme of our model for a row of whiskers is plotted in Figure 1E, and the model [based on a previous model of Hill et al. (2008)] is fully described in Appendix A [including its reference parameter set (see Table A2)]. Specifically, each whisker is connected to springs and dampers in the plane of the skin and in the plane of the plate, and we ignore passive forces applied by muscle stretching (see Rhoades and Pflanzer, 2003). The direction along the row is denoted by the x-direction (posterior–anterior), and the direction from the plate to the skin is the y-direction (medial–lateral). There are four main differences between our model and that of Hill et al. (2008).
In our model, each whisker is coupled by springs and dampers to the tissue, and not to the adjacent whiskers. This modification was made because the resistance of tissue to movements of one whisker does not disappear if the adjacent whiskers are removed. Furthermore, our simulations of the model of Hill et al. (2008) revealed that when every whisker in a row is coupled to its neighbors, contraction of one intrinsic muscle causes motion of all the whiskers in a row with similar protraction angles, a finding that contradicts the experimental results of Herfst and Brecht (2008).
Each follicle is coupled at its bottom to the tissue by springs (with spring constants kπ,y) and dampers (with spring constants ζπ,y) in the y-direction. This coupling is needed to prevent slow oscillations of the whole row. In our model, coupling in the y-direction occurs near the plate and not near the skin because the resistance of the tissue near the plate, closer to the bone, to movements in the medial–lateral direction is considerably larger than the resistance of the tissue near the skin, adjacent to the outside air. As a result of attachment to the plate, the vertical spring and damper constant, kπ,y and ζπ,y, respectively, are larger than their horizontal counterparts, kπ,x and ζπ,x.
The most anterior whisker is coupled by one intrinsic muscle only to its posterior neighbor (Dörfl, 1982). In our model, in agreement with this anatomical finding (Dörfl, 1982) and in contrast to the model of Hill et al. (2008), there is no intrinsic muscle connected anteriorly to the most anterior whisker.
An intrinsic muscle is attached to its posterior whisker at the skin, and to its anterior whisker somewhere along the follicle, closer to the plate than to the skin. In our model, the distance along the follicle from this attachment point to the skin is denoted by ld, where ld is a parameter of the system, and not necessarily equal to −2C as in Hill et al. (2008), where C is the distance along the whisker from the whisker's center of mass to the skin. C < 0 because the whisker's center of mass is below the skin.
Here, we consider a row of NW whiskers and denote the index of whiskers in a row by i, the position of the whisker's center of mass with respect to its position at rest (steady state with no muscle contraction) by x⃗i(t) = (xi(t), yi(t)), and the angle between the whisker and the tangent to the pad at the point of whisker implant (on the left side of the x-axis) by θi(t) (Fig. 1C). The value of θi(t) at rest is θ0i, and the protraction angle is θi(t) = θi(t) − θ0i. For simplicity, we assume that all the resting angles are equal: θ0i = θ0. The mass and moment of inertia of each whisker are M and I, respectively. The equations of motion for x⃗i(t) and θi(t) are governed by the forces and the torques applied on the ith whisker. Each whisker is affected by the total force (F⃗), and the torque in the z-direction (Nz) applied by the springs (superscript s), dampers (superscript d), and muscles (superscript m): The minus sign in Equation 10 appears because the angle θi(t) is defined relative to the negative x-axis.
Quantifying temporal and spatial whisking dynamics.
Suppose that in response to j consecutive spikes fired by an MN, the protraction angle of the ith whisker is θij(t). Then the linearity index of the response to j spikes, LIj, is defined as The integration times, Tj, should be large enough to enable the whisker to get back to rest. For analysis of data from the facial motor nerve (FMN) stimulation and from the numerical simulations, we use Tj = 200 ms. For analysis of data from vMN stimulations, we use Tj = (j − 1)Δt + 40 ms. The linearity index LIj is 1 if the summation is linear, and is >1 or <1 if the summation is supralinear or sublinear, respectively.
In response to stimulation of one MN, or several MNs simultaneously, a whisker moves with a profile θi(t). If the whisker only protracts, θi(t) > 0. If there are time intervals during which the whisker retracts, θi(t) < 0 during these intervals. The maximal protracting angle is denoted by θimax and the minimal angle, i.e., the retraction angle with the maximal absolute value, is denoted by θimin. The “peak amplitude” θipeak is defined to be θimax if the whisker does not retract at all or if θimax > |θimin|, and is defined to be θimin otherwise.
For comparison of the protraction of adjacent whiskers, we use R, the “peak amplitude ratio,” which is defined as the ratio between the peak amplitudes of the anterior and the posterior whiskers, θanteriorpeak and θposteriorpeak, respectively:
Results
Temporal summation of muscle force and whisker protraction
FMN stimulation mainly affects intrinsic muscles
To quantify temporal summation properties of whisker protraction in response to MN firing for protraction angles that are not small, we induced electrical stimulation to the upper buccal motor branch of the FMN. The nerve was stimulated every 200 ms either by one pulse or by bursts of two and three consecutive pulses with equally spaced interpulse intervals Δt (see Materials and Methods). The time courses of the whisker protraction angles θ(t) of nine whiskers were measured; negative values of θ(t) indicate retraction.
We tracked xσ(t), the x-coordinate (i.e., the coordinate along the posterior–anterior direction) of the point where the whisker emerges from the skin (Fig. 1E and Appendix A), and used this coordinate to quantify pad movement. Sample trajectories of θ(t) and xσ(t) are presented in Figure 2A. During each electrical stimulation pulse of the FMN, many MN axons were stimulated at once. In principle, these axons can innervate both intrinsic and extrinsic muscles. We assessed which muscle type dominates whisker protraction in our experiments by using previous measurements of Hill et al. (2008) (their Fig. 7G). Those authors plotted pad translation versus whisker angle at peak movement during stimulation of three muscles types: intrinsic muscles (INT), the protracting extrinsic muscle m. nasalis (NA), and retracting extrinsic muscles (RET), m. nasolabialis and m. maxillolabialis (see Fig. 2B). Values of xσ(t) as a function of θ(t) for our nine whiskers were plotted under the same conditions as the data from Hill et al. (2008). The trajectories of our xσ(t)s cluster around those for the INT of Hill et al. (2008) (Fig. 2B) (p = 0.6992, Wilcoxon test), and are distinct from their trajectory for NA (p = 0.0039, Wilcoxon test). The mean slope of our FMN stimulations (0.037 ± 0.005 mm/deg) was similar to the slope of their (Hill et al., 2008) INT (0.033 mm/deg). These results indicate that whisker movements resulting from FMN stimulations are governed mostly by contraction of intrinsic muscles.
Temporal summation of responses to consecutive pulse stimulations of the FMN
The protraction angles θ(t) of one whisker in response to two and three pulses with Δt = 17 ms, each averaged over 10 trials, are depicted in Figure 3, A and B, respectively. We also plotted the trial-average response to one pulse, and the computed response to two and three pulses, under the assumption of linear temporal summation. The protraction in response to two pulses is larger than that expected by the linear calculation (Fig. 3A). This supralinear summation is even more apparent in the response to three consecutive pulses (Fig. 3B).
We quantified the temporal summation by calculating the linearity indices LI2 and LI3 (Eq. 11) for two and three pulses, respectively, each averaged over nine whiskers for Δt ≥ 5 ms. For 1 ms ≤ Δt ≤ 4 ms, the average was calculated over three whiskers. The whisker-averaged values of LI2 and LI3 versus Δt are plotted in Figure 3C; the corresponding values for the single whisker (see Fig. 3A,B) are plotted in the inset. For Δt = 1 ms, both LI2 and LI3 are smaller than 1, reflecting sublinear temporal summation. The two linearity indices increase as Δt increases above that value, and become larger than 1 for Δt between 3 and 6 ms. LI2 and LI3 remain larger than 1 below ΔtL, which is the minimal interval for linear summation (∼35 ms). This reflects supralinear temporal summation. Above ΔtL, LI2 and LI3 are essentially 1, and the temporal summation is linear. The two linearity indices (LI2 and LI3) increase with Δt for 1 ms < Δt < ΔtM ∼ 17 ms, and decrease for Δt > ΔtM. Interestingly, ΔtM is similar to the Δt value for which protraction in response to a single twitch is maximal (Fig. 3A). The maximal values of LI3 over Δt are plotted in Figure 3D as a function of the corresponding maximal values of LI2 for the nine whiskers. The
Temporal summation of responses to consecutive spike of a single MN
Nerve stimulations excite many axon terminals and many muscle cells. Therefore, we examined whether temporal summation is also supralinear when only a single motoneuron MN fires. This was achieved by our reanalysis of the results of a previous study (Herfst and Brecht, 2008) in which single MNs were stimulated by juxtacellular injections of depolarizing step currents (see Materials and Methods). Those step currents evoked one or more spikes in an MN, which resulted in whisker rotation. The ISIs varied from trial to trial, and could not be directly controlled by the experimenter. Examples of whisker trajectories resulting from one MN spike and from two consecutive MN spikes are shown (Fig. 4A,B). The supralinear temporal summation is demonstrated by comparing the whisker response to two spikes with the computed putative trajectory in response to those two spikes under the assumption of linear temporal summation. Plotting the linearity index LI2 as a function of Δt for four motoneurons (Fig. 4C) shows that temporal summation is supralinear for 13 ms < Δt < 30 ms, and linear for larger Δt. Supralinear summation was evident for each of the neurons studied, including one whose responses were classified before as sublinear and linear [Herfst and Brecht (2008), their Fig. 6C,D] (red symbols in Fig. 4C here); this neuron showed supralinear responses like the other neurons once the oscillations that were superimposed on the retraction phase of all its trajectories were removed. The similarity of the supralinear behavior of this whisker to the other studied whiskers (Fig. 4C) suggests that the oscillations superimposed on its trajectories did not have muscular origin; these oscillations might have resulted from a mechanical resonance induced by the label attached to the whisker.
Modeling temporal summation of whisker protraction responses for small protraction angles
The dependence of the linearity index LIj(Δt) on Δt and the number of spikes j was examined by analyzing our model of the vibrissal response to firing of a single MN (for detailed model description, see Material and Methods and Appendix A). Our model describes muscle contraction (Eqs. 1–8) and whisker-follicle translation and rotation (Eqs. 9, 10, A1–A30). We start by considering small protraction angles θ(t) and small whisker translations x⃗(t), for which the dynamical equations of the biomechanical model of the plant can be linearized to the first order of the dynamical variables (Appendix A, Eqs. A2, A3). Under these conditions, we show in Appendix B that if only one intrinsic muscle contracts in response to j consecutive spikes with a temporal force profile Fm,j(t) (where m stands for muscle and j for in response to j spikes), or if the two intrinsic muscles on the two sides of a whisker contract with the same temporal profile Fm,j(t), the linearity index LIj depends only on the muscle properties, not on the properties of the biomechanical motor plant, i.e., Therefore, the dependence of LIj on Δt can be accounted for solely by the muscle model (Eqs. 1–8).
We calculate, using Equation 11, the protraction angle and LI2 by simulating the plant model (with five whiskers in a row), while activating the muscle between the second and the third whisker and computing θ3(t). Responses of the model to one spike (Fig. 4A) and to two spikes (Fig. 4B) with Δt = 19 ms fired by a single MN are depicted for the reference parameter set and A = 1.33. The dependence of the linearity index LI2 on Δt is shown in Figure 4C. The model responses fit the experimental recordings of whisker responses (Fig. 4) very well. We use the model and the fact that LIj depends only on muscle dynamics (Eq. 13) to explain the temporal summation properties of whisker responses to consecutive MN spiking. Note that since the protraction angles and the whisker translations are small, the muscle contraction length ratio z = lm(t)/l0m(t) (Eq. 6) obeys |z − 1| < zh (see Eq. 7), and hence the force–length function FL(z) equals 1. Therefore, Fm(t) = FC([Ca2+]i) (Eq. 8).
The temporal summation properties of the muscle force response are the result of two competing processes. The process r(t), which represents the fraction of open RyRs, exhibits sublinear temporal summation properties (Eq. 1) if Δt is small in comparison with τr. The dependence of the force term FC([Ca2+]i) on [Ca2+]i is supralinear (Eq. 5), unless Δt is considerably larger than τc. For very small Δt, the process r(t) is almost unaffected by spikes that follow the first spike, and LIj = 1/j. This accounts for the muscle refractoriness. As Δt increases, the additional value of r(t) in response to the second and third spikes becomes more substantial, and the temporal summation properties of the r(t) process become more linear. Concomitantly, the nonlinear dependence of FC([Ca2+]i) on [Ca2+]i leads to supralinear temporal summation. Those two factors cause the indices LIj to increase with Δt, until they reach their maximal values (well above 1) at approximately ΔtM. Above this time interval (ΔtM), as Δt increases much beyond the sequestration time constant τc, the level of [Ca2+]i depends only on the time elapsed from the latest spike (Eqs. 3, 4). As a result, the muscle force in response to the second and third spikes becomes similar to the one developed in response to the first spike, and LIj becomes closer to 1.
Modeling temporal summation for protraction angles that are not small
When protraction angles are not small, Equation 13 is no longer strictly valid. Moreover, under such conditions, the linearity index is significantly affected by the force–length function, FL(z), and by the biomechanics of the plant. Therefore, for protraction angles that are not small (θ > 5°), we study temporal summation by simulating the full model of the vibrissa plant. Results from experiments in which the facial nerve was stimulated were compared to those of simulation of our model of a row of five whiskers, when all five intrinsic muscles are activated simultaneously by one “pulse.” Traces of computed angles θ3(t) of the intermediate whisker in response to one, two, and three consecutive pulses with Δt = 15 ms, and the linearity indices LIj versus Δt for the simulated and experimental results, are shown (Fig. 5). The corresponding experimental results are from a single tracked whisker out of five in B row. The computed results fit the experimental observations well. A good fit of modeling and experimental results for protraction angles θ > 20° is achieved when our model includes the force–length function, FL(z). Without this function, our model predicts larger protraction angles in response to three consecutive spikes than the angles observed experimentally (supplemental Fig. S1, available at www.jneurosci.org as supplemental material).
Protraction of two adjacent whiskers in response to contraction of their connecting muscle
Simulation results for small protraction angles
The viscoelastic properties of the tissue are modeled by coupling the whiskers to the pad with springs and dampers (Fig. 6A). Therefore, contraction of a single intrinsic muscle may cause movement only of two whiskers, the posterior and the anterior to that muscle. A reduced system composed of one intrinsic muscle, its adjacent whiskers, and the springs and dampers coupled to these whiskers are therefore the basic whisking motor unit (WMU).
We characterize the response of the WMU to MN spikes by studying the protraction dynamics of the two whiskers next to the activated muscle. We start with a small muscle force (and, therefore, small protraction angles), and use the reference parameter set and ld = 3 mm. An example of plant morphology in which whiskers point posteriorly during rest with resting angle θ0 = 70° is depicted (Fig. 6A). Time traces of the protraction angles θ (Fig. 6B) and the coordinates of the whisker center of mass (x, y) in response to a single small muscle twitch (Fig. 6C) are plotted as functions of time for the two adjacent (posterior and anterior) whiskers. Both whiskers protract, albeit the anterior whisker more. The contraction component of the muscle in the x-direction causes the whiskers to attract each other, and therefore the x-coordinate of the centers of mass of the posterior and the anterior whiskers are positive and negative, respectively (Fig. 6C, top panel). Since the muscle is attached to the posterior whisker near the skin and wraps around the anterior whisker more medially, i.e., closer to the plate, contraction of the muscle causes the coordinate y of the posterior and the anterior whiskers to be negative and positive, respectively (Fig. 6C, bottom panel). In the other example of plant morphology presented here, the whiskers point anteriorly during rest with resting angle θ0 = 95° (Fig. 6D). Time traces of the protraction angle θ for this example are shown (Fig. 6E). As in the previous example, in response to a single, small muscle twitch, the two whiskers protract, but in this case the posterior whisker protracts more. The movements of x and y are similar to those for θ0 = 70° (data not shown).
The two examples in Figure 6A–E demonstrate general behavior observed in all our simulations: the peak amplitude of the posterior whisker, θposteriorpeak, increases with θ0, whereas the peak protraction amplitude of the anterior whisker θanteriorpeak, decreases with θ0 (Fig. 6F). The two amplitudes are equal for θ0 = θ0,eq (∼88° for the parameters of Fig. 6). When θ0 < θ0 = θ0,pr (∼50° for the parameters of Fig. 6), the posterior whisker retracts, whereas the anterior whisker protracts (Fig. 6G). When θ0 > θ0 = θ0,ar (∼108° for the parameters of Fig. 6), the anterior whisker retracts, whereas the posterior whisker protracts (Fig. 6H). The sudden jumps in the curves of the peak amplitudes of the posterior and anterior whisker, for θ0 = 50° and θ0 = 108°, respectively (Fig. 6F), occur when the maximal whisker movement changes from positive (protraction) to negative (retraction), and θpeak switches signs (see Materials and Methods, Quantifying temporal and spatial whisking dynamics). The peak amplitude ratio, R = θanteriorpeak/θposteriorpeak (Eq. 12) decreases with θ0 for all θ0 values except for angle θ0,pr, where R switches from −∞ to ∞ as θposteriorpeak switches from negative to positive values (Fig. 7A, inset).
ld is the distance measured along the follicle from the skin to the attachment point where the muscle cell slings around the anterior whisker (Fig. 1E). If ld is larger, the torque applied on the anterior whisker is also larger, and thus, the protraction of this whisker is expected to be larger relative to the torque on the posterior whisker, and therefore R should be larger. As long as the posterior whisker protracts, this expectation is confirmed in our simulations (Fig. 7A), where R decreases with θ0 and increases with ld. In addition, increasing ld increases θ0,pr, and θ0,ar, and therefore, increases the θ0 range in which the posterior whisker retracts (Fig. 7A, right side) and decreases the θ0 range in which the anterior whisker retracts.
The vertical spring constant kπ,y and vertical damping constant ζπ,y define the viscoelastic properties of the connective tissue along the y-axis (medial–lateral). We study how varying these constants, while keeping the ratio ζπ,y/kπ,y fixed, affects the protraction ratio R (Fig. 7B). For kπ,y = 0 and ζπ,y = 0, R increases slightly with θ0 and the value of R is always smaller than 1, i.e., the posterior whisker protracts more than the anterior whisker. Moreover, under these conditions, the whisker coordinates θ, x and y exhibit slow oscillations that decay slowly with time (data not shown). These oscillations disappear if kπ,y and ζπ,y are positive. Elevating kπ,y and ζπ,y first causes R to decrease with increasing θ0, and then increases the steepness of the dependence of R on θ0. As long as the posterior whisker protracts and θ0 < 90°, R increases with kπ,y and ζπ,y. θ0,pr and θ0,ar also increase with kπ,y and ζπ,y.
The horizontal spring constant kπ,x and horizontal damping constant ζπ,x define the viscoelastic properties of the tissue along the x-axis at the bottom of the follicle (Dörfl, 1982). For our reference parameter set, we choose these constants to be equal to the corresponding constants at the skin: kπ,x = kσ,x and ζπ,x = ζσ,x. Here, we examine the effect of increasing kπ,x and ζπ,x, while keeping kσ,x, ζσ,x, and the ratio ζπ,x /kπ,x fixed. When the posterior whisker protracts, R decreases as the elastic constants kπ,x and ζπ,x increase near the plate for all the values of θ0 we test (Fig. 7C), and R increases as the elastic constants kσ,x and ζσ,x increase near the skin (data not shown).
Our study of the effects of the dynamical parameters of whiskers on protraction dynamics reveals that the maximal amplitudes θposteriormax and θanteriormax depend only weakly on the center of mass C (Fig. 8A). Furthermore, varying the whisker mass, M, and the whisker moment of inertia, I, from 0 to three times their reference values also affects those amplitudes weakly (Fig. 8B,C). Specifically, the maximal amplitudes have a defined limit when M and I approach 0, and this limit is close to their value for the reference parameter set (see Table A2). This means that the peak amplitudes of whisker movement are determined primarily by interplay between the muscle forces and viscoelastic properties of the tissue, and only marginally by the dynamic properties of the whisker.
Analytical theory for small protraction angle, based on steady-state approximation
We find numerically that, in response to a single muscle contraction, the ratio R between the maximal protraction angles of the anterior and posterior whiskers decreases along with θ0 and along with kπ,x and ζπ,x. The reason for this behavior can be examined by developing a reduced model of the WMU (Fig. 6A,D) under certain limits. First, take the small θ limit, and assume that the protraction angles of both whiskers are small. Second, assume that the spring constants of the tissue in the y-direction, kπ,y, are much larger than those in the x-direction, kπ,x and kσ,x, because it is difficult to push the tissue toward the bone. Third, note that the dependency of R on the center of mass of the whisker, C, is very weak (Fig. 8A). Therefore, for simplicity, we arbitrarily perform the calculation for a whisker with a center of mass located in its bottom (near the plate), i.e., −C = lf, where lf is the length of the follicle. Fourth, since R increases with ld, examine the most extreme case for which ld = lf. Fifth, look at the steady-state condition, by assuming that the muscle force, Fm, is constant in time. This assumption enables us to simplify the dynamics considerably, and the simplification is accurate when the muscle twitch is slow enough. Under such conditions, we can calculate the ratio between the protraction angles of the anterior and the posterior whiskers, θanterior and θposterior, respectively (see Appendix C). We find (Eq. C14) that if kπ,x = kσ,x = kx, where s is the distance between two adjacent follicles. R > 1 if θ0 < θ0,eq (where θ0,eq = 90°) and R < 1 otherwise. The posterior whisker retracts if the denominator of Equation 14 is negative, i.e., θ0 < θ0,pr = arccos(s/lf) (see also Eq. C16); θ0,pr = 60° for our reference parameter set (see Table A2). The anterior whisker retracts if the numerator of Equation 14 is negative, i.e., θ0 < θ0,ar = arccos(−s/lf) (see also Eq. C17); θ0,ar = 120° for our reference parameter set. These results are similar to those obtained from simulations of the full system dynamics (Fig. 6F).
These analytically obtained results can be explained intuitively as follows. The whisker movement can be visualized as the relative movement between the skin translation, xσ, and the plate translation, xπ (see also Eq. C5). If θ0 < 90° (Fig. 6A), the muscle force and torque applied on the posterior whisker cause the bottom point of the posterior follicle to move anteriorly (xposteriorπ > 0) (Eq. C11). The more the point xposteriorπ moves anteriorly, the less the point where the whisker emerges from the skin, xposteriorσ, moves anteriorly (and it may even move posteriorly) (see Eq. C9). Therefore, if θ0 < 90°, the protraction of the posterior whisker is small, and it may even retract if θ0 is small enough. In contrast, If θ0 > 90° (Fig. 6D), the point xposteriorπ moves posteriorly (xposteriorπ < 0) (Eq. C11), the point xposteriorσ moves anteriorly, and the posterior whisker always protracts. The point where the anterior whisker emerges from the skin, xanteriorσ, does not move because no net force is applied on it, and no torque is applied on the anterior whisker (since ld = −C = lf). The force equation for xanteriorπ (Eqs. C6, C7) shows that the coordinate xanteriorπ moves posteriorly, and therefore the anterior whisker protracts, unless the top of the posterior follicle is more anterior than the bottom of the anterior follicle, which can only happen only if θ0 is large.
If kπ,x ≠ kσ,x, we find (Eq. C14) that This means that when the two whiskers protracts and R is positive, R decreases with the ratio kπ,x/kσ,x. Similarly, the angle θ0,pr decreases with kπ,x/kσ,x (Eq. C16). These results are consistent with those obtained by simulating the full model dynamics (Fig. 7C). Interestingly, the angle θ0,ar computed analytically is independent of kπ,x/kσ,x (Eq. C17). In simulations, we observe a weak dependency of θ0,ar on kπ,x/kσ,x.
This analytical calculation above was derived for ld = lf. We can explain, however, why the protraction of the anterior whisker is larger in comparison to that of the posterior whisker as ld increases. This stems from the fact that increasing ld increases the torque on the anterior whisker in the direction that causes it to protract.
Beyond small protraction angles
When the protraction angles of whiskers are not small, the θ traces of two whiskers following muscle contraction do not necessarily depend linearly on the amplitude, A, of the muscle force. Here, this dependency is studied via numerical simulations of the contraction of one muscle. The protraction profile of two whiskers surrounding one contracting muscle is shown for θ0 = 70° and A = 40 (Fig. 9A). The temporal profile of the anterior whisker looks similar to that obtained with smaller A (A = 1.33) (Fig. 6B), except for the larger amplitude. However, the profile of the posterior whisker exhibits protraction followed by retraction. In this case, in contrast to that with small angles (see above), the maximal protraction angle of the posterior whisker θposteriormax is larger than that of the anterior whisker θanteriormax. The angles θanteriormax and θposteriormax and the maximal retraction angles of the anterior and posterior whiskers, θanteriormin and θposteriormin, are shown as functions of A (Fig. 9B). Beyond the small angle regime (say, of θ < 3°), the increase of the angle θanteriormax with A becomes more moderate, until θanteriormax saturates for large A. The angle θanteriormin remains 0, i.e., the anterior whisker does not show any retraction. In contrast, the increase of θposteriormax with A is accelerated for large values of A, and θposteriormax eventually exceeds θanteriormax. Furthermore, if A is not small, the posterior whisker exhibits a retraction phase after the protraction phase. For θ0 = 95°, the dominance of θposteriormax over θanteriormax, which is already observed for small angles (Fig. 6E), also extends for larger angles (A = 40) (Fig. 9C). In this case, no prominent retraction of any whisker is observed. The dependencies of θanteriormax and θposteriormax on A are shown (Fig. 9D). As in the case of θ0 = 70°, the increase of those maximal angles with A for moderate forces is much smaller than for small forces. Video demonstrations of the model's response to contractions of a single intrinsic muscle for four values of θ0 (45°, 70°, 95°, and 120°) are available at www.jneurosci.org as supplemental material.
Stimulation of single motoneurons: a test of theoretical predictions
Based on our modeling and theoretical work, we predict that, for small protraction angles, the protraction amplitude ratio, R, between the anterior and posterior whisker can be larger than 1 when the whiskers point posteriorly at rest, and decreases with decreasing θ0. We examined whether these predictions are consistent with the experimental results described by Herfst and Brecht (2008). In that study, there were three video recordings in which the movements of two adjacent whiskers could be tracked, and therefore, θ0 could be determined. Data from one recording is presented (Fig. 10A). Stimulation of a single MN elicited a single spike, which caused an intrinsic muscle contraction. The anterior whisker, C3, had a larger protraction angle than the posterior whisker, C2. The R values of the three MNs are also shown (Fig. 10B). This small, and thus limited, data sample is consistent with our model.
Spatiotemporal summation of forces applied by two muscles
Typical whisking bouts usually involve contraction of several (or even all) intrinsic muscles in a coordinated manner during the protraction phase. We examined the spatial summation of forces applied by the two intrinsic muscles that are attached to the same whiskers (Fig. 11A), assuming the same force amplitude (A) for both. Note that in contrast to conditions described in Figures 3⇑–5, each muscle is stimulated here only once, and the two stimulating pulses are applied to two separate muscles. We define Δt = ta − tp, where ta and tp are the times of the stimulation pulses given to the anterior and the posterior muscle, respectively. The protraction angle of the whisker at time t following a pulse stimulus to the posterior muscle with force A at time Δt is denoted by θpost,muscle(t; A, Δt), and its maximal value over t is θpost,musclemax(A, Δt). θant,muscle(t; A, Δt) and θant,musclemax(A, Δt) are defined similarly for the anterior muscle. If the posterior muscle is stimulated at time 0 and the anterior muscle is stimulated at time Δt, each with a force amplitude A, the protraction angle of the whisker is θpost+ant,muscle(t; A, Δt) and its maximal value over t is θpost+ant,musclemax(A, Δt). In this analysis, we compare four quantities: θpost,musclemax(A, 0); θant,musclemax(A, 0); the maximum of the linear sum of the whisker responses,
The four functions defined above are plotted in Figure 11B for θ0 = 70° and Δt = 0 (simultaneous contraction of the two muscles). In this case, the spatial summation is linear for small protraction angles, i.e.,
We define the “spatiotemporal linearity index,” S(A, Δt), to be the following: The ratio S(A, Δt) is plotted versus the resting angle θ0, and Δt for A = 10 (Fig. 11D), and S(A, 0) as a function of θ0 is plotted in the inset. The value of S(A, 0) increases with θ0 for small θ0, and reaches a peak larger than 2 for θ0 = 73°. For larger θ0 values, S(A, 0) drops steeply to approximately 1.2 at θ0 = 95°, and then increases again. Therefore, for large force amplitudes, the spatial summation (Δt = 0) is supralinear and the magnitude of this supralinearity depends strongly on θ0. Similar behavior is observed for small Δt values (|Δt| below ∼20 ms). Furthermore, for all θ0 values, nonlinear temporal summation is maximal for positive Δt.
Finally, we investigate the peak protraction angle θpeak of a whisker whose two adjacent muscles are simultaneously stimulated, and compare this value with θpeak of its neighboring whiskers, each stimulated by only one muscle. When the muscle forces are weak (Fig. 11E, same parameter values as in Fig. 6F), spatial summation is linear. Therefore, since the effect of the anterior muscle (to which this whisker is posterior) increases with θ0 and the effect of the posterior whisker decreases with θ0, θpeak depends on θ0 only weakly. The peak protraction angle θpeak is almost independent of θ0 for large A values, as shown in Figure 11F for A = 10, despite the fact that the angles θpeak of the anterior, and (especially) the posterior whiskers to these two contracting muscles depend strongly on θ0. We conclude that when several intrinsic muscles contract simultaneously and with similar forces, the dependence of the maximal protraction angle θpeak on θ0 for whiskers that are stimulated by two muscles is weak for both small and large muscle forces.
Discussion
We characterized the transformation from MN spikes to vibrissae movements. The temporal summation of consecutive MN spikes is nonlinear if the ISI is on the order of 10 ms. It is sublinear for briefer ISIs, supralinear for moderate ISIs, and linear for more prolonged ISIs. We explain these properties by an interplay between two processes that affect force generation in muscle cells: one process that causes saturation for brief ISIs, and one that enhances muscle responses for moderate ISIs. Stimulation of a single intrinsic muscle leads to movements of its adjacent whiskers. The amplitude of the posterior whisker increases with the resting angle θ0, and the amplitude of the anterior whisker decreases with θ0. For small muscle forces and θ0 less than ∼90°, the anterior whisker protracts more than the posterior one. When a whisker is pulled by its two adjacent muscles with similar forces, the protraction amplitude depends only weakly on the resting angle for both small and large muscle forces. Counterintuitively, stimulation of intrinsic muscles may cause whisker retraction: the posterior whisker retracts for small θ0, and the anterior whisker retracts for large θ0. When several muscles contract, the spatial summation of whisker movement is linear for small muscle forces and supralinear for larger forces.
Supralinear temporal and spatial summation
Nonlinear temporal summation of consecutive contractions of skeletal muscles is well known, and can be either supralinear or sublinear (for examples, see Cooper, 1930; Burke, 1967; Duchateau and Hainaut, 1986; Rhoades and Pflanzer, 2003). Previously, linearity of summation was quantified as the ratio between the peak amplitude in response to the second stimulus and the corresponding computed peak amplitude based on response to one stimulus and assuming linearity (Herfst and Brecht, 2008). We propose to use the linearity index LI (Eq. 11) for the following reasons. First, LI is >1 if, and only if, the response is supralinear. Second, for the vibrissa system and small forces, LI depends only on muscle contraction, not on the mechanics of the pad. Therefore, this linearity index can be used to describe physiological processes in the muscle. Using this measure, we found experimentally that supralinearity increases with Δt for small interpulse intervals Δt, and reaches a peak at approximately Δt = 17 ms. This value of Δt is small with respect to that for skeletal muscles, and may reflect the fast contracting properties of the intrinsic whisker muscles (Jin et al., 2004). The dependency of LI on Δt (Figs. 3⇑–5) can be explained by two competing processes: saturation of [Ca2+]i at low ISIs, and based on Bobet and Stein (1998), nonlinear dependency of the force on [Ca2+]i. This dependency can be enhanced by other types of nonlinearities. For example, Ca2+-induced Ca2+-release channels on the SR may generate supralinear responses to a fast trains of action potentials (Otazu et al., 2001).
The small protraction angle approximation
Our model, like other models of its type (for examples, see Bobet and Stein, 1998; Ding et al., 1998; Hill et al., 2008), is based on a linear representation of the viscoelastic properties of the tissue by springs and dampers, and therefore, is strictly valid only for small θ. We performed two sets of experiments to cover the two regimes of small and moderate θ: single MN stimulations yielded θ on the order of a few degrees, and facial nerve stimulations yielded moderate protraction angles (∼10°–20°). As in models of other systems (e.g., Golomb and Hansel, 2000), our modeling work yields testable predictions that fit experimental observation even for moderate θ values (Fig. 5).
For small muscle forces and small θ (on the order of a few degrees), we find that the anterior whisker protracts more than the posterior one for small θ0 (Eq. 15). Beyond these angles, we find other dynamical effects. When θ is already >20°, the additional θ generated by muscle force can be large in comparison with the corresponding additional θ that generated when θ is small (supplemental Fig. S1, available at www.jneurosci.org as supplemental material). However, this effect is compensated for by the reduction of muscle force due to the force–length curve (Fig. 5). Furthermore, for large θ, protraction of the posterior whisker increases relative to the anterior one (Fig. 9).
Mechanical implications of the model
Although the structure of our model builds on and is similar to that of Hill et al. (2008), differences between the models have significant consequences. For example, the whiskers in our model are coupled to the pad, not to each other. This generates the concept of a WMU, and enables two whiskers to move while others are stationary, which is consistent with experimental observations. Using our model, we show that the muscle attachment position, ld, and the elastic properties of the tissue in the lateral–medial direction, kπ,y and ζπ,y, have a strong impact on whisker movements. Interestingly, the dynamical properties of the whisker (the center of mass, C, the mass, M, and the moment of inertia, I) have only minor effects on whisker trajectory (Fig. 8). This implies that the system is not sensitive to the location of a “virtual” pivot along the follicle.
Model predictions and experimental data
Unlike existing models of sensation-targeted motor control (such as oculomotor control), which describe control of eye movements at levels of details that are no finer than population firing rates, the model presented here allows us to compare predictions related to single motor spikes with experimental data. Our model predicts supralinear temporal summation of interspike intervals with Δt < 30 ms. Except for a single example, this prediction is in good agreement with experimental data (Herfst and Brecht, 2008) (Fig. 4C). In that example, an obscuring whisker oscillation appeared to be superimposed on the whisker trajectory. Filtering out the superimposed oscillation uncovered the same supralinear summation that characterized the other examples.
Our findings suggest that when an intrinsic muscle contracts, its two adjacent whiskers usually move. While this is true for all resting angles, except θ0,pr (see Fig. 6F), this could be hard to measure, especially when R is significantly different from 1. So far, whenever we tracked a moving whisker and both its adjacent whiskers, at least one of the adjacent whiskers moved as well (see Fig. 10).
We predict that in response to contraction of an intrinsic muscle, the posterior whisker retracts if θ0 < θ0,pr (∼60°) and the anterior whisker retracts if θ0 > θ0,ar (∼120°). In trained, head-restrained rats, two adjacent whiskers could move simultaneously in opposite directions: the anterior whisker protracted while the posterior whisker retracted (Sachdev et al., 2002), albeit both movements were small (<5°). The angles θ0 in that study were 33° < θ0 < 76° [see Sachdev et al. (2002), their Fig. 5], which is consistent with divergent whisker movement being generated by contraction of one intrinsic muscle.
Another prediction of our model involves the most anterior whisker in a row. Since this whisker is attached to only one intrinsic muscle, the model predicts that if only intrinsic muscles are activated, and with similar forces, the most anterior whisker will protract less than the others in the row. However, extrinsic muscles may compensate for this effect, and allow larger protraction amplitudes of those anterior whiskers. In fact, protracting extrinsic muscles attached to the most anterior whiskers in rows A and B were recently found (S. Haidarliu, personal communication). Herein, we treat θ0 as a constant parameter. Tonic activation of extrinsic muscles may cause movement of the whole mystacial pad, and therefore, modify θ0, which in turn would affect the whisker protraction response to the phasic activation of intrinsic muscles.
The basic unit for whisking control
Our model shows that the vibrissal system cannot exclusively control single whiskers; the minimal unit for control is a pair of whiskers that are connected by an intrinsic muscle, which we call a WMU. During free whisking in the absence of object contact, many intrinsic muscles are activated in a coordinated manner and many whiskers move simultaneously in a coherent manner (Towal and Hartmann, 2006; Mitchinson et al., 2007). Fractionated movements of an individual or a small number of whiskers are also occasionally seen in behaving rats, particularly during object contact (Knutsen et al., 2006; Mehta et al., 2007). We show that interactions between adjacent WMUs are linear for small angles (<∼4°), and become nonlinear with large angles (Fig. 11, Appendix B). Control of individual WMUs can be considered independent only up to the spatial summation that occurs at linked whiskers. Specifically, with small angles, the motion equation for a given whisker row can be obtained by superposition of the equations of all the WMUs that compose that row.
Inverse model of whisking
The strong forward correlation between MN spikes and whisker trajectories indicates that the system is deterministic to a large extent, with every MN spike affecting whisker motion (Simony et al., 2008). Therefore, the detailed model presented here can, in principle, be inversed to compute putative patterns of MNs spikes that could generate specific high-resolution trajectories of whisker movements, which are observed experimentally. Such an inverse model can also be used as a first approximation for a neuronal inverse model (Kawato et al., 1987; Wolpert and Kawato, 1998) used by the rat vibrissal system for whisking control. A spike-based inverse model derived from the model presented here can complement existing rate-based descriptions of inverse models (e.g., Ghasia et al., 2008; Lisberger, 2009).
Our proposed model can serve as a basic building block in models describing one or more motor–sensory closed loops of the vibrissal system. Its direct translation from MN spike patterns to whisker trajectory allows efficient analysis of both closed loops controlling whisking and those controlling interactions with objects.
Appendix A: Model of the Biomechanical Plant for Rat Whiskers
Our model of the biomechanical plant is based on that of Hill et al. (2008). In our model, there are NW whiskers in a row, and these whiskers are considered to be rigid bodies (Fig. 1E). The x-axis of the system is defined as being in the posterior–anterior direction. The y-axis is defined as being in the medial–lateral direction, from the plate to the skin, with y = 0 representing the line of the center of masses. The z-axis is perpendicular to the x and y axes. The length of a row in the mystacial pad is w, and the distance between the centers of mass of two adjacent whiskers is s. The follicle length is lf, and the position of the center of mass along the whisker is C (C < 0). The mass of the whisker (follicle + hair) is M, and the moment of inertia with respect to the z-axis is I. For simplicity, we consider these parameters to be uniform for all the whiskers. The position X⃗icm = (Xicm, Yicm) of the center of mass of the ith whisker at rest, i.e., when all the muscles are relaxed and at steady state, is When a whisker moves, the position of its center of mass x⃗icm(t) = (xicm(t), yicm(t)) is x⃗icm(t) = X⃗icm + x⃗i(t), where x⃗i(t) = (xi(t), yi(t)) is a vector of the deviations of the centers of mass from their resting values. The resting angle of the whisker with respect to the posterior side of the posterior–anterior direction is θ0i, the protraction angle is θi(t), and the angle between the whisker and the posterior–anterior direction is θi(t) = θ0i + θi(t). For simplicity, and in agreement with Hill et al. (2008), we assume that all the resting angles are equal: θ0i = θ0. The state of the ith whisker, i = 1 … NW, is defined by six dynamical variables: xi, ẋi, yi, ẏi, θi, and θ̇i. The equations of motion are as follows: where F⃗total,is, F⃗total,id, and F⃗total,im are the total forces applied on the ith whisker by the springs (s), dampers (d), and muscles (m), respectively, and Nz,is, Nz,id, and Nz,im are the total torques applied on that whisker by the springs, dampers, and muscles, respectively. A minus sign appears in Equation A3 because the angle θi is defined as an angle from the negative x-axis. For clarity, a list of the superscripts defining the parameters and variables in the model is given in Table A1, and a list of the reference model parameters is given in Table A2.
Follicle edges and attachment points of springs and dampers
The edge of the follicle located at the skin is denoted by “σ,” and the edge of the follicle located at the plate is denoted by “π.” The coordinates of these points are as follows, for 1 ≤ i ≤ NW: A whisker is connected to tissue by springs and dampers. In the skin plane, the “anchor points” of the springs and dampers are assumed to be spread horizontally between two adjacent whiskers. Similar anchor points are assumed to be in the plate plane. Moreover, there is another anchor point for each whisker below the bottom of the follicle, toward the bone. Each anchor point is denoted by two superscripts. The first superscript is “σ” for skin or “π” for plate. The second superscript is “p” for a posterior anchor point, “a” for an anterior anchor point, or “b” for the bone. The coordinates of the anchor points are as follows:
Forces and torques of springs
The absolute value of the force applied by each spring is where k is the spring constant, and for the spring l⃗ is the vector length, |l⃗| is the length, and L0 is the length for which the spring does not produce any force. For simplicity, we assume that L0 is equal to the length of the spring at rest (without any muscle forces or movements). We also assume that the spring is always linear, and that the direction of the spring force is along the spring. The vector length of a spring, l⃗iα,n, is defined as being the distance from the spring's anchor point to the contact point of the spring with the follicle: where the index α denotes the plane of the contact point (σ for skin and π for plate), and the index n denotes the direction of the spring (p, a, and b for posterior, anterior, and bone, respectively). From Equations A11 and A12, the force applied by each spring is We assume that the horizontal spring constants are independent of the direction (posterior or anterior), and define the spring constant in the x-direction as follows: kσ,x ≡ kσ,a = kσ,p and kπ,x ≡ kπ,a = kπ,p. The spring constant in the y-direction, kπ,y ≡ kπ,b, can be different from kπ,x. The force applied on the ith whisker by all the springs at a certain point α = σ (skin) or π (plate) is as follows: The total force applied by the springs on the ith whisker is The torque applied by the springs on the whisker is Substituting Equations A4 and A5 into Equation A16, we obtain
Forces and torques of dampers
The absolute value of the force applied by each damper is From Equations A18 and A12, the force applied by each damper is As with the spring constant, we assume that the horizontal damping constants are independent of the direction (posterior or anterior), and define ζσ,x ≡ ζσ,a = ζσ,p, ζπ,x ≡ ζπ,a = ζπ,p, and ζπ,y ≡ ζπ,b. The parameters k and ζ (Table A2) were chosen to ensure over-damped dynamics (Hill et al., 2008). The force applied on the ith whisker by all the dampers at a certain point is The total force applied by the dampers on the ith whisker is Using Equations A4 and A5, the torque applied by the dampers on the whisker is
Forces and torques of muscles
The intrinsic muscles are small muscles that associate solely with vibrissal follicles, without any bony attachment (Dörfl, 1982). Each intrinsic muscle, except for the most posterior, connects two adjacent follicles of the same row. These muscles form a “sling” that surrounds the lower third of the anterior follicle and connects to the superior part of the posterior follicle near the skin. The follicles within a row, except for the most anterior, e.g., B4, C7, are connected to two intrinsic muscles, and the most anterior follicles are connected by only one muscle. Each of the most posterior follicles (and, particularly, the straddlers: α, β, γ, and δ) is connected by one or two intrinsic muscles to one or two more anterior whiskers, and by another intrinsic muscle, the edges of which blend with the fibers of extrinsic muscles (m. levator labii superioris and m. maxillolabialis).
Here, we present a model of a row of whiskers with intrinsic muscles being the only muscles considered. Each intrinsic muscle, except the most posterior (i = 0), connects two whiskers. The ith muscle, 1 ≤ i ≤ NW − 1, is attached to the more posterior follicle (i) at the skin plane, and it wraps around the more anterior follicle (i + 1) at a distance ld along the follicle from the skin, with 0 < ld ≤ lf. The attachment point of an intrinsic muscle to its anterior follicle is denoted by “μ.” Since our model includes only intrinsic muscles, we do not incorporate the coupling of the straddlers to the extrinsic muscle. Instead, based on Hill et al. (2008), we use an intrinsic muscle that connects the most posterior whisker to an anchor point x⃗0μ on the pad, as follows: The most anterior follicle (i = NW) is coupled only to its posterior neighbor (open boundary condition).
The forces developed by the muscles are The total force applied by the muscles on the ith whisker, 1 ≤ i ≤ NW − 1, is For the most anterior whisker, F⃗total,NWm = −F⃗NW−1m. The total torque applied by the muscles on the ith whisker, 1 ≤ i ≤ NW − 1, is By substituting Equations A4 and A23 in Equation A28, we obtain for that whisker (1 ≤ i ≤ NW − 1) For the most anterior whisker,
Appendix B: For Small Protraction Angles, the Linearity Index Depends Only on the Muscle Activity
Assume the 6NW-dimensional state vector V⃗ includes all the dynamical variables in the system where ( … )T denotes a column vector. Assume also that the NW-dimensional vector f⃗ includes the absolute value of the forces produced by all the muscles Using Equations A25 and A26, the equations of motion, Equations A2 and A3, can be written in a vector form (see Strogatz, 1994), as follows: The vector function Φ⃗ includes the autonomous terms in Equations A2 and A3, which do not depend on the forces of the muscles. The matrix M is a 6NW × NW matrix that represents the dependence of the direction of the vector forces on the state of the vibrissae. At rest, V⃗ = 0, f⃗ = 0, and Φ⃗ = 0. For small protraction angles, all components of the vector V⃗ are small, and therefore, based on continuity arguments, Φ⃗(V⃗) is also small. For V⃗ to be small, f⃗(t) should be small, because the matrix M(V⃗) is not small for V⃗ = 0. Therefore, we linearize the right-hand side of Equation B3 with respect to V⃗ by keeping the first-order terms in the expansion of Φ⃗, but only the zeroth-order term in the expansion of M. Deviations of M(V⃗) from M(0) introduce only high-order terms. For small V⃗, Equation B3 becomes where the elements of the 6NW × 6NW matrix A are Without any muscle forces, the rest state is stable because of damping. Therefore, all the eigenvalues of the matrix A are negative. Suppose that the system is at rest at t = 0, i.e., V⃗(t = 0) = 0. Then, the solution of Equation B4 is The linearity index LI,ij (Eq. 11) is calculated from the protraction profiles θij(t) and θi1(t) of the ith whisker in response to j spikes and one spike, respectively. Thus, we compute the integral Jij of the function θij(t) over t, as follows: where f⃗j denotes the forces developed by the muscles in response to j spikes. The vector u⃗iT is used to select the component of V⃗ representing θi. Its kth component ui,kT equals 1 if k = 6(i − 1) + 5 and 0 otherwise. After changing the order of integration and making the change of variables from t to q = t − t′, we obtain where B = ∫0∞ dq exp(Aq)M(0) is a 6NW × NW matrix. The matrix B is finite because all the eigenvalues of A are negative.
Suppose that only one muscle, with an index k, contracts when it receives j spikes from its MN, and develops force Fkm,j. The resulting Jij is Therefore, Thus, the linearity index depends only on the biophysics of the muscle, not on the mechanical properties of the mystacial pad. This result is general, since it does not depend on the details of the biomechanical model of the motor plant.
For our model of the motor plant (Fig. 1E), the ith whisker (except the most anterior one) is affected by two intrinsic muscles, i − 1 and i. Hence, and Equation B10 is valid if the two muscles adjacent to the whisker undergo the same temporal response to spikes. Similarly, Equation B10 always applies for the most anterior whisker, which is affected by only one muscle, k = NW − 1.
Appendix C: Steady-State Whisker Mechanics for Small Protraction Angles
Here, we consider a system of two whiskers: i = 1 (posterior) and i = 2 (anterior), coupled by a muscle (Fig. C1). For this system, we compute the protraction angle during steady state (θ̇i = 0), when the muscle applies a small, time-independent force Fm. The angle ϕ denotes the angle between the muscle and the posterior direction of the x-axis. From Figure A1, In steady state, deviations of the edges of follicles from their resting value for Fm = 0 are denoted by xiσ and xiπ, where i = {posterior,anterior}. The calculation is performed for the following conditions: (1) −C = ld = lf; (2) kπ,y ≫ kπ,x, kσ,x, and therefore, movements in the y-direction are negligible; and (3) θi ≪ 1 for i = 1, 2.
At steady state, Equations A2 and A3 become Since, from Figure C1, for small θi, The muscle force on the anterior whisker is applied at the plate. We seek a solution with xanteriorσ = 0, which entails Fpoint,anteriors,σ = 0. Since ld = lf (the center of mass of the whisker is located in the plate), there are no muscle or spring torques, and Equation C3 is satisfied for the anterior whisker (i = 2). Thus, for the anterior whisker, Equation C2 becomes therefore, For the posterior whisker, the force balance equation, Equation C2, is Hence, The torque balance equation, Equation C3, for the posterior whisker is Using Equations C8 and C10, we obtain Thus, the bottom point of the posterior follicle (near the plate), xposteriorπ, moves to the right (anteriorly) if θ0 < π/2 and to the left (posteriorly) if θ0 > π/2. The movement of xposteriorπ to the right causes the protraction angle θposterior to be small or even negative for small θ0, as we show below.
Substituting Equations C7, C9, and C11 in Equation C5 yields Using Equations 12, C1, C12, and C13, we obtain Therefore, R decreases with θ0 and with the ratio kπ,x/kσ,x. The value R = 1 is obtained for θ0 = θ0,eq, where The angle θposterior is negative, i.e., the posterior whisker retracts, if Similarly, θanterior is negative, i.e., the anterior whisker retracts, if If kσ,x = kπ,x, then θ0,eq = π/2. For our reference parameter values, s = 2 mm and lf = 4 mm, we obtain θ0,pr = 60° and θ0,ar = 120°.
Footnotes
The research was supported by European Union Grant BIOTACT (ICT-215910) to M.B., E.A., and D.G., by Grant 2007121 from the United States–Israel Binational Science Foundation to E.A. and D.G., by Israeli Science Foundation Grants 311/04 to D.G. and 959/06 to E.A., by the KAMEA program, Ministry of Absorption, Israel, to K.B., and by the Minerva Foundation funded by the Federal German Ministry for Education and Research to E.A. We thank Sebastian Haidarliu for helpful discussions, Naama Rubin for animation programming, Dan Hill and David Kleinfeld for carefully reading the manuscript, and Barbara Schick for editing the manuscript.
- Correspondence should be addressed to either of the following at the above addresses: David Golomb, golomb{at}bgu.ac.il; or Erez Simony, erez.simony{at}weizmann.ac.il