Planar Covariation of Hindlimb and Forelimb Elevation Angles during Terrestrial and Aquatic Locomotion of Dogs

The rich repertoire of locomotor behaviors in quadrupedal animals requires flexible inter-limb and inter-segmental coordination. Here we studied the kinematic coordination of different gaits (walk, trot, gallop, and swim) of six dogs (Canis lupus familiaris) and, in particular, the planar covariation of limb segment elevation angles. The results showed significant variations in the relative duration of rearward limb movement, amplitude of angular motion, and inter-limb coordination, with gait patterns ranging from a lateral sequence of footfalls during walking to a diagonal sequence in swimming. Despite these differences, the planar law of inter-segmental coordination was maintained across different gaits in both forelimbs and hindlimbs. Notably, phase relationships and orientation of the covariation plane were highly limb specific, consistent with the functional differences in their neural control. Factor analysis of published muscle activity data also demonstrated differences in the characteristic timing of basic activation patterns of the forelimbs and hindlimbs. Overall, the results demonstrate that the planar covariation of inter-segmental coordination has emerged for both fore- and hindlimbs and all gaits, although in a limb-specific manner.


Introduction
Quadrupeds have the ability to generate adaptive coordination and exhibit versatile gait patterns (walk, trot, pace, bound, gallop, etc.) in order to move at different speeds or under different environmental conditions. All gaits are either symmetrical or asymmetrical, and involve a lateral-sequence or a diagonal-sequence [1]. Several modes of legged terrestrial locomotion can be simplified in terms of two general mechanisms, a pendulum and a spring-mass system, which are utilized either separately or in combination [2][3][4][5]. In addition, for both terrestrial and aquatic locomotion, midline stabilization and maneuverability can be increased by controlling side-to-side, mutually opposing forces [6]. Dynamic changes in the gear ratio of muscles can also enhance the performance of skeletal muscles by maintaining them at the shortening velocities that maximize their power or efficiency in trotting and galloping [7].
The inter-segmental limb coordination typically shows adaptive behavior during different gaits [16,35]. In particular, a planar covariation of the temporal changes of limb segment elevation angles has been demonstrated during different gaits in a few animal species, including humans [14], Rhesus monkeys [9] and cats [36,37]. Studying how the coordination patterns change in different locomotion conditions can lead to a better understanding about how the central pattern generators (CPG) control the timing of inter-limb coordination. In this respect, aquatic locomotion may represent a distinctive condition for revealing the intrinsic properties of the inter-and intra-limb coordination, without the constraint of the terrain substrate for the limbs during stance. Evolutionary conservation of ancestral neural networks [38][39][40][41] involves both biomechanical and neurophysiological aspects of quadruped limb coordination. The inter-limb phase and the inter-segmental coordination pattern can be controlled by symmetrically organized unit burst generators for each joint, limb segment, or groups of muscles [42,43] and may emerge from the coupling of neural oscillators with limb mechanical oscillators [35,37]. Therefore, investigating both inter-limb and inter-segmental phase patterns may characterize phase relationships between neural oscillators and CPG organization in different gaits [44]. The application of principal component analysis and neural networks to myoelectric signal analysis may capture the general structure of neural control for generating rhythmical motor patterns [45][46][47][48][49][50]. Indeed, the dynamic behavior of the musculo-skeletal system can be modelled through a linear combination of a small number of basic muscle activation patterns that reflect the shaping function of CPG [51]. Muscle activations tend to intervene during limited time epochs [52] and the biomechanical correlates of each activation pattern have been described [53][54][55]. It is worth noting that the characteristics timing of muscle activation is gait dependent [56,57], although their inter-limb dependence received less attention [58].
In this study, we aimed at comparing the inter-segmental coordination pattern during various forms of quadrupedal locomotion. To this end, we recorded the kinematics of different gaits (walk, trot, gallop, and swim) of six dogs (Canis lupus familiaris), all of the same or closely related breeds (Golden or Labrador Retrievers) and similar size (weight *35 kg, height *57 cm at the withers). An important feature of our study is that the behavior of all the animals was observed outdoors in a naturalistic setting. Because we used a completely non-invasive (markerless) field-recording, some control of accuracy was inevitably sacrificed. However, we separately verified the results obtained with the field-recording with those obtained with a high-performance 3D motion-capture system, and found good agreement. We examined the inter-limb coupling and the inter-segmental coordination in forelimbs (FL) and hindlimbs (HL) and, in particular, the planar covariation of limb segment elevation angles using the model of the tri-segmental limb [9,14,22,37]. Finally, to investigate a characteristic timing of muscle activation and its potential link to the observed differences in the inter-segmental phase pattern of FL and HL, we used available published data on muscle activity during canine locomotion [59] and decomposed them into basic activation patterns using a factorization algorithm [60].

Ethics statement
No special permission is required in Italy for non-invasive observation of animals (here, the dogs) outside laboratory settings in behavioral studies like the present one (Italian Normative 26/2014). All recordings were made in collaboration with members of the "SICS, Scuola Italiana Cani Salvataggio Tirreno" during normal training procedures of dogs. Permission for video recording of dog behavior was obtained from dog owners who were briefed as to what type of gait would be encouraged for the dog and gave consent before the test could commence.

