A system for tracking whisker kinematics and whisker shape in three dimensions.

Quantification of behaviour is essential for biology. Since the whisker system is a popular model, it is important to have methods for measuring whisker movements from behaving animals. Here, we developed a high-speed imaging system that measures whisker movements simultaneously from two vantage points. We developed a whisker tracker algorithm that automatically reconstructs 3D whisker information directly from the 'stereo' video data. The tracker is controlled via a Graphical User Interface that also allows user-friendly curation. The algorithm tracks whiskers, by fitting a 3D Bezier curve to the basal section of each target whisker. By using prior knowledge of natural whisker motion and natural whisker shape to constrain the fits and by minimising the number of fitted parameters, the algorithm is able to track multiple whiskers in parallel with low error rate. We used the output of the tracker to produce a 3D description of each tracked whisker, including its 3D orientation and 3D shape, as well as bending-related mechanical force. In conclusion, we present a non-invasive, automatic system to track whiskers in 3D from high-speed video, creating the opportunity for comprehensive 3D analysis of sensorimotor behaviour and its neural basis.

Introduction

Substantial progress towards the long-standing ambition of measuring “total movements made by the intact animal” [1] is coming from the application of machine vision methods to video recordings of behaving animals [2]. Since the whisker system is a major experimental model and since the whiskers are readily imageable [3,4], the whisker system is ideally suited to this endeavour. Tracking the whiskers of mice/rats has already deepened our understanding of active sensation and refined our capacity to relate behaviour to neural mechanisms [5-14]. Our aim here was to develop a method to track whisker movements and whisker shape in 3D in behaving mice at millisecond temporal resolution.

Both the movement of whiskers and the mechanical forces of whisker-object contact are 3D. During each whisking cycle, the whisker follicles translate with respect to the head and each whisker rotates in 3D. Although only the horizontal component of this movement is typically measured, whiskers also move vertically [15] and rotate around their longitudinal axes (‘roll’) [5]. The mechanical forces of whisker-object contact that are the drivers of neural activity are also 3D. When a mouse ‘whisks’ against an object, the whiskers bend. Again, although only the horizontal component of bending is typically measured, bending can occur in all directions [16]. In the trigeminal ganglion, all directions of deflection are represented [17-20], indicating that 3D bending information is both encoded and transmitted to the brain.

Starting with the first “cinematographic” study of whisking by Welker in 1964, there is a 50 year history of increasingly sophisticated efforts to measure whisker movement from behaving animals [4]. Most studies have measured whisker movement only in the horizontal plane, using either linear, Charge Coupled Device (CCD) arrays [21,22] or high-speed imaging [6,23-25]. However, horizontal plane imaging provides direct measurement of only one of the 3 angles that define 3D whisker orientation. Moreover, estimates of whisker bending moment obtained by imaging apparent curvature of a whisker in the horizontal plane [24,26,27] can be contaminated by roll [5]. This is significant, since bending moment is a primary driver of contact-related mechanotransduction [11,13,14]. High-speed cameras sufficient to form the basis of a 3D whisker imaging system have long been available: the main bottleneck to achieving 3D whisker tracking has been the computational complexity of the 3D reconstruction problem. A few studies have measured aspects of 3D whisker movement in vivo [5,15,16,28] and ex vivo [29], but no automatic approach has so far been developed that measures both 3D whisker orientation and 3D whisker shape from high-speed video of behaving animals. Here, we obtained 3D information noninvasively from head-fixed mice using a high-speed imaging system consisting of two cameras, positioned to minimise occlusion, and developed computational methods to reconstruct 3D whisker state. We used the system to track up to 8 whiskers in parallel, and to obtain a 3D description of each whisker, encompassing both its 3D orientation and 3D shape.

Results

3D imaging of whisking behaviour

To obtain a video data set with which to develop 3D whisker tracking, we trained head-fixed mice to detect objects with their whiskers (n = 6). On each trial, a vertical pole was presented in either an anterior location out of reach of the whiskers (‘no-go trial’) or a posterior location within reach (‘go trial’). Mice learned to perform the task accurately (81±17%, mean ± SD over mice) and performed 135±22 trials per daily session. When the pole moved up at the onset of a trial, mice would typically commence exploratory whisking. On go trials, one or more whiskers typically contacted the pole; on no-go trials, there was no contact. In this way, we obtained a varied data set, which included episodes of whisking both with and without contact (average behavioural session ~0.5M video frames).

We recorded high-speed video of mice using a system of two high-speed cameras (Fig 1). One camera imaged whisking in the horizontal plane. The other camera imaged in an off-vertical plane, the orientation of which (25⁰ off coronal, 10⁰ off horizontal) was optimised to minimise occlusion of whiskers against the background of the body of the mouse (Fig 1). For brevity, we refer to this second image plane as ‘vertical’.

Fig 1

Experimental set-up for 3D imaging.

A) Schematic showing the camera angles and 3D head-centred xyz coordinate frame. B) Horizontal and vertical views, with corresponding 2D coordinate frames.

Reconstructing whiskers in 3D from 2D views

Using the two-camera set-up, we imaged mice (1000 frames/s) as they performed the pole detection task. This resulted in a time series of image pairs (horizontal and vertical views): we refer to each such image pair as a ‘frame’. Our next aim was to develop an automated algorithm to track multiple whiskers in 3D. For two reasons, we focussed on the basal segment of the whisker shafts. First, during whisker-object contact, whiskers bend and the associated mechanical forces/moments drive mechanoreceptors located in the follicle. Thus, changes in whisker shape at the base of the whisker are intimately related to neural activity in the ascending whisker pathway [11,14]. Second, tracking only the basal segment (not the whole whisker) reduces the number of parameters required to describe the shape of a whisker. To accurately describe the shape across its entire length requires at least a 5th degree polynomial [27] which, in 3D, has 15 parameters. However, the basal segment of a whisker is well-approximated by a quadratic curve [24,29] which, in 3D, has 9 parameters.

The basic computational problem is to reconstruct the 3D coordinates of multiple whiskers from two 2D images. 3D whisker reconstruction is particularly challenging, since whiskers are similar to each other and their images can overlap. Also whiskers lack intrinsic features which can be matched across the two images. Knutsen et al overcame this by marking the shaft of one whisker with spots of dye [5]. The two images of each dye spot could then be matched and, after camera calibration, their 3D locations reconstructed. This is an elegant solution but there is the possibility, particularly for thin whiskers such as those of a mouse, that application of dye perturbs whisker mechanical properties. In a more recent, non-invasive approach, Huet et al achieved 3D tracking of a single whisker; tracking automatically in one image plane and tracking manually in a second, orthogonal plane [16]. Since there was only one whisker, coordinates describing the whisker in one plane could be uniquely matched to those in the other plane. However, manual tracking is impractical for large data sets [6,7,9,11,14]. In sum, no method currently exists that achieves automatic and non-invasive tracking of multiple whiskers in 3D. Our innovation here was to achieve this by applying constraints expressing ‘prior knowledge’ of natural whisking to the 3D reconstruction problem.

