Morphological peculiarities of a harbour seal (Phoca vitulina) whisker revealed by normal skeletonisation

Of all mammalian vibrissae, those of certain species of pinnipeds are exceptional. Researchers believe that their curious undulating form evolved for hydrodynamic detection. Our understanding of how these whiskers work depends on a geometrical model that captures the crucial pertinent features of the natural vibrissae including its tapering and curvature. It should also account for the form of the whisker when it flexes under external loading. We introduce and study a normal skeleton of a two-dimensional projection of a harbour seal whisker. The normal skeleton is a complete shape descriptor that involves reduction to the centreline equipped with a thickness function of the orthogonal cross-section. The contours of the whisker shape are extracted from a 2D greyscale scan. Our analysis reveals correspondence between the undulations of the width and oscillations of the centreline curvature as functions of arc length. We discuss two possible explanations for that remarkable feature: one based on consideration of growth and the other of plastic deformation. For the latter we employ a mechanical model to demonstrate appearance of curvature oscillations caused by extensive deflection of the undulating whisker due to external loading.


Introduction
Among all mammals, the whiskers of certain phocids, especially of harbour seals (Phoca vitulina), figure 1, have probably the most intricate shape: close inspection reveals that they have an undulating thickness [1][2][3]. Those undulations are superimposed on a tapered and intrinsically curved whisker shaft (figure 2). Pinnipeds use their vibrissae for fine-scale tactile discrimination and hydrodynamic detection [4][5][6][7]. It was suggested that the peculiar bumpy shape of harbour seal whiskers is evolved to reduce hydrodynamic drag, to suppress self-generated noise during swimming and to improve sensitivity of detecting vortex trails left by prey [1,8,9]. This has inspired researchers to investigate high-sensitivity flow sensors shaped in a similar form, as well as aero-and hydropropulsion flow structures [10][11][12]. Our understanding of how these whiskers work depends on a geometrical model that captures the crucial pertinent features of the natural vibrissae, including its tapering and curvature. It should also account for the form of the whisker when it flexes under external loading.
An ad hoc geometric model was proposed that catches conspicuous features of harbour seal whiskers [3,13,14]. However, it does not provide a rigorous specification of how the model parameters can be measured for real vibrissae. In particular, it does not specify the location of the centreline, which is assumed to be a straight line. We need a well-defined and accurate skeletonisation method serving both the proper geometrical description and mechanical characterisation of whiskers. In this paper we set out a new approach, one that is consistent with conventional mechanical engineering models of slender elastic structures. In presenting that approach we draw attention to interesting morphological aspects of the  harbour seal vibrissae that have not been discussed to date.
The development of our model evolved from extensive analysis of the curved and tapered forms of mystacial vibrissae belonging to a wide variety of species, of which pinnipeds are a subset [15]. Inspection of 2D greyscale scans of harbour seal and grey seal (Halichoerus grypus) vibrissae revealed their characteristic slender curved profile confined by distinctly wavy boundaries. By applying a standard image processing technique we find the 2D contours of the whiskers. It is convenient to characterise the whisker shape by reducing it to a centreline (skeleton) accompanied with a function that describes the variation of width as we move along this centreline. The process of computing that pair is called skeletonisation [16]. The 2D shape of the whisker can be reconstructed on the basis of that skeleton and the width function.
From a variety of the previously proposed and used geometrical skeletonisation techniques, we recall here just three: the widely used medial axis [17] and its derivatives, the symmetry chord axis [18] and the symmetry axis defined by the process-inferring symmetry analysis or PISA [19]. In all these cases we have a set of maximal inscribed disks, whose diameters clearly tend to overestimate the perceived thickness. Note that though the symmetry chord axis may be defined without appealing to the inscribed disks and the width can be measured along the symmetry chord, the latter is not necessarily orthogonal to the centreline. As mentioned by Leyton, '. . . different symmetry analyses serve different purposes and are therefore appropriate for different circumstances' [19]. What we need for the whisker analysis is to require that the centreline bisects the normal cross-section everywhere. We call the centreline that satisfies this condition a normal skeleton. The normal skeleton is a more appropriate tool for revealing both the centreline and width of an extended shape of a whisker than other symmetry axes. Furthermore, the normal skeleton is . Normal centreline r(s) (black) for a shape bounded by two edges r ± (s ± ) (blue and red), cross sections shown in green; t and t ± are the unit tangent vectors to the centreline and the boundaries, resp., these vectors make angles θ, θ ± , resp., with the x-axis. The tangent vector r + is drawn in magenta and r − in light blue. The yellow quadrangle is spanned by r ± . consistent with one-dimensional models for slender elastic bodies, in particular for the Cosserat rod [20].
In this note we set out a procedure for constructing a normal skeleton for a 2D projection of a whisker. In section 2 we introduce the normal skeleton. In section 3 we apply this model to a natural harbour seal whisker which is representative of the form that we observe in this species. In doing so, we find new features, which are addressed in sections 4 and 5. We finish this note with concluding remarks.