Animals and protocols
We used six healthy domestic dogs (Canis lupus familiaris): three Golden Retrievers and three Labrador Retrievers (weight 35±4 kg [mean±SD], height at the withers 0.57±0.4 cm, see Table 1). All recordings were made in collaboration with members of the "SICS", that provided the dogs. SICS had previously trained the selected dogs for companionship and rescue exercises, but it did not train them to use one or another specific gait. Consequently, all dogs were very compliant to the instructions during the experiments on ground and in water. In particular, we were able to instruct the dogs to perform each locomotion task in the required direction, i.e., roughly orthogonal to the direction of the recording video camera. For the terrestrial locomotion session, the dogs walked, trotted and galloped in a large horizontal open space at their preferred speed. Data collection took place outdoors, in areas where there was ample space for dogs to move. For the swimming session, the recordings were performed during a dog show, in which the dogs were swimming in a special pool (*10 m × 5 m) with transparent walls that allowed the audience to see their movements under water. For galloping and swimming, the dogs ran and swam alone towards the handlers at the opposite end of the pathway (rescue dogs are trained to reach as fast as they can a person that is yelling for help). For the other gaits, the handlers walked or ran aside the dog without any leash restriction. For walk, the dogs followed the walking handlers after a vocal command specific for this gait. For trot, the command was similar, but the handlers encouraged the dog to maintain a trot gait. Specific gait in dogs is typically determined by speed [30] and its categorization (Table 1) was confirmed a posteriori: all dogs displayed a lateral sequence walk, a trot with synchronized diagonal limb movements and a transverse gallop with a forelimb-initiated aerial phase [59]. a The height is measured at the withers during standing, and corresponds to the distance between the ridge of the scapular blades and the ground.

Data recordings
The recordings were made using a Fujifilm Camera (FinePix SL1000, at 60 Hz) in order to obtain two-dimensional coordinates of selected landmarks ( Fig 1A). The dogs moved in a direction roughly perpendicular to the optical axis of the recording camera to minimize errors in 2-D kinematic analysis [61]. The distance between the camera and the dog was about 8-10 m during terrestrial recordings and *5 m for aquatic locomotion, allowing us to record about 2-7 strides in each trial depending on the stride length (more strides were recorded during swimming due to slower motion). The recordings of locomotion were performed in both directions in order to analyze the kinematics of both the right and left side. For symmetrical gaits (walk, trot, swim), the kinematic data of the left and right limbs were pooled together because both sides have similar locomotion characteristics and functions. For the asymmetrical gait (gallop), we analyzed separately the gaits in which the limb was acting as either the trailing limb (i.e. the first limb to touch down) or the leading limb (i.e. the second limb to touch down).
For two dogs, we failed to record the kinematics of the leading or trailing limb during gallop, as it did not face the camera.

Data analysis
The types of gait were manually distinguished. We identified successful sequences of stepsthose when the gait occurred in the dog's sagittal plane steadily and on a straight path (starts and stops excluded) roughly perpendicular to the optical axis of the camera. On average, for each terrestrial gait we analyzed 14±7 successful strides per animal (247 strides total), and for swimming we analyzed 9±2 successful strides (54 strides total, Table 1). While the number of recorded strides differed somewhat between gaits/animals (Table 1), the inter-stride kinematic data for dog locomotion are known to be repeatable (e.g., [62]). Moreover, it is worth noting that a comparison between FL and HL segment planar covariation was performed at the same speed and on about the same number of strides. Finally, we provided information about the orientation of the covariation plane for individual strides (see Results).
Once we obtained the video, the reconstruction was performed using the software Tracker (v.4.85), a video analysis and modeling tool built on the Open Source Physics Java framework. The anatomic landmarks of the ipsilateral side (with respect to the camera) tracked in the reconstruction were: hip joint, knee joint, ankle joint, metatarsophalangeal joint (MP), endpoint (digital tip) of the hindlimb (HEP), scapular fulcrum, shoulder joint, elbow joint, wrist joint, and endpoint of the forelimb (FEP). These landmarks were used for the kinematic analysis and assessment of the inter-segmental coordination in the fore-and hindlimbs. Tracking of scapula during swimming was sometimes problematic (since the scapula was in the close vicinity of the surface of water) Nevertheless, in all analyzed strides, its waveform showed similar back-and-forth angular oscillations. To assess the diagonality of gaits (the sequence of foreand hind-footfalls), the endpoints of the contralateral forelimb and hindlimb were also tracked. The trunk length of each dog (scapula-hip, Fig 1A) was measured prior to the experiments and was used as a metric scale to convert the 2-D video coordinates into real-world 2-D coordinates. Since the trunk length changes during locomotion (especially during galloping), we computed its mean length across each trial and used it for scaling under the assumption that it corresponds to that during quiet standing. While there might be a small discrepancy between the mean trunk length during locomotion and standing, this may only slightly influence the estimated speed ( Fig 1C). Importantly, the estimates of angular movements of the limb segments and the principal component analysis are not affected by trunk length changes.
We first low-pass filtered (10 Hz, fourth-order dual-pass Butterworth) the kinematic data, and then we used a model-based algorithm that optimizes the locations of joint centers by constraining changes in the limb segment lengths (similar to the concept of a template-skeleton for accurate tracking with markerless motion capture systems, [63]. As a template, we used the average segment lengths calculated over all frames of the trial. The average accuracy of kinematic reconstruction, assessed as the mean coefficient of variation of the limb segment lengths during recorded locomotion, was 0.0012 ± 0.0011% (pooling together the data for all segments and gaits). We will describe further checks on the accuracy of kinematic analysis in the section "Potential inaccuracies of the markerless procedure for segment motion reconstruction.".

General gait parameters
The gait cycle for each limb was defined as the time-interval between two successive maxima of the horizontal motion of the distal point of the respective limb (relative to the most proximal marker of the limb). Speed was calculated from the distance the dog (hip landmark) covered during the stride. Duration of rearward (relative to the body) limb movements (that is the stance phase during terrestrial locomotion and the power stroke during swimming) was calculated for each limb. Limb endpoint excursion was determined separately for fore-aft and up-down (relative to the body) movements. To compare different dogs, the calculated values were normalized to hindlimb length (L). To compare the waveforms of limb segment motion across gaits (walk, trot, gallop and swim), we time-normalized (stretched) the data of all gaits to the same relative stance duration as that of walking. To this end, all gaits were scaled to 100 points originally and then scaled/normalized again to the same stance duration as walking (58% gait cycle) using the 'interp1.m' function of Matlab (so that all strides were again finally scaled to 100 points).

Inter-limb coupling
To evaluate the inter-limb coupling, the phase lag (PL) between limbs was determined using methods previously described [1,[64][65][66]. In brief, the relative timing of limb cycle onset was expressed as a percentage of the gait cycle: where Δt i is the interval of time between the cycle onset of the hindlimb ipsilateral relative to the recording camera and the cycle onset of the i th limb (i = 2,3,4), and T is the cycle duration of the hindlimb. According to this method, lateral sequence footfalls (ipsilateral fore/hind limb cycle onsets at similar instances) are determined at a value of 0%, while diagonal sequence footfalls (contralateral fore/hind limb cycle onsets at similar instances) are determined at a value of 50%. Intermediate values (*25%) correspond to no limb pairing. For galloping, PL was calculated with respect to the limb ipsilateral to the recording camera and acting as the trailing limb.

The tri-segmented limb
The model of the dog limbs consisted in an interconnected chain of rigid segments (Fig 1A), following the tri-segmental limb scheme of mammals [67]. The analyzed segments were: trunk (scapula spine-hip), thigh (hip-knee), shank (knee-ankle), foot (ankle-MP), scapula (scapula spine-shoulder), upper arm (shoulder-elbow), lower arm (elbow-wrist) and hand (wrist-FEP, including toes). The elevation angle of each segment in the sagittal plane corresponds to the angle between the projected segment and the vertical ( Fig 1B). The angles are positive when the distal marker is located anterior to the proximal marker. Movements of the toes in the hindlimb contribute relatively little to the kinematics of forward progression and this distal part of the hindlimb is relatively short. If we omit these most distal segments, the hindlimb is functionally tri-segmented [67]. Therefore, from this functional perspective, we used the serially homologous forelimb segments for comparing fore-and hindlimbs: upper arm-thigh, lower armshank, and hand-foot ( Fig 1B). Nevertheless, the scapula also undergoes significant rotations in the sagittal plane in most mammalian groups during locomotion. Therefore, in a complementary analysis of the inter-segmental coordination we included a proximal segment (scapula) in the serial tri-segmental forelimb model [67]: scapula-upper arm-lower arm ( Fig 1B).

Potential inaccuracies of the markerless procedure for segment motion reconstruction
Even though the Tracker software allowed a reliable reconstruction of the recorded body landmarks during both terrestrial and aquatic locomotion (see Data analysis above), we compared the results of the analysis of inter-segmental coordination obtained from the video camera with those obtained with a high-performance 3D motion-capture system. To this end, in a separate set of experiments, we recorded walking in one dog at a natural speed over-ground using simultaneously two different systems: a) the same Fujifilm camera used for the main study, and b) the SIMI Motion system with 3 cameras (Unterschleissheim, Germany, at 100 Hz). The latter system monitored infrared reflective markers (diameter 2.4 cm) firmly attached to the appropriate anatomical landmarks ( Fig 1A). We interpolated the data of SIMI Motion at the same rate as the data acquired with the Fujifilm camera. In a previous study [62], it has been shown that the kinematic data of the sagittal motion of canine hindlimbs during walking obtained with a 2D system correlate well with those obtained with a 3D system, and the data obtained with the 2D system are repeatable. In our analysis, we found that the waveforms of limb segment elevation angles estimated independently by the two recording systems were very similar: the mean RMS difference between the angular waveforms obtained by the two systems was 3±2°(range 1-7°), with the highest difference (7°) for the hand segment and the lowest difference (1°) for the thigh and scapula segments. Nevertheless, the 7°error for the hand segment corresponds to only *4.5% error (since the range of motion of the hand segment was *150°, see Results). Furthermore, an average correlation coefficient between angular waveforms estimated by the two systems was 0.984 (obtained by pooling all segments and all steps together, range 0.985-0.992). These waveforms were also very similar to those obtained during walking in the main experimental session of this dog (when we recorded all three overground gaits using the Fujifilm camera), with an average correlation coefficient of 0.981. Thus, the adequacy of markerless 2D recordings was supported by the low RMS differences and high correlations calculated between the angular waveforms obtained by the two systems and in different experimental sessions. It is also worth noting that, since the main focus of the study was on the intersegmental coordination estimated by angular motion covariation (principal component analysis, see below), high correlations (*0.99) resulted in almost identical characteristics of the covariation plane orientation obtained using the two systems.
As a further test of the adequacy of the kinematic analyses employed for dog locomotion, we added a random noise (±1.5 cm) to the time-series of the coordinates of all reconstructed anatomic landmarks in one dog, and verified how this noise affected the parameters of the inter-segmental coordination. We found that the average correlation coefficient between the original and the noisy angular waveforms was 0.99 (all elevation angles and all gaits being pooled together). Moreover, the parameters of the inter-segmental coordination (u 3 and PV 3 , see below) were very similar: the index of planarity (PV 3 ) increased only by *0.3% and the azimuth and elevation angles of u 3 (orientation of the covariation plane) changed only by 4°on average.

Inter-segmental coordination
The kinematic data were time interpolated over individual gait cycles to fit a normalized 100-point time base. The time course of the elevation angle of each limb segment was expanded into a Fourier series using the fft routine of Matlab. Phase and percent of variance of the first harmonic were computed. The inter-segmental coordination of the elevation angles of hindlimbs and forelimbs segments (thigh, shank, foot, and upper arm, lower arm, hand, respectively) was evaluated in position space as previously described using the principal component analysis [14,35,68]. In humans and primates, the temporal changes of the elevation angles covary during walking [9,14,16,68]. When these angles are plotted in three dimensions (3D), they describe a path that is least-squares fitted to a plane over each gait cycle. Here, we computed the covariance matrix of the ensemble of time-varying elevation angles (after subtraction of their mean value) over each gait cycle. The first two eigenvectors identify the best-fitting plane of angular covariation. The third eigenvector (u 3 ) is normal to the plane, and defines its orientation. The planarity of the trajectories was quantified by the percentage of total variation (PV 3 ) accounted for by the third eigenvector (for ideal planarity, PV 3 = 0% and the 3 rd eigenvalue = 0). For each dog and gait, the parameters of planar covariation (u 3 and PV 3 ) were averaged across strides. The u 3 vector was averaged across animals using spherical statistics on directional data, the direction cosines of the mean vector ( x, y, z) being defined as [69]: where The inter-subject variability in u 3 and PV 3 was also assessed. Variability in PV 3 was assessed as SD across animals. To assess variability of u 3 , we calculated the confidence cones centered on the u 3 mean direction. Briefly, for each gait we calculated the points corresponding to the projection of the normal to the plane for each dog onto the unit sphere, the axes of which are the direction cosines with the semi-axis of the thigh, shank, and foot. The 2D distribution of these points for each gait in the plane orthogonal to the u 3 mean direction was quantified using the appropriate scaling factor for the 95% confidence ellipse depending on whether the data had a Fisher or Kent distribution [69,70]. The resulting 95% confidence cones, based on these confidence ellipses [70], were drawn in the 3D space defined by the elevation angles to characterize spatial distribution of the normal to the covariation plane in each gait and for each limb. We also computed the area of these ellipses to quantify the variability of u 3 .

Basic activation patterns of EMG profiles during walk, trot and gallop
To characterize a gait-specific timing of muscle activation and its potential link to the differences in the planar covariation of FL and HL segment motion, we computed the basic activation patterns derived from averaged EMG profiles of 23 extrinsic FL and HL muscles using previously published data from 12 mixed-breed dogs while they walked, trotted and galloped on a level treadmill [59]. Published graphs were scanned, digitized manually taking the touchdown of each limb as time reference, time-interpolated to fit a normalized 100-points time base, and low-pass filtered using a sliding window of ±7 points (in order to avoid the boundary effects, the signal waveform was replicated and concatenated prior to averaging). These data included average EMG recordings from the following ipsilateral muscles: HL (12 muscles The hypothesis that muscle activity profiles (m i ) during different gaits are generated by the nervous system through a linear combination of a small number of basic activation patterns (p j ) [54,56] is quantitatively formulated in terms of the equation: where w ij are weighting coefficients. We applied a non-negative matrix factorization of the data using the algorithm described by Lee and Seung [71] that constrains the basic patterns and weights to be nonnegative (the rationale was that our data consisted of the non-negative values of rectified EMG activity). To determine the number of significant basic patterns p j , we performed an iterative reconstruction of the EMG profiles using j = 1,. . .,6 patterns, until the cumulative variance accounted for by these patterns was closest to 90%, that is, the residual error accounted for~10% of data variance.

Statistics
Descriptive statistics included the calculation of the mean and SD. For each subject, the parameters were averaged across cycles for the subsequent statistical analysis. Shapiro-Wilk test was used to verify the normality distribution of data. One-way (effect of gait) or two-ways (effects of gait and limb) repeated measures ANOVA were used to evaluate differences in general gait parameters and kinematics. Since for galloping we analyzed separately the gaits in which the limb (ipsilateral to the video camera) was acting as the trailing limb or the leading limb, we recorded the trailing limb in 5 dogs and the leading limb in 4 dogs. Missing data for the ANOVA were replaced by the unweighted mean value estimated from all other dogs. If ANOVA resulted in a significant effect for gait (and/or limb when assessed), then post-hoc tests and multiple comparisons analysis were performed by means of Tukey HSD test. In the Results, we report only selected important differences related to general gait parameters and kinematics. One-tailed t-test was used to verify whether the phase lag between the contralateral limbs significantly differed from 50% (symmetrical gaits). Statistics on correlation coefficients was performed on the normally distributed, Z-transformed values. A Matlab Toolbox for circular statistics [72] was used to characterize phases of the first Fourier harmonics of angular waveform and provided a parametric Watson-Williams test to compare them. Statistical analysis of spherical data [69] was used to characterize the mean orientation of the normal to the covariation plane (see above) and its variability across steps. Reported results are considered statistically significant for p<0.05.

General gait parameters
Gait speeds and cycle durations are plotted in Fig 1C and 1D. Swimming was the slowest gait, but its typical speed (1.44±1.15 m/s) was relatively high taking into account water resistance (Fig 1C), and the cycle duration was not significantly different from that of walking (Fig 1D, p = 0.45, Tukey HSD test). A cycle can be divided in two phases: backward (relative to the body) limb endpoint excursion, and forward excursion. During swimming, these two phases can be labeled as a virtual 'stance' and 'swing' phase respectively, and rearward leg movements are considered as power strokes and forward movements as return movements. During walking, the relative duration of rearward limb movement was longer than 50% of gait cycle, while for all the other gaits it was less than 50% ( Fig 1E). There was also asymmetry in this parameter between FL and HL for swimming (symmetrical gait) and for the trailing limb in galloping (p<0.00002 and p = 0.05, respectively, Fig 1E). Horizontal limb excursions were significantly higher for FL than for HL in all gaits (p<0.03, Tukey HSD tests) except for comparing the leading FL and HL in galloping (p = 0.11) (Fig 2B). Vertical limb excursions were also significantly higher for FL than for HL in all gaits (p<0.04). Moreover, it is worth noting that they exceeded 1 L for the forelimb in swimming, likely in relation to the need to provide a sufficient vertical stroke to avoid sinking and to keep the head outside the water.

Inter-limb coupling
The fore-aft movement of the limb endpoints was used to determine the gait patterns (Fig 3A  and 3B). In general, the dogs maintained a 1:1 frequency relationship between the limbs, despite some inter-stride variability in the onset of footfalls. Fig 3B shows the phase lag (PL) between the limbs determined as the relative timing of the limb cycle onset (relative to HL ipsilateral to the recording camera), expressed as a percentage of the gait cycle. For the pairs of contralateral limbs, PLs were on average 49.8±1.9% for FLs and 50.2±1.3% for HLs for walking, trotting and swimming, as one would expect for symmetrical gaits, whereas they were significantly smaller than 50% for galloping (18.3±4.5% for FLs and 26±23.1% for HLs, p<0.001 for both limbs, one-tailed t-tests). For the overall pattern during symmetrical gaits, the results showed significant variations in the inter-limb coordination: the canine gait pattern ranged from a lateral sequence of footfalls during walking to a diagonal sequence in trotting and swimming (Fig 3B). Fig 4A shows the average joint and elevation angle profiles (across all dogs) plotted as a function of the normalized gait cycle. In general, the angular waveforms differed across gaits, in particular due to differences in the relative duration of rearward and forward limb movements. In addition, there were differences in the amplitude of angular motion between HL and FL. For instance, the range of motion of the most distal joints differed by as much as *2-3 times (ankle vs. wrist), depending on the gait (Fig 4B). For the elevation angles, there was a similar tendency for the most distal segments (foot vs. hand). The smallest range of motion among different elevation angles was observed for the scapula (Fig 4B).

Joint and segment angular motion
There were also similar features in the kinematic patterns across gaits (Fig 4A). These general features of angular waveforms were also similar to those of walk and trot obtained in the previous studies [30,32,73,74]. The angular waveforms of trot, gallop and swim correlated much better with those of walking if they were time-normalized (stretched) to the same relative stance duration as that of walking (see Methods, Table 2). After normalization, stance duration of the gaits different from walking increased, while the swing duration decreased (according to Fig 1E) to match the same (58% gait cycle) relative stance duration. On average (all segments and gaits being pooled together), the correlation coefficient was 0.54 for non-normalized waveforms and 0.91 for normalized ones. Also, there was a similar temporal sequence of minima in the elevation angles across gaits for each limb (marked schematically by green arrows in Fig  4A). On the other hand, the waveforms were substantially different for HL and FL in all gaits (Fig 4A). Fig 5 illustrates the characteristics of the first harmonic and the temporal sequence of the minima in the elevation angles. The first harmonic accounted for a considerable proportion of data variance for all segments and gaits (on average, 87.8%, Fig 5A). Thus, its relative phase captures basic phase relationships of the inter-segmental coordination. Fig 5B (bottom panels) shows that there was a phase-lead of the thigh segment with respect to the other hindlimb segments in all gaits (on average, the phase was 37°±13°for thigh, 57°±9°for shank, and 45°±11°f or foot), and a phase-lead of the lower arm segment with respect to the other forelimb Inter-Segmental Coordination in Dog Locomotion segments of FL (on average, 44°±17°for upper arm, 38°±13°for lower arm, and 53°±14°for hand). The parametric Watson-Williams test confirmed that the phases of paired segmentsfoot vs thigh, shank vs foot, upper arm vs lower arm and hand vs upper arm-were all different (p<0.02) except for foot vs shank of the leading HL in galloping (p = 0.46). Thus, the temporal sequence of the phase of the first harmonic was maintained in all gaits, namely: 'thigh-footshank' for HL and 'lower arm-upper arm-hand' for FL, consistent with a temporal sequence of minima in the elevation angles (Figs 4A and 5B).  computing the percentage of variance accounted for by the third eigenvector (PV 3 ) of the data covariance matrix: the closer PV 3 is to 0, the smaller the deviation from planarity. The results demonstrated that the planar covariation was generally maintained for all gaits (PV 3 = 0.66-3.49%, Fig 6D), though PV 3 depended on gait (F(4, 20) = 10.204, p = .0001, RM ANOVA). A post-hoc Tukey test showed that PV 3 was significantly lower in swimming compared to trotting and leading and trailing limbs in galloping (p = 0.007, p = 0.0001 and p = 0.001, respectively, Fig 6D, and tended also to be lower than walking though not-significantly, p = 0.1), consistent with better planarity of the 3D gait loop in the absence of foot-support interactions in humans [75].

Planar covariation of limb segment elevation angles
The best-fitting planes of the corresponding loop trajectories are illustrated in Fig 6A. The third eigenvector (u 3 ) of the covariance matrix is the normal to the plane, and thus characterizes its orientation. Fig 6C shows the direction cosines of the plane normal for all individual strides and dogs as a function of locomotion speed, while Fig 6E presents the direction cosines of the plane normal averaged across dogs for each gait and limb. Despite some inter-stride and inter-subject variability ( Fig 6C) and differences in the waveforms of angular motion (Fig 4), the direction cosines of the covariation plane tended to be similar across gaits (Fig 6E). To describe the plane orientation, we also analyzed the 95% confidence cones that characterize the spatial distribution of the normal to the covariation plane between dogs (see Methods). The area of the confidence ellipse represents the inter-subject variability (for HL 0.056, 0.080, 0.095, 0.063 and 0.274 and for FL 0.252, 0.199, 0.157, 0.277 and 0.113 for walk, trot, leading and trailing limbs in gallop and swim, respectively,). There was an overlap of confidence cones across different gaits for each limb (Fig 6B, left panel-HL, right panel-FL), consistent with similarities of the normalized angular waveforms ( Table 2, right columns).
While the covariation plane orientation tended to be similar for different gaits, it differed systematically between HL and FL (Fig 6B, 6C and 6E). Indeed, the confidence cones of the plane normal did not overlap for HL and FL (Fig 6B), and the mean u 3 vectors were different: for HL, the elevation angle of u 3 was 35°, 33°, 25°, 27°and 29°, and the azimuth angle was -180°, -174°, -172°, -161°and -160°in walk, trot, leading and trailing limbs in gallop and swim, respectively, for FL, the elevation angle of u 3 was 13°, 6°, 6°, 11°and 15°, and the azimuth angle was 156°, 149°, 149°, 147°and 168°, respectively. The analysis of the inter-segmental coordination for the forelimb (Fig 6) was performed using the forelimb segments (upper arm, lower arm and hand) serially homologous to those of the hindlimbs. Since the scapula also undergoes appreciable rotations in the sagittal plane (Fig 4B), in a complementary analysis we verified Inter-Segmental Coordination in Dog Locomotion whether the observed difference in plane orientation between FL and HL ( Fig 6A) holds also in an alternative serial tri-segmental forelimb model [67]: scapula-upper arm-lower arm (Fig 7, even though during swimming the reconstruction of scapula is less precise due to vicinity to the water surface). The covariation plane orientation of this alternative FL model (Fig 7B) was more similar to that of the previous FL model than to HL (cf. Figs 6A and 7). Decomposition of published EMG data during canine locomotion into basic temporal patterns We examined the patterns of muscle activation from a large set of muscle recordings published by [59], in which the average EMG activity from 12 HL and 11 FL muscles was determined for a standard step cycle during walk, trot and gallop (with respect to touchdown of each limb).
Muscle activity tends to occur in bursts with specific timings during a gait cycle (Fig 8A), consistent with 'drive pulse' rhythmic elements or primitives in the spinal circuitry of animals [76,77]. Although the activation patterns in the present data set appear muscle specific, there are clearly preferred phases of activation during the gait cycle. To characterize a gait-specific timing of muscle activation, EMG data were decomposed into basic activation patterns ( Fig  8B) using non-negative matrix factorization (see Methods). We were especially interested in looking for a temporal correspondence of hypothetical pulsatile burst generators for HL and FL for each gait.
characterizing spatial distribution of the normal to the covariation plane between dogs in each gait and for each limb (left panel-HL, right panel-FL) in the 3D space defined by the elevation angles (see Methods). The foot (and hand for FL) semiaxis is positive, and the shank and thigh (lower and upper arm for FL) semiaxes are negative. Note overlapping of confidence cones in different gaits. C: direction cosines (u 3t , u 3s , u 3s for HL and u 3ua , u 3la , u 3h for FL) of the normal to the covariation plane (u 3 vector) for all individual strides and all dogs versus locomotion speed. D: percent of variance (+SD) accounted for by u 3 (PV 3 ). Asterisks denote significant differences. The results showed that four basic temporal patterns accounted for *90% of the total waveform variance across the recorded HL and FL muscles for all gaits (range 88-95%). These basic patterns are illustrated in Fig 8C, designated in chronological order of the main peak relative to the gait cycle. For gallop, the patterns for leading and trailing limbs were superimposed, since pairing limbs in the asymmetrical gait is uncertain. Overall, the timing of basic patterns differed between walk, trot and gallop, as it happens for human walking and running in relation with corresponding changes in the relative stance/swing duration [54]. However, it is worth stressing that, for each gait, there are different phases of muscle activity with respect to touchdown of HL and FL (Fig 8C), in accordance with the previously described differences in the planar covariation of limb segment motion (Fig 6). Fig 8D summarizes a schematic representation of the canine motor program for walking as a sequence of activation pulses [52,55,76].

Discussion
We examined the kinematics of HL and FL movements during quadrupedal locomotion in dogs. A novel finding was that all limbs exhibited the planar constraint of the inter-segmental coordination, even though FL and HL substantially differ in musculoskeletal anatomy and show opposing phase relationships and limb-specific covariation plane orientation (Figs 4-6). The analysis of published muscle activity patterns also suggested a characteristic limb-specific phase control of muscle activation for each gait (Fig 8). Below we discuss the findings in the context of limb-specific control of the inter-segmental coordination in different gaits.

Kinematics of terrestrial and aquatic locomotion
The rationale of our study was to compare various forms of quadrupedal locomotion. Swimming is a relatively frequent behavior among mammals, normally performed when they actively cross short stretches of water. It is valuable to place the locomotion in water in context with that of terrestrial creatures, swimmers, and fliers, and so contribute to the emerging integrative view of biolocomotion [4,[78][79][80]. For instance, locomotion in water in some terrestrial insects, which occasionally swim, shows a stereotyped tripod pattern similar to that used by many insects for running on land [81]. Estimations of the forward thrust from 3D recordings of insect leg movements showed that thrust was mainly produced by the front legs (and to a lesser extent by the middle legs), while the hind legs contributed drag [81]. Nevertheless, several arthropods exhibit different modes of locomotion in/on water and on land. Typically, during the power stroke, the front limbs of swimming animals are extended and moved rapidly, whereas in the return stroke they move slower and are held closer to the body [4]. The hind legs contribute only little to the overall thrust, suggesting that they mainly serve for stabilization and steering. The data on swimming mammals are still somewhat limited [82,83].
Aquatic locomotion represents an interesting example of body movements without the ground constraint due to a common support for the limbs during stance. Our results showed that the sequence of footfalls was basically reversed in walking and swimming: HL-FL-HL contr -FL contr in walking and FL-HL-FL contr -HL contr in swimming (Fig 3). Overall, the canine gait pattern was characterized by a lateral sequence of footfalls during walking and a diagonal sequence in swimming (Fig 3). Early tetrapods from the Devonian period exhibited diagonal-sequence locomotion [84], suggesting that this type of quadrupedalism is a very old locomotor trait [85], has been preserved for millions of years and may involve both biomechanical and behavioral advantages [86].
For swimming, we observed a high planarity of inter-segmental coordination (PV 3 < 1%, Fig 6). The limb strokes for swimming had significantly larger amplitudes than those of walking and trotting (Fig 2). Some variations in the covariation plane orientation across gaits or  [59]. In dark gray-HL, from bottom to top (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12) schematic of motor modules. Simulated example of muscle activity profiles as weighted sum of basic patterns (p j ): m i = ∑ j p j (t)Áw ij + residual. The outputs of the first (green), second (blue), third (red) and forth (violet) modules are summed together to generate overall muscle activation (black envelope). C: basic activation patterns of HL (black) and FL (gray) muscles during walk, trot and gallop. EMG data were decomposed into basic activation patterns using nonnegative matrix factorization. For gallop, the patterns for the leading (solid lines) and trailing (dotted lines) limbs were superimposed. The four basic patterns account for *90% of variance in each limb and gait and are each characterized by a relatively narrow peak of activation (Gaussian-like) at a particular phase of the cycle (with respect to touchdown of each limb). Components are designated in chronological order. D: the timing of the main peaks of 4 basic patterns across gaits and limbs. Zero timing for HL and FL refers to touchdown of HL and FL, respectively. E: locomotion motor program as a sequence of activation pulses [52,54]. The schematic diagram of activation pulse timings corresponds to those of walk. TD = touchdown. Note different phases of basic muscle activity in HL and FL for each gait. steps (Fig 6) could be explained in part by the effect of speed, also observed in humans [68]. Nevertheless, it is remarkable that, despite substantial changes in the duration and amplitude of angular motion (Fig 4) and variable inter-limb phase patterns (Fig 3), the intra-limb coordination pattern was fairly well conserved in different gaits (Fig 6). This apparent discrepancy may reflect different degrees of flexibility of phase connections of neural oscillators that control inter-limb and inter-segment coordination in different gaits. Similarities in the inter-segmental phase pattern across gaits may also explain a relatively rapid dog's adaptation to aquatic locomotion (as soon as a dog is plunged into water it can swim), since a similar locomotor program may underlie phase relationships. In contrast, albeit arm and legs are coordinated [87], humans display quite different styles of coordination during swimming and they typically need to learn how to do it.

Inter-segmental coordination in forelimbs and hindlimbs
Our results extend to canine locomotion previous observations in humans, monkeys and cats that the temporal changes of limb segments with respect to the direction of gravity co-vary according to a law of planar co-variation during a variety of tasks [9,16,35,37,[88][89][90][91].
The spatial orientation of the covariation plane of the hindlimb of dogs is very different from that of human lower limbs (cf. Fig 6 with Fig 2 in [92]. Even though biomechanics contributes to the emergence of the planar covariation of limb segment motion [92], the orientation of the covariation plane reflects basic phase relationships between segment motions, and adaptation to different walking conditions [68,93]. For instance, the orientation of the planar covariation shows systematic changes with walking speed [68] and in different gaits (running, staircase and uphill walking, obstacle avoidance, etc.) in adults [16,89,94,95], while the orientation is relatively fixed in different gaits in toddlers at the beginning of independent walking [93], suggesting its neural control rather than simple biomechanical consequence. Moreover, the planar covariance is markedly different in human newborns [77] and in adults with a spinal cord lesion [92]. In prosthetic gait of transfemoral amputees, the best-fitting plane rotates around the long axis of the gait loop with increasing walking speed even more for the sound limb of expert amputees than in control subjects [94]. The authors of this latter study suggested that their results reveal a centrally commanded compensation strategy. It is also worth stressing that the orientation of the covariation plane is different between dogs and humans (Fig 9), further suggesting different phase relationships between segment movements. Finally, the modelling studies confirmed changes in the covariation plane orientation of simple mechanical oscillators with appropriate phase shifts [17]. In particular, the difference in the orientation between dogs and humans (Fig 9) might be attributable to a digitigrade (dogs) vs. plantigrade (humans) gait [9,16]. In fact, the correlation between the foot and shank elevation angles was much lower in dogs (Fig 4A) than that in humans. The current results therefore suggest that the strategy by which the central nervous system achieves inter-segmental coordination and adapts its spatiotemporal structure in dogs (and non-human primates) differs somewhat from the kinematic principles that operate in human gait control, consistent with early studies on variances in the distal limb segment control in humans and dogs [22].
Interestingly, we observed a limb-specific inter-segmental coordination pattern (Figs 4-7) likely due to a different parametric tuning in the phase-relationship of inter-segmental coordination between FL and HL. Despite differences in inter-limb coupling (the frequent use of diagonal-couplet interlimb timing in primates), the organization of intra-limb coordination during walking in dogs shows a number of similarities with primates [9], suggesting that common principles may operate during stepping among a wide range of mammals. In particular, the orientation of the covariation plane is also different between FL and HL in Rhesus monkey gait (cf. our Fig 6A with Fig 6A in [9]). The general features of angular waveforms ( Fig 4A) were also strikingly similar to those obtained in previous studies of dog kinematics during walk and trot [32,73,74], suggesting that they are characteristic for each joint in dogs of similar morphology [73].
The key anatomical inter-limb difference consists in the orientation of the forelimb and hindlimb; the elbow is facing posteriorly and the knee joint anteriorly, and the scapula is either held by the clavicle or, in most mammalian groups, completely freed from any connection with the trunk [67]. This distinctive functional orientation of limb segments may impose specific parametric tuning in the phase-relationships and specific control of FL and HL segments. The theory of coupled oscillators is a useful tool for studying synchronized, periodic dynamics of physical and biological systems [44,96], based on specific phase resetting and the stabilization of an isolated cycle as the attractor of animal dynamics [10,97]. For instance, the fact that a temporal sequence of minima in the elevation angles around the stance-to-swing transition was different for each limb (Figs 4A and 5) may indicate that the so-called 'leading' segment [98] is also different for the HL and FL control: proximal segment (thigh) for HL and intermediate segment (lower arm) for FL. Barliya et at. [17] introduced a mathematical model that represented the rotations of the elevation angles in terms of simple oscillators with appropriate phase shifts between them. Again, such analysis showed a clear phase-lead of the first harmonic of the thigh and lower arm segments for HL and FL, respectively (Fig 5B), consistent with different covariation plane orientation (Fig 6) and different 'leading' segments [98] for HL and FL.
Hypothetical limb-specific control of inter-segmental phase patterns Differences in the intra-limb coordination strategy (Fig 6) must also be taken into consideration when studying the coupling of neural oscillators with limb mechanical oscillators [35,37]. If one accepts that the inter-segmental coordination pattern can be controlled by symmetrically organized unit burst generators for each joint or synergistic sets of muscles [17,42], the observed findings (Figs 4A, 5B, 6A and 6B) may point to the differences in the organization for cervical and lumbosacral central pattern generators. In particular, it has been proposed that a multi-layered organization of mammalian locomotor CPG includes a rhythm-generating layer Covariation of thigh, shank, and foot elevation angles is shown during walking at a natural speed. Gait loops of each species are represented with respect to the best-fitting planes (grids). Note different covariation plane orientation in human and dog locomotion, attributable to different phase relationships between limb segment oscillations. Data from human were adapted from [92]. and a pattern-generating layer [55,[99][100][101], and the locomotor output can be represented as a characteristic timing of muscle activation [52]. Furthermore, the descending cortical signal interacts with the interneuronal networks in the spinal cord to ensure the appropriate limb movement control and the basic locomotor rhythm [102][103][104]. Subpopulations of motor cortical neurones, active sequentially during the step cycle, may regulate the activity of small groups of synergistic muscles during quadrupedal locomotion, and these synergies, identified by a cluster analysis, are defined by periods of muscle activity that are coextensive [103].
An interesting approach to capture the essence of biological input/output transformations consists in analyzing the spatiotemporal characteristics of the spinal motor output and modelling the locomotor control system using artificial neural networks [45,46,51,55]. The results of this approach demonstrate the ability of neural networks to model the transformation between a kinematic movement plan and the necessary muscle activations [48], even though they do not directly specify the actual structures involved. The dynamic behavior of the musculo-skeletal system can be modelled through a linear combination of a small number of basic muscle activation patterns linked to specific kinematic events [48,51]. Both supraspinal input and sensory information have a major role in determining the timing of the motor activation patterns during steady state locomotion [99,105].
Our analysis of available published data for canine locomotion [59] showed that multi-muscle activity patterns of both HL and FL (Fig 8A) can be decomposed into a small set of 4 basic temporal patterns that account for~90% of total variance, and are each characterized by a relatively narrow peak of activation at a particular phase of the cycle (Fig 8C and 8D). Their timing is not invariant but differs between different gaits in conjunction with changes in the relative stance/swing duration (Fig 1E), as it does for human walking and running [54]. Since we were mainly interested in the temporal structure of the motor patterns, the decomposition analysis was applied to the normalized EMG profiles (they were scaled to the maximum amplitude observed during each gait [59]). Further investigations are needed to elucidate muscle synergies (weights of basic activation patterns, Fig 8B) common to each gait mode. However, it is worth stressing that, for each gait, there are different phases of muscle activity with respect to touchdown of HL and FL (Fig 8C).
The previous works also revealed some limb-specific features in the organization and coupling between FL and HL controllers [20,[106][107][108][109]. Both descending and ascending connections between cervical and lumbosacral CPGs, as well as an intrinsic rhythmogenesis capacity of the thoracic spinal network, have been described in quadrupedal animals [110]. Nevertheless, there are essential differences in their neurotransmitter systems [108] and in the interlimb influences [20]. For instance, forelimb movements in the dog may facilitate or even trigger hindlimb stepping while the opposite influences are much weaker [20]. Inter-limb coordination may also reflect supraspinal control; hindlimb-related neurons in motor cortex respond to changes in forelimb movements during locomotion and vice versa, although the percentage of neurons from the HL and FL area that are modulated by motion of the other pair of limbs is very different [111]. In sum, the observed limb-specific inter-segmental phase pattern (Figs 4-8) may be in accordance with the specific biomechanical function of FL and HL.
Finally, it is worth noting that the intra-limb coordination pattern is relatively conserved (Figs 5 and 6) while muscle activations tend to intervene during different time epochs in different gaits (Fig 8), consistent with the idea that EMG commands are subservient to the kinematic reference. Therefore, the coupling of neural oscillators with limb mechanical oscillators may be more complex than it was previously thought, since the rhythm-generating layer or 'time-keeping function' of the CPG for locomotion [51,55,99] appear to be kinematically-driven. Therefore, although it is often assumed that CPGs control patterns of muscle activity, another hypothesis is that they control patterns of limb segment motion instead [17,37].
Supporting Information S1 Dataset. A compressed archive (.zip) contains two folders and a text file with a detailed description of the data organization. S1_Dataset.zip contains two folders: 1) Files in the 'Dogs Data' folder are the kinematic data in Matlab format (.mat) for each subject and each gait; 2) Files in the 'comparison' folder are the data used for the validation of the markerless motion capture system. (ZIP)