## 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 Ca^{2+} 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. 1*A*,*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, Δ*t*s, 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. 1*C*). 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. 1*A*), 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. 1*B*). 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 *t _{k}*,

*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. 1

*D*).

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 Ca

^{2+}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 Ca^{2+}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*), Ca^{2+}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 Ca^{2+}concentration, [Ca^{2+}]_{i}, and to have a time scale τ_{c}. Therefore, [Ca^{2+}]_{i}evolves according to the following equation: where*r*_{0}is a constant. The solution of Equation 2 with the initial condition [Ca^{2+}]_{i}(*t*=*t*) =_{k}*C*is where_{k}*C*= 0 for_{k}*k*= 1 and for*k*≥ 1, where Δ*t*=_{k}*t*_{k}_{+1}−*t*._{k}The muscle force (

*F*) depends on the level of [Ca^{2+}]_{i}(Rhoades and Pflanzer, 2003), and this dependency is denoted by*F*([Ca_{C}^{2+}]_{i}). The function*F*([Ca_{C}^{2+}]_{i}) is nonlinear because several (*n*) Ca^{2+}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*A*_{0}= 1 mg · mm/ms^{2}. 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),

*F*(_{L}*z*), where*l*^{m}(*t*) is the instantaneous muscle length, and*l*_{0}^{m}is the muscle length at rest. Following Hill et al. (2008),*F*(_{L}*z*) is modeled as where 0 <*z*_{h}<*z*_{l}< 1. Note that*F*(_{L}*z*) = 1 for small protraction angles (Fig. 1*D*). Following Hill et al. (2008), we ignore the velocity dependence of the muscle force.The total muscle force

*F*^{m}(*t*) generated by consecutive MN spikes is a product of*F*([Ca_{C}^{2+}]_{i}) and*F*(_{L}*z*): We use the following parameter set, referred to as the “reference parameter set,” throughout the paper unless otherwise stated:*r*_{0}= 1.9, τ_{r}= 5 ms, τ_{c}= 6 ms,*z*_{h}= 0.1,*z*_{l}= 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 1*E*, 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

*l*, where_{d}*l*is a parameter of the system, and not necessarily equal to −2_{d}*C*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 *N*_{W} 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*) = (

*x*(

_{i}*t*),

*y*(

_{i}*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. 1

*C*). 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

*i*th whisker. Each whisker is affected by the total force (

*F⃗*), and the torque in the

*z*-direction (

*N*

_{z}) 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 *i*th whisker is θ_{i}^{j}(*t*). Then the linearity index of the response to *j* spikes, *L*_{I}^{j}, is defined as
The integration times, *T _{j}*, 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

*T*= 200 ms. For analysis of data from vMN stimulations, we use

_{j}*T*= (

_{j}*j*− 1)Δ

*t +*40 ms. The linearity index

*L*

_{I}

^{j}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 θ_{i}^{max} and the minimal angle, i.e., the retraction angle with the maximal absolute value, is denoted by θ_{i}^{min}. The “peak amplitude” θ_{i}^{peak} is defined to be θ_{i}^{max} if the whisker does not retract at all or if θ_{i}^{max} > |θ_{i}^{min}|, and is defined to be θ_{i}^{min} 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, θ_{anterior}^{peak} and θ_{posterior}^{peak}, 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. 1*E* and Appendix A), and used this coordinate to quantify pad movement. Sample trajectories of θ(*t*) and *x*^{σ}(*t*) are presented in Figure 2*A*. 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. 7*G*). 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. 2*B*). 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. 2*B*) (*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. 3*A*). This supralinear summation is even more apparent in the response to three consecutive pulses (Fig. 3*B*).