The normal skeleton equations
We start with a short description of the normal skeleton for a 2D domain. In the next section this formalism will be applied to description of natural whiskers.
Let r(s) = (x(s), y(s)) be a curve parameterised by arclength s. The tangent to that curve makes an angle θ(s) to the x-axis of the Cartesian coordinates in the plane. Then we can write for the unit tangent vector t(s) := r (s) = (cos θ, sin θ), here and in the following the prime denotes the derivative with respect to s (figure 3). We also define the unit normal n(s) = (−sin θ, cos θ).
We require that each point of r(s) to be at the same distance w(s) from two points at the boundaries, i.e. the points r(s) ± w(s)n(s) to belong to the boundary curves r ± (s ± ) each parameterised by its arclength s ± . In other words, we wish the curve r(s) to be a normal centreline for a domain bounded by the given curves which are assumed to be defined by their Whewell equations, i.e. the tangential angles θ ± (s ± ) are known continuous functions of the corresponding arclengths (along with coordinates of the initial points r ± (0)). For the boundary unit tangent vectors ds ± . It may be shown that the normal centreline is defined by the following sixth-order system of autonomous ODEs: Equations (1)-(5) have six dependent variables: half-width w(s), the tangential angle of the normal centreline θ(s), two arclengths s ± (s) that specify the positions at the boundaries of the incident points with the centreline normal line, and the coordinates of the centreline x(s), y(s).
Equations (1)-(3) become singular when θ + (s + ) = θ − (s − )mod π. In these points the chord is orthogonal to both tangents and θ = θ − = θ + . We can use these points to set initial conditions. Let the chords connecting points at the edges be