Our whisker tracker describes each target whisker as a 3D Bezier curve (Fig 2). This is a parametric curve segment b(s) = (x(s),y(s),z(s)), where 0≤s≤1 parameterises location along the curve. In our case, s = 0 marked the end closest to the whisker base and s = 1 marked the end furthest from the base. A Bezier curve is defined by 2 or more ‘control points’ (Fig 2), the number of which controls the complexity of the curve. We used quadratic Bezier curves, each of which has 3 control points since, as noted above, this is the lowest-degree curve that is accurate for our purposes.

Fig 2

Description of whiskers by quadratic 3D Bezier curves.

Left: schematic of a 3D Bezier curve representing a whisker (blue line), defined by its three control points cp,cp and cp (blue dots). Middle, right: projection of the 3D Bezier curve, and its control points, onto horizontal and vertical image planes.

Whisker tracking algorithm

The essence of our algorithm was to track one or more target whiskers by fitting 3D Bezier curves to the image data. The core principle was, for each frame, to tune the control points of the Bezier curves so that their projections onto the two image planes matched as closely as possible the images of the basal segments of the target whiskers. The degree of match was quantified by the following cost function (Methods, Eq 4):

Here E and E (defined in Eqs 5 and 6) measured the average image intensity (the line integral) along a given Bezier curve projected into the horizontal and vertical image plane respectively; R1 and R2 (defined in Eqs 7 and 8) were regularising terms. Since whiskers imaged as black (low pixel intensity) and background as white (high pixel intensity), E and E were low when a Bezier curve b(s) coincided with a whisker (Fig 2), and at a local minimum with respect to local changes to the positions of the control points. In contrast, E and E were high if they coincided with a background region. Thus, a target whisker was tracked by minimising the cost function with respect to the control point positions of the associated Bezier curve.

It was possible to track multiple whiskers in parallel by taking advantage of prior knowledge of natural whisker motion. First, since whiskers move and deform smoothly over time, the location of a whisker in a given frame was (with sparse exceptions, see below) predictable from that in the previous frames. Thus, by seeding the control point positions of a Bezier curve representing a given whisker using corresponding positions from previous frames (Methods), it was possible to maintain accurate ‘locking’ between each Bezier curve and its target whisker. Second, we used prior knowledge to constrain the cost function. As detailed in Methods, a ‘temporal contiguity’ constraint (R1) penalised discontinuous, temporal changes in Bezier curves (Methods, Eq 7) and a ‘shape complexity’ constraint (R2) penalised unnaturally complex curve shapes (Methods, Eq 8). The full whisker tracking pipeline (Fig 3) is detailed in Methods.

Fig 3

Whisker tracking pipeline.

Left: Initialisation of control points for a given target whisker (see Methods for details). Initial values for control points in horizontal (top, white circles) and vertical views (bottom, white circles). White dotted lines in vertical view represent the range of z values consistent with each of the (x,y) points in horizontal view. Middle: Estimation of snout contour (yellow). Right: Fitting of 3D Bezier curves to image data. Projections of the 3D Bezier curve for one whisker (blue lines) and of its control points (blue dots) are shown in horizontal (top) and vertical (bottom) views. Yellow dots indicate intersections between snout contour and extrapolated Bezier curves.

Tracking multiple whiskers in 3D

To test the algorithm, we applied it to the image data from our task. We found that we were able to track several whiskers at the same time (Fig 4; S1 Movie). Fig 4C shows a 12 ms sequence of whisker-pole contact from a mouse where 8 whiskers were intact and the others had been trimmed to the level of the fur (Fig 4A and 4B). The algorithm successfully tracked changes in both orientation and shape of the 8 whiskers. Different types of motion were tracked: some whiskers bent against the pole whilst others slipped past it (Fig 4D). The outcome of the tracker was a sequence of 3D curve segments, each representing the basal segment of a given whisker in a given frame (Fig 4E).

Fig 4

Tracking multiple whiskers in 3D.

A-B) 8 whiskers were tracked in a 3.5 s video sequence (1000 frames/s). C) A sequence of 12 frames showing Bezier curves for all tracked whiskers, projected into horizontal and vertical views, taken from the example video (S1 Movie). Whiskers are colour coded as in panel A. D) Tracking solutions for 2 whiskers (colour coded as in panel A) across 12 frames projected onto horizontal and vertical views. E) 3D tracking solutions for 8 whiskers across a sequence of 30 frames, including the sequence of panel D.

To assess tracking accuracy, we tracked a randomly selected set of 100 trials (50 go, 50 no-go) where a mouse was performing the task with 3 whiskers, all others trimmed to the level of the fur. This dataset comprised 350,000 frames. During ‘free whisking’, changes in whisker position/shape were entirely due to whisking motion, and such changes were smooth as a function of time, so that the ‘temporal contiguity’ and ‘shape complexity’ constraints of the cost function (Methods, Eqs 6 and 7) were accurate. Such errors as did occur were mainly due to (1) occlusion and (2) whisker overlap. (1) On occasion, a whisker was occluded against either the mouse’s body (ear or cheek) or the experimental apparatus (pole). However, by optimising the view angles of the cameras (see above) and by minimising the image footprint of the apparatus, we minimised these effects. On the no-go trials, we detected 0.11 occlusion events/whisker per 1000 frames. Since occlusion was rare, such events were dealt with by skipping affected frames and restarting tracking afterwards. (2) On occasion, whiskers overlapped each other in either horizontal or vertical view. Because the tracker applies prior knowledge of natural whisker shape/location, our algorithm was relatively robust to such events. For example, the tracked video sequence illustrated in Fig 4 includes overlap between whiskers C1 and D2. However, errors did sometimes occur. The most common overlap error was when a small whisker overlapped a large one to such a degree that the Bezier curve tracking the small whisker locked onto the large whisker. Such events were minimised by trimming all non-target whiskers to the level of the fur. Also, since our software includes capability for manual curation, we were able to rectify these errors when they did occur, by using GUI tools to nudge the control points of the errant Bezier curve back onto its target whisker. The incidence of overlap errors increased with the number of intact whiskers. It depended also on their location: there was typically less overlap within a row of whiskers than across an arc. On no-go trials of the test data set, there were 0.01 overlap errors/whisker per 1000 frames.

Tracking was more challenging when videos included not only whisking motion but also whisker-object contact. During contact, changes in whisker position/shape from frame to frame were most often smooth and gradual but, occasionally, a whisker slipped off the pole at high speed (‘slip event’), generating discontinuous whisker change between adjacent frames (Fig 4D and 4E). During such slips, tracking errors sometimes occurred, since the tracker’s routine for estimating the location of a whisker based on previous frames assumes smooth motion. On go-trials, high-speed slips occurred in a small fraction of video frames (0.23 slips/whisker per 1000 frames). As above, using the software’s GUI curation tools, we were able to correct errors resulting from slips.

Measuring 3D whisker orientation and 3D whisker shape

Having tracked one or more whiskers in a set of videos, the next step was to use the tracking data to estimate 3D whisker kinematics and 3D whisker shape. To this end, we measured the 3D orientation of each whisker [5,16,29] in terms of its azimuth (θ), elevation (φ) and roll (ζ) (Fig 5; Methods; S2 Movie).

Fig 5

Description of a whisker in terms of 3D kinematic and 3D shape parameters.

