Tactile Modulation of Whisking via the Brainstem Loop: Statechart Modeling and Experimental Validation

Rats repeatedly sweep their facial whiskers back and forth in order to explore their environment. Such explorative whisking appears to be driven by central pattern generators (CPGs) that operate independently of direct sensory feedback. Nevertheless, whisking can be modulated by sensory feedback, and it has been hypothesized that some of this modulation already occurs within the brainstem. However, the interaction between sensory feedback and CPG activity is poorly understood. Using the visual language of statecharts, a dynamic, bottom-up computerized model of the brainstem loop of the whisking system was built in order to investigate the interaction between sensory feedback and CPG activity during whisking behavior. As a benchmark, we used a previously quantified closed-loop phenomenon of the whisking system, touched-induced pump (TIP), which is thought to be mediated by the brainstem loop. First, we showed that TIPs depend on sensory feedback, by comparing TIP occurrence in intact rats with that in rats whose sensory nerve was experimentally cut. We then inspected several possible feedback mechanisms of TIPs using our model. The model ruled out all hypothesized mechanisms but one, which adequately simulated the corresponding motion observed in the rat. Results of the simulations suggest that TIPs are generated via sensory feedback that activates extrinsic retractor muscles in the mystacial pad. The model further predicted that in addition to the touching whisker, all whiskers found on the same side of the snout should exhibit a TIP. We present experimental results that confirm the predicted movements in behaving rats, establishing the validity of the hypothesized interaction between sensory feedback and CPG activity we suggest here for the generation of TIPs in the whisking system.


Introduction
Many behaviors in which animals move their sensors in a repetitive, stereotyped pattern, such as whisking, sniffing, tasting, and looking, are based on an active, often periodic sampling of the environment by the sensor organs [1][2][3][4][5][6][7][8][9]. Periodic sampling is a class of closed-loop motor-sensory interactions that is hypothesized to be driven by the interaction of central pattern generator (CPG) activity and sensory feedback. However, the nature of the modulation of CPG activity by sensory signals, and the dependency of motor patterns on sensor-environment interactions are not yet understood. Additionally, it is unclear where in the brain these interactions occur. In this study, we used the well-studied brainstem loop of the vibrissal system of rats and a wellcharacterized vibrissal movement, ''touch-induced pump'' (TIP), to explore the interactions between CPG activity and sensory feedback.
A whisking CPG(s) is thought to dominate the control of ''freeair exploratory whisking'', a perceptual behavior in which the vibrissae are swept back and forth rhythmically [1,2,4,6-8] and in a coordinated manner [6,10,11]. The existence of a CPG was demonstrated by observations that showed that whisking in adult rats continues in the absence of peripheral sensory inputs [4,12,13] and of descending (e.g. cortical) control mechanisms [14,15]. Such a CPG was shown to reside within the brainstem [16,17]. The existence of whisking motor modulation based on whisker-derived sensory feedback was demonstrated behaviorally [11,18]; however, no explicit mechanism has been suggested so far to account for this modulation. We herein present such a mechanism for the induction of TIPs.
TIP is a whisking behavior that occurs upon whisker-object contact, and consists of a rapid retraction and protraction of the whiskers, during which contact with the object is preserved [4,11]. The TIP was chosen as a benchmark for theoretical modeling of sensory modulation of CPG-driven behavior, as it is sufficiently quantified and hypothesized to convey sensory feedback via brainstem circuitry.
We combined experimental observations with a dynamic, bottom-up computerized model, which we built using the visual language of statecharts [19] to investigate the role of the CPG, sensory feedback, and the interactions between them in controlling touch-induced vibrissae motion. Several biological systems have proven to be highly suitable for statechart-based modeling [20][21][22][23]. Here we present the first statechart model of a neural system. Using statecharts, we modeled the brainstem loop of the whisking system, including the whisker-follicle complex, the primary afferents, the trigeminal sensory nuclei, and the facial nucleus. The main advantage of such a statechart model is its object-based modeling approach, which allows a simple and unbiased analysis of the studied system (see Modeling approach and Model behavior in statecharts, under Materials and Methods). In this model, periodic whisking is assumed to be induced by a tri-phasic CPG (see Materials and Methods). As our interest was in low-level sensory modulations of CPG activity, top-down modulations of the CPG and of the brainstem loop are not modeled.
Our model provides a possible mechanism for TIP generation, solely dependent upon sensory modulation of CPG activity at the level of the brainstem. Furthermore, the model provided several predictions regarding touch-induced vibrissae movements, for which we herein provide empirical validation. In a wider context, our results can be generalized to help in the understanding of CPG-sensory feedback interactions in any adaptive, periodic behavior.

Materials and Methods
The model Modeling approach. This study was performed using a computerized model implemented in the visual language of statecharts [19], using the Rhapsody tool [24]. This language provides the means required to describe the complex behavior of a system in an understandable and unbiased form [25]. Using a bottom-up approach, the behavior of the system is defined through the behavior of its various components. The relatively simple behaviors of individual components are described separately. These are then assembled in a dynamic fashion, without regard to how they will work in the assembly, to give rise to more complex, higher-level behaviors. Thus we get non-predefined, dynamic, complex behaviors of the system, which can be analyzed.
The modeling and analysis tools: Statecharts, Rhapsody, Matlab and Visual Studio. Using the object-oriented variant of the visual language of statecharts [19] described in [24], a statebased hierarchical and concurrent transition diagram (statechart) was generated to graphically describe the discrete behavior of each component of the modeled brainstem loop, via the use of (1) states, (2) events that cause transitions between states, and (3) actions that generate events (and are transmitted from one component to another). In statecharts, states can be nested inside other states, creating sub-states, which enables description at multiple levels. States may also be divided into orthogonal states, thus modeling concurrency, and allowing each component to reside simultaneously in several independent states.
Using the Rhapsody tool [24], the behavior defined by the various statecharts was executed collectively, giving rise to a dynamic simulation that allowed interactions between the different objects (via actions), responded to events and changed objects' states and parameters during run time.
During model execution, a server-client interface was used in order to allow the model (written in Java code) to frequently use a Matlab function that transformed motor neurons' firing rate to whisker motion (see Muscle forces & whisker motion in File S1).
Analysis of the dynamics of the resulting integrated behavior was done using a Visual Studio application that provided a 2D animation to help visualize vibrissae motion.
Model specifications. The model consisted of six types of elements that participate in composing the brainstem loop: 1. Neurons: (1) Primary afferents (SN1s), which directly innervate the whisker follicles and convey raw sensory input to the brain. SN1 cell bodies are located in the trigeminal ganglion (TG).
(2) Secondary afferents (SN2s), which receive synaptic input from SN1s and project to the motor neurons. SN2 cell bodies are located in the trigeminal nuclei (TN).
(3) Efferent motor neurons (MNs), which receive synaptic inputs from SN2 cells and project to mystacial pad muscles. MN cell bodies are located in the lateral facial nucleus (FN).
2. Mystacial pad muscles, which were innervated by MNs and responded to the motor output by contraction, which moved the whiskers attached to them. This group of elements was divided into three subgroups, according to both anatomical and functional properties of the different types of mystacial pad muscles in the rat: (1) Intrinsic muscles (Int_muscle), which were innervated by Int_MNs and whose contraction resulted in an active forward movement of the whiskers attached to them. (2) Extrinsic protracting muscles (ExtP_muscle), which were innervated by ExtP_MNs and whose contraction resulted in an active forward movement of the whiskers attached to them, and in a forward pad translation. (3) Extrinsic retracting muscles (ExtR_muscle), which were innervated by ExtR_MNs and whose contraction resulted in an active backward movement of the whiskers attached to them, and in a backward pad translation.
5. An obstacle, which, if present in the whiskers' sweeping range, induced whisker-obstacle contacts. 6. A manager, which acted as an external environment object that passed information between objects and supported technical issues. This component did not simulate any biological component directly.
Model assumptions. Five rows of whiskers located at one side of the snout were modeled. Two rows were implemented as four-whisker rows, representing rows A-B in the rat, and the other three as seven-whisker rows, representing rows C-E ( Figure 1A).
In both types of rows, each whisker was attached to several muscles [31] ( Figure 1B): (1) Two intrinsic muscles: one rostral and one caudal. In all rows (A-E), the most rostral whisker was attached to only one intrinsic muscle, caudally, as observed in the rat [31].
(2) One superficial extrinsic protracting muscle. In rows A-B, the most rostral whisker was attached to an additional deep extrinsic protracting muscle -a pseudo-intrinsic muscle [31]. (3) Two extrinsic retracting muscles: a superficial muscle and a deep muscle.
The extrinsic muscles in the model and their equivalents in the rat can be found in Table 1. A combined, tri-phasic activation of the above muscles simulated free-air whisking motion ( Figure 2) [30]: First, extrinsic protracting muscles were activated, pulling the pad forward and initiating whisker protraction (Phase 1). Second, the intrinsic and pseudo-intrinsic muscles were activated, further moving the vibrissae forward (Phase 2). Third, relaxation of protracting muscles (both extrinsic and intrinsic) occurred, while extrinsic retracting muscles were activated, pulling the pad and the whiskers backward (Phase 3). Phases 1 and 2 gave rise to a forward whisker Each whisker (black) is attached to two intrinsic muscles (blue), two extrinsic retracting muscles (red), and one extrinsic protracting muscle (green). In rows A-B, the rostral intrinsic muscle of the most rostral whisker is replaced by a pseudo-intrinsic muscle (light green), whereas in rows C-E this muscle is removed. Intrinsic muscles and whiskers are indexed in an increasing order, starting at 1, from caudal to rostral. doi:10.1371/journal.pone.0079831.g001 motion (protraction), which was followed by a backward motion (retraction), evoked by phase 3.
The timed activation of the different muscles was controlled by the coordinated firing patterns of the different CPGs, exciting the different types of MNs at the appropriate times. Each CPG stimulated all the MNs of a certain type, which together innervated all the modeled whiskers ( Figure 2B). Thus, the activation of a CPG affected the motion of all the ipsilateral whiskers in a same manner. In this model, each CPG was part of an open loop and thus its behavior was not regulated. Free-air whisking motion induced by the behavior of the CPGs is further described in the ''Model behavior in Statecharts'' section.
CPG-induced MNs activity, and thus whisker motion, could be modulated by another pre-synaptic source, SN2s, which took part in a feedback loop. Separated sensorimotor feedback loops were implemented for the different whiskers, to allow each whisker to affect its own motion ( Figure 2B). Feedback loops implementation is described in more detail in Figure 3 and in File S1.
Although each whisker was innervated by a separate pool of neurons, sensory feedback derived from one whisker could affect the motion of other whiskers, via muscles that connected the whiskers. The extrinsic muscles (either protracting or retracting) were assembled into groups that affected the motion of several rows of whiskers, as described in Table 1 [32]. For example, when MNs in row C activated the superficial extrinsic protracting muscle of row C (part of the PMI complex in Table 1), all three superficial extrinsic protracting muscles in rows C-E contracted, moving all whiskers in rows C-E. In addition, sensory feedback that activated the intrinsic muscles attached to the innervated whisker, affected the movement of that whisker together with its two neighboring whiskers, both the caudal and the rostral whisker (if they existed). The whiskers within a row were assumed to move within a single plane, without rotating about their own axis. For simplicity, the whiskers in the model were considered rigid.
Model behavior in statecharts. To illustrate a model simulation, we describe here the dynamic process of generating free-air whisking motion. In the model, whisking motion resulted from the dynamically combined behavior of four types of elements: a CPG, a MN, a muscle, and a whisker. We start by describing the behavior of each of these four elements separately, and then describe how their combined behaviors gave rise to whisking motion. Similar descriptions of other model elements can be found in File S1 (see Model behavior in statecharts) and in Figure S2.
The CPG statechart, whose transition diagram is shown in Figure 4A, describes the behavior of a CPG. The CPG has two states, and at any moment it can be in exactly one state: (1) ''Activate'' -a stimulation period, during which the CPG stimulates all MNs of the corresponding type simultaneously, and at a constant rate, and (2) ''Relax'' -a ''silent'' period, during which no stimuli are sent to the MNs. The CPG can move between states following the occurrence of a triggering event specified next to the transition arrow. In the model, the CPG stays several tens of milliseconds in each state, and continuously switches between the two states upon the occurrence of a timeout (tm) trigger.
The MN statechart has four states ( Figure 4B): (1) ''Rest'' -in which the MN is inactive, (2) ''Generate action potential'' (generateAP) -in which the MN prepares to evoke an action potential (AP) and to transmit the signal to its post-synaptic muscle/s, (3) ''Absolute refractory period'' (ARP) -in which the MN is inactive, but in comparison to the ''Rest'' state, it cannot respond to any stimulus it receives. Thus, when being in this state the cell can never evoke an AP, (4) ''Relative refractory period'' (RRP) -which is very similar to the ''Rest'' state, during which the MN is inactive, but from which it can become activated. The difference between this state and the ''Rest'' state is that here the threshold is higher; i.e., the stimulus should be stronger.
Since statecharts are fully executable, one can capture the dynamics of the modeled behaviors by executing the model. When the simulation is initiated, many components are created for each element, with each component receiving its own copy of its parent's initialized statechart and initialized set of parameters, which together describe the component's independent behavior. The number of components in the model and a full list of model parameters are specified in Figure S1 and in File S1. The simultaneous execution of the statecharts of all the model's components resulted in a combined dynamic behavior of the system, as follow: Upon model execution, three copies of the CPG statechart were generated, one of each type (a CPG parameter): extrinsic protractor (ExtP_CPG), intrinsic (Int_CPG), and extrinsic retractor (ExtR_CPG). Each CPG started in the ''Activate'' state and its currentCycle parameter was assigned to zero. After a few milliseconds (cycleDuration), the CPG exited the ''Activate'' state and reached a condition connector, from which it could move to either the ''Activate'' state or the ''Relax'' state, depending on the fulfillment of a specific condition: If the value of currentCycle was lower than the constant cyclesNum, another parameter of the CPG, the CPG transmitted a stimulation event to all of its MNs (of the corresponding type) and re-entered the ''Activate'' state (and incremented its currentCycle), otherwise it moved to the ''Relax'' state. The condition held for several iterations, keeping the CPG in the ''Activate'' state for several tens of milliseconds. When the condition was no longer valid, the CPG moved to the ''Relax'' state. After a tm of several tens of milliseconds (silenceDuration) the CPG returned to the ''Activate'' state (and reset its currentCycle).
While in the ''Activate'' state, each MN innervated by the CPG received successive stimulation events. Upon the first stimulation, the MN exited the initial ''Rest'' state and entered the ''GenerateAP'' state in which it prepared to fire. After a few milliseconds (t relay ), the MN exited the state, transmitted a stimulation event to its post-synaptic target/s, and entered the ''ARP'' state. This period also lasted several milliseconds (ARP), during which the MN could not respond to any stimulus. Next, the cell entered the ''RRP'' state and stayed there for a few tens of milliseconds (RRP). If stimulated during the RRP period, it moved to the ''GenerateAP'' state, and otherwise, it returned to ''Rest''. Many MNs of a certain type innervated a single muscle of the corresponding type. Upon the first stimulation event transmitted to the muscle, the muscle exited the initial ''Rest'' state and entered the ''Contract'' state. The stimulus affected several muscle parameters, including its cytoplasmic Ca +2 concentration (see Muscle's Ca(t) in File S1). The muscle stayed in the ''Contract'' state until its Ca(t) decreased below a certain threshold (Ca 0 ). If while in the ''Contract'' state the muscle received another stimulus, it re-entered the ''Contract'' state with the updated parameter values.
Every time the muscle entered the ''Contract'' state it transmitted a move event to its whisker(s). Every whisker can receive a move event from several muscles, and upon the first stimulation event, the whisker exited its initial ''Rest'' state and entered the ''Move'' state, in which it stayed until its angle and angular velocity returned to their resting values (i.e., h t ð Þ~h 0 and _ h h t ð Þ~_ h h 0 ). The angle and angular velocity were constantly updated each time a muscle that was attached to the whisker transmitted a move event. The muscles' force magnitude (A) and the resulting whisker motion were calculated using a Matlab function (see Muscle forces & whisker motion in File S1).
The three types of CPGs were activated during different time periods along the whisk cycle (see CPGs' active period in File S1), The three phases of whisking are generated by a tri-phasic activation of three different groups of muscles: extrinsic protractor (Phase 1, green), intrinsic (Phase 2, blue), and extrinsic retractor (Phase 3, red) muscles. Phases 1 and 2 result in a forward whisker motion (protraction); phase 3 yields a backward motion (retraction). Note that in four-whisker rows the rostral intrinsic muscle of the most rostral whisker is replaced by a pseudo-intrinsic muscle, and that in seven-whisker rows (not shown) it does not exist. The CPG-induced tri-phasic activation may be modulated by sensory feedback. doi:10.1371/journal.pone.0079831.g002 The complexity of the dynamic behavior of the system described here makes it very suitable for statechart modeling. This complexity stems from highly concurrent and time-intensive changes of the system, due to constant interactions of the (hundreds of) system's components with each other. Such complex behavior can be analyzed when modeled by statecharts, using the kinds of execution and analysis techniques offered by the Rhapsody tool [24]. The resulting dynamic behavior of the system can be quantified via the emerging whisker motion. Comparing the simulated whisker motion to that observed in the rat allows testing the model.
Model configurations. In this study, three hypothesized whisking-inducing mechanisms were inspected, differing from each other in the pre-synaptic sources that innervated the MNs during free-air whisking ( Figure 3B-D). Each type of MNs was innervated by: (1) a CPG of the corresponding type, (2) by SN2_Ws, or (3) by both a CPG of the corresponding type and SN2_Ws. Free-air whisking described in the last section, ''Model behavior in statcharts'', was induced by mechanism (1).
In all three configurations, the cumulative stimulus that a given MN received had to cross a threshold in order for the MN to successfully fire. In model configuration (1) or (2), stimuli by the CPG or SN2s alone (respectively) allowed the MNs to cross threshold and fire. In contrast, in model configuration (3), only a mutual stimulation by both SN2s and the CPG allowed the MNs to cross the threshold and successfully fire. The tuning of synaptic strengths required for each configuration is assumed to be accomplished during development by behaviorally-controlled local learning rules [33][34][35][36].
Independent of the hypothesized whisking-inducing mechanism, this study also examined three hypothesized neural mechanisms for TIP-generation ( Figure 3E

Animals
Whisker movements obtained by the model were compared to those measured in head-fixed rats by Deutsch et al. [11], and to whisker movements measured here in freely-moving animals. In the following sections we describe the methods used in this study for collecting data from behaving rats.
The whisking patterns of albino Wistar rats (n = 3) were measured. All whiskers were trimmed, except for one row of whiskers (C) on each side of the snout. This configuration was chosen in order to simplify the tracking of the whiskers in postprocessing (see below). Trimmed whiskers were clipped close (,1 mm) to the skin during Dormitor anesthesia (0.05 ml/100 g, S.C.) the day prior to behavioral recording. Recordings were performed prior to and following transectioning of the infraorbital branch of the trigeminal maxillary nerve (IoN).

Ethics statement
Animal maintenance and all experimental procedures were conducted in accordance with the guidelines of the National Institutes of Health (USA) and The Weizmann Institute of Science. The protocol was approved by the Council for Experiments on Animals of the Weizmann Institute of Science (Application Number: 01260212-2). Surgeries were performed under medetomidine hydrocholoride and ketamine anesthesia. Efforts were made to minimize suffering, including analgesiac injections post-operatively (Rimadyl). After the experiment was completed, animals were sacrificed using an overdose of barbituates (Pentobarbital).

Transectioning of the IoN
Before surgery, rats were removed from their cage and anesthetized (medetomidine hydrocholoride, 0.05 ml/100 g, S.C.; ketamine 79 mg/kg, S.C.). The rats were then mounted in a stereotaxic device. An incision was made caudal to the whisking pad. The IoN was exposed where it emerges from the eye socket. A 2 mm section of the nerve was cut out in order to prevent nerve regeneration. The wound was then sutured. Post-operatively, rats were given antibiotics (penicillin and streptomycin; Pen-Strep, 2 ml/kg, sc), analgesiac (Rimadyl, 5 mg/kg, SC in 1 ml saline), and ad libitum food and water.

Experimental apparatus
Behavioral experiments was performed in a darkened, quiet room. The behavioral apparatus consisted of a holding cage (25 cm width, 35 cm length, 29.5 cm height), with a small door (6.9 cm height, 6 cm width), through which the rats could emerge into the experimental area ( Figure 5). Both the holding cage and the experimental area were held approximately 15 cm above the surface of a table. The experimental area consisted of a glass plate with 1-2 objects (metal poles or plexiglass cubes and cylinders) placed on the plexiglass. The location of the objects was changed between experiments, and the two obstacles were far enough from each other to prevent simultaneous touch of both obstacles. The glass plate was back-lit by an IR-lamp (880 nm wavelength, 23 cm623 cm, Metaphase, USA). The experimental area was filmed from above by a high-speed, high-resolution camera (128061024 pix, 500 fps, CL60062, Optronics, East Musckogee, OK). An in-house program, written by Dr. Enrico Segre, triggered the high-speed camera whenever the rat emerged from the holding cage into the experimental area. Video recording stopped when the rat returned to the holding cage.

Behavioral task
An experimental session consisted of recording an animal's whisking behavior during an emergence task, over a period of 30-120 minutes. Preceding a behavioral session, the animal was removed from its cage and placed in the experimental holding cage for 15 minutes. During the acclimation period, the door in the holding cage was blocked. The experimental period began with unblocking the door in the holding cage to allow the animal to leave the holding cage and explore the experimental area at will. The experimental period varied, depending on the animal's behavior and the amount of recorded video. The experimenter's interference and contact with the animal was minimized during the experimental session and there was no interference or contact while the animal was in the experimental area.

Video analysis
Video analysis was performed using the BIOTACT Whisker Tracking Tool [http://bwtt.sourceforge.net ; 37]. This tracker provides the head position and head angle in each video frame through calculations based on the location of the rat's nose tip and the center of its snout. Additionally, the tracker identifies the rats' whiskers and provides the base angle of each of the trackeridentified whiskers with respect to the mid-saggital plane of the rat. The number of tracked whiskers varied between sessions due to differences in the whiskers' length and changes in focus owing to rat's movement.

TIP analysis
Whisker movements were inspected both qualitatively (by eye) and quantitatively, by a computerized whisker tracker, in order to identify whiskers that ''pumped'' in response to touch. The identity of the touching whisker(s) varied between the different TIP events.
(1) Qualitative inspection: TIP events and the identity of the ''TIPing'' whiskers were identified manually, by three independent observers. Since most whiskers were trimmed, individual whiskers were easily identified. Only TIPs that were detected by at least two observers were counted. The observers were in agreement on the occurrence of TIPs and (2) Quantitative inspection: All manually counted TIPs were also analyzed using the Matlab-based WhiskerTracker image processing software [37]. This allowed a better inspection of individual whiskers' movement for determining the occurrence of TIPs.

Results
(1) Sensory feedback is required for TIPs generation in naturally behaving rats Sensory feedback was indicated to be necessary for TIP generation in head-fixed rats [11]. Here we verify in freely moving rats that TIPs indeed depend on sensory feedback, by examining TIP occurrence in behaving rats whose infraorbital branch of the trigeminal maxillary nerve (IoN) was cut bilaterally. We postulated that in IoN-transected rats, where transmission of sensory information from the vibrissae to the brain was eliminated, TIPs would be absent if sensory feedback is indeed required.
Rats in which all whiskers were trimmed besides row C (bilaterally) were filmed using a high-speed video camera. Twelve movies with durations of 0.67-4 minutes were acquired before (seven) and after (five) bilateral IoN transection. In contrast to a TIP occurrence of 27% (34 TIPs out of 126 touch events) in response to whisker-obstacle contacts in intact rats, only 3% (2/61) of touch events yielded TIPs in IoN-cut rats (p = 0.028, one-tailed binomial test). Thus, bilateral IoN transection brings TIP occurrence rate to its spontaneous level (7% [11], p = 0.19, binomial test). Besides the significant reduction in TIP probability, IoN transection did not affect significantly characteristics of whisking in air as measured via cycle duration, protraction duration, whisking set-point or whisking amplitude (Dependent ttest for paired samples, n = 28; p = 0.161, p = 0.869, p = 0.125, p = 0.053, respectively). These results confirm that sensory feedback is required for TIP generation in naturally behaving rats while not significantly affecting parameters of whisking in air.
(2) Establishing model credibility Before inspecting possible control mechanisms that could modulate whisking motion in response to touch, we first established a model that mimicked the motion to be regulated, i.e., free-air whisking. In the rat, whisking is assumed to result from tri-phasic activation of mystacial pad musculature [30]. We implemented all vibrissal muscles that participate in generating free-air whisking [32]: extrinsic protracting, intrinsic, and extrinsic retracting muscles (see Model specifications under Materials and Methods). The sequential contraction of these three types of muscles, timed based on Fisher et al. [25], resulted in a periodic motion of all mystacial pad whiskers back and forth in a sweeping motion (Figure 2). Simulated movements obtained by the model were qualitatively and quantitatively compared with corresponding motions observed in behaving rats in order to fine-tune the mimicry of exploratory whisking motion.
Accurate mimicry of free-air whisking motion. During free-air whisking, rats repeatedly move their whiskers along the rostral-caudal axis in order to scan their immediate environment [1,4,38]. The large vibrissae on each side of the rat's snout are arranged in a grid of five rows and several (4-7) arcs ( Figure 1A). Although each whisker has some capability for independent movement, the whiskers on each side of the snout generally move together [10,11]. Thus, in principle, the modeling of an individual whisker would be sufficient to describe the motion trajectories of all ipsilateral whiskers. Yet, we modeled the entire array of ipsilateral whiskers, in order to later on investigate the effect of contact made by one whisker on the motion trajectories of other whiskers located in different rows and arcs. We start by describing the motion in a single row.
Our model generated whisking trajectories that can be mapped to the low-frequency, low-amplitude, and square-wave range of the in-vivo repertoire. The movement of the five rows of whiskers found on one side of the snout was simulated. Two rows were implemented as four-whisker rows, equivalent to rows A-B in the rat, and the other three as seven-whisker rows, equivalent to rows C-E (see Model assumptions under Materials and Methods). For simplicity, a synchronized movement within the whisker array was assumed [10]. We present here whisker motion of two representative rows, since, in the absence of contacts, similar motion was obtained in rows with a similar number of whiskers. Figure 6B,C displays the angle of whiskers found in two rows as a function of time. The graph shows a coordinated movement of the whiskers while moving back and forth. In the model, a whisk cycle lasted about 150 msec, and was composed of protraction and retraction phases that lasted about 80 and 70 msec, respectively, thus mimicking biological whisking at the low end of the frequency spectrum [11]. Interestingly, the squared shape of single whisk cycles generated in these low frequencies changed to a sinusoidal shape when the whisking frequency was increased (not shown). Whisking amplitude was set to be ,10u, at the lower end of the amplitude range observed in awake rats, in order to allow accurate calculation of whisking motion (see Muscle forces & whisker motion in File S1). From a resting angle of 70u, all whiskers reached a very similar maximal protracting angle (79-81u) except for the rostral-most whiskers in the seven-whisker rows, which only reached ,74.5u. This reduction in protraction amplitude resulted from the model's assumption (based on anatomical data) that in a seven-whisker row, the rostral-most whisker is not attached to any protracting muscle that is active during the second phase of whisking motion (see Model assumptions under Materials and Methods). This model prediction could not be tested in our rats due to insufficient whisker length for reliable tracking.
Consistency with experimental data. The results described above were obtained using a neural configuration in which whisking was generated solely by the activity of a CPG. Similar results were obtained when whisking was controlled by both a CPG and sensory feedback (data not shown). In this latter configuration, CPG activity could not activate the MNs alone. Rather, the MNs required simultaneous input from both the CPG and the SN2s (see Model configurations under Materials and Methods; Note that independently of these mechanisms, three hypothesized TIP-inducing mechanisms are examined). In contrast to these two configurations, a third configuration, in which whisking was hypothesized to be generated by the activation of MNs by SN2s alone, could not produce the desired whisking pattern (data not shown). These results are consistent with experimental data that show that a CPG is necessary and sufficient to induce free-air whisking.

(3) Model predictions
After establishing model's credibility, we examined possible control mechanisms responsible for TIPs. Deutsch et al. [11] demonstrated that, in head-fixed rats, TIPs occur due to sensory feedback which modulates whisking. They suggested that this sensory feedback occurs during the second phase of the whisk cycle (activation of the intrinsic and pseudo-intrinsic muscles), but not during the other phases. Briefly, they rejected the possibility of feedback during the first phase of whisking (activation of extrinsic protracting muscles), since this could only generate TIPs early in the whisk cycle, whereas TIPs typically occur later, near the peak of a whisking cycle. They also ruled out the possibility of feedback that modulates the backward motion generated during the third phase of whisking (activations of the extrinsic retracting muscles), since whisker velocity measured at the beginning of TIP retraction was much smaller than that at the beginning of end-of-whisk retraction, indicating that the negative velocity in the TIP is unlikely to be the beginning of the ''end-of-whisk'' retraction. In addition, Deutsch et al. [11] proposed that TIP-inducing sensory feedback is relayed along an anatomically-short neural loop, which would allow for the short latency of TIPs in response to touch.
Guided by these assumptions, we aimed at examining all possible sensory feedback mechanisms that could be mediated by the anatomically shortest neural loop that connects the whiskers to the mystacial muscles, the brainstem loop, in order to induce TIPs during the second phase of the whisk cycle. These mechanisms included: (1) ''E-P'', excitation of protracting MNs (ExtP_MNs and/or Int_MNs), which would produce a short-term, enlarged protraction, (2) ''I-P'', inhibition of the already-active second phase protracting components -Int_MNs or their pre-synaptic CPG, which would result in two discrete, non-overlapping phases of whisker protraction, or (3) ''E-R'', excitation of retracting MNs (ExtR_MNs), which would interrupt the ongoing whisker protraction, and would possibly result in a reversal of movement direction. The E-P mechanism was less probable, since it would only allow producing TIPs early in protraction [11 ( figure S4A)], whereas most TIPs occur in the late part of protraction, as mentioned. In contrast, the I-P and E-R mechanisms would allow generation of TIPs during the entire second phase of the whiskcycle. These two remaining TIP-inducing mechanisms were examined here, where in the I-P mechanism we further distinguished between the two possible inhibition routes: direct inhibition of Int_MNs by trigeminal-ganglion-to-facial-nuclei inhibitory projections, ''direct I-P'' [42], and indirect inhibition of Int_MNs via Int_CPG inhibition, ''indirect I-P''. Accurate mimicry of TIPs and further characterization of this motion by model predictions. Frequently, when a rat's whisker touches an obstacle it does not simply keep on protracting throughout the whole protraction phase, but, rather, ''pumps'' shortly (tens of milliseconds) after the beginning of touch. During the pump, the whisker briefly protracts while bending, and then a retraction of about 1u occurs, followed by another protraction ( Figure 7A). During retraction, the whisker usually remains touching the object without detaching from it, by gradually straightening up [11].
A similar movement pattern was obtained by the different TIPinducing model configurations. Figure 7B displays the angles of four whiskers found in a single row, as a function of time (similar results were obtained for seven-whisker rows, data not shown). In the figure, the touching whisker started to retract 17 msec after the beginning of touch, moving 0.6u backwards, while the other whiskers in its row retracted by about 1u. This first TIP was followed by a second TIP (see Discussion). In this example, TIPs were induced by the E-R mechanism.
As in free-air whisking, we inspected the accuracy of the simulated TIP movements obtained by the different mechanisms, by comparing the values of several quantitative characteristics of a TIP movement to the matching values measured in the rat. Our analysis only refers to the first TIP observed in each whisk cycle, since in the rat only one TIP per cycle is usually observed (see TIP occurrence, below). The results are summarized in Table 2.
TIP delay. In Figure 7B, the 17 msec delay in TIP onset time, relative to touch onset time, was obtained for whiskerobstacle contacts at a radial distance of 40% of the touching whisker's length (measured from the base of the whisker). This delay was obtained by all three model configurations for the specified radial distance, and its value was consistent with the mean delay observed in head-fixed rats (18.165.8 msec, mean6SD), in response to contact at a similar distance [11].
However, in the model, this delay varied for different radial distances of contact. All three model configurations predicted dependency between the delay in TIP onset time and the radial distance of contact: as the contact point became closer to the base of the whisker, the delay decreased (Figure 8, black). Analysis of experimental measurements from head-fixed rats revealed that such an increase in TIP delay was significant (Figure 8, red, p = 0.001). The model prediction was based on the finding that the delay in pressure cells firing (SN1_P) increases as the radial distance of contact increases [43].
TIP delays in the model were then fine-tuned to match the delays observed in head-fixed rats, which varied between 12-30 msec. This range is in agreement with experimental data, where 95% of TIPs in the rat start 12-30 msec after contact onset time.
TIP amplitude of the touching whisker. In contrast to the high compatibility of the delays in TIPs between all three model configurations and the rat, TIP retraction amplitude of the touching whisker obtained by only two mechanisms matched the values measured in head-fixed rats (about 1u, in response to contacts at a radial distance of 40%): E-R or direct I-P resulted in a TIP retraction amplitude of 0.5-0.75u or 2-2.5u, respectively, whereas indirect I-P yielded a too-large retraction amplitude of 5u. The difference in retraction amplitudes between the E-R vs. I-P (either direct or indirect) configurations stemmed from the different magnitudes of forces applied by the different types of muscles involved in TIP induction. In the first modeling stage, in which model's credibility was established, force magnitudes of the modeled muscles were set to certain values that yielded an accurate mimicry of free-air whisking amplitude. Accordingly, intrinsic muscles were set to apply larger force than the extrinsic retracting muscles, in agreement with experimental data [13 (figure 2B),30 ( figure 7C,E)]. Hence, I-P (either direct or indirect) resulted in a larger whisker retraction than E-R. The difference in Figure 7. Comparison of TIPs between the rat and the model. (A) Example of whisking trajectory against an object of a single whisker (C2) in a head-fixed rat (filtered at 80 Hz). Whisker-object contact is indicated by bold. This example is a zoom-in of Figure 1D in [11]. (B) Simulated TIPs in a single row as predicted by the ''excitation of retractor MNs'' configuration. While only whisker 4 touches an obstacle, all row's whiskers pump. At touch onset time, the touching whisker's angle is 74.7u, and the radial distance of contact is 40% of the whisker's length. TIP delay, relative to touch onset time: 17 msec; TIP retraction amplitude: ,1u for all the non-touching whiskers, 0.6u for the touching whisker. doi:10.1371/journal.pone.0079831.g007 retraction amplitudes between the direct vs. indirect I-P configurations stemmed from different number of intrinsic muscles that were involved. In the direct I-P, sensory feedback derived from the touching whisker resulted in the relaxation of only two intrinsic muscles that were attached to it, whereas in the indirect I-P, sensory feedback led to the relaxation of all ipsilateral intrinsic muscles, amplifying the effect (see Model assumptions under Materials and Methods). We emphasize that neural response times to whisker detachment were similar in all three model configurations, rejecting the possibility that larger TIP retraction amplitudes were obtained due to longer lags in response to whisker detachment.
Spatial spread of TIP across whiskers. The different model configurations also differed in the predicted identity of the whiskers that would ''pump'' in response to the contact of a single whisker with an object. The direct I-P configuration predicted that in addition to the touching whisker, both of its caudal and rostral neighbors (if they existed) would also ''pump'' ( Figure 9A). This resulted from the assumption that sensory feedback derived from the touching whisker leads to the relaxation of the two intrinsic muscles attached to it, as each muscle is also attached to one of the touching whisker's neighbors (see Model assumptions under Materials and Methods). In contrast, the other two model configurations predicted that all whiskers found on the same side of the snout would ''pump'' ( Figure 9B,C). In the E-R configuration, the prediction resulted from the assumption that sensory feedback derived from the touching whisker excites the superficial extrinsic retracting muscles, which are attached to all ipsilateral whiskers. Alternatively, in the indirect I-P configuration, the prediction resulted from the assumption that the sensory feedback derived from the touching whisker inhibits the Int_CPG, and thus leads to the relaxation of all ipsilateral intrinsic muscles (see Model assumptions under Materials and Methods).
Following model predictions, we examined the identity of the ''pumping'' whiskers in behaving rats. Whisking rats in which all whiskers were trimmed except for row C (bilaterally, 5-7 whiskers on each side) were filmed using a high-speed video camera. Whisker movements were inspected both qualitatively, by three independent observers, and quantitatively, by a computerized whisker tracker, to identify TIP events and the participating ''pumping'' whiskers (see Materials and Methods). No rats with fully intact whisker pads were inspected, as whisker tracking and identification was not feasible with full whisker pads by any of the inspection strategies.
In 53% of the detected TIP events, all the traceable whiskers in the row ''pumped'' (18 of 34 studied TIPs). Out of the remaining 16 TIP events, in 11 cases 1-2 whiskers slowed down but did not exhibit detectable retraction; the other 3-5 whiskers exhibited clear TIPs. In the other five cases, 1-2 whiskers were either out of focus or too short to observe a TIP; the other 3-5 whiskers exhibited clear TIPs. The number and identity of the touching and ''pumping'' whisker(s), as well as the number of traceable whiskers, varied between the different TIP events, and are summarized in Table 3 (whisker identity is not specified). An example of a TIP observed in all tracked whiskers, both touching and non-touching, is displayed in Figure 10. All whiskers found on the same side of the snout as the touching whisker All inspected whiskers found on the same side of the snout as the touching whisker** The presented values, of both the simulations and the rat, are of TIPs induced by contacts at a radial distance of 40% from the base of the touching whisker. Values of simulated TIPs that were induced by contacts at a wide range of radial distances (20-70%) are displayed in brackets. For each TIP-inducing mechanism, values were obtained while using either one of two whisking-evoking mechanisms: generation of whisking by CPG alone or by both the CPG and sensory feedback. *,10u for all non-touching whiskers in rows A-E, except for the most rostral whiskers in rows C-E, which retracted by ,5u. **Observed in both head-fixed and behaving rats. In head-fixed rats, TIP was observed in the two additional non-touching whiskers that where examined [11]; in behaving rats TIP was observed in all the whiskers found on the same row as the touching whisker. doi:10.1371/journal.pone.0079831.t002 Figure 8. Dependency of TIP delay on radial distance of contact in head-fixed rats. Both in the model (black) and in the rat (red graph and table, p = 0.001), the mean TIP delay elongated as the radial distance from the base of the whisker increased. The radial distance was divided by the length of the touching whisker. TIP delay is relative to onset time of whisker-object contact. doi:10.1371/journal.pone.0079831.g008  Our results are consistent with those of Deutsch et al. [11], who tracked the simultaneous movements of three whiskers (C1, C2, and D1) bilaterally in head-fixed rats. Deutsch et al. found that in addition to TIPs in the touching whisker (C2), the other two ipsilateral whiskers usually also ''pump''.
TIP amplitude of the non-touching whiskers. Not only did the identity of the non-touching ''pumping'' whiskers differed between the different model configurations, but the predicted TIP retraction amplitude of these whiskers differed as well. Following contact of a single whisker with an object at a radial distance of 40%, the E-R configuration predicted retraction amplitude of 1-2u, direct I-P predicted 3-3.5u, and the indirect I-P configuration predicted 10u. These differences resulted from similar considerations as those for the touching whiskers (see TIP amplitude of the touching whisker under Results).
Comparison of the predicted values with those measured in head-fixed rats pointed at one feedback mechanism that best fit experimental data. Deutsch et al. [11] measured a retraction amplitude of about 2u of a non-touching whisker, C1, in response to contacts at a radial distance of 40%, in accordance with the E-R configuration.

Discussion
This study examines the roles of CPG, sensor-derived sensory feedback, and CPG-sensory feedback interactions in shaping periodic sensor movements. As a case study, we examined a wellcharacterized vibrissae movement that occurs following vibrissaobject contact (touch-induced pump, TIP), and which is assumed to result from the modulation of whisking-inducing CPG activity by vibrissa-derived sensory feedback. Using a dynamic, objectbased computerized model built with the language of statecharts, and following Deutsch et al. [11], we examined three possible feedback mechanisms that may generate TIPs, and which pass along the brainstem loop, one excitatory and two inhibitory. Comparison of simulated and real TIP movements ruled out all hypothesized mechanisms but one: excitation of retractor motor neurons (E-R), which consistently and accurately simulated the corresponding motion observed in head-fixed rats [11] (Table 2). Moreover, several model predictions that further characterized TIP motion were empirically validated in head-fixed and freely-moving rats (Figures 8 and 10, respectively), further supporting the proposed mechanism. Previous electrophysiological studies provide additional support, as only excitatory neuronal feedback has been shown to occur in the brainstem loop in vivo [44]. These results can be generalized to help understand the nature of CPGsensory feedback interactions in other adaptive periodic behaviors.

Statecharts modeling
Computational models provide helpful means for understanding how behavior (e.g., whisker motion) emerges from the activity of a network of neurons, by allowing the experimenter to reproduce complex dynamic activities which can then be analyzed and explained. While both object-based and equation-based models accommodate this, the former have one major advantage in biological modeling [45]: Object-based models allow both qualitative (states) and quantitative description of biological behaviors, offering an alternative if precise quantitative biological relationships are unknown, if too many variables are involved, or if they change over time, depending on certain events. Among object-based modeling methods, the visual language of statecharts [19] provides a major practical advantage. Using the Rhapsody tool [24], statecharts is accessible through very intuitive graphical interfaces, allowing experimentalists to easily define behaviors for improving their understanding of complex biological phenomena.

The E-R brainstem loop
The E-R mechanism is similar to the flexor (pain) reflex, where a bi-synaptic excitatory feedback activates antagonist muscles to oppose the activity of the agonist muscle [46]. The flexor reflex is usually referred to as a protection reflex whose role is to pull the body away from damaging stimuli. TIPs, in contrast, lead to a more frequent object palpation, probably for a better estimation of object's features. An interesting possibility that is raised by our results is that flexor 'reflex' loops may also serve perceptual functions, in cases where stimulus intensity is not painful.
Alternatively, TIPs can be considered as 'fixational whisker movements', as they result in the whisker centering, or 'fixating', on a point of interest. This can be compared with similar processes in other modalities, such as fixational eye movements in primates [47]. In both cases the scanning movements of the sensory organs involve large and small movements, where the large movements shift the scanning from one point of interest to another and the small ones explore a given point of interest. It may thus be interesting to consider an E-R like mechanism for participating in the control of fixational eye movements, and a similar principle participating in the control of sniffing [31,48] or tasting of novel food [49].

Model predictions
Our model provided several predictions about vibrissae movements: 1. In free-air whisking simulations, the model predicted that the most rostral whiskers in seven-whisker rows (corresponding to rows C-E in the rat) would exhibit smaller protraction amplitude (,74.5u) than all other whiskers (79-81u) ( Figure 6B,C), since no muscles were attached rostral to these whiskers [31]. 2. In TIP simulations, when TIP was evoked by E-R, the model predicted a difference in TIP retraction amplitude between two groups of whiskers: those in rows A-B, and those in rows C-E. Whiskers in one group (either A-B or C-E) were hypothesized to retract by a larger extent than those in the other group, when the touching whisker was among them. This difference stemmed from the activation of deep extrinsic retracting muscles, which were separated into two different groups that were attached to either rows A-B or C-E (see Model assumptions under Materials and Methods, Table 1).
Testing these predictions is beyond the resolution of our current technology. Future inspection of these predictions by better means would allow us to further test the validity of the proposed model.

Beyond the brainstem loop: Top-down modulations
Although our brainstem-loop model mimics within-cycle dynamics of whisker motion very well, it is limited in its ability to capture longer timescale, across-cycles dynamics. Comparison of simulated TIPs with those measured in head-fixed rats demonstrates two such dynamics. First, in head-fixed rats, TIPs induced by contacts at a radial distance of 40% occur at about 25% of the whisk cycles, such that a whisker usually ''pumps'' once in a whisk cycle, and more rarely twice [11]. In contrast, in our model, TIPs occurred every cycle, for all radial distances of contact, and repeated one to two times per cycle, depending on the radial distance of contact (once for radial distance .60%, twice otherwise). Second, in head-fixed rats, contact duration is doubled during whisks that exhibit TIPs in comparison to whisks with contact but no TIPs [11]. In contrast, in our model, TIP occurrence did not affect contact duration (data not shown).
We postulate that these differences stem from top-down regulation, which may modulate the magnitude and timing of brainstem-derived responses significantly [50][51][52][53]. More specifically, TIP probability in the rat could result from top-down regulation based on the following findings: (1) In the rat, the occurrence of TIPs was associated with a more protracted set-point [4,11], a parameter that is believed to be controlled by higher brain areas than the brainstem -i.e., cortical [1,54], thalamic [55] or midbrain (e.g., superior colliculus; [56]) areas; (2) TIPs are temporally clustered, occurring in successive whisks rather than distributed randomly upon touch events [11]. These findings suggest top-down gating that switches between facilitating and suppressing the excitability of sensory neurons (in the model these would be touch cells SN2_C and SN2_P), resulting in promoted and obstructed TIP periods, respectively. Further assuming that these facilitating and suppressing processes are slow processes that gradually develop over several whisk cycles, then such top-down gating could also explain the occurrence of more than one TIP per cycle: in the first facilitating cycles, the touch cells would receive relatively weak stimuli that would allow them to fire only while at ''Rest'', during which their threshold is relatively low (see Figure 4B). These initial, weak stimuli would not allow the touch cells to cross the higher threshold while at ''relative refractory period'' state, possibly resulting in one TIP per cycle. In subsequent cycles, the stronger stimuli could allow these cells to also fire upon excitation from the ''RRP'' state, resulting in two TIPs per cycle, as observed by Deutsch et al. [11] and as shown in Figure 7A. The longer contact periods in the rat could result from top down regulation that resets the activity of the tri-phasic CPG in response to touch, which would elongate contact periods [16]. Thus, we speculate that integrating higherloops into our brainstem-loop model, as well as closing the CPGs' loop, would allow for the capturing of such slow modulations.
Our bottom-up modeling approach enables us to gradually add higher loops in a simple and natural manner. In this way, increasingly complicated behaviors of the system emerge step by step, thus enabling the explanation of each incremental change. As observed in this study, extending our model may not only explain familiar behaviors, but can also reveal new, unfamiliar behaviors of the studied system. We hope that as our model evolves we will obtain more insights that will further promote our understanding of the way the vibrissal system works. Figure S1 The number of elements that compose a single whisker's loop. Each whisker is innervated by a separated pool of primary afferents (SN1s), secondary afferents (SN2s) and motor efferents (MNs), which contain tens of neurons of several subtypes, as indicated in the scheme. For example, a single whisker is directly innervated by 162 SN1s which include 28 detach (D), 73 whisking (W), 33 contact (C), and 28 pressure (P) cells. Each type of SN1s innervates the corresponding type of SN2s, where a single SN2 is innervated by randomly chosen 50% SN1s of the corresponding type. Depending on the TIP-inducing configuration, different types of SN2s innervate different types of MNs (as indicated in figure 3B-G in the paper), with each MN innervated by all SN2s of the matched type. The ''E-R'' TIP-inducing configuration is displayed here. Each type of MNs innervates the corresponding type of muscle/s attached to the whisker (as indicated in figures 1B, 2B in the paper). In addition to this closed loop, all whiskers' MNs are innervated by the model CPGs, with each CPG innervating all MNs of the corresponding type. Note that no connections exist between sub-populations of neurons of a certain type (e.g., between whisking (SN1_W) and pressure (SN1_P) cells). * Two (instead of one) extrinsic protractor muscles are attached to the most rostral whisker in rows A-B. (TIF) Figure S2 The statechart of the Obstacle element. The behavior of the Obstacle is described in File S1. (TIF)

Supporting Information
File S1 Supporting information.