## Abstract

During exploratory behavior, rats brush and tap their whiskers against objects, and the mechanical signals so generated constitute the primary sensory variables upon which these animals base their vibrissotactile perception of the world. To date, however, we lack a general dynamic model of the vibrissa that includes the effects of inertia, damping, and collisions. We simulated vibrissal dynamics to compute the time-varying forces and bending moment at the vibrissa base during both noncontact (free-air) whisking and whisking against an object (collision). Results show the following: (1) during noncontact whisking, mechanical signals contain components at both the whisking frequency and also twice the whisking frequency (the latter could code whisking speed); (2) when rats whisk rhythmically against an object, the intrinsic dynamics of the vibrissa can be as large as many of the mechanical effects of the collision, however, the axial force could still generate responses that reliably indicate collision based on thresholding; and (3) whisking velocity will have only a small effect on the transient response generated during a whisker–object collision. Instead, the transient response will depend in large part on how the rat chooses to decelerate its vibrissae after the collision. The model allows experimentalists to estimate error bounds on quasi-static descriptions of vibrissal shape, and its predictions can be used to bound realistic expectations from neurons that code vibrissal sensing. We discuss the implications of these results under the assumption that primary sensory neurons of the trigeminal ganglion are sensitive to various combinations of mechanical signals.

## Introduction

During tactile exploratory behavior, rats tap and brush their large macrovibrissae against objects (Vincent, 1912; Welker, 1964). When a vibrissa touches an object, reaction forces and bending moments are generated at the vibrissa base (Kaneko et al., 1998; Solomon and Hartmann, 2006; Birdwell et al., 2007; O'Connor et al., 2010; Quist and Hartmann, 2012). These mechanical signals are transduced to electrical activity by mechanoreceptors within a densely innervated follicle (Ebara et al., 2002; Lottem and Azouz, 2011), and the activity from many mechanoreceptors is then integrated by primary sensory neurons of the trigeminal ganglion, Vg (Szwed et al., 2003, 2006; Jones et al., 2004; Leiser and Moxon, 2007; Khatri et al., 2009; Lottem and Azouz, 2011).

The forces and bending moments at the vibrissal base thus constitute the primary physical variables to which Vg neurons respond, and form the basis for vibrissotactile perception by the rat of its world. The difficulty, however, is that it is not possible to directly observe or measure forces and bending moments at the vibrissal base. Mechanical modeling is necessary to calculate these variables.

Previous studies have developed quasi-static models that relate the bending of the vibrissa to the forces and moments at the vibrissal base (Solomon and Hartmann, 2006, 2011; Birdwell et al., 2007; Quist and Hartmann, 2012). Other studies have developed models of vibrissal resonance (Hartmann et al., 2003; Neimark et al., 2003; Boubenec et al., 2013; Yan et al., 2013). To date, however, we lack a general dynamic model of the vibrissa that includes the effects of inertia, damping, and collisions.

Incorporating collisions into a dynamic model is a particularly challenging problem, from both an analytic and a computational perspective. The last few decades have seen an increase in the number of approaches to incorporate impacts in dynamical models (Brogliato, 2000). These developments have been driven largely by animation and gaming applications (Baraff, 1989; Trinkle et al., 1997; Stewart, 2000), which focus on results looking correct. Fortunately, there has also been a recent surge of interest in principled and physically meaningful ways of solving the simulation problem (Marsden and West, 2001; Pandolfi et al., 2002; Fetecau et al., 2003; Lew et al., 2004), with applications in the realms of modeling and control (Pekarek and Marsden, 2008).

In the present work, we exploited a subset of the latter group of approaches to develop a variational integrator-based model of vibrissa dynamics. The model enables us to quantify the time-varying forces and moments at the vibrissa base during both noncontact (“free-air”) whisking and as the vibrissae tap against an object (“collisions”). With the assumption that Vg neurons are sensitive to various combinations of forces and moment at the vibrissal base, the model generates a set of predictions that bound realistic expectations for the input to these neurons and how they may code features of the environment during vibrissotactile exploratory behavior.

## Materials and Methods

The model of the vibrissa developed in the present article applies to the mechanics of any macrovibrissae on a rat of either sex.

To improve readability, the first five subsections in Materials and Methods provide an overview of all the methods. These five subsections will provide sufficient information for most readers. The next three sections, titled “Details of the simulation method,” “Governing equations and their implementation within the algorithm” and “Linearization of the system for parameter optimization,” provide more technical details about the model.

The simulation in the present work makes use of a framework for scalable variational integration (Marsden and West, 2001; Fetecau et al., 2003; Johnson and Murphey, 2009). Equations for collision resolution were developed by Seghete and Murphey (2010).

#### Schematic of the vibrissa and the definition of variables

All simulations were constrained to two dimensions. For simplicity of exposition in the present work, the vibrissa was set to have 0 intrinsic curvature. That is, the undeflected vibrissa was straight.

The shape of the vibrissa was represented as a series of 13 rigid links, as shown in Figure 1*A*. The choice to use 13 links was somewhat arbitrary, reflecting a tradeoff between resolution and computation time. Each link was a conical frustum whose mass was modeled as a point mass at the center of mass of the link.

The joints between links were defined as nodes, which are shown as black dots in Figure 1*A*. Each node contains a torsional spring (defined by a spring constant, *k*, where τ = −*k*θ) and a torsional damper (defined by a damping constant, *b*, where τ = −*b*ω*)*. The first node was defined to lie at the origin (0,0). The position of each subsequent node *n* was then defined in terms of a rotation and translation from the previous node *n* − 1, as schematized in Figure 1*B*. The first link, which connects the first and second nodes, was considered to represent the vibrissa inside the follicle.

In all simulations of the present study, the base of the first link was rigidly controlled. This implies that the spring constant at the first node is infinite and the damper at the first node is 0.

#### Model parameters: masses, springs, and dampers

In all simulations, we initialized our model parameters measured from the D1 vibrissa of Hartmann et al., 2003: a base radius of 91 μm, a tip radius of 1.75 μm, an arc length of 49 mm, and a total mass of 0.498 mg. The vibrissa had linear taper. Each link had a length of 3.78 mm (49 mm/13 links).

To determine how to distribute the stiffness and damping across the nodes of the model, we optimized these two parameters so that the frequency response of the model matched the experimentally measured frequency response of the D1 vibrissa shown in the study by Hartmann et al., 2003 (their Figure 7).

The optimization assumed a constant density of 1.14 mg/mm^{3}, and this parameter was not allowed to vary. The mass of each link was set equal to the density multiplied by the volume of the link.

The optimization allowed Young's modulus to vary linearly along the length of the whisker (Quist and Hartmann, 2012). The stiffness at each node was set equal to Young's modulus multiplied by the area moment of inertia of the associated link, divided by the link length.

To determine the damping in the optimization routine, we imagine a pendulum extending from each node. The length of the pendulum is the distance from that node to the center of mass of the remainder (distal) portion of the whisker. For each run of the optimization, this pendulum has some (known) mass and some spring constant (*k*, which is also being optimized). We represent the damping for the pendulum as the variable ζ, and this is the damping variable we optimize over. The damping parameter *b* for that node is then calculated as *b* = 2ζ*L _{p}*(

*M*)

_{p}k_{p}^{[1/2]}, where

*L*and

_{p}*M*are the length and mass of the pendulum, and

_{p}*k*is the spring constant at that node (and also the spring constant of the pendulum).

_{p}We then linearized the vibrissa model about its resting shape and simulated shaking the vibrissa to generate a frequency–response curve. The model parameters were optimized by fitting the frequency response of the linearized system to the resonance curve from Hartmann et al. (2003) using TREP for the linearization (Johnson and Murphey, 2011) and the derivative-free optimization routines that are part of the SciPy python software package. The cost function for the optimization consisted of the L^{2} norm of the distance between the theoretical frequency response and the experimentally determined frequency response. We included terms penalizing deviations from our estimates for spring constants. The density of each link was not allowed to vary.

The final model parameters are shown in Table 1.

#### Simulating vibrissal movement and computing forces and moments at the vibrissal base

Vibrissa movement was simulated by rotating and/or translating the anchor (Fig. 1*A*, link 1). To compute reaction forces and moments at the vibrissa base, link 1 was permitted 3 degrees of freedom in a global coordinate frame: rotation about the *z*-axis (*R _{z}*), translation along the

*x*-axis (

*T*), and translation along the

_{x}*y*-axis (

*T*). Reaction forces and moments at the vibrissal base were then computed by finding the constraints that would be necessary to completely eliminate these translations and rotations. When

_{y}*R*,

_{z}*T*, and

_{x}*T*were each constrained, they yielded the bending moment (

_{y}*M*), the axial force (

*F*), and the transverse force (

_{x}*F*), respectively. These forces and moments represent those that would be exerted by the vibrissa on the follicle during active behavior.

_{y}#### Simulation method: overview

We simulate the vibrissa using numerical methods that are tailored for generalized coordinates (Fig. 1), and that excel at handling contact and impacts. The generalized coordinates are *q* = (*y*_{0}, θ_{1}…θ* _{n}*), where

*y*

_{0}is the linear translation (in the

*y*-direction) of the base of the whisker and θ

*is the angle between the two links that connect at node*

_{i}*i*.

We assume the vibrissa is very stiff in its axial direction. As will be explained in Discussion, we did not choose to use a finite-element software package in part because the axial stiffness would lead to a multiscale problem. By using a Lagrangian approach to modeling the vibrissa, we can describe the dynamics in the coordinate system in Figure 1 and avoid numerically simulating the fast dynamics associated with the axial stiffness. Forward dynamics are then computed for our Lagrangian system using variational integrators, as developed by Marsden and West (2001). This type of integrator provides discrete time simulations of constrained mechanical systems involving impact (Fetecau et al., 2003). The resulting discrete-time equations are of the following form:
where *L _{d}* is a discrete approximation of the continuous-time Lagrangian (kinetic energy minus potential energy) and

*q*is the configuration at time

_{k}*k*.

Equation 1 expresses conservation of discrete momentum. For comparison, it may be helpful to note that in continuous time momentum is defined as the derivative of the Lagrangian with respect to velocity. Note that energy is not explicitly represented in Equation 1. Instead, the energy and momentum conserving characteristics of Equation 1 are a consequence of the integration method (Marsden and West, 2001).

Equation 1 is an implicitly defined root-finding problem, in contrast with the class of explicit integration techniques such as Euler integration. We can solve for *q _{k}*

_{+ 1}(in term 2 on the left-hand side) given the initial configurations

*q*and

_{k}*q*

_{k}_{− 1}(in terms 1 and 2). Additional terms are required in this equation to add damping and geometric constraints to the vibrissa model. Geometric constraints allow us to drive the vibrissa in any desired way (e.g., 8 Hz motion) and to ensure that the vibrissa does not pass through an object (e.g., rotation into a peg).

Briefly, the simulation consists of solving Equation 1 at every time step *k*. First, we compute the first term, based on configuration states *q _{k}*

_{− 1}and

*q*. Second, we vary the value of

_{k}*q*

_{k}_{+ 1}within the second term until the sum of the first and second terms equals 0, as required by Equation 1. In practice, we solve for

*q*

_{k}_{+ 1}until the norm of the left-hand side of Equation 1 is below a threshold (e.g., close to machine precision). With an acceptable solution for

*q*

_{k}_{+ 1}, we increment

*k*and repeat the process. More details of the simulation algorithm are provided in the subsection titled “Governing equations and their implementation within the algorithm.”

#### Adding collisions: overview

The process described above allows us to compute vibrissal dynamics during noncontact whisking, but does not include collisions.

The calculation to resolve a collision requires three steps (Seghete and Murphey, 2010). First, we find the time when the collision occurred. Second, we find the energy of the collision (*E _{c}*) as well as the constraint force required to prevent the vibrissa from penetrating the object surface. Third, we compute the configuration of the system for the time (

*k*+

*1*). Details about each of these three steps are provided within the subsection titled “Details of the stimulation method.”

The effect of a collision is parameterized by the coefficient of restitution. The coefficient of restitution is a value between 0 and 1 that indicates the amount of kinetic energy lost by the system due to the collision. A coefficient of restitution equal to 0 is associated with a perfectly inelastic collision. Kinetic energy is not conserved. A coefficient of restitution equal to 1 is associated with a perfectly elastic collision. Kinetic energy is the same before and after the collision.

As an example, suppose a ball is dropped from rest from a particular height onto a flat surface. The coefficient of restitution will determine how high the ball bounces back. In this example, the configuration variable for the ball is the height (considered a variable translation in the *y*-direction). A perfectly elastic collision means that the ball returns to the original height from which it was dropped. A perfectly inelastic collision means that the ball stops at the surface and does not bounce at all.

#### Details of the simulation method

We reiterate that the first five subsections of Materials and Methods provided an overview of the mechanical modeling methods with sufficient information for many readers. The present section now provides details of the algorithms used.

##### General algorithm.

The general algorithm for computing the free dynamics of the vibrissa is outlined in Figure 2. The algorithm uses conservation of discrete momentum to generate the equations needed for a variational integrator. These equations are then solved by a Newton–Raphson iterative root finder to obtain the configuration *q _{k}*

_{+ 1}at time

*t*

_{k}_{+ 1}. The index of the time step is then incremented, the equations updated, and the configuration at the next time step found through another run of the iterative root finder.

The algorithm presented in Figure 2 assumes, apart from the two initial configurations, that the user has a way of calculating *h*(*q*) (the constraints on the system, usually known in closed form), *f*_{d}^{±} (the discrete external forcing, which includes damping terms), and also *L _{d}* (the discrete Lagrangian in Eq. 1). For simple systems, these functions can be evaluated and implemented analytically. As the number of links increases, however, the expressions become rather unwieldy, and a more careful approach is needed (Johnson and Murphey, 2009).

The algorithm consists of a main for loop, which iterates through the time steps by incrementing the index *k*. For each time step, a root-finding problem is initialized, where the 0 crossing of a function, *g _{k}*(

*q*, λ), is to be found. This is done using a Newton–Raphson iterative algorithm [the (repeat … until) loop] in which the initial guesses for

*q*and λ [the Lagrange multipliers that enforce the constraint

*h*(

*q*)] are refined with the aid of the derivative

*Dg*(

*q*, λ) in Equation 1. The process is then repeated until the value of

*g*is deemed sufficiently close to 0.

If no collisions occur, the general algorithm shown in Figure 2 is sufficient to solve for the motion of the system. On the other hand, if a collision occurs, it must be resolved using a second algorithm, which is described in detail in the following section.

##### Resolving a collision.

The calculation to resolve a collision in many ways follows the steps required to compute the general dynamics of the system. However, the following two additional steps are required: first, the precise time of the collision must be determined; and second, the change in energy due to the collision must be computed based on the coefficient of restitution. Below, we briefly describe the methods of collision resolution, as developed by Seghete and Murphey (2010). The algorithm is shown in Figure 3.

To determine the time of the collision, a functional form of the surface is used to determine the distance of each point to the surface. When the vibrissa passes into the surface, the sign of this distance becomes negative. This sign change allows the use of a Newton–Raphson root solver algorithm to determine the exact time at which the 0 crossing occurred. The configuration at the precise moment of impact is thus known. Analogously, the time of detachment from a surface is indicated by a sign change in the Lagrange multiplier λ associated with that surface.

The second step requires us to compute the energy removed from the system. We must first choose a value for the coefficient of restitution. The *E _{c}* is determined by the constraint force required to prevent the vibrissa from penetrating the object surface. If the collision is elastic, then this step directly computes the subsequent

*q*value. If the collision is perfectly inelastic, then the collision energy is computed to ensure that the vibrissa lies directly on the object surface. At the same time, the configuration for the vibrissa in this state is computed for time

_{k}*t*. If the collision is neither perfectly elastic nor perfectly inelastic, the following three steps take place: first, the energy of the collision is computed as though the collision were perfectly inelastic; second, the energy removed by the collision is calculated by multiplying

_{k}*E*by the coefficient of restitution; and third, this new energy of the collision,

_{c}*E*

_{c}

^{*}, is held fixed while the configuration of the system at time

*t*is recomputed.

_{k}The final step of the collision resolution involves computing the configuration of the system for the time *k* + 1. This uses the configuration at timestep *k* and also the midpoint average of the configuration at timestep *k*, and the configuration at the time of the collision. The simulation then returns to the general algorithm and continues as previously described in Figure 2.

#### Governing equations and their implementation within the algorithm

This subsection is intended for those readers interested in the details of the algorithms schematized in Figures 2 and 3.

We begin with a discussion of continuous Lagrangian mechanics. We next show how discrete Lagrangian mechanics differs from continuous Lagrangian mechanics. Finally, we describe the implementation of the discrete Lagrangian equations within the algorithm framework. In general, this section follows the explanation in Johnson and Murphey (2009). Please see that article for further information on the derivation and implementation of the following equations.

##### Continuous Lagrangian mechanics.

The Lagrangian, *L*(*q*,*q̇*), of a system is defined as the kinetic energy of the system, *T*(*q*,*q̇*), minus its potential energy, *V*(*q*), as follows:
The integral of the Lagrangian along a trajectory from *t*_{0} to *t _{f}* is defined as the action,

*S*, as follows: The “principle of least action” states that the system will follow a trajectory that extremizes the action (Marsden and Ratiu, 1999). The variational principle is used to extremize the action integral, producing the following Euler–Lagrange equations: The Euler–Lagrange equations can be solved using standard integration techniques to produce the trajectory of the system.

If rat vibrissae moved only in free air, and did not collide with objects, we could use these standard integration techniques to solve for the dynamics of the system. The present study, however, is primarily interested in how vibrissae will respond when they collide with objects. Equation 4 cannot be used to solve for interactions that involve collisions because it assumes smooth trajectories. We can avoid this restriction by turning to discrete Lagrangian mechanics, as described in the next section.

##### Discrete Lagrangian mechanics.

To incorporate the effect of collisions, our simulations use discrete Lagrangian mechanics. In this formulation, we define an *L _{d}*, which approximates the action integral of Equation 3 as the action sum, as follows:
where
The action sum in Equation 5 is also extremized according to a variational principle, producing the discrete Euler–Lagrange equations, as follows:
Equation 7 is the same as Equation 1, but here we have explicitly included notation to indicate the time step for

*L*. This additional notation is necessary when solving for the effect of a collision.

_{d}Next, we want to add damping to the whisker and add the ability to drive the whisker at its base. To do this, we need to introduce external forcing through the Lagrange–d'Alembert principle (a slightly modified version of the least-action principle) and also impose trajectory constraints (e.g., the bottom node must move at 8 Hz). Any external force (e.g., damping or friction) must be converted from its continuous form to a discrete form that is consistent with the discrete mechanics framework. The equation to approximate a continuous force with a discrete force is as follows:
where *f*(*q*,*q̇*,*t*) is the expression for the continuous external forcing, including damping. The expressions for *f*_{d}^{+} and *f*_{d}^{−}, just like *L _{d}*, are required by algorithms 1 and 2.

We then add constraints, *h*(*q _{k}*,

*t*), and external forcing,

_{k}*f*, to the left side of Equation 7 to produce the constrained, forced discrete Euler–Lagrange equations: where λ

_{d}*are Lagrange multipliers that must be found to enforce*

_{k}*h*(

*q*,

_{k}*t*).

_{k}The constraints *h*(*q _{k}*,

*t*) depended on how the whisker was being driven. For example, if we rotated the base of the whisker in a sinusoidal whisking motion, then the constraints would be

_{k}*h*(

*q*,

_{k}*t*) = [

_{k}*y*

_{0}− constant, θ

_{1}−

*A*sin(ω

*t*+ φ)], where

_{k}*t*is the time,

_{k}*y*

_{0}is the base translation in the

*y*-axis, θ

_{1}is base rotation about the

*z*-axis, ω is the angular frequency of whisking, and φ is the phase offset.

Damping is implemented using external forcing. Any other external force (e.g., friction) would also be applied through the external forcing terms. For example, to add linear dampers, *f* would include a velocity-dependent term:
where, in most cases, *C* is a diagonal matrix of damping coefficients.

##### Constrained, forced system (general set of equations).

Now that we have defined the basic constrained, forced discrete Euler–Lagrange equations, we can solve for the motion of the system.

It is important to note that when solving the discrete Euler–Lagrange equation, including damping and constraints (Eq. 9), there are always three discrete time steps involved. The first two time steps are used to solve for the configuration at the third time step. Thus, Equation 9 uses discrete steps *t _{k}*

_{− 1}and

*t*to compute

_{k}*t*

_{k}_{+ 1}.

Suppose then, for example, that we want to solve Equation 9 for a forced, constrained system with *N* configuration variables. In other words, *q* is a vector of length *N*. Also, assume the system has *M* constraints, *h*(*q _{k}*,

*t*). The

_{k}*h*(

*q*,

_{k}*t*) are functions that must equal 0 when the constraint is satisfied. There must be the same number of Lagrange multipliers (λ

_{k}*) as there are constraints*

_{k}*h*(

*q*,

_{k}*t*).

_{k}For the constrained system with no collisions, we therefore have (*N* + *M*) unknowns (*q _{k}*

_{+1}is

*N*-dimensional and λ

*is*

_{k}*M*-dimensional), and the following (

*N*+

*M*) equations:

##### Constrained system with a collision (step 1 of 3), find the time of collision.

Collisions are added by introducing an object (e.g., a surface). The object boundary is defined as a function of the system state, Φ(*q*). When Φ(*q*) > 0, the configuration is considered valid and the system (e.g., our model whisker as represented by Eq. 9) is not in contact with the object. When Φ(*q*) < 0, the configuration is considered invalid, and the system has penetrated the object. Finally, when Φ(*q*) = 0, the system exists on the object boundary.

With this framework in mind, Φ is evaluated at each new configuration *q* to determine whether a collision has occurred (i.e., Φ(*q*) < 0).

When a collision has occurred, the first step is to determine the time of the collision, *t _{A}*, at the configuration

*q*. This adds one new unknown to the general set of equations. In addition, we have to step back two time steps to

_{A}*k*− 2 to determine the time of collision. A schematic representation of solving a collision is shown in Figure 4. For the constrained system with collisions, we therefore have (

*N*+

*M*+1) unknowns (

*q*is

_{A}*N*-dimensional, λ is

*M*-dimensional, and

*t*adds one dimension), and the following (

_{A}*N*+

*M*+ 1) equations:

##### Constrained system with a collision (step 2 of 3), find the first postcollision configuration.

The next step is to determine the new state of the system resulting from the collision. This new state is directly related to the type of collision that occurred: elastic (no energy is removed due to the collision, coefficient of restitution = 1); inelastic (enough energy is removed so that the system exists on the object boundary, coefficient of restitution = 0); or something in between (a coefficient of restitution somewhere between 0 and 1). We now consider each of these types of collisions in turn.

##### Case 1: perfectly elastic collision.

For this type of collision, no energy is lost due to the collision. This does not mean the system itself is not losing energy (e.g., due to damping at each node). Instead, this means that the instantaneous collision event, itself, does not remove energy. Because no energy is lost at the instant of collision, only the constraint force required to keep the system from penetrating the object must be computed. The constraint force is defined as λ_{C}. Thus, there is 1 additional unknown. For the constrained system with collisions, we therefore have (*N* + *M* +1) unknowns (*q _{k}* is

*N*-dimensional, λ

*is*

_{A}*M*-dimensional, and λ

*adds one dimension), and the following (*

_{C}*N*+

*M*+ 1) equations:

##### Case 2: perfectly inelastic collision.

For a perfectly inelastic collision, enough energy is removed from the system so that all momentum normal to the object surface is eliminated. This amount of energy is given the name *E _{C}*

_{.}As with the elastic collision, we also must solve for the constraint force required to keep the system from penetrating the object. The constraint force is again defined as λ

*.*

_{C}Thus, there are two additional unknowns, for a total of (*N* + *M* + 2). The configuration *q _{k}* is

*N*-dimensional, λ

*is*

_{A}*M*-dimensional, and λ

*and*

_{C}*E*both add one dimension. Correspondingly, there are (

_{C}*N*+

*M*+ 2) equations:

##### Case 3: collision with a coefficient of restitution between 0 and 1.

Resolving a collision that has an intermediate coefficient of restitution (C.O.R.) requires solving the set of equations twice.

The first time the equations are solved, the system is treated as undergoing a perfectly inelastic collision. This yields the maximum energy that could be lost by the system, *E _{C}*, due to the collision.

The coefficient of restitution (C.O.R.) is applied to *E _{C}* via multiplication as follows:
Finally, the set of Equations 19–22 is solved as though the system had a perfectly elastic collision with the new term

*E*

_{C}

^{*}added as a constant to Equation 21.

##### Constrained system with a collision (step 3 of 3), find the second postcollision configuration.

Finally, we must find the second postconfiguration of the system after the collision has occurred at time step *t _{k}*

_{+ 1}.

##### Cases 1 or 3: perfectly elastic collision or collision with an intermediate coefficient of restitution.

If the collision was perfectly elastic or had a coefficient of restitution that resulted in a Φ(*q*) > 0, then the general set of equations is used. We therefore have (*N* + *M*) unknowns (*q _{k}*

_{+1}is

*N*-dimensional and λ

*is*

_{k}*M*-dimensional) and the following (

*N*+

*M*) equations:

##### Case 2: perfectly inelastic collision.

If the collision were perfectly inelastic, then the same Equations 24 and 25 are used, but now with the addition of a new constraint, Φ(*q _{k}*

_{+ 1}) = 0. We therefore have (

*N*+

*M*+ 1) unknowns (

*q*

_{k}_{+1}is

*N*-dimensional, λ

*is*

_{k}*M*-dimensional, and λ

_{Φ}adds one dimension) and the following (

*N*+

*M*+ 1) equations:

#### Linearization of the system for parameter optimization

In a separate calculation, we linearized the model for noncontact whisking to determine the optimal parameters to match the frequency–response curve of a real whisker (Hartmann et al., 2003, their Figure 7, whisker marked “D1”). The linearization procedures used standard techniques.

The linearized Lagrangian was defined as the kinetic energy minus the potential energy, as follows:
where *n* is the link number, *q*^{0} represents the translation of the base of the vibrissa along the *y*-axis, *L _{S}* is the length of a single link, and

*q*is the

^{i}*i*th element of the vector. The superscript indexing should not be confused with the exponent operator.

This was then used to derive the Euler–Lagrange equations that govern the system dynamics, as follows:
With *x* = [*q*, *q̇*]* ^{T}*, this system of differential equations can be rewritten as follows:
where

*x*now represents the full state of the linear system under investigation. We add external forcing through a vector

*u*(the velocity at the base of the whisker) and compute the value of

*y*that is, the deviation at the tip of the whisker, as follows: where, in this case,

*B*= [1, 0, 0,…]

*and*

^{T}*C*= [1,

*L*,

_{S}*L*,…]. The values of these matrices give us a state space representation of our linear system. Converting the state space system into a transfer function representation allows us to calculate the frequency response and compare it to experimental data.

_{S}## Results

Results are divided into three main sections. First, the model is validated against a quasi-static elastic beam model that accurately matches experimental data (Solomon and Hartmann, 2006, 2008; Birdwell et al., 2007; Quist and Hartmann, 2012). The frequency–response curve of the model is then compared with the experimentally determined frequency–response curve that was used to optimize the model parameters. Second, the model is used to predict the forces and moments that will occur during noncontact whisking, that is, as the rat whisks in free air. Third, the model is used to predict the forces and moments that will occur at the vibrissal base as the rat whisks against objects.

Each section of the Results that makes predictions concludes with a paragraph that directly addresses behavioral and physiological relevance. For these paragraphs, it is important to keep in mind that each Vg neuron is almost certainly sensitive to a combination of mechanical signals. There is no evidence, for example, that distinct classes of afferents are tuned to axial and transverse forces. Vg neurons sensitive to axial force might also be sensitive to transverse force, and/or bending moment.

### Slow deflections of the dynamic vibrissa model match quasi-static results

We first validated the dynamic model in the limiting case of extremely slow motion. In this case, we expect results to converge with a quasi-static solution. Quasi-static solutions for beam bending depend only on beam stiffness and applied loads; they do not include inertial or damping effects.

Using parameters representative of a large, caudal vibrissa (see Materials and Methods), we simulated the very slow rotation of a vibrissa against a peg placed either at 10% out along its length, or at 95% along its length. We then compared results from the dynamic model during this slow collision with results from a well validated quasi-static model (Solomon and Hartmann, 2006, 2008; Birdwell et al., 2007; Quist and Hartmann, 2012). The quasi-static model used the same number of links (13) as the dynamic model. Figure 5 shows results for two radial distances along the arc length (*s* = 10% and *s* = 95%). Both the moment at the base (Fig. 5*A*,*B*) and all force components (Fig. 5*C*,*D*) matched to within 10^{−5} or better between the two models.

Some consideration of the signs and shapes of the curves is useful. In these simulations, the vibrissa base is placed at the origin of a standard Cartesian coordinate system, and the vibrissa is assumed to rotate counterclockwise against a peg in the first quadrant. As the vibrissa deflects against the peg, the bending moment increases in magnitude, but its sign is negative because the force applied by the object is in the negative *y*-direction. If the peg had been placed in the fourth quadrant and the vibrissa had rotated clockwise against it, the moment would have been positive in sign.

Moment is equal to the cross-product of the following two vectors: (1) the vector connecting the vibrissa base to the point of object contact; and (2) the applied force vector. The magnitude of the cross product is equal to the product of the magnitude of each vector multiplied by the sine of the angle between the two vectors. As the vibrissa rotates through an increasing angle against the peg, the sine of the angle between the two vectors in the cross-product begins to decrease, while the magnitude of the applied force vector increases. This is why the moment trace follows the transverse force closely, but not exactly.

Transverse force and axial force are the components of the applied force projected into an *x*/*y*-coordinate system that rotates with the vibrissal base. A negative axial force means the vibrissa is in compression (i.e., being pushed into the follicle), while a positive axial force means the vibrissa is in tension (i.e., being pulled out of the follicle). A negative transverse force means the vibrissa is being pushed in the negative *y*-direction, while a positive transverse force means the vibrissa is being pushed in the positive *y*-direction.

In Figure 5, *C* and *D*, the transverse force and the axial force are both negative. This means the applied force is pushing the vibrissa in the negative *y*-direction and also generating a compressive axial force (i.e., pushing the vibrissa into the follicle).

### Frequency response of the vibrissa model optimized to match experimental data

As described in Materials and Methods, we used the frequency–response curve from experiments on a D1 vibrissa (Hartmann et al., 2003) to tune the parameters of our dynamic model. The base of the vibrissa was translated sinusoidally through a range of frequencies (1–150 Hz) with an amplitude that scaled inversely with the angular frequency squared (constant acceleration).

The top and bottom rows of Figure 6 reveal that a good match was obtained between model and experiment using the optimization algorithm described in Materials and Methods. The subplots in Figure 6*A* illustrate the displacement amplitudes of the vibrissal base and tip. Both simulation results and the experimental data show a first mode resonance peak near 40 Hz and a second mode resonance peak near 120 Hz.

The subplots of Figure 6*B* illustrate the magnification ratios for simulated and experimental data. Magnification ratios were found by dividing the tip amplitude by the base amplitude. As noted in Hartmann et al. (2003), and explained in Boubenec et al. (2013), the magnification curve of the real vibrissa displays a tendency to increase with frequency. This trend was captured by the model (Fig. 6*B*).

Notice that the base of the whisker was not simulated to oscillate with exactly the same amplitude as in the experiment (Fig. 6*A*, compare the amplitudes of the two base displacements). Nonetheless, the magnification ratio matches well between the experiment and the model (Fig. 6*B*). This demonstrates that the model yields the correct dynamics independent of exact excitation characteristics.

### Noncontact whisking: frequency content of moment and transverse force reflect the driving frequency, and the frequency content of axial force reflects twice the driving frequency

We next use the model to make predictions about the mechanics at the vibrissal base during whisking behavior. Specifically, we quantified the forces and moments that would be present at the base of a vibrissa as the rat whisks in free air, without making contact with an object (“noncontact” whisking). The vibrissa was started from rest and rotated at its base at 8 Hz through an amplitude of 20°, assuming a rigid constraint for the first link. This meant that the first link was constrained to follow exactly the 8 Hz driving motion. The motion of the second link and all subsequent links was free.

Figure 7 illustrates the results of a time and frequency analysis of the following three variables of interest at the vibrissa base: moment, transverse force, and axial force.

As expected, moment and transverse force both followed the 8 Hz driving frequency, and the traces were very similar (Fig. 7*A*,*B*, gray trace). Initially, both traces exhibit a high-frequency transient response near the first mode resonance frequency (50 Hz). This transient quickly damps out to leave only the 8 Hz rhythm in the steady state. Frequency analysis of the moment (Fig. 7*C*) and transverse force (Fig. 7*D*, gray trace) reveals the strong 8 Hz driving frequency and a resonance response near 50 Hz more than three orders of magnitude weaker than the driving frequency response.

In contrast, axial force exhibits a periodicity at twice the driving frequency (Fig. 7*B*, black trace) and is approximately six times smaller than the transverse force. Axial forces were greatest at the mid-phase of the whisk, peaking when the moment and transverse force crossed 0. These features of the axial force signal are clear when plotted in the frequency domain (Fig. 7*D*, black trace). Note that the axial force is positive because it is the force required to hold the vibrissa in the follicle (a tensile force, not a compressive force).

Three points are worth further exposition. First, none of the plots in Figure 7 would change significantly if the vibrissa were completely rigid. Only a small fraction of the moment in Figure 7*A* is due to the bending of the whisker; instead, the trace primarily reflects the moment required to accelerate the vibrissa. This means that the signal during noncontact whisking is extremely close to the product of the total moment of inertia and the angular acceleration of the vibrissa. Second, inertial effects are small. Both the moment and the transverse force are almost completely in phase with the position of the vibrissa. If additional mass were added, the vibrissa would have more difficulty following the 8 Hz driving frequency. Third, the axial force in Figure 7 is exclusively positive because the resting shape of the vibrissa is modeled as being straight. If the vibrissa had intrinsic curvature, the axial force would exhibit small negative “dips” during the trajectory. Although outside the scope of the present study, a preliminary analysis of the effects of intrinsic curvature suggests that its effects on dynamics are likely to be relatively small.

Five key points relevant to behavior emerge from this simulation. First, as suggested by Knutsen et al. (2008), the vibrissa can be well approximated as a rigid body up until the time of collision. Because inertial effects are small, the rat could regulate the position of the vibrissa base using a feedback law that is essentially an inverse kinematics calculation. Second, precisely because inertial effects are small, Vg responses to noncontact whisking are expected to be small compared with those generated when the vibrissa makes contact with an object. This prediction is consistent with physiological recordings from Vg neurons showing relatively small responses during free-air whisking behavior (Zucker and Welker, 1969; Leiser and Moxon, 2007; Khatri et al., 2009). The Discussion section suggests experiments that could determine the extent to which these “whisking” responses are due to inertial effects or to muscles squeezing down on the follicle. Third, during noncontact whisking, neurons sensitive to bending moment or transverse forces may respond near the whisking frequency. This could permit them to provide a signal related to the angular position or (spatial) phase of the vibrissa. Fourth, during noncontact whisking, neurons sensitive to axial forces are likely to respond twice during a whisk (Khatri et al., 2009, their Fig. 2, neuron 3). This would permit them to provide a signal related to the speed (absolute value of the velocity) of the whisk. Finally, a key characteristic of all mechanical signals during noncontact whisking is that they are very small. We therefore expect Vg neurons to respond largely in a stochastic manner to noncontact whisking, consistent with the high degree of variability observed experimentally (Leiser and Moxon, 2007; Khatri et al., 2009).

### Whisking against an object: energy losses due to the inelasticity of the collision

Colloquially, the term collision generally connotes a violent event, such as a vehicular crash. In the present work, however, a collision refers to very gentle tapping and pushing of the vibrissae against an object, with magnitudes that would typically occur during real whisking behavior. An important parameter choice in simulating collisions is the stiffness with which the vibrissa is held at its base. For all simulations, we chose to model the vibrissa as though it is held perfectly rigidly in the follicle. This means that all results represent the upper bound of signal strength expected for the behaving rat.

To begin to examine collision dynamics, we first evaluated the extent to which a vibrissa will lose energy due to inelasticity of the collision. The coefficient of restitution (see Materials and Methods) defines how much energy is removed at the instant of the collision. When the coefficient of restitution is 1, the collision is defined to be “perfectly elastic,” and none of the energy of the collision event is removed. When the coefficient of restitution is 0, the collision is defined as “perfectly inelastic,” and the energy of the collision event is entirely removed. The energy of the collision event is the amount required to ensure that the velocity of the entire vibrissa normal to the object is 0.

It is important to note that energy losses due to damping do not contribute to energy losses at the instant of collision. Whether a collision is elastic or inelastic is independent of damping.

The coefficient of restitution associated with a collision of a vibrissa against an object is not known, and is not practical to measure experimentally. Therefore, we deliberately bracketed the range from perfectly elastic to perfectly inelastic. We simulated both perfectly elastic impacts and perfectly inelastic impacts for the vibrissa rotating at 705°/s against a peg at both a proximal (arc length *s* = 40% of total whisker length) and distal (arc length *s* = 85%) location. 705°/s was chosen because it is the average speed of an 8 Hz whisk through a ∼45° amplitude. Rotation was stopped immediately upon impact with the peg to eliminate any effects of bending as the vibrissa further rotated against the peg. Figure 8, *A* and *B*, illustrates that the energetic difference between purely elastic and purely inelastic collisions was extremely small. The energy difference between the two conditions was only 0.19 nJ for the proximal collision, and 0.034 nJ for the distal collision. These values represent the difference in the total energy of the system at the time steps immediately before and immediately after the collision. Much larger energy losses are observed in the time steps that follow the collision. These losses are the result of damping.

To gain more physical insight into the effects of the collision, we quantified the motion of the colliding link normal to the object surface. Results are shown in Figure 8, *C* and *D*. In the case of the perfectly inelastic collision, the distance to the surface was, by definition, identically 0 (Fig. 8*C*,*D*, black trace). The link collides with the object and remains in contact. In the case of the perfectly elastic collision (Fig. 8*C*,*D*, gray trace), the link bounces off the object, but the amplitude is extremely small. The amplitude of the bounce for the proximal collision peaked near 20 pm, and near 400 pm for the distal collision. Subsequent bounces for the elastic collision were even smaller and are therefore not shown. In general, the amplitude of the bounce increases with the radial distance of the collision and decreases with damping. If damping were larger, the curves in Figure 8, *A* and *B*, would have steeper energy losses, and the peak amplitude of the bounce during an elastic collision (Fig. 8*C*,*D*) would be smaller.

Behavioral relevance: After a vibrissa collides with an object, it tends to resonate (Hartmann et al., 2003; Neimark et al., 2003; Ritt et al., 2008; Wolfe et al., 2008; Boubenec et al., 2013). Two independent factors contribute to why these resonant oscillations ultimately disappear. The first is the energy loss due to impact itself, occurring on a nearly instantaneous time scale. The second is the loss due to damping, which occurs over a much longer time scale (tens of milliseconds) as the vibrissa vibrates. The present simulations lend qualitative insight into the mechanical behavior one expects due to a collision: vibrissae will tend not to “bounce” after colliding with an object.

This finding has two important implications. First, when the rat explores a complex surface, vibrissae will tend to make surface contact only once. This fact may help ensure that tactile signals accurately reflect the object's spatial features. Second, the results imply that it is reasonable to model the collision of a vibrissa with an object as perfectly inelastic (coefficient of restitution equal to 0), but that it is important to pay careful attention to damping parameters, as these will almost certainly have a large influence on tactile exploration. In particular, multiple collisions with objects are more likely to be driven by the environment than by self-motion, such as the “stick-slip-ring” events that characterize texture exploration (Ritt et al., 2008; Wolfe et al., 2008). The precise timing and mechanical details of stick-slip-ring events will depend on a complex interaction among texture properties, intrinsic whisker dynamics, and active control. The damping characteristics of the vibrissa will almost certainly affect the degree to which the “ring” that occurs after a vibrissa “slips” off a ridge influences the precise timing of the following “stick.”

In line with the finding that energy losses due to impacts are small, all subsequent simulations in the present study are run assuming perfectly inelastic collisions.

### The velocity of the vibrissa at the instant when it collides with an object has little effect on the magnitude of the signals generated at the vibrissa base

Having determined that collisions are well approximated as perfectly inelastic, the next question that naturally arises is as follows: how does the velocity of the vibrissa at the time of the collision affect the mechanical signals at the vibrissa base? We performed simulations to answer this question.

The two graphs in the top row of Figure 9, *A* and *B*, are identical to each other. They illustrate how the vibrissa was rotated against the peg. To model a behaviorally realistic whisking trajectory, we varied the frequency of a half-sine wave to change the velocity at the time of collision. A half-sine trajectory is more realistic than a ramp profile because rats do not whisk at constant velocity or suddenly decelerate their vibrissae to 0. The frequency of the base of the vibrissa was defined to be a sine-wave with a frequency of 4, 8, or 16 Hz, and the base of the vibrissa was always rotated from −14° to 6°. The vibrissal velocities at the time of collision were 232, 465, and ∼955°/s. These velocities were obtained by ensuring that the whisker always hit after the whisker had gone 14° through its cycle.

The remaining plots of Figure 9 illustrate the bending moment and rate of change of bending moment at the vibrissal base during a proximal (40%) and a distal (85%) collision. Results are not shown for transverse and axial force because results were very similar to those illustrated for bending moment.

Figure 9, *C* and *D*, shows results for a proximal collision. The negative moments observed before collision correspond to the initial acceleration of the vibrissa toward the object (Fig. 8, illustrating noncontact whisking). The large negative slope that occurs immediately after the collision results from the vibrissa bending against the object. The higher-frequency sine wave (16 Hz) results in a steeper slope than the lower frequencies (8 and 4 Hz). The steeper slope is not a result of the collision per se, but occurs because the whisker is bending more rapidly.

Note that all signals in Figure 9*C* change on a time scale of ∼20–40 ms and can therefore be thought of as signals to which slowly adapting (SA) ganglion neurons might be sensitive. By contrast, rapidly adapting (RA) ganglion neurons would tend to be more sensitive to the rate of change of mechanical signals, for example, the rate of change in bending moment shown in Figure 9*D*. In Figure 9*D*, it is clear that the bending following the collision causes a sharp rate of change of moment. The sharpest peak is for the collision that occurs at the highest frequency, again because the postcollision bending is more rapid.

Figure 9, *E* and *F*, shows results for a distal collision. In this case, the mechanical signals associated with the collision are approximately the same in amplitude as those associated with noncontact whisking (Fig. 7*A*,*B*). Only for the very highest frequency collision (velocity at time contact ∼955°/s) is a small negative “divot” observed as a result of the collision.

We can draw two broad conclusions from Figure 9. First, dynamic effects scale with velocity at the time of collision. As the velocity at the time of impact increases, so do the tiny vibrations associated with the collision. Second, however, these dynamic effects are extremely small compared with the effects of bending and deceleration. These much larger mechanical signals will depend on how far and how fast the rat chooses to push its vibrissae against the object after initial contact.

To demonstrate this effect more clearly, we ran the same simulations as in Figure 9 but varied the velocity of the collision by changing the frequency of the half-sine wave in 2 Hz increments. The results of this simulation are shown in Figure 10. For comparison, we also modeled noncontact whisking by running the identical simulations without any object present.

For noncontact whisking, we defined “peak moment” as the maximum of the absolute value of the moment generated during noncontact whisking. The noncontact whisking trace of Figure 10 shows that a greater peak moment is required to accelerate the vibrissa to higher velocities during noncontact whisking. This result follows from straightforward physics because the vibrissae is behaving essentially as a rigid body (compare to the results shown in Fig. 7).

For whisking against a peg, we defined peak moment as the magnitude of the largest negative peak after the collision occurred. The negative peak was chosen because the positive peak was extremely small for proximal collisions (Fig. 9*C*, compare positive and negative peaks). The traces labeled *s* = 40% and *s* = 85% in Figure 10 illustrate that when a collision occurs, the peak moment after the collision depends only weakly on the velocity of the collision. For example, if the vibrissa collides with an object near its tip (*s* = 85%) the peak moment is nearly constant, on the same order of magnitude as noncontact whisking. The slight upward trend in the last two data points reflects the small oscillations that occur after the collision, which were seen earlier in Figure 9. When the vibrissal–object collision occurs closer to the vibrissal base (*s* = 40%), the overall magnitude of the peak moment is larger because the vibrissa bends more, but the dependence of the peak bending moment on the velocity at the time of collision is small.

Behavioral relevance: Numerous behavioral and physiological studies have focused on whisking velocity as a key parameter that likely to influence the responses of Vg neurons (Gibson and Welker, 1983; Pinto et al., 2000; Shoykhet et al., 2000; Szwed et al., 2003; Jones et al., 2004; Leiser and Moxon, 2007; Towal and Hartmann, 2008; Khatri et al., 2009). Collectively, these studies have mixed two behavioral conditions of interest. The first behavioral condition involves the velocity of noncontact (free-air) whisking. As shown in Figure 7, our model predicts that the axial force could generate weak Vg responses related to the speed of noncontact whisking, while the bending moment could generate weak Vg responses related to the position of the whisker. The second behavioral condition involves the instantaneous whisking velocity at the time of collision with an object. Figure 10 illustrates that whisking velocity at the time of collision will have only a small influence on the mechanical signals that will be received by Vg neurons. Instead, the subsequent active control of the rat over vibrissal deceleration will be of paramount importance to Vg signals.

### Periodic protraction and retraction against a peg

Together, the results of Figures 9 and 10 illustrate that the signals generated at the vibrissal base by the impulse due to the collision in and of itself are very small. Instead, it is the bending and deceleration associated with the collision that cause the largest magnitude signal at the vibrissal base. Figure 9 further shows that these bending and deceleration effects can extend over a duration of ∼25–40 ms. This means that a collision occurring during any given whisk may affect the dynamics of the subsequent whisk.

To examine these effects, we simulated periodic 8 Hz whisking against a point object (a peg). The vibrissa was simulated to rotate 3°, 6°, or 9° against the object. As before, two radial distances were selected, representing a proximal collision (*s* = 40%) and a distal collision (*s* = 85%).

The top row of Figure 11 illustrates the bending moment, transverse and axial forces associated with the proximal collision, in which the object is placed 40% out along the vibrissa length. When these data are compared with those in Figure 7, *A* and *B*, we see that the mechanical signals associated with the collision are large relative to these same signals during noncontact whisking. For example, Figure 7*A* shows that the bending moment during noncontact whisking ranges between approximately −0.06 and 0.06 mN/mm. By comparison, Figure 11*A* shows that the bending moment associated with collision is nearly 25 times greater (−1.4 mN/mm) even for the low-amplitude rotation into the object of 3°. Similarly, the transverse and axial forces during noncontact whisking are in the range of −2.8 to 2.8 and 0 to 0.5 μN, respectively (Fig. 7*B*). During the collision, these forces are a minimum of −70.8 and −7.6 μN (Fig. 11*B*,*C*). These results mean that it would be relatively easy for the rat to detect a proximal collision, simply by thresholding the signal at a particular value.

The effect of vibrissal dynamics becomes increasingly important, however, when the collision occurs at the distal portion of the vibrissa (Fig. 11, bottom row). In this case, the magnitudes of the bending moment and the transverse force associated with the collision are in the same range as those associated with noncontact whisking (Fig. 7*A*,*B*). In addition, the vibrations generated by the first collision with the peg interfere with the signal generated by the second collision; vibrations from the second collision interfere with the third, and so on. Collision detection based on thresholding either of these two mechanical signals is likely to be difficult. In contrast, the axial force signal associated with the distal collision stays well above the signal associated with noncontact whisking (Fig. 11*F*).

Behavioral relevance: The only way that a rat can detect that a vibrissa has touched an object requires change in force and/or moment at the vibrissal base. These signals are, physically, what constitutes a “touch.” The results of Figure 11 illustrate that if the rat taps on an object near the vibrissal tip, the changes in bending moment and transverse force associated with the touch are comparable in magnitude to those associated with noncontact whisking. It therefore seems unlikely that a simple threshold of the responses of Vg neurons sensitive to bending moment or transverse force could reliably indicate that the vibrissa had touched an object. In contrast, reliable collision detection could occur if second-order neurons (e.g., in the trigeminal nuclei) performed a threshold of the responses of Vg neurons sensitive to axial force. These results suggest a particularly important role for a compressive (negative) axial force in the detection of distal collisions. It is important to note that these effects are likely to be more significant for the more caudal whiskers, which are larger and have more mass. In the case of smaller, more rostral whiskers, the dynamic response from one whisk may not interfere with the next.

### Proximal vibrissa collisions are well approximated with a quasi-static model

One striking feature of Figure 11 is how relatively smooth the signals appear for a proximal collision (*s* = 40%) compared with a distal collision (*s* = 85%). Although the vibrations are of comparable magnitude for proximal and distal collisions, the bending is much larger for the proximal collision, making the vibrations appear much smaller by comparison.

Figure 12 illustrates the relative effects of quasi-static bending and vibrissal dynamics. Results are shown only for bending moment because similar results were found for both transverse and axial force. Figure 12*A* plots the bending moment at the base of the vibrissa for quasi-static and dynamic models for proximal (*s* = 40%) and distal (*s* = 85%) collisions. For both collision locations, the quasi-static model clearly accounts for a considerable amount of the signal.

To begin to quantify this effect, we computed the total energy present in both the quasi-static and dynamic models. For the quasi-static model, the total energy of the system is exactly equal to the potential energy stored in the torsional springs. For the dynamic model, the total energy is the sum of the kinetic energy and potential energy. Figure 12*B* illustrates that for the proximal collision, the total energy matches extremely well between the two models. In the distal collision, however, it is clear that kinematic energy makes a significant contribution to the total energy.

Further quantification of the difference between quasi-static and dynamic results is shown in Figure 12*C*. In Figure 12*C*, values on the *y*-axis represent the percentage difference in total energy between dynamic and quasi-static signals, normalized to the energy in the dynamic signal. If the quasi-static model accounts for all the energy in the signal, the *y*-value will be 0. If the quasi-static model accounts for none of the energy in the signal, then the *y*-value will be 100%. We simulated collisions for the full range of radial distances for three values of the push angle θ* _{p}* (3°, 6°, and 9°; as in Fig. 11).

Two key trends can be observed in Figure 12*C*. First, the quasi-static model captures less of the signal as the collision occurs more distally along the vibrissa. This result is evidenced in the positive slope of each trace. Second, if we imagine drawing a vertical line at a single radial distance, across curves, we see that the quasi-static model can explain a greater fraction of the overall signal as θ* _{p}* becomes larger. This result makes good intuitive sense, as the bending of the vibrissa is by definition larger for larger values of θ

*.*

_{p}Behavioral relevance: Figure 12*C* illustrates that very proximal collisions (up to ∼30% out along the vibrissal length) can be fairly well characterized by quasi-static bending, even if the rat does not rotate the vibrissa very far against the object. At 80% out along the vibrissal length, however, dynamic effects account for nearly 50% of the total energy for a rotation of 3°. As the rat increasingly rotates its vibrissae against the object, the quasi-static approximation again becomes valid. Figure 12*C* suggests that Vg neurons classified as RA and SA will be differentially active, depending on both the radial location of contact as well as the value of θ* _{p}*. RA neurons are likely to be sensitive to the higher-frequency dynamic effects that characterize distal collisions and small pushing angles. SA neurons are more likely to be sensitive to proximal collisions and larger pushing angles.

Figure 12*C* may be particularly useful to researchers studying rat whisking behavior. Suppose, for example, that high-speed video is used to quantify the shape of the vibrissa as the rat whisks against an object. The shape of the vibrissa at all points in time can be fit to a quasi-static bending model to calculate forces and moments at the base (Kaneko et al., 1998; Solomon and Hartmann, 2006, 2011; Birdwell et al., 2007; O'Connor et al., 2010; Quist and Hartmann, 2012). Figure 12*C* allows the experimentalist to quantify how well that quasi-static model represents the true mechanical signals at the base for any radial distance, and for increasing angles through which the vibrissa has pushed against the object.

## Discussion

Neurons in the vibrissal–trigeminal system must respond to forces and bending moments at the vibrissal base. Because these primary stimulus variables are not directly observable, however, to date it has not been possible to characterize neural responses in terms of these variables. The present study has developed a dynamic model that we anticipate will allow us to relate neural responses directly to the forces and moments that characterize whisking behavior.

### Modeling approach and caveats

Standard techniques for dynamic modeling (e.g., FEM software) are poorly suited to model vibrissal dynamics because vibrissae are disproportionately stiff along their long axis, vibrissae are constrained when in contact, and vibrissae routinely experience impact. Combined, these characteristics tend to lead to “nonphysical” outcomes in standard FEM packages such as ANSYS (www.ansys.com). These nonphysical outcomes can be solved by introducing numerical stabilization, but this stabilization distorts the mechanical meaning of the simulation. Previous work has described vibrissal dynamics based on resonance models (Hartmann et al., 2003; Neimark et al., 2003; Yan et al., 2013), or as the linear superposition of quasistatic bending and resonance (Boubenec et al., 2013). These models have yielded excellent matches to experimental data under some conditions (e.g., Boubenec et al., 2013) but are not physically guaranteed to generalize. The present model generalizes to large deflections, describes both non-contact whisking as well as collisions and impacts, and is easily extended from two to three dimensions.

Techniques from the emerging field of discrete variational mechanics (Marsden and West, 2001; Fetecau et al., 2003) allow us to simulate the vibrissae with constraints and impacts, and do not require arbitrary parametric choices for numerical stabilization (Hairer et al., 2006). The vibrissa can be modeled as rigid along its length while being otherwise flexible, and Lagrangian mechanics guarantees that the computed behavior will be physical. Thus, variational integrators provide flexibility to model the vibrissae in a way that includes contact mechanics without relying on ad hoc numerical methods that may undermine interpretation.

Two primary caveats pertain to the interpretation of these simulation results. First, the springs, masses, and dampers at each node of the model were fit to the parameters of a single vibrissa; in the future, a full system identification should be performed. Second, we did not model the intrinsic curvature of the vibrissa. This can be added to future models and is likely to add asymmetries to mechanical signals, as noted in the text describing Figure 7.

### Implications for rat exploratory behavior

We interpret the results of the present study to suggest that rodents may use a position–control strategy to detect and localize objects, and then shift to a force–control strategy during object exploration and feature extraction. We elaborate this idea as it pertains to three stages of rat exploratory behavior.

#### Stage 1: noncontact whisking

During noncontact whisking, the rat scans the environment to detect and localize objects. There are no significant forces other than the rat's own muscles acting on the vibrissae. To detect an object, the rat must sense the mechanical transients generated by a collision. To localize the object in head-centered coordinates, the rat must determine the position of its whiskers relative to its head at the instant of collision.

The present work provides evidence that the rat could directly control the position of its whiskers right up until a collision. Specifically, Figure 7 shows that the first 60–70% of the whisker behaves as a rigid body (Knutsen et al., 2008) and that the bending moment closely follows the whisk cycle. Because inertia is low, the rat could regulate the position of the vibrissa base using a feedback law that is essentially an inverse kinematics calculation. This type of direct position control would not be possible if dynamic effects were larger.

Several studies have provided evidence that cortical structures send command signals to control the set point and position or phase of whisker motion (Mehta et al., 2007; Curtis and Kleinfeld, 2009). One problem, however, is that there are no proprioceptors in the whisker muscles (Ebara et al., 2002), so there is no obvious “handshake” signal to ensure that the vibrissae are actually moving in the manner dictated by central structures. So how does the rat know where its whiskers are relative to its head, at any given time during the whisk cycle? One possible answer is that no tactile handshake is required, and the rat relies entirely on the approximate position signal obtained by efference copy. A second possibility is that tactile feedback is provided by skin stretch around the whiskers. A third alternative is suggested by Figure 7. Figure 7 shows that during noncontact whisking all mechanical signals are extremely small. However, both *M* and *F _{y}* reliably reflect the driving frequency, and

*F*is maximized twice during the whisk. When averaged over a sufficient number of cycles or neurons, these signals could represent a handshake for central motor commands.

_{x}Notably, only weak Vg responses to noncontact whisking are observed in awake animals, and these responses are not strongly correlated with the kinematics of the overall whisking cycle (Leiser and Moxon, 2007; Khatri et al., 2009). We hypothesize that these Vg responses reflect the tiny forces and moments generated by the inertia of the whisker. This would explain the inexact correlations with overall whisk cycle kinematics: because the effects are so small, the relationship would be observed only if the mechanical signals were calculated for each whisker individually. If our hypothesis is correct, then Vg responses to noncontact whisking should vary systematically as a small mass is added at different locations along the vibrissal length. If, alternatively, the responses are largely due to muscles squeezing on the outside of the follicle, or to skin stretch, then they should be largely invariant to the location or magnitude of the mass.

#### Stage 2: collision detection

Although a handshake signal during noncontact whisking may be useful, it should not reduce sensitivity to the mechanical transients that indicate a collision. Figure 11 shows that during a distal collision, the bending moment and transverse force signals are on the same order of magnitude as those generated by free-air whisking. Thus, a simple threshold of these signals will not reliably indicate a collision. Vg neurons sensitive to these signals would likely be observed as whisking neurons, equally sensitive to whisking and light touch (Szwed et al., 2003; Leiser and Moxon, 2007). In contrast, the axial force signal maintains a particularly high signal-to-noise ratio during object contact and would permit reliable collision detection based on thresholding (Stüttgen et al., 2008).

Figure 8 shows further that whisker–object collisions tend to be inelastic. Thus, only a single contact with the object is likely to occur, ensuring high temporal precision and thereby accurate localization. In addition, Figure 10 shows that whisking velocity at the time of collision will have almost no influence on the mechanical signals at the vibrissal base. The rat can therefore expect every collision of a given vibrissa to generate a similar mechanical signature, regardless of whisking speed. This result may seem counterintuitive because it is easy to confuse the effects of the collision itself (which occur instantaneously) with the bending and deceleration that immediately follow. The mechanical signals generated by the bending and deceleration following a collision will, indeed, be large (Fig. 9).

Together, these results lend significant support to previous studies suggesting that the rat could determine the horizontal location of an object in head-centered coordinates (θ_{impact}) based on a temporal or phase code (Knutsen et al., 2008, 2009; Curtis and Kleinfeld, 2009; Kleinfeld and Deschênes, 2011). The essential idea is that there is a linear relationship between θ_{impact} and the time difference between whisk onset and object collision. The rat could monitor the time difference; the constant of proportionality is the speed of the whisk, and the intercept is the angular position of the vibrissa at the start of the whisk.

The mechanical results of the present work support this coding scheme in many ways, as follows: (1) low inertia enables direct position control; (2) a weak *F _{x}* signal is correlated with whisking speed during noncontact whisking and retains a high signal-to-noise ratio at the instant of collision; (3) large mechanical transients will immediately follow a collision; (4) vibrissae will tend to experience temporally precise, inelastic impacts; and (5) the rat can expect collisions to have a similar mechanical signature regardless of whisking velocity.

An alternative (or complementary) possibility, first suggested by Knutsen and Ahissar (2009), is that the “roll” of the whisker about its own axis could generate collision-related signals in different regions of the follicle, thereby providing a cue to the horizontal angle of contact.

#### Stage 3: object exploration

In contrast to noncontact whisking, significant external forces act on the vibrissae as a rat whisks against an object. Quasi-static models have suggested that these contact forces could permit the location of the vibrissal–object contact point to be determined (Solomon and Hartmann, 2006, 2010, 2011; Stüttgen et al., 2008; Williams and Kramer, 2010; Bagdasarian et al., 2013; Hires et al., 2013; Pammer et al., 2013).

The results shown in Figure 11, however, would seem to pose a problem for this coding scheme—when a rat whisks against an object, the mechanical signals generated during one whisk will interfere with those generated by the next. With the exception of axial force, dynamic effects take at least 20 ms to damp out.

Figure 9 helps to illustrate what we predict to be the key for the rat to solve this problem: the primary determinant of the postcollision mechanical signals at the vibrissal base is the active control by the rat of the deceleration profile (Mitchinson et al., 2007; Grant et al., 2009; Deutsch et al., 2012) We therefore specifically predict that during feature extraction, rats will “press in” their vibrissae against the object surface for durations between 20–60 ms, so as to damp the dynamic response and monitor forces and moments at the vibrissal base in the quasi-static regime. These extended contact durations have already been observed in some studies (Deutsch et al., 2012) and would also help to ensure that many vibrissae are all in contact with the object surface simultaneously.

Overall, our mechanical results point to the idea that the timing and sequence with which the whiskers make contact with an object are critical for object localization (Ahissar and Arieli, 2001; Knutsen and Ahissar, 2009) but may be largely irrelevant to feature extraction, which is likely to depend on contact forces.

## Footnotes

This work was supported by National Science Foundation (NSF) Awards IOB-0446391, IOS-0818414, and EFRI-0938007 to M.J.Z.H.; and CMMI-0951688 to T.D.M. B.W.Q. was partially sponsored by an NSF graduate research fellowship as well as a National Institutes of Health Ruth L. Kirschstein National Research Service Award (F31 NS065630). L.A.H. was partially sponsored by the DoD, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a.

The authors declare no conflicting financial interests.

- Correspondence should be addressed Mitra J. Z. Hartmann, Biomedical Engineering and Mechanical Engineering Departments, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208. hartmann{at}northwestern.edu