A) Azimuth (θ), elevation (φ) and roll (ζ) angles. These angles are defined with respect to the tangent to the Bezier curve b(s) describing the whisker, at s = 0. Azimuth describes rotation about the vertical (dorso-ventral) axis through s = 0; elevation describes rotation about the horizontal (anterior-posterior) axis through s = 0; roll describes rotation about the x′ axis, defined in panel B. B) Left. Whisker-centric coordinate frame with origin at s = 0 (Eqs 9–11). The x′ axis is tangent to b(s) at s = 0; the y′ axis is the direction in which b(s) curves; the z′ axis is orthogonal to the x′−y′ plane. Middle. Components of moment in the whisker-centric coordinate frame. Right. 2D and 3D whisker curvature (Eqs 13–15). rh and r denote the radii of the circles that best fit the projection of b(s) into the horizontal and vertical image planes respectively (at a given point s); r denotes the radius of the circle that best fits b(s) itself.

To illustrate the method, we first estimated 3D whisker angles during free whisking, and illustrate results for a mouse with whiskers C1-C3 intact (Fig 6). Azimuthal whisker angle was highly correlated across whiskers (whiskers C1-C3, Pearson correlation coefficients ρ = 0.98–0.99), as was elevation (ρ = 0.94–0.99) (Fig 6A). Elevation was highly anti-correlated with azimuth (ρ = -0.89–0.96; Fig 6B left). Roll angle correlated with azimuth/elevation but, consistent with [5], the degree of correlation was whisker-dependent (ρ = 0.13–0.80; Fig 6B middle, right). Whisker-object contact perturbed these angle relationships and partially decorrelated the azimuth-elevation relationship (ρ = -0.49–0.79).

Fig 6

3D whisker kinematics during free whisking.

A) Changes in 3D angles for whiskers C1, C2 and C3 during a 3.5 s episode of free whisking. B) Relationships between angles.

An important feature of video-based whisker tracking is that it allows non-invasive measurement not only of kinematics but also of mechanical forces/moments acting on the whiskers [4,16,24,26]. The physical basis for this is that the shape of a whisker contains information about these mechanical factors. In particular, when a whisker quasi-statically bends against an object, there is a linear relationship between whisker curvature and bending moment. On this basis, studies combining video-based tracking with measurement of neural activity have discovered that bending moment is the best single predictor of primary afferent activity during whisker-object contact [4,11,13,14]. However, it should be noted that, in other situations (whisking in the absence of object contact and dynamic situations such as texture exploration and collisions), curvature does not have the same clean, mechanical interpretation [30-32].

Previous studies have sought to measure the bending moment related to whisker-object contact by imaging in the horizontal plane [6,7,9-11,14,24,25]. However, there are limitations of the planar imaging approach. First, it senses only the component of bending moment in the horizontal image plane and necessarily misses any out-of-plane bending. Second, since whiskers roll during the whisking cycle, the shape of a whisker, as projected in the horizontal plane can change purely due to roll even during free whisking (Fig 7). A benefit of 3D imaging is that it addresses these limitations. First, 3D imaging enables bending in any direction to be measured. Second, it permits the intrinsic shape of a whisker to be teased apart from both its position and angular orientation. We measured the intrinsic shape, which we term κ3 : κ3(s) expresses the 3D curvature at each point s along the whisker shaft and has the key property of being invariant to curve location and to curve orientation (Methods; Eq 13).

Fig 7

Tracking and estimating 3D curvature for a rigid test object (panels A-B), whiskers of a behaving mouse (panel C) and an ex vivo whisker (panels D-F). A) Tracking the edge of a coverslip. The coverslip was mounted, like a lollipop, on a rod; the rod was oriented in the mediolateral direction and rotated around its axis. Red lines indicate tracking results (30 frames, 10 millisecond intervals, 1000 frames/s). B) Top: Azimuth angle for two trials (black and grey traces). Bottom shows measured curvature: horizontal curvatures (dotted lines), κ3 (solid lines) and true curvature (orange). C) Horizontal and 3D curvatures during free whisking (same trial as Fig 6). Solid lines represent κ3 and dotted lines indicate horizontal curvatures for C1-3 (colours coded as in Fig 6). Fluctuations in vertical curvature were similar to those in horizontal curvature (|ρ|>0.49). D) Variation in κ3 for a stationary ex vivo whisker (C3) as a function of roll angle. E) Azimuth angle for ex vivo trials with simulated whisking at different speeds. F) κ3 as a function of whisking phase.

First, to verify that our system accurately measures 3D curvature we tested whether it could correctly recover the shape of a rotating, rigid object. To this end, we simulated whisking (motion with both forward-backward and rolling components) of a rigid disk, whose edge had curvature comparable to that of a typical mouse whisker (a 13 mm diameter, glass coverslip). We imaged it as described above and tracked its edge (Fig 7A). Due to the rolling motion, curvature in the horizontal and vertical planes, κ and κ, oscillated strongly (Fig 7B). Despite this, our algorithm correctly recovered that 3D curvature κ3 was constant—fluctuations (SD) in κ3 were 5.6% that of κ and 8% that of κ—and matched the true curvature of the object (Fig 7B).

However, when we measured κ3(s) of whiskers from behaving mice, we found substantial fluctuations, even during free whisking (Fig 7C): SD of κ3 was 40–91% that of κ and 51–68% that of κ. This suggests that, in contrast to rat where the proximal segment of the whisker shaft undergoes rigid motion during whisking [5], during natural whisking, mouse whiskers do not behave as rigid objects. To determine whether this lack of rigidity occurs in the absence of whisker movement (‘static’) or whether it is dependent on movement (‘dynamic’), we first positioned a (stationary) ex vivo mouse whisker (C3) in the apparatus at a series of roll angles, imaged it, tracked it and measured κ3. Varying roll angle changed κ3 by up to 6% (Fig 7D). Next, we ‘whisked’ the whisker backwards and forwards in the horizontal plane at different velocities (Fig 7E) and measured κ3 as a function of whisking phase (Fig 7F). We found that κ3 changed by up to 12% over a cycle. This change in κ3 increased linearly with whisking velocity (Fig 7F; R2 = 0.95). Linear extrapolation to average whisking velocity observed during free whisking in mice performing our task (0.78 ⁰/ms) implied a κ3 change of 66% (C3). In comparison, changes in κ3 during free whisking were 56–150% (whiskers C1-3). Overall, these data indicate that whisker motion during free whisking in mouse is non-rigid and predominantly a dynamic effect.

Estimation of κ3 allowed construction of a simple proxy to the magnitude of the bending moment, which we term Δκ3 (Methods). This quantity is a generalisation of the corresponding measurement from horizontal plane imaging, referred to here as Δκ. Since the mice were whisking against a vertical pole, the predominant changes in curvature were in the horizontal plane and, thus, one might expect minimal benefit from the 3D approach. Even here, however, we found Δκ to be markedly contaminated by roll. A simple instance of this effect is shown in Fig 8A (S3 Movie). Here, the whisker initially curved downwards (roll angle -90⁰): whisker-pole contact from time 0–45 ms rolled the whisker in the caudal direction (roll angle -180⁰) with only minimal change in 3D curvature but with substantial effect on the curvature projected in the horizontal plane. Thus, Δκ increased by 0.05 mm-1, introducing a marked mismatch between Δκ and Δκ3 (Fig 8A bottom). A more typical and complex instance of the effect of roll angle is shown in Fig 8B (S4 Movie). Here there were large fluctuations in Δκ (e.g., at times 310–370 ms) which almost entirely reflected changes in roll angle in the absence of change to 3D curvature. On average, Δκ explained 44% of the variation in Δκ3 (touch periods from 47 trials). In this way, 3D imaging permits more accurate measurement of mechanical forces acting on whiskers.

Fig 8

Comparison of 2D and 3D curvature as mouse whisks against a pole (whisker C2): curvatures (upper panel)-, 3D kinematics (middle panel) and curvature change (bottom panel, Δ, Δ, and Δ). A) Contact episode where both movement and bending of the whisker were largely restricted to the horizontal plane. In this case, Δκ3 and Δκ were highly correlated. Grey shading indicates periods of whisker-pole contact. See S3 Movie. B) Example with same whisker as panel A for contact episode with significant vertical component of whisker motion. See S4 Movie.

Discussion

In order to obtain a comprehensive description of 3D whisker movements and 3D whisker-object interactions, we imaged mouse whisking behaviour using a high-speed ‘stereo’ imaging system. We developed software, first, to fit a 3D curve segment to each of one or more target whiskers and, second, to extract 3D kinematic and shape parameters from them. The new method allows both the 3D orientation (azimuth, elevation and roll) of a whisker and its intrinsic 3D curvature to be measured at millisecond frame rate, during both free whisking and whisker-object contact.

The vast majority of previous work on automatic whisker tracking has focussed on imaging in the horizontal plane [6,7,9,14,22-24,27,30]. However, as detailed in the Introduction, single-plane imaging necessarily captures only a fraction of the full 3D kinematic and 3D shape parameters that characterise a whisker. The advance here is a system able to extract a full 3D description of both whisker kinematics and whisker shape. The tracking algorithm is automatic (after initialisation) and achieves accurate performance under our typical conditions. The system also includes GUI error correction tools, which allow accuracy to be further improved, and extend the applicability of the system to more challenging situations.

Our work builds on previous advances in 3D whisker tracking in behaving rodents. Methods based on linear CCD arrays and markers provide simple and effective means to measure some aspects of 3D whisker kinematics, but also have limitations. CCD arrays [15] image whiskers at a single point along the shaft: this method confounds whisker translation with whisker rotation, cannot measure whisker shape and cannot recover roll angle. Also, CCD study of 3D whisking has been limited to one whisker at a time. The approach of marking a whisker with spots of dye [5,28] has been used to recover all 3 orientation angles but has not been used to measure whisker shape changes during object touch. Also, although application of dye marker dots was reported not to disturb rat whisking, mouse whiskers are thinner, less stiff and hence more liable to perturbation. Our method addresses these limitations by leveraging the extra information available from video, and achieves automatic, non-invasive multi-whisker tracking, as well as the ability to recover not only three-angle kinematic information but also information about whisker shape changes during object contact. Previous 3D tracking [16] required manual tracking for the vertical view. Our algorithm runs on a standard workstation and the system only requires addition of a second camera to an existing 2D imaging set-up.

The method presented here has some limitations. First, although the error rate was low, the experimental conditions were designed to avoid both extensive object-whisker occlusion (by using a stimulus object with small foot-print) and extensive whisker-whisker overlap (by trimming non-target whiskers). Error rates are likely to be higher in the presence of experimental apparatus where there is more occlusion or if no whisker trimming is carried out. Use of additional cameras provides a potential way to reduce error rates further (and can be incorporated through adding terms analogous to Eqs 5 and 6 to the cost function of Eq 4). Second, our tracker describes the basal segment of the whisker, not the full whisker, since our focus is primarily on elucidating the fundamental mechanical events that drive neural activity in the whisker system. To track a whisker across its entire length, a quadratic curve description is insufficient. One direction for future work is to investigate whether the tracker can be extended to fit higher degree Bezier curves without excessive loss of robustness. Third, our method has been developed for head-fixed mice. This has the advantage that it simplifies the tracking problem. However, it would also be useful to extend the approach to freely moving animals, perhaps by combining whisker tracking with head-body tracking [33-41].

In conclusion, although whisker movements are well-known to be 3D, previous, automatic methods for tracking whiskers from high-speed video were limited to the horizontal plane. Here, we obtained 3D information using a two-camera, high-speed imaging system and developed computational methods to reconstruct 3D whisker state from the video data. Our method permits measurement of both 3D whisker kinematics and whisker shape changes at millisecond frame rate from awake, behaving mice. The method can be combined with cellular resolution neural activity measurement and thus has potential to advances our understanding of sensorimotor behaviour and its neural basis.

Methods

All parameters and variables used in this section are summarised in Table 1. All computer code was written in MATLAB (The MathWorks Inc., Natick, MA) and run on a standard workstation (Core i7, 16GB RAM).

Table 1

Parameters and variables summary.

Axes
x	Anterior-posterior axis, with positive posterior
y	Medio-lateral axis, with positive medial
z	Dorsal-ventral axis, with positive dorsal
v	Medio-lateral axis from the vertical view, with positive medial
w	Dorsal-ventral axis from the vertical view, with positive ventral
x′	Axis of the whisker centred coordinated frame described as tangent to the whisker at s = 0. Positive indicates toward the tip of the whisker.
y′	Axis of the whisker centred coordinated frame defined by the second derivative of the Bezier curve at s = 0. Positive indicates s>0
z′	Axis of the whisker centred coordinated frame defined as the cross product between x′ and y′
Calibration
pH	Point in the horizontal plane with coefficients (x,y)T
H	Matrix related to mapping between p3D and pH. [1 0 0; 0 1 0]
p3D	Point with coefficients (x,y,z)T
h	Vector related to the mapping between p3D and pH.
pV	Point in the horizontal plane with coefficients (v,w)T
V	Matrix related to mapping between p3D and pV. Values of V are fitted during the calibration procedure
v	Vector related to the mapping between p3D and pV.Values of v are fitted during the calibration procedure
Bezier curves and fitting process
b(s)	Bezier curve evaluated at s where 0≤s≤1
cpi	Control point i = 0,1,2 with coordinates (x,y,z)T
j	whisker
f	frame
Ev(j)	Term from objective function (Eq 4) related to vertical image
Eh(f)	Term from objective function related to horizontal image
Ih(x,y)	Intensity at the point (x,y) in the horizontal image
Iv(v,w)	Intensity at the point (v,w) in the horizontal image
R1	Regularising term from the objective function related to temporal contiguity
R2	Regularising term from the objective function related to shape complexity
σ1	Selectable parameter that weights the first regularising factor R1
σ2	Selectable parameter that weights the second regularising factor R2
q	Vector related to Eq 8.
Extraction of kinematic parameters
i′	Unit vector that point in the direction of x′
j′	Unit vector that point in the direction of y′
k′	Unit vector that point in the direction of z′
ζ	Rotation angle of the whisker respect to x′
θ	Azimuth angle defined as the angle between the x axis and the projection of tangent at the base of the whisker in the horizontal plane
φ	Elevation angle defined as the angle between the -z axis and the projection of tangent at the base of the whisker in the vertical image
κ3D(s)	3D Curvature evaluated at s
κh(s)	Curvature of the projection of the Bezier curve in the horizontal plane evaluated at s
κv(s)	Curvature of the projection of the Bezier curve in the vertical image evaluated at s