We quantified the temporal summation by calculating the linearity indices *L*_{I}^{2} and *L*_{I}^{3} (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 *L*_{I}^{2} and *L*_{I}^{3} versus Δ*t* are plotted in Figure 3*C*; the corresponding values for the single whisker (see Fig. 3*A*,*B*) are plotted in the inset. For Δ*t* = 1 ms, both *L*_{I}^{2} and *L*_{I}^{3} 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. *L*_{I}^{2} and *L*_{I}^{3} remain larger than 1 below Δ*t*_{L}, which is the minimal interval for linear summation (∼35 ms). This reflects supralinear temporal summation. Above Δ*t*_{L}, *L*_{I}^{2} and *L*_{I}^{3} are essentially 1, and the temporal summation is linear. The two linearity indices (*L*_{I}^{2} and *L*_{I}^{3}) increase with Δ*t* for 1 ms < Δ*t* < Δ*t*_{M} ∼ 17 ms, and decrease for Δ*t* > Δ*t*_{M}. Interestingly, Δ*t*_{M} is similar to the Δ*t* value for which protraction in response to a single twitch is maximal (Fig. 3*A*). The maximal values of *L*_{I}^{3} over Δ*t* are plotted in Figure 3*D* as a function of the corresponding maximal values of *L*_{I}^{2} for the nine whiskers. The *L*_{I}^{3} increases with *L*_{I}^{3} (linear regression, slope = 1.005, *R*^{2} = 0.6652), and the *L*_{I}^{3} is slightly, but not significantly, larger than *L*_{I}^{2} (*p* = 0.0547, Wilcoxon test).

#### 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. 4*A*,*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 *L*_{I}^{2} as a function of Δ*t* for four motoneurons (Fig. 4*C*) 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. 6*C*,*D*] (red symbols in Fig. 4*C* 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. 4*C*) 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 *L*_{I}^{j}(Δ*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 *F*^{m,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 *F*^{m,j}(*t*), the linearity index *L*_{I}^{j} depends only on the muscle properties, not on the properties of the biomechanical motor plant, i.e.,
Therefore, the dependence of *L*_{I}^{j} on Δ*t* can be accounted for solely by the muscle model (Eqs. 1–8).

We calculate, using Equation 11, the protraction angle and *L*_{I}^{2} 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. 4*A*) and to two spikes (Fig. 4*B*) 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 *L*_{I}^{2} on Δ*t* is shown in Figure 4*C*. The model responses fit the experimental recordings of whisker responses (Fig. 4) very well. We use the model and the fact that *L*_{I}^{j} 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* = *l*^{m}(*t*)/*l*_{0}^{m(t)} (Eq. 6) obeys |*z* − 1| < *z*_{h} (see Eq. 7), and hence the force–length function *F _{L}*(

*z*) equals 1. Therefore,

*F*(

^{m}*t*) =

*F*([Ca

_{C}^{2+}]

_{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 *F _{C}*([Ca

^{2+}]

_{i}) on [Ca

^{2+}]

_{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

*L*

_{I}

^{j}= 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

*F*([Ca

_{C}^{2+}]

_{i}) on [Ca

^{2+}]

_{i}leads to supralinear temporal summation. Those two factors cause the indices

*L*

_{I}

^{j}to increase with Δ

*t*, until they reach their maximal values (well above 1) at approximately Δ

*t*

_{M}. Above this time interval (Δ

*t*

_{M}), as Δ

*t*increases much beyond the sequestration time constant τ

_{c}, the level of [Ca

^{2+}]

_{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

*L*

_{I}

^{j}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, *F _{L}*(

*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

*L*

_{I}

^{j}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,

*F*(

_{L}*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. 6*A*). 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 *l _{d}* = 3 mm. An example of plant morphology in which whiskers point posteriorly during rest with resting angle θ

_{0}= 70° is depicted (Fig. 6

*A*). Time traces of the protraction angles θ (Fig. 6

*B*) and the coordinates of the whisker center of mass (

*x*,

*y*) in response to a single small muscle twitch (Fig. 6

*C*) 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. 6

*C*, 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. 6

*C*, bottom panel). In the other example of plant morphology presented here, the whiskers point anteriorly during rest with resting angle θ

_{0}= 95° (Fig. 6

*D*). Time traces of the protraction angle θ for this example are shown (Fig. 6

*E*). 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 6*A–E* demonstrate general behavior observed in all our simulations: the peak amplitude of the posterior whisker, θ_{posterior}^{peak}, increases with θ_{0}, whereas the peak protraction amplitude of the anterior whisker θ_{anterior}^{peak}, decreases with θ_{0} (Fig. 6*F*). 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. 6*G*). When θ_{0} > θ_{0} = θ_{0,ar} (∼108° for the parameters of Fig. 6), the anterior whisker retracts, whereas the posterior whisker protracts (Fig. 6*H*). The sudden jumps in the curves of the peak amplitudes of the posterior and anterior whisker, for θ_{0} = 50° and θ_{0} = 108°, respectively (Fig. 6*F*), 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* = θ_{anterior}^{peak}/θ_{posterior}^{peak} (Eq. 12) decreases with θ_{0} for all θ_{0} values except for angle θ_{0,pr}, where *R* switches from −∞ to ∞ as θ_{posterior}^{peak} switches from negative to positive values (Fig. 7*A*, inset).

*l _{d}* is the distance measured along the follicle from the skin to the attachment point where the muscle cell slings around the anterior whisker (Fig. 1

*E*). If

*l*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

_{d}*R*should be larger. As long as the posterior whisker protracts, this expectation is confirmed in our simulations (Fig. 7

*A*), where

*R*decreases with θ

_{0}and increases with

*l*. In addition, increasing

_{d}*l*increases θ

_{d}_{0,pr}, and θ

_{0,ar}, and therefore, increases the θ

_{0}range in which the posterior whisker retracts (Fig. 7

*A*, 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. 7*B*). 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. 7*C*), 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 θ_{posterior}^{max} and θ_{anterior}^{max} depend only weakly on the center of mass *C* (Fig. 8*A*). 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. 8*B*,*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. 6*A*,*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. 8*A*). 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* = *l _{f}*, where

*l*is the length of the follicle. Fourth, since

_{f}*R*increases with

*l*, examine the most extreme case for which

_{d}*l*=

_{d}*l*. Fifth, look at the steady-state condition, by assuming that the muscle force,

_{f}*F*

^{m}, 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}=

*k*

^{x}, 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*/

*l*) (see also Eq. C16); θ

_{f}_{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*/

*l*) (see also Eq. C17); θ

_{f}_{0,ar}= 120° for our reference parameter set. These results are similar to those obtained from simulations of the full system dynamics (Fig. 6

*F*).

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. 6*A*), the muscle force and torque applied on the posterior whisker cause the bottom point of the posterior follicle to move anteriorly (*x*_{posterior}^{π} > 0) (Eq. C11). The more the point *x*_{posterior}^{π} moves anteriorly, the less the point where the whisker emerges from the skin, *x*_{posterior}^{σ}, 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. 6*D*), the point *x*_{posterior}^{π} moves posteriorly (*x*_{posterior}^{π} < 0) (Eq. C11), the point *x*_{posterior}^{σ} moves anteriorly, and the posterior whisker always protracts. The point where the anterior whisker emerges from the skin, *x*_{anterior}^{σ}, does not move because no net force is applied on it, and no torque is applied on the anterior whisker (since *l _{d}* = −

*C*=

*l*). The force equation for

_{f}*x*

_{anterior}

^{π}(Eqs. C6, C7) shows that the coordinate

*x*

_{anterior}

^{π}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. 7*C*). 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 *l _{d}* =

*l*. We can explain, however, why the protraction of the anterior whisker is larger in comparison to that of the posterior whisker as

_{f}*l*increases. This stems from the fact that increasing

_{d}*l*increases the torque on the anterior whisker in the direction that causes it to protract.

_{d}#### 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. 9*A*). The temporal profile of the anterior whisker looks similar to that obtained with smaller *A* (*A* = 1.33) (Fig. 6*B*), 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 θ_{posterior}^{max} is larger than that of the anterior whisker θ_{anterior}^{max}. The angles θ_{anterior}^{max} and θ_{posterior}^{max} and the maximal retraction angles of the anterior and posterior whiskers, θ_{anterior}^{min} and θ_{posterior}^{min}, are shown as functions of *A* (Fig. 9*B*). Beyond the small angle regime (say, of θ < 3°), the increase of the angle θ_{anterior}^{max} with *A* becomes more moderate, until θ_{anterior}^{max} saturates for large *A*. The angle θ_{anterior}^{min} remains 0, i.e., the anterior whisker does not show any retraction. In contrast, the increase of θ_{posterior}^{max} with *A* is accelerated for large values of *A*, and θ_{posterior}^{max} eventually exceeds θ_{anterior}^{max}. Furthermore, if *A* is not small, the posterior whisker exhibits a retraction phase after the protraction phase. For θ_{0} = 95°, the dominance of θ_{posterior}^{max} over θ_{anterior}^{max}, which is already observed for small angles (Fig. 6*E*), also extends for larger angles (*A* = 40) (Fig. 9*C*). In this case, no prominent retraction of any whisker is observed. The dependencies of θ_{anterior}^{max} and θ_{posterior}^{max} on *A* are shown (Fig. 9*D*). 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. 10*A*). 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. 10*B*). 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. 11*A*), 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* = *t*_{a} − *t*_{p}, where *t*_{a} and *t*_{p} 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,muscle}^{max}(*A*, Δ*t*). θ_{ant,muscle}(*t*; *A*, Δ*t*) and θ_{ant,muscle}^{max}(*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,muscle}^{max}(*A*, Δ*t*). In this analysis, we compare four quantities: θ_{post,muscle}^{max}(*A*, 0); θ_{ant,muscle}^{max}(*A*, 0); the maximum of the linear sum of the whisker responses, _{post,muscle}(*t*; *A*, 0) + θ_{ant,muscle}(*t*; *A*, Δ*t*)); and θ_{post+ant,muscle}^{max}(*A*, Δ*t*).

The four functions defined above are plotted in Figure 11*B* 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., _{post,muscle}(*t*; *A*, 0) + θ_{ant,muscle}(*t*; *A*, 0)) < 3°. This result is expected because the biomechanical model of the plant can be linearized for small protraction angles, with the vector of muscle forces being an inhomogeneous term in a linear equation (Appendix B). For larger protraction angles, the response to simultaneous innervations of the two muscles is larger than the sum of the responses to each muscle separately, i.e., θ_{post+ant,muscle}^{max}(*A*, 0) > _{post,muscle}(*t*; *A*, 0) + θ_{ant,muscle}(*t*; *A*, *0*)). Interestingly, even if in response to the contraction of the anterior muscle only, the whisker protracts and then retracts (Fig. 9*A*,*B*), stimulating the anterior muscle together with simultaneous stimulation of the posterior muscle enhances the protraction angle supralinearly without any whisker retraction. To examine the spatiotemporal summation properties of the system, we compute the four functions versus Δ*t* for *A* = 10 (Fig. 11*C*). The summation is supralinear for small values of |Δ*t*| (−25 ms ≤ Δ*t* ≤ 19 ms) and is linear for larger |Δ*t*| values. Interestingly, maximal nonlinearity is obtained for positive Δ*t* (Δ*t* = 4 ms for the parameters of Fig. 11*C*), i.e., when the posterior muscle is stimulated before the anterior one.

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. 11*D*), 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. 11*E*, same parameter values as in Fig. 6*F*), 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 11*F* 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 *L*_{I} (Eq. 11) for the following reasons. First, *L*_{I} is >1 if, and only if, the response is supralinear. Second, for the vibrissa system and small forces, *L*_{I} 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 *L*_{I} on Δ*t* (Figs. 3⇑–5) can be explained by two competing processes: saturation of [Ca^{2+}]_{i} at low ISIs, and based on Bobet and Stein (1998), nonlinear dependency of the force on [Ca^{2+}]_{i}. This dependency can be enhanced by other types of nonlinearities. For example, Ca^{2+}-induced Ca^{2+}-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, *l _{d}*, 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. 4*C*). 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. 6*F*), 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 *N*_{W} whiskers in a row, and these whiskers are considered to be rigid bodies (Fig. 1*E*). 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 *l _{f}*, 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⃗*

_{i}

^{cm}= (

*X*

_{i}

^{cm},

*Y*

_{i}

^{cm}) of the center of mass of the

*i*th 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⃗*

_{i}

^{cm}(

*t*) = (

*x*

_{i}

^{cm}(

*t*),

*y*

_{i}

^{cm}(

*t*)) is

*x⃗*

_{i}

^{cm}(

*t*) =

*X⃗*

_{i}

^{cm}+

*x⃗*(

_{i}*t*), where

*x⃗*(

_{i}*t*) = (

*x*(

_{i}*t*),

*y*(

_{i}*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

*i*th whisker,

*i*= 1 …

*N*

_{W}, is defined by six dynamical variables:

*x*,

_{i}*ẋ*,

_{i}*y*,

_{i}*ẏ*, θ

_{i}_{i}, and θ̇

*. The equations of motion are as follows: where*

_{i}*F⃗*

_{total,i}

^{s},

*F⃗*

_{total,i}

^{d}, and

*F⃗*

_{total,i}

^{m}are the total forces applied on the

*i*th whisker by the springs (s), dampers (d), and muscles (m), respectively, and

*N*

_{z,i}

^{s},

*N*

_{z,i}

^{d}, and

*N*

_{z,i}

^{m}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 θ

*is defined as an angle from the negative*

_{i}*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* ≤ *N*_{W}:
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 *L*_{0} is the length for which the spring does not produce any force. For simplicity, we assume that *L*_{0} 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 *i*th whisker by all the springs at a certain point α = σ (skin) or π (plate) is as follows:
The total force applied by the springs on the *i*th 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 *i*th whisker by all the dampers at a certain point is
The total force applied by the dampers on the *i*th 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 *i*th muscle, 1 ≤ *i* ≤ *N*_{W} − 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 *l _{d}* along the follicle from the skin, with 0 <

*l*≤

_{d}*l*. 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

_{f}*x⃗*

_{0}

^{μ}on the pad, as follows: The most anterior follicle (

*i*=

*N*

_{W}) 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 *i*th whisker, 1 ≤ *i* ≤ *N*_{W} − 1, is
For the most anterior whisker, *F⃗*_{total,NW}^{m} = −*F⃗*_{NW−1}^{m}. The total torque applied by the muscles on the *i*th whisker, 1 ≤ *i* ≤ *N*_{W} − 1, is
By substituting Equations A4 and A23 in Equation A28, we obtain for that whisker (1 ≤ *i* ≤ *N*_{W} − 1)
For the most anterior whisker,

## Appendix B: For Small Protraction Angles, the Linearity Index Depends Only on the Muscle Activity

Assume the 6*N*_{W}-dimensional state vector *V⃗* includes all the dynamical variables in the system
where ( … )* ^{T}* denotes a column vector. Assume also that the

*N*

_{W}-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 6

*N*

_{W}×

*N*

_{W}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 6

*N*

_{W}× 6

*N*

_{W}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

*L*

_{I,i}

^{j}(Eq. 11) is calculated from the protraction profiles θ

_{i}

^{j}(

*t*) and θ

_{i}

^{1}(

*t*) of the

*i*th whisker in response to

*j*spikes and one spike, respectively. Thus, we compute the integral

*J*

_{i}

^{j}of the function θ

_{i}

^{j}(

*t*) over

*t*, as follows: where

*f⃗*denotes the forces developed by the muscles in response to

^{j}*j*spikes. The vector

*u⃗*

_{i}

^{T}is used to select the component of

*V⃗*representing θ

*. Its*

_{i}*k*th component

*u*

_{i,k}

^{T}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 6

*N*

_{W}×

*N*

_{W}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 *F*_{k}^{m,j}. The resulting *J*_{i}^{j} 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. 1*E*), the *i*th 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* = *N*_{W} − 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

*F*

^{m}. 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

*F*

^{m}= 0 are denoted by

*x*

_{i}

^{σ}and

*x*

_{i}

^{π}, where

*i*= {posterior,anterior}. The calculation is performed for the following conditions: (1) −

*C*=

*l*=

_{d}*l*; (2)

_{f}*k*

^{π,y}≫

*k*

^{π,x},

*k*

^{σ,x}, and therefore, movements in the

*y*-direction are negligible; and (3) θ

*≪ 1 for*

_{i}*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

*x*

_{anterior}

^{σ}= 0, which entails

*F*

_{point,anterior}

^{s,σ}= 0. Since

*l*=

_{d}*l*(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 (

_{f}*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),

*x*

_{posterior}

^{π}, moves to the right (anteriorly) if θ

_{0}< π/2 and to the left (posteriorly) if θ

_{0}> π/2. The movement of

*x*

_{posterior}

^{π}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 *l _{f}* = 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