Normal skeleton of a harbour seal whisker
We have previously scanned several mystacial vibrissae of semi-aquatic mammals [15]. Here we explain how we compute the normal centreline for a harbour seal taken from that database [results for other whiskers can be found in supplementary material (https://stacks.iop.org/BB/17/034001/mmedia)].
The whisker was plucked from the mystacial pad, then placed on the bed of an Epson V600 scanner (Epson, Tokyo) and scanned at 12 800 dpi for a pixel resolution of approximately 2 microns (figure 2). The image segmentation was performed to obtain a set of coordinates of points at the two edges of the 2D projection [15]. We have N − = 26 163 points belonging to the leading edge, and N + = 24 482 points at the trailing edge. We first approximate each of the two edges by a curve of class at least piecewise s ± arc length (0 corresponds to the tip and L ± to the base). We compute the tangential angles θ ± as functions of their arclengths s ± (figure 4) so that the edges are described by the Whewell equations.
The minimal chords connecting the edges are found and the intervals of the normal centrelines between the chords are computed by integrating the governing equations (1)- (5). Then the entire centreline is assembled for the full length L = 35.41 mm (note that to find the length of any (slender) shape, we need to define its centreline first).
The tangential angle of the obtained normal centreline is presented in figure 5, together with its approximation by a parabola. The half-width w(s), together with the approximating polynomial,  is shown in figure 6, where we observe pronounced oscillations for the thick part (ca 2/3 of the total length, see the red dot in figure 6). It is interesting to note that the polynomial has an inflection point that coincides with the emergence of high amplitude oscillations. Close inspection shows that the similar oscillations also occur in the θ function (figure 5) within the same interval.
Indeed, analysis reveals that the oscillations of the curvature θ (s) and of the width correlate. Furthermore, peaks of curvature tend to arise at troughs. To investigate that, we prepared two functions: the half-width w(s) and the smoothed curvature κ(s), computed by application of the LOWESS smoothing of θ (s) (Maple procedure Lowess with bandwidth 0.1). We then centralised them to keep only variations relative to the slowly varying polynomial backbone expressed as low-order polynomials (for the curvature we take the first-order). So, we consider the functions Δw To characterise coordinated variation of the width and the centreline shape, we compute where c ι,δ (s 1 , s 2 ) is the Pearson product-moment correlation coefficient of (1) Δw(s) at the interval [s 1 , s 2 ] and (2) Δκ(s) at the shifted interval [s 1 + δ, s 2 + δ], σ ι + p s 1 < s 2 L − p, σ 1 = 0, σ 2 = s . Parameter p is chosen equal to the half of the period of the width oscillations. We finally find C ι (s 1 , s 1 + ).
The correlation functions C ι ( ) serve as quantitive indicators of maximal lengths of intervals at which  They confirm that variations of the whisker width and of its centreline curvature are coordinated. The above observations have not been reported before. Further experimental studies should be carried out to determine their causes. Nevertheless, from a mechanics perspective, we can make a couple of remarks that we believe may be helpful.
On the one hand, assuming that the whisker geometry we observe is intrinsic, the curvature oscillations must have been caused by a specific growth mechanism: one that coordinates the centreline extension, its bending and the width variation (see section 4). On the other hand, we would expect such a form of the centreline if an initially straight rod with undulating thickness were bent by an external force. We note that the curvature sine-like waves are almost absent in the part of the vibrissa that was grown first, i.e. near the tip (s < s = 12.17 mm, see figure 6), when it was relatively less slender, i.e. when the width increases faster than a linear function. The high curvature undulations for s > s (the younger part of the whisker closer to the base) may have come about because this part is more prone to deformations with amplitudes exceeding the elastic limit. In section 5, we apply a simple beam model to investigate this hypothesis.

Whisker growth
The normal skeleton not only delivers a convenient means to encode a whisker shape, but it may also serve as a formal tool to describe the kinematics of its growth. The growth occurs in the follicle and the vibrissal shafts are made of dead cells. Thus, it can be considered as an intrinsic accretive growth [21]. Although we surmise that the harbour seal vibrissae (similarly to other pinnipeds) grow non-linearly [22][23][24][25], the From this perspective, the tangent vectors t and r ± represent the rates and directions of the growth (the growth velocity vectors) in the middle point and at each of the boundaries, resp., at the 'time' moment s. It may be shown that the tips of all three vectors lie on a straight line, therefore we come to the linear gradient of growth across the normal chord (see the yellow quadrangle in figure 3, spanned by the growth vectors of the points at the normal chord). The growth gradient manifests itself as the curvature of the centreline. Variations of this gradient with 'time' (i.e. along the centreline) encode non-uniformity of the centreline curvature [26][27][28]. Since such non-uniformity is observed, we conjecture that it is caused by corresponding oscillations of the growth gradient whilst the whisker is growing. The cross-sectional growth gradient that controls the centreline curvature should be steeper near the troughs.

Deflection of a straight model whisker
Alternatively, the appearance of the local maxima of the centreline curvature associated with the trough points may be caused by plastic deformations of the whisker material when the whisker is deflected under action of external forces, including contact forces with environmental objects and hydrodynamic drag. In this section we demonstrate, on a simple model, that a deformed whisker with undulating profile bends non-uniformly so that high strain and stress are localised near the points of local thinning.
We consider an intrinsically straight elastic rod of length L = 35.41 mm with thickness 2w(s) mm, where s arclength, s ∈ [0, L], (figure 6). We assume that the normal cross-section is an ellipse with major axis w(s) and minor axis b(s) = q(2w(s) − w(s)), wherew(s) is the sixth-order polynomial linear fit (defined in the caption of figure 6) and q < 1 is a positive dimensionless parameter that characterises an average noncircularity of the cross-section. This model is consistent with that described in [14]. The second moment of area with respect to the axis z, normal to the plane of the whisker, equals I(s) = π 4 w 3 (s)b(s) and the bending rigidity of the rod B(s) = E(s)I(s), where E is the Young modulus. Following [7], we approximate the latter as a linear function E(s) = E 0 + E 1 s/L with E 0 = 2.0 GPa and E 1 = 4.0 GPa.
The planar shape of the centreline is governed by x (s) = cos ϑ(s), where ϑ is the tangential angle, M(s) is the bending moment and N equals a normal force applied to the rod's centreline at some point in the xy-plane [29], e.g. if applied at the tip, then this represents a standard cantilever. We numerically integrate the above system with the initial conditions x(s 0 ) = x 0 , y(s 0 ) = y 0 , M(s 0 ) = 0, ϑ(s 0 ) = ϑ 0 over the interval [s 0 , L], where s 0 is the point of application of the force N. For certainty, we choose s 0 = s . Figure 9 shows a solution to equations (6)-(9), for the case q = 1/2, N = −0.1 N. We observe an oscillatory behaviour of the curvature M(s)/B(s) with its minima (maxima) at the same arclength coordinates as maxima (minima) of the w(s). This depends on neither the value of the parameter q nor the position s 0 . The graphs of the curvature multiplied by the half-width w(s) and, in addition, by the Young modulus, better represent the strain and traction near the whisker surface (see black and magenta curves, resp., in figure 9). The whisker material must yield if and where the strain exceeds of the elastic threshold. Note that the graphs are, on average, less inclined to the horizontal axis compared to the curvature (green), i.e. the traction extreme values are more uniform along the length of the whiskers, which suggests that the elastic limit may be violated at many isolated locations at the same time.
Note that the assumption here of an intrinsically straight rod, does not detract in any way from our conclusions. Rather, an intrinsically curved rod would serve just as well in demonstrating that elastic-plastic deformations could be the cause of these oscillations in the curvature.

Concluding remarks
We additionally tested a whisker of another individual harbour seal and found essentially the same pattern of oscillations of the curvature function as that reported above (see supplementary material). We also computed the normal skeletons for whiskers of two other species: grey seal and California sea lion (Zalophus californianus). The shape of the first also has undulations though not so pronounced as that of the harbour seal. The oscillations of the centreline curvature look random and their periodicity is not observed. The whiskers of the second species lack periodic variations of width; correspondingly, the centreline curvature does not show any signs of periodicity.
In a previous study, involving a wide variety of mammalian whiskers, we adopted a different procedure for the computation of an approximate normal skeleton [15]. That procedure was designed on the assumption that both the curvature and width are some parameterised functions of arclength. We showed that linear functions worked quite well in most cases. However, the linear functions are obviously not able to characterise oscillations neither of the centreline curvature nor of the whisker thickness. It is straightforward to use other functions that potentially provide better accuracy, though the penalty is likely an increase in the number of fitting coefficients. Here we have employed a novel approach, one that does not require imposing any conditions on either the possible shape of the centreline or the width function. The normal skeleton delivers an efficient formalism for an accurate shape representation, for mechanical modelling of a whisker as a slender structure and for growth description. The normal skeleton concept can be generalised to describe a 3D shape. We believe that because of the aforementioned features, the new shape descriptor can be helpful in the design of artificial whisker sensors, and other slender structures, inspired by the intricate geometry of harbour seal whiskers [8].