Ethics statement

This research has been approved by the University of Manchester Animal Welfare and Ethical Review Board and by the United Kingdom Home Office.

Behavioural apparatus

Mice (C57; males; 6 weeks at time of implant) were implanted with a titanium head-bar as detailed in [11]. After surgery, mice were left to recover for at least 5 days before starting water restriction (1.5 ml water/day). Training began 7–10 days after the start of water restriction.

Mice were trained and imaged in a dark, sound-proofed enclosure using apparatus adapted from [11]. A head-fixed mouse was placed inside a perspex tube, from which its head emerged at one end. The stimulus object was a 0.2 mm diameter, vertical carbon fibre pole which could be translated parallel to the anterior-posterior (AP) or medio-lateral (ML) axes of the mouse by a pair of linear stepper motors and rotated in the horizontal plane to ‘go’ or ‘no-go’ locations by a rotatory stepper motor. To allow vertical movement of the pole into and out of range of the whiskers, the apparatus was mounted on a pneumatic linear slide, powered by compressed air. The apparatus was controlled from MATLAB via a real-time processor. Mouse response was monitored by a lick port located anterior to the mouth. Licks were detected as described in [6]. Each lick port consisted of a metal tube connected to a water reservoir via a computer-controlled solenoid valve. Lick port position was monitored using an infrared camera and adjusted using a micromanipulator.

Behavioural task

Head-fixed mice were trained to detect the presence of a metal pole using their whiskers, using behavioural procedures similar to [9]. On each trial, the pole was presented either within reach of the whiskers (‘go trial’) or out of reach (‘no-go trial’). At the start of each trial, the computer triggered the pole to move up (travel time ~100 ms). The pole stayed up for 1 s, before moving down. On go trials, the correct response was for the mouse to lick a lick port. Correct responses were rewarded by a drop of water (~10 μl). Incorrect responses on go trials (not licking) were punished by timeout (3–5 s). On no-go trials, the correct response was to refrain from licking and incorrect responses (licking) were punished by timeout and tone (frequency 12 kHz).

High-speed stereo whisker imaging

Whiskers were imaged, based on the methods of [11], except that, to provide 3D information, two cameras were used. The whiskers were imaged using two high-speed cameras (Mikrotron LTR2, Unterschleissheim, Germany; 1000 frames/s, 0.4 ms exposure time) via telecentric lenses (Edmunds Optics 55–349, Barrington, NJ) as illustrated in Fig 1. Illumination for each camera was provided by a high-power infrared LED array (940 nm; Roithner LED 940-66-60, Vienna, Austria) via diffuser and condensing lens. The imaging planes of the two cameras were horizontal (spanning AP and ML axes) and vertical respectively. The field of views were typically 480 x 480 pixels, with pixel width 0.047 mm.

The two cameras were synchronised by triggering data acquisition off the computer-generated TTL pulse that initiated a trial. Typically, imaging data were acquired in an interval starting 0.5 s before pole onset, ending 1.8 s after pole onset. To provide an independent check that data files from the two cameras came from corresponding trials, an IR LED was positioned in the corner of the field of view of each camera and, starting at pole onset on each trial, flashed a binary sequence that encoded the trial number. Onset of this LED signal also served to verify camera synchrony.

Coordinate frame and calibration

To describe the location of whiskers in 3D, we used a left-handed Cartesian coordinate frame, fixed with respect to the head of the animal (Fig 1). The axes were x (AP, with positive x posterior); y (ML, with positive y medial) and z (DV, with positive z dorsal). In standard anatomical convention, the x−y (AP-ML) plane was horizontal; the x−z (AP-DV) plane sagittal and the y−z (ML-DV) plane coronal. The column vector p3D = (x,y,z) denotes a point with coefficients along the x, y and z axes respectively. Throughout, we denote vectors by lower-case bold (e.g., p), scalars by lower-case italic (e.g., s) and matrices by upper-case bold (e.g., M).

Analogously to stereoscopic vision, our whisker tracker reconstructs 3D whisker location/orientation from images obtained from two viewpoints—horizontal and vertical. Pixel locations in the horizontal image were defined using the x and y axes of the 3D frame (Fig 1). Thus, the location of a point pH in the horizontal image was described by a column vector (x,y). Pixel location pV in the vertical image was described by a column vector (v,w), defined with respect to axes v and w (Fig 1). Due to the orientation of the vertical-view camera detailed above, the v,w coordinate frame was rotated and translated with respect to the x,y,z frame. The relation between the x,y,z and v,w coordinate frames was determined as follows.

With a telecentric lens, only light rays parallel to the optical axis pass through the camera aperture and contribute to the image. In contrast to standard lenses, a telecentric one provides equal magnification at a wide range of object-lens distances and forms an orthogonal projection, facilitating the reconstruction of 3D objects [42]. Since the imaging was done with telecentric lenses, the mappings from 3D to the two image planes were described as orthogonal projections. The projections of a 3D point p3D onto a 2D point pH in the horizontal image plane and a 2D point pV in the vertical plane were:

H and V were 2x3 matrices; h and v, 2-element column vectors. The 3D coordinate frame and the horizontal image frame had common x and y axes and a common x,y origin: hence, H = [1 0 0; 0 1 0] and h = [0;0].

To determine the mapping from 3D to the vertical image plane (V and v), we performed the following calibration procedure. Using stepper motors, we moved an object with 2 protruding pins on a 3D path through the region of the behavioural set-up where the target whiskers were typically located, and recorded a sequence of 100 corresponding images on each camera. The z location of each pin was known in each image. We then tracked the tips of the 2 pins in each horizontal and vertical image frame to obtain a time series consisting of the (x,y,z) and (v,w) coordinates of each pin tip. Using Eq 2, V and v were then estimated by linear regression. The variance of the regression residuals was low (<0.1% of total variance for data of Fig 6).

Bezier curve framework for whisker tracking

Our aim was to develop ‘whisker tracker’ software to track the orientation and shape of one or more target whiskers. The whisker tracker described each target whisker segment as a Bezier curve, since these have convenient mathematical properties (Fig 2). A Bezier curve is a parametric curve segment b(s) = (x(s),y(s),z(s)), where 0≤s≤1 parameterises location along the curve segment: in our case, s = 0 marked the end closest to the whisker base and s = 1 marked the end furthest from the base. The shape, orientation and position of a Bezier curve are determined by its ‘control points’, the number of which determines the complexity of the curve. We used quadratic Bezier curves, which have 3 control points cp where i = 0,1,2, each with coordinates (x,y,z) These control points were termed “proximal” (cp0), “middle” (cp1) and “distal” (cp2) according to their distance from the whisker base. cp0 defined the location of the basal end of the whisker segment, cp2 the distal end and cp1 the shape. In terms of these control point parameters, a quadratic Bezier curve b(s) was expressed as:

Whisker tracking pipeline

Overview

Whisker tracking was operated via a Graphical User Interface (GUI). The GUI allowed a user to load a pair of corresponding videos (horizontal view and vertical view). The first step was to calibrate, as detailed above. Next, to initialise tracking, target whiskers were specified (automatically or manually) by defining approximate locations for Bezier control points. After initialisation, each video frame was processed automatically, in turn. Tracking direction could be set to be either forwards or backwards in time. First, the contour of the snout was located in both horizontal and vertical views. Second (except in the first frame), initial estimates for the Bezier control points were calculated by linear extrapolation from their locations in the previous frames. Third, each Bezier curve was fitted to the image data by adjusting its control points to minimise the cost function defined below (Eq 4). Provided the quality of fit for a given Bezier curve met a minimum threshold, tracking of that curve proceeded automatically to the next frame.

Manual initialisation (Fig 3 left panel)

When tracking a video for the first time, the first step was to specify the target whiskers. In the first frame of the video, the user employed a graphical user interface (GUI) to select sets of control points, specifying one or more target whiskers. For each target whisker, the user defined approximate locations for control points specifying a curve segment corresponding to the basal segment of the whisker, by making computer-mouse clicks within the video images. For each target whisker, the user first specified (x,y) coordinates for the 3 control points in the horizontal view. Since the imaging geometry was described by linear equations (Eqs 1 and 2), each such point corresponded to a line in the vertical view. In the vertical view, the user specified a point (v,w) along each vertical line where it intersected the target whisker. From these (x,y,v,w) data, the z coordinates of the control point estimates were calculated from the calibration equation (Eq 2). Once initial values for all 3 control points of a given target whisker were specified, refined estimates were calculated by the fitting procedure described below. In order to obtain a reference value for the length of each curve, the arc length of each Bezier curve was calculated. To obtain a reference value for the distance of the proximal control point from the snout, a second order polynomial was fitted to the Bezier curve and extrapolated to find its intersection with the snout contour.

Automatic initialisation

Typically, an experiment will result in many videos taken of the same mouse under identical experimental conditions. Once one video was tracked using the manual initialisation procedure described above, other videos could then be initialised automatically through a template-matching approach. To initialise tracking of a new video, the user selected a previously tracked video (via the GUI). For each target whisker from this file, a sample of Bezier curve ‘templates’ were extracted (typically, the solution in every fifth video frame) and goodness of fit of each sample curve to the first frame of the new video was calculated (using the cost function, Eq 4). The lowest cost template was then selected. The template was refined by optimising the fit with respect to translations along both x and z axes (within the range ±5 pixels).

Snout contour detection (Fig 3 middle panel)

First, to isolate the contour of the snout in a given video frame, fine structure such as the hairs of the fur and the whiskers were removed by median filtering (5 x 5 pixels, 0.23 x 0.23 mm) of the images followed by smoothing with a Gaussian filter (SD = 12 pixels, 0.56 mm). Next, the spatial gradient of each filtered image was calculated in a direction approximately normal to the snout contour. This gradient was small, except at the edge of the snout where it had a large peak. In the horizontal image, the snout contour was estimated as a function of the x coordinate by minimising the gradient with respect to y. In the vertical image, the snout contour was estimated as a function of the w coordinate by minimising the gradient with respect to v.

Bezier curve fitting (Fig 3 right panel)

To achieve 3D tracking, we fitted 3D Bezier curves to the horizontal and vertical view image data by varying the locations of their control points so as to minimise the following cost function. Control points for each target whisker were optimised independently:

Here E(f) is the cost (or mismatch) between the image data of frame f and the Bezier curve b(f,s) = (x(f,s),y(f,s),z(f,s)), defined by control points cp0(f), cp1(f) and cp2(f). E(f) and E(f) quantified how well b(f,s) described, respectively, the horizontal and vertical image data of frame f. R1(f) and R2(f) were regularising terms (defined below). In the following, to keep down clutter in the notation, dependence on frame and whisker is omitted except where necessary for clarity.

E and E were defined as line integrals over the projection of b(s) in the horizontal/vertical images respectively:

Here: I(x,y) was the intensity at point (x,y) in the horizontal view image, calculated by linear interpolation between pixel values and I(v,w) the analogous quantity for the vertical view image; (x(s),y(s)) was the projection of b(s) in the horizontal view image and (v(s),w(s)) its projection in the vertical view image (Eqs 1 and 2).

Except at occasional stick-slip events, whiskers move smoothly and, when imaged at 1000 frames/s, changes in location and shape from frame to frame were usually small, particularly for the basal segment. The regularising term R1 formalised this prior knowledge of natural whisking behaviour (‘temporal contiguity’):

Here: cpi(f) was the location of control point i of the whisker in frame f and was its location estimated by linear extrapolation based on its location in the previous two frames; σ1 was a variable gain that the user could set from the GUI.

Additional regularisation was necessary to address degeneracy that could arise when tracking near-straight whiskers. Since a line segment is fully described by the location of its two ends, a straight whisker is fully defined by its proximal and distal control points–in this case, the middle control point is ill-defined. We found, in such situations, that the middle control point tended to migrate towards the whisker base and to generate high-curvature, unnatural shapes when extrapolating the curve to the snout contour (see above). To address this, we used a second regularising term which penalised deviations of the middle control point away from the midpoint between the proximal and distal control points:

Here q(f) = cp2(f)−cp0(f) and σ was a user-adjustable gain. R2 measured deviation of the component of the middle control point cp1(f) (relative to cp0(f)) in the direction of cp2(f) (relative to cp0(f)) away from the midpoint of the cp0(f)−cp2(f) line.

Nonlinear cost functions can be difficult to minimise due to local minima. However, in the present case, due to the smooth motion of whiskers referred to above, we expected control point solutions usually to be close to their values in the previous frame. Not only was it therefore effective to use a local search strategy, where the initial value for a given control point was set by extrapolating its values from the previous two frames (), but this also made it possible to track multiple whiskers independently. The cost function (Eq 4) was minimised (using MATLAB function ‘fminunc’) with respect to components of the control points. To counter-act possible drift of b(s) along the whisker shaft over time, or change in the arc-length of b(s) over time, we minimised the cost function with respect to components of cp0 and cp2 normal to b(s) at s = 0 and s = 1 respectively. This procedure also had the advantage of reducing the number of free parameters from 9 to 7. Furthermore, after convergence in a given frame, both the arc length of b(s) and the distance of cp0 to the snout were normalised to equal their reference values set in the first frame (see above), whilst preserving curve shape.

Error correction

As noted above, tracking of each target whisker proceeded automatically to the next frame, so long as the cost E (Eq 4) remained less than a user-defined threshold (adjustable via the GUI). Should the threshold be exceeded, for example when a Bezier curve was ‘left behind’ by rapid, discontinuous motion of its target whisker during a slip event, tracking of that whisker ceased. To correct such an error, the GUI had tools allowing the user to nudge control points back onto the target whisker, and to restart automatic tracking.

Extracting 3D kinematics of the tracked whiskers

The next step was to use the tracking data to estimate 3D whisker kinematics and 3D whisker shape. Since whiskers bend during whisker-object contact, and since this contact-induced whisker bending is a fundamental driver of neural activity (see Introduction), it was important to develop a general procedure for describing 3D whisker motion, applicable to non-rigid whisker movement We separated changes to the orientation of a quadratic curve from changes to its shape in the following manner.

Formally, we described whisker orientation by the following ‘whisker-centred’ Cartesian coordinate frame x′y′z′, with origin at s = 0 [16]. In contrast to the head-centred coordinate frame xyz, the x′y′z′ frame is time-dependent; rotating and translating along with its target whisker. The x′-axis is aligned to the longitudinal axis of the whisker (tangent to b(s) at s = 0). The y′-axis is orthogonal to the x′-axis, such that the x′−y′ plane is that within which b(s) curves. The z′-axis is orthogonal to both x′ and y′ axes. Let i′,j′ and k′ be unit vectors that point in the direction of the x′, y′ and z′ axes respectively:

Here i′×j′ denotes the cross product of vectors i′ and j′. The orientation of a whisker was then described by the 3D angle of the x′y′z′ coordinate frame with respect to the xyz coordinate frame. We translated the frames to have a common origin and then calculated the 3D rotation matrix that rotates the xyz frame to the x′y′z′ frame [43]. This rotation can be described as the net effect of an ordered sequence of three elemental rotations with angles θ (azimuth), φ (elevation) and ζ (roll), and was expressed as a matrix (θ,φ,ζ). Azimuth describes rotation in the horizontal (x−y) plane, about an axis parallel to the z axis through the whisker base; elevation describes rotation in the vertical (x−z) plane, about an axis parallel to the y- axis; roll describes rotation around the axis of the whisker shaft (Fig 5AB). We determined the angles θ,φ,ζ for a given whisker at a given time point by minimising the error function:

Here i, j and are column unit vectors parallel to the x,y and z axes and the summation is over all matrix elements.

Extracting 3D shape and bending moment of the tracked whiskers

Having described the orientation of a whisker, the next task was to describe its shape. By ‘shape’, we intend those geometric properties of a curve that are invariant to its location and orientation. As noted above, we described whiskers by quadratic curve segments, which curve entirely within a plane (geometric torsion τ(s) = 0). The intrinsic shape of a quadratic curve is fully described by a curvature function κ3(s) [44]:

Here |a| denotes the magnitude (2-norm) of vector a. κ3(s) has units of 1/distance and is the reciprocal of the radius of the circle that best fits the curve at point s. We computed planar curvatures as:

Here x(s),y(s),z(s) are the components of b(s) in the x,y,z coordinate frame. Note, as detailed in Results, that, in contrast to κ3(s), κ(s) and κ(s) are not invariant measures of geometric shape; they depend also on curve orientation.

In whisker-centric coordinates, bending corresponds to changes in shape of b(s) in the x′−y′ or x′−z′ planes (with, respectively, component m defined in the direction of the positive z′ axis and m defined in the directions of the positive y′ axis) (Fig 5;[45]). Since b(s) is a quadratic curve, it has zero torsion and its curvature is entirely confined to the x′−y′ plane: κ3(s) is the curvature in this plane; the only non-zero component of bending moment is m. Applying the standard relation between bending moment about a given axis and curvature in the plane normal to that axis [26,29], it follows that m(s) is proportional to: where κ3(s) is the curvature when the whisker is free from contact and in its resting state. All results presented here were evaluated at s = 0.

Tracking example of 8 whiskers.

Left: horizontal view of the whiskers and tracking superimposed. Colours are shown as in Fig 4A. Axes are shown as in Fig 1B. Middle: Vertical view of the whiskers and tracking superimposed. Right: Bezier curves in the three dimensional space. Axes are shown as in Fig 1A.

(MP4)

Click here for additional data file.

Whisker tracking and variables.

Top: Horizontal and vertical views with tracking superimposed. Colours are shown as in Fig 4A. Bottom: 3D kinematic and 3D shape parameters: Horizontal angle (Azimuth), vertical angle (Elevation), horizontal and vertical curvature, κ3 and roll for each tracked whisker.

(MP4)

Click here for additional data file.

Whisker tracking example of movement restricted to the horizontal plane (Fig 8A).

Top from left to right: Horizontal and vertical views with tracking of C2 superimposed. Bezier curve in the whisker centred coordinate frame isolating roll angle (Fig 5A). Comparison of whisker shape over time: Dashed line represents the whisker shape at t = 0 ms and solid line represents whisker shape of current frame. Bezier curve was rotated using azimuth, elevation and roll angle to be captured in the two dimensional plane. Bottom from left to right: Horizontal angle (Azimuth), Vertical angle (Elevation), Roll, horizontal and 3D curvature over time. Colours are shown as in Fig 8A.

(MP4)

Click here for additional data file.

Whisker tracking example of movement with significant vertical components (Fig 8B).

Top from left to right: Horizontal and vertical views with tracking of C2 superimposed. Bezier curve in the whisker centred coordinate frame isolating roll angle (Fig 5A). Comparison of whisker shape over time: Dashed line represents the whisker shape at t = 0 ms and solid line represents whisker shape of current frame. Bezier curve was rotated using azimuth, elevation and roll angle to be captured in the two dimensional plane. Bottom from left to right: Horizontal angle (Azimuth), Vertical angle (Elevation), Roll, horizontal and 3D curvature over time. Colours are shown as in Fig 8A

(MP4)

Click here for additional data file.

30 Oct 2019

Dear Dr Petersen,

Thank you very much for submitting your manuscript, 'A system for tracking whisker kinematics and whisker shape in three dimensions', to PLOS Computational Biology. As with all papers submitted to the journal, yours was fully evaluated by the PLOS Computational Biology editorial team, and in this case, by independent peer reviewers. The reviewers appreciated the attention to an important topic but identified some aspects of the manuscript that should be improved.

We would therefore like to ask you to modify the manuscript according to the review recommendations before we can consider your manuscript for acceptance. Your revisions should address the specific points made by each reviewer and we encourage you to respond to particular issues Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.raised.

In addition, when you are ready to resubmit, please be prepared to provide the following:

(1) A detailed list of your responses to the review comments and the changes you have made in the manuscript. We require a file of this nature before your manuscript is passed back to the editors.

(2) A copy of your manuscript with the changes highlighted (encouraged). We encourage authors, if possible to show clearly where changes have been made to their manuscript e.g. by highlighting text.

(3) A striking still image to accompany your article (optional). If the image is judged to be suitable by the editors, it may be featured on our website and might be chosen as the issue image for that month. These square, high-quality images should be accompanied by a short caption. Please note as well that there should be no copyright restrictions on the use of the image, so that it can be published under the Open-Access license and be subject only to appropriate attribution.

Before you resubmit your manuscript, please consult our Submission Checklist to ensure your manuscript is formatted correctly for PLOS Computational Biology: http://www.ploscompbiol.org/static/checklist.action. Some key points to remember are:

- Figures uploaded separately as TIFF or EPS files (if you wish, your figures may remain in your main manuscript file in addition).

- Supporting Information uploaded as separate files, titled 'Dataset', 'Figure', 'Table', 'Text', 'Protocol', 'Audio', or 'Video'.

- Funding information in the 'Financial Disclosure' box in the online system.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com  PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

We hope to receive your revised manuscript within the next 30 days. If you anticipate any delay in its return, we ask that you let us know the expected resubmission date by email at ploscompbiol@plos.org.

If you have any questions or concerns while you make these revisions, please let us know.

Sincerely,

Robyn A. Grant

Guest Editor

PLOS Computational Biology

Samuel Gershman

Deputy Editor

PLOS Computational Biology

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Thank you for submitting your manuscript.

Your manuscript has been seen by reviewers, and I am pleased to report that, subject to satisfactory minor revision, it is considered potentially acceptable for publication.

Comments from the reviewers will help to make the paper more accessible and applicable to researchers outside of your direct field of research, and also to justify some of your statements and practices.

I look forward to receiving your next version.

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: General comments:

Seems like a very nice method, nice write up too. Would like to hear a little more about applicability, maybe, to make it as useful as possible to other experimenters.

Specific comments:

The authors use feed-forward tracking rather than full-sequence (forward-backwards) tracking. Is there any comment worth making here, trade-off of algorithmic or computational complexity vs performance?

Fig 4 - nice result. Is it cherry-picked, given the discussion in lines 170-178?

l208 "In particular, there is a linear relationship between whisker curvature and bending moment." - is this true in general? E.g. for a tapered whisker model? For large bends? Citation?

Why did the authors choose three control points? Was this because this is the minimum number required to measure a bending moment? Or some other reason? l105 "is well-approximated by" is subjective, if practical.

Fig 7 I got confused about there being only one "horizontal curvature" trace per trial in the lower panel. Worth saying "kv not shown for clarity" (that's right, is it?)?

Fig 7D is tough to read - is that a polar plot with a false radial origin? Would this be easier to read as a Cartesian plot, perhaps? The actual data, rather than the axes/annotations, is rather cramped.

Fig 7F - I can see why this is plotted polar, because it's periodic, but it's still quite tough to read, I guess because it's a bit cramped, again. Moving things around a bit so that the actual plot lines are a lot larger might help, I think. If 7D is intended to marry up with this one, and that's why it's polar, then same comment there - could be reorganised to use the available space better.

l260 "In this way, 3D imaging permits more accurate measurement of mechanical forces acting on whiskers." Certainly seems that way, good stuff.

l293 "Use of additional cameras provides a potential way to reduce error rates further." Can the authors make a statement on the applicability of the current method to 3+ cameras?

Perhaps a little more discussion on applicability? Computational load? Equipment requirements? Obstacles to running this online?

Reviewer #2: review is uploaded as an attachement

Reviewer #3: The study of Petersen et al. describes a novel combination of known computer-vision technologies, and an impressive auto-tracker, to the specific problem of tracking whisker shape in 3D. While potentially suitable for publication as a Methods paper in PLoS CB, there are a few concerns.

1. As it currently stands, the manuscript gives the impression that the 3D merging works because Bezier curves were used. In fact, as indicated in Methods, the 3D merging approach is viable because the authors performed a point-to-point calibration between the cameras using two pins mounted on a motor.

It will likely be helpful to readers to explain the basic problem more clearly: in order to 3D merge, it is necessary either to ensure that several points along the whisker length can be directly matched between the two cameras, or to perform pixel-to-pixel calibration between cameras. The first approach was used by Knutsen et al., 2008, who placed dye at known points along the whisker length. The second approach was used by Huet et al., 2015 and is used in the present work. This information would help readers better compare technical approaches between studies, and would highlight the novelty of the present work, which lies in automating the approach for multiple whiskers and in making the code publicly available.

2. Because this is a Methods paper, it would be helpful for the introduction to be more clear that the approach is currently limited to the head fixed animal, to a single row of whiskers, and for very specific locations of the cameras.

3. As currently written the manuscript risks misleading the field about the relationship between whisker curvature and the mechanical signals at the whisker base. Specifically, the manuscript states: “… there is a linear relationship between whisker curvature and bending moment,” and that several studies have computed bending moment based on curvature.

But the linear relationship between curvature and bending moment holds only in quasistatic conditions. During non-contact (“free-air”) whisking and during collisions, the quasistatic approximation does not hold. To compute mechanical signals at the whisker base under these conditions requires a dynamic model of the whisker. Notably, computing whisker dynamics requires (among many other things) knowledge of the entire whisker shape all the way to the tip, not just the proximal portion.

The manuscript should add a paragraph explicitly warning the reader that even though the curvature can be computed during non-contact whisking, during texture exploration, and immediately after collision with an object, this curvature is not a good proxy for any particular mechanical signal at the whisker base. The curvature is a good proxy for bending moment only during pure bending, as the whisker deflects against an object, after vibrations from the initial collision have damped out.

Note that curvature might still correlate with a variety of neural signals, but it is not uniquely related to the mechanical signals at the whisker base (except during pure bending).

Minor:

The sentence: “Moreover, estimates of whisker-object bending force (‘bending moment’) obtained by imaging apparent curvature of a whisker…” makes it sound as though force and bending moment are synonymous. They are very different. Also the idea of “bending force” is not clear.

The sentence: “This is significant, since bending moment is the primary driver of contact-related mechanotransduction (23–25).” should be changed to “one of the primary drivers” as the axial (longitudinal) force may also play a significant role (Stuttgen and Schwartz). In addition, the transverse forces and twisting moment could also be important – no study has yet shown otherwise.

The sentence: “First, during whisker-object contact, whiskers bend and the associated mechanical force/moment drives mechanoreceptors” should be changed to forces and moments to indicate that they are distinct quantities and both plural.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: None

Reviewer #3: No: The code is available on github, but it is not clear where to find the data underlying figures and results on figshare.

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Ben Mitchinson

Reviewer #2: No

Reviewer #3: No

Submitted filename: PCOMPBIOL-D-19-01570-review.docx

Click here for additional data file.

26 Nov 2019

Submitted filename: PCOMPBIOL_reviewers_comments_final.docx

Click here for additional data file.

5 Dec 2019

Dear Dr Petersen,

We are pleased to inform you that your manuscript 'A system for tracking whisker kinematics and whisker shape in three dimensions' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. Please be aware that it may take several days for you to receive this email; during this time no action is required by you. Once you have received these formatting requests, please note that your manuscript will not be scheduled for publication until you have made the required changes.

In the meantime, please log into Editorial Manager at https://www.editorialmanager.com/pcompbiol/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production and billing process.

One of the goals of PLOS is to make science accessible to educators and the public. PLOS staff issue occasional press releases and make early versions of PLOS Computational Biology articles available to science writers and journalists. PLOS staff also collaborate with Communication and Public Information Offices and would be happy to work with the relevant people at your institution or funding agency. If your institution or funding agency is interested in promoting your findings, please ask them to coordinate their releases with PLOS (contact ploscompbiol@plos.org).

Thank you again for supporting Open Access publishing. We look forward to publishing your paper in PLOS Computational Biology.

Sincerely,

Robyn A. Grant

Guest Editor

PLOS Computational Biology

Samuel Gershman

Deputy Editor

PLOS Computational Biology

15 Jan 2020

PCOMPBIOL-D-19-01570R1

A system for tracking whisker kinematics and whisker shape in three dimensions

Dear Dr Petersen,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Matt Lyles

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol