Appendix 1

Coupling length

The interpretation of coupling length starts with selecting four tracks, two pes and two corresponding manus tracks. Since alternative pairs of manus tracks could have been selected to pair with the given pes tracks, the subsequent analysis begins with one such choice of pairings, which will then be generalized to explore other possible pairings of manus and pes tracks.

Next, we construct two couplers, one for the pair of pes locators and another for the pair of manus locators. In general, a coupler is a line segment representing the pairing of two locators (xl, yl) and (xr, yr), and its midpoint (gamma ) is simply given by:

$$ gamma = left( {frac{{x_{l} + x_{r} }}{2},frac{{y_{l} + y_{r} }}{2} } right) $$


The midpoint (gamma ) of a coupler approximates the center of a limb girdle, therefore the coupling length CL is the instantaneous separation between the pes and manus couplers, specifically:

$$ CL = sqrt { left( {x_{gamma ,pes} – x_{gamma ,manus} } right)^{2} + left( {y_{gamma ,pes} – y_{gamma ,manus} } right)^{2} } $$


where ({x}_{gamma ,pes}) corresponds to the x coordinate of the pes coupler, and so forth. Figure 5 shows the pes and manus couplers as line segments between corresponding pes and manus locators, and the line segment CL is shown connecting the midpoints of the two couplers.

Coupling length is subject to measurement uncertainty since it is a computation based on uncertain locator positions. The positional uncertainty of a locator reflects the underlying track measurement error when in support, and that uncertainty is necessarily greater when the limb is stepping between tracks. When the locator is in mid-step, its position is calculated using linear interpolation between the coordinates of the previous and next track centers, using a conventional interpolation variable α (0 ≤ α ≤ 1). Initially (at α = 0), the locator is still at the previous track, and is just finishing support. For intermediate values of α (0 < α < 1) the locator is in transition and its position is interpolated between that of the previous and next tracks. When α reaches 1 the locator is again in support at the next track. Note that α is a relative phase, i.e., relative to when the locator begins transition, which depends on LP and DF and the particular limb. We will need to also compute the absolute phase (i.e., relative to PL) at which to start the interpolation for that locator, which is specific to the limb, the limb phase, and the duty factor. Given that, we can then compute the position of each of the four limbs even if one limb is in transition, at any point within the overall step cycle.

Since we consider only walking gaits (DF > 0.5), both hindlimbs are simultaneously in support during two intervals of each step cycle. We measure CL only at the end of each such interval (when either pes PL or PR just finishes its support phase and begins its next step). This choice of timing will be shown to permit CL to converge to GL. At those two moments within a cycle, we know the position of the pes coupler ({gamma }_{pes}) since both pes locators are in support, placed on their respective tracks. In the following, we describe only the case for PL since the reasoning is similar for PR. Since DF > 0.5 we also know that one manus locator is also in support, hence we know the position of the track on which it stands, and finally, we can also interpolate the position of the other manus locator. With both manus locators known, we can then solve for the position γM of the manus coupler and thus compute CL = (({gamma }_{manus}-{gamma }_{pes})) according to Eq. 4 for any combination of limb phase and duty factor DF ≥ 0.5.

Coupling length for the walk gait

The equation for coupling length will be developed first for the familiar case of a G2 walk (LP = 0.25), then will be generalized to other gaits. Figure 13 shows phase diagrams for the walk for three duty factors (DF = 0.5, 0.75, and 0.8). As in Fig. 2, the blue and red intervals indicate support and transition, respectively. In each case, a step cycle begins with the left pes PL just beginning support, followed by the left manus ML at t = 0.25, the right pes PR at 0.5, and the right manus MR at t = 0.75.

Coupling length is defined twice per step cycle (when either the left or right pes completes support and begins transition), but for purposes of development we consider when the left pes PL just finishes support and begins to step. The state of the four locators is examined at this instant for three examples of duty factor (Fig. 13). For the very short duty factor of 0.5, note that while three locators are in support (indicated by green dots), MR is in transition (red dot), midway between the previous and next track at that instant. Note that for longer duty factors (wherein PL remains in support longer), MR will complete more of its transition and converge on the location of the next track at that moment of measurement. For sufficiently long duty factors (i.e., DF ≥ 0.75), all four locators are in support at that moment.

Fig. 13
figure 13

The step cycle for a G2 walk for three different duty factors (DF), with blue indicating support and red indicating transition states. We consider the case when PL finishes support, and examine the state of the other three locators. For DF = 0.5 (top graph), at that moment, indicated by the colored dots, we see that only three of the four locators are in ground contact (green), while MR is in transition (red) and its location must be interpolated. For DF ≥ 0.75 (lower two graphs) MR has completed its step and all four limbs are in ground contact. The dashed black line indicates the relative phase between MR and PL

The left pes PL begins support at time t = 0.0 and it remains in support until t = DF whereupon it begins transition. At that moment, the right manus MR is still in transition for any duty factor less than 0.75, and its position would have to be interpolated as a fraction α of the distance between the last position for MR and its next position. For the case of a G2 walk (LP = 0.25), the interpolation value αwalk depends only on the duty factor DF, and is given by:

$$ alpha_{walk} = frac{0.25}{{1.0 – DF}} $$


The numerator of (5) is the relative phase between MR and PL, indicated by the dashed black line in Fig. 13 highlighting the fraction of a step cycle between the end of support for PL and the end of support for MR. The denominator of (5) is the duration in which the limb is in transition (i.e., it is the relative complement of the duty factor). The quotient αwalk is therefore the fraction of the way that MR has traveled from the previous track to the next track. For DF = 0.5, MR is half the distance between the previous and next tracks, and for any value DF ≥ 0.75 the interpolation would be complete, and MR would be at the next track. What are the consequences of interpolation on the coupling length computation CL? Since PL and PR are in support, the pes coupler ({gamma }_{pes}) is fixed, and ML is also in support at this moment, so that as MR converges on the next track as DF increases, CL converges on GA.

Figure 14 shows this diagrammatically for DF = 0.5 (upper trackway diagram) and DF = 0.75 (lower). The gait in both cases is a G2 walk (LP = 0.25). Both diagrams show the configuration of the four locators at the instant the left hindlimb is about to initiate another step. For DF = 0.5 the right manus is in transition midway between two manus tracks while for DF = 0.75 its transition phase is complete, all locators correspond to track positions, and CL converges to the conventional GL computation (Fig. 1). For values of DF > 0.75, CL would remain as shown in the lower diagram of Fig. 14.

Fig. 14
figure 14

An idealized trackway with two simulations wherein the left pes is about to finish the support phase. In the upper diagram the duty factor is only 0.5 and therefore the right manus MR is still in transition. In the lower diagram, the duty factor is 0.75 and now MR has just completed its step. Notice the difference in CL for the two cases

Generalization to other gaits

The above discussion, while specific to the G2 walk, also applies to the other lateral sequence gaits such as G1 and G3. Recall that in the case of the walk, depending upon the duty factor, the right manus might still be in transition when we need to compute the location of the manus coupler. We interpolate the position of MR between the previous and next right manus tracks using an interpolation value αwalk. While this interpolation value depended only upon the duty factor in Eq. 5. The generalization of (5) to apply to other gaits will depend on both LP and DF.

To illustrate, Fig. 15 diagrams the G1 gait (lateral sequence lateral couplet, LP = 0.125), for three duty factors (DF = 0.5, 0.6, and 0.9). Recall that the numerator in (5) was graphically indicated by the dashed black line in Fig. 13. The generalization of the numerator, which will be called the α-phase and given the variable ({phi }_{alpha }) below (in Table 3 and Eq. 6), is graphically indicated by the dashed black line in Fig. 15.

The generalization of the interpolation value α is now an equation with two unknowns, the α-phase ϕα (which is itself a function of LP and DF), and the duty factor DF:

$$ alpha left( {phi_{alpha } , DF} right) = fleft( x right) = left{ { begin{array}{*{20}c} {1 ;left| {if phi_{alpha } ge 1.0 – DF} right.} \ {frac{{phi_{alpha } }}{1.0 – DF};left| {otherwise} right.} \ end{array} } right. $$


Fig. 15
figure 15

The step cycle for a G1 ‘pacing walk’ for three different duty factors (DF). For DF = 0.5 and 0.6 the right manus is in transition and must have its position interpolated. For DF ≥ 0.625 all four limbs are in ground contact. As in Fig. 13, the dashed black line indicates α-phase ({phi }_{alpha }) the relative phase between MR and PL

Table 3 shows the values of α-phase and the minimum duty factor DFMin required to achieve four-limb support, for each of the eight gaits considered. This generalization extends to diagonal sequence gaits G5G7, wherein the three-limb support is provided by PL, PR, and MR and ML in transition (and the α-phase is then measured between ML and the pes PL).

Table 3 The α-phase ϕα and the minimum duty factor DFMin required to achieve four-limb support for each given gait

Coupling length CL can now be measured as the distance between the pes and manus couplers for any walking gait. In each case, the two pes locators are in support and one or the other manus locator may require linear interpolation, by:

$$ gamma_{manus} = left( {frac{{x_{1} + x_{2, prev} + alpha left( {x_{2, next } – x_{2, prev} } right)}}{2}, frac{{y_{1} + y_{2, prev} + alpha left( {y_{2, next } – y_{2, prev} } right)}}{2} } right) $$


where (psi 1) is the locator in support, ({(x}_{psi 1}, {y}_{psi 1})) is the position of that locator, and (psi 2) refers to the other locator which is in transition, the position of which is being interpolated.

This formula reduces to the original GA formula when α is 0 or 1, i.e., when the duty factor is sufficiently long to permit all four limbs to be in simultaneous ground contact for that given gait. Coupling length is therefore a discretely sampled, periodic value within a trackway that can be measured twice per gait cycle (each time a pes support cycle ends).

Given that the S1 trackway consists of 10 complete step cycles, each trial resulted in 20 coupling length measurements (twice per cycle for 10 cycles). Median CL was then computed for each trial for the 20 coupling lengths. Figure 16 plots the results sorted by median CL length. No coupling lengths were excluded, resulting in some improbably short-coupled (CL = 0.77 m) and some improbably long-coupled (CL = 4.7 m) trackmakers for the size of the tracks in this trackway. Each computed value of CL also has an associated uncertainty which reflects error propagation from the uncertainty in the track locations. Despite the excellent preservation of this trackway, the small uncertainties in track position result in a large degree of overlap between median coupling length values in neighboring solutions. While trials with shorter CL have a greater relative uncertainty than those with longer values of CL, the median values by themselves provide very little information. Not finding any distinguishing behavior of median CL across any of the 40 distinct gaits, we turned to examining how CL varied locally along the length of the trackway for each trial condition, hence we developed a quality-of-fit measure.

The next step is to create a metric by which we can evaluate the variation in coupling length along the trackway for each trial, in order to find a ‘best’ solution.

Fig. 16
figure 16

Median Coupling lengths by trial. Note the color code indicating gait (G0G7)

Appendix 2

Evaluating quality of fit

We considered how CL varies along the trackway for each condition. We began by focusing on the G2 walk gait (specifically G2-1-60 and G2-2-60) and plotting their coupling length values versus cycle along the length of the trackway. Figure 17 shows the results for G2-1-60 and G2-2-60, respectively, where the configuration G2-2 has one additional stride length separation between the pes pair and the manus pair. Each of the following plots has a one-meter range on the ordinate centered on the median coupling length value for that trial. The observed variations in coupling length values along each trial arise from two potential sources. The first is noise, which is inherent in any biological system. The second is a characteristic perturbation to CL in response to the track pattern perturbation. It’s this second type of characteristic variation that we will exploit, provided it can be distinguished from the noise.

With reference to Fig. 17, there are clear differences between the responses of the G2-1-60 and the longer-coupled G2-2-60 to this trackway. Since we hold constant both the gait and the duty factor, the difference between the two trials is essentially the ‘body length’ (by adding one stride length to the choice of manus tracks to associate with the given pes tracks). We seek to establish parameters derived from the coupling lengths in each trial that can be used to quantify the characteristic variations between trials. These parameters will form a “feature space” (a set of possible combinations of parameter values) that can be used to determine the fitness of each simulation trial as a solution for the trackway. However, as with many complex systems, the characteristic variations are not perfectly orthogonal to the noise. A feature space cannot be created that represents only the characteristic variations within trials. Some amount of noise will be included in each parameter and will have to be accounted for in the analysis using uncertainty and error propagation techniques.

Fig. 17
figure 17

Two trials for which coupling length is computed along the BEB-500-S1 trackway for a G2 walk (LP = 0.25, DF = 0.6). Upper: shows the G2-1-60 gait with track separation of 1 (refer to Fig. 6 for the track configuration corresponding to G2-1). Lower: the same trackway and gait but for G2-2-60 (i.e., a track separation of 2)

What interpretation should be given to the variations in measured coupling length as it is applied along a trackway? This computation results from the application of fixed specific values for the free parameters to the positions of actual track locations. The irregularity in computed CL therefore results from the irregularity in the trackway, and even a trackway as subtly irregular as that in Fig. 7 could result in substantial variation in CL. This variation is difficult to explain in biological terms, as that would suggest the axial length of the trackmaker varied quite substantially during its passage. If we reject that conclusion, then it comes down to either having applied the ‘wrong’ choice of free parameters or having incorrectly assumed that the trackmaker walked with a constant gait as it created that trackway. To resolve this, we seek a metric on which to judge the suitability of a given gait as a model for the trackmaker’s movements along an actual trackway. We do not seek to solve for the ‘right’ choice of gait parameters; rather we seek to rule out ‘wrong’ choices, based on a fitness function which we will develop.

Persistent residuals

Consider a fictional trackway where two simulation trials, a and b, produce coupling length plots with the same median value but rather different profiles (Fig. 18). We wish to quantify the salience of these differences.

Fig. 18
figure 18

In a, coupling length increases from a median value at two separate points along the trackway, while in b, this shift occurs once, but for a longer duration. Note that the area under the two curves is identical

A residual analysis is a common approach to distinguish these cases, wherein the difference between each sample and the expected value is calculated. Each difference is squared to make all values positive and to weigh relatively larger differences more strongly than smaller ones. If we specify the consistent baseline value as the expectation value, we can calculate this residual sum of the squares (RSS) as:

$$ RSS = mathop sum limits_{i = 1}^{N} left( {x_{i} – overline{x}} right)^{2} $$


where (overline{x}) is the baseline value. For the two plots in Fig. 18, however, the RSS results are identical, hence a standard residual analysis would conclude that the two trials are equally efficacious solutions for that trackway.

We can improve the residual analysis by considering the physical restrictions on the trackmaker. In trial a, the two deviations are isolated and short in duration (consisting of one sample above the baseline in each case) while in trial b there is one deviation for twice the duration. Given we are computing coupling length along a trackway, one single step of either a pes or manus that deviates substantially from the median step length could produce an isolated deviation in CL from the baseline. A slightly longer (or shorter) than usual pes step at the same moment as a slightly shorter (or longer) manus step could also combine to create a short-isolated deviation in the CL computation. But the story is different for trial b in Fig. 18, where CL deviates from the baseline for two successive measurements. For the trackmaker to sustain a deviation from baseline for two or more successive measurements (one step cycle or more) requires some change in gait parameters, if one is to believe the trackmaker cannot substantially change its axial length during that protracted interval. Taking this into account requires altering the way the residuals are calculated. Our concern is not that a residual exists, but that a residual persists for a longer interval than would reflect an isolated deviation of the trackmaker. We can define this notion of persistence by replacing the square of each residual in the RSS calculation (Eq. 8) with the product of nearest neighbor residuals, i.e.,

$$ PRSS = mathop sum limits_{i = 1}^{N – 1} left| {x_{i} – overline{x}} right|left| {x_{i + 1} – overline{x}} right| $$


Applied to our examples in Fig. 18, trial a has a value of PRSS = 0 while trial b has a PRSS = 0.023. The persistent residual effectively removes noise variations caused by single-sample deviations. With it we can conclude that trial a is a more efficacious solution for the trackway. In the example above we specified a reasonable, but somewhat arbitrary, value for the expectation value (underline{x}). For more general simulation trials, PRSS is better expressed as:

$$ PRSS = mathop sum limits_{i = 1}^{N – 1} left| {CL_{i} – overline{CL} } right|left| {CL_{i + 1} – overline{CL} } right| $$


where (overline{CL}) is the median length across the given simulation trial. Applying this PRSS formulation (Eq. 10) to the BEB-500 S1 trackway yields the following plot of the persistent residual sum of squares for coupling length (Fig. 19).

Fig. 19
figure 19

PRSS values for the 40 trials of CL computations along BEB-500-S1

It is apparent that there is a scaling issue in these results. We did not observe this in our fictional example because the two trials shared the same expectation (baseline) values. Median coupling values generally differ in actual simulation trials, however. Normalizing the residuals by their median coupling values eliminates the problem:

$$ PRSS = mathop sum limits_{i = 1}^{N – 1} left| {frac{{CL_{i} }}{{overline{CL} }} – 1} right|left| {frac{{CL_{i + 1} }}{{overline{CL} }} – 1} right| $$


There is another subtle issue caused by boundary conditions at the beginning and end of the trackway. The number of samples, N, in a trial depends on the gait for a given trial. For example, manus locators in a long-coupled trial (such as G2-2) are further ahead in a trackway than in a short-coupled trial (G2-1) at the same moment in a simulation, and therefore will complete fewer samples. This relatively penalizes those simulation trials with more samples, which is resolved by simply dividing the PRSS by the number of residuals in the summation:

$$ PRSS = frac{1}{N – 1 }mathop sum limits_{i = 1}^{N – 1} left| {frac{{CL_{i} }}{{overline{CL} }} – 1} right|left| {frac{{CL_{i + 1} }}{{overline{CL} }} – 1} right| $$


PRSS can now be calculated by Eq. 12 for every simulation trial on a trackway and used to compare one trial to another. For BEB-500 S1 this yields the graph in Fig. 20.

Fig. 20
figure 20

PRSS values normalized to account for expected increases in CL across the different trial conditions

We can now use PRSS values to compare the solution efficacy of the two trials a and b in terms of standard deviations:

$$ Delta PRSS = frac{{|PRSS_{a} – PRSS_{b} |}}{{sqrt {sigma_{a}^{2} } + sigma_{b}^{2} }}, $$


where ({sigma }_{(a, b)}) are the uncertainties for the a and b trials.

While pairwise comparisons are useful, a global comparison across all trials would be preferable and allow the solution efficacy of every trial to be ranked at once. We therefore generalize the previous pairwise comparison formula (Eq. 13) to:

$$ Delta PRSS = frac{{|PRSS_{i} – PRSS_{min } |}}{{sqrt {sigma_{i}^{2} + sigma_{min }^{2} } }} $$


Comparisons are made for each trial relative to the smallest persistent residual sum-of-squares value, PRSSmin, from the collection of solutions. A value of 0 indicates that there is no observed difference between a particular trial and the ‘best’ solution while values greater than 0 indicate a less efficacious solution than PRSSmin. This makes ΔPRSS suitable as a fitness parameter for persistent residuals (Fig. 21).

Fig. 21
figure 21


Consider another fictional trackway where this time two simulation trials a and b produce the two coupling length CL plots in Fig. 22. Each trial has two persistent deviation regions, one large and one small. While the two persistent deviation regions are both positive in trial a, in trial b they swing from positive to negative.

Fig. 22
figure 22

In a, coupling length deviates from a median value at two separate points along the trackway with the same sign, while in b the deviations have opposite sign. Note that the area under the two curves is identical

The value of PRSS for these two trials is 0.025 and thus fails to distinguish the differences in the sign of the residuals, which represents a greater range of CL values for trial b than trial a. Another parameter is therefore needed to express the overall range of coupling length values within a trial. We define this parameter, which we will call swing, as the difference between the minimum and maximum coupling length values:

$$ swing = frac{{left| {CL_{max } – CL_{min } } right|}}{{CL_{median} }} $$


It is normalized by the median coupling length value in order to avoid the scaling issue that arose earlier in the PRSS formulation when comparing trials with different median coupling lengths. For trials a and b in Fig. 22 the calculated swing values are 0.15 and 0.2, respectively. The smaller swing for trial a suggests that a is a better solution than trial b (i.e., the trackmaker using gait a would have had a lesser transient deviation from steady state than if using gait b).

But the swing formulation is not yet complete. We previously established the importance of persistence in the formulation of the PRSS parameter. Persistence pertains to swing as well. While deviation for the course of one sample can be attributed to noise, deviation for more samples likely indicates more than noise as the underlying cause. To illustrate this, suppose another fictional trackway produced the following two coupling length plots for two trials a and b (Fig. 23).

Fig. 23
figure 23

The results of the two trials are identical but for the third sample, which is above baseline in b

In this example, the only difference between trials a and b is the third coupling length sample, which, in Fig. 23a has zero deviation from the baseline value, while in Fig. 23b that sample has a non-zero deviation from the baseline value. Using the formulation for swing outlined above, these two example trials have the same swing value despite the greater persistence of a deviation in trial b than in trial a. This is efficiently resolved by introducing a forward-moving average, which is generally defined in the form

$$ overline{x}_{i} = frac{{x_{i} + x_{i + 1} }}{2} $$


Next, we need to incorporate the uncertainties of the samples into this forward-moving average, resulting in a weighted version:

$$ overline{CL}_{i} = frac{{frac{{CL_{i} }}{{sigma_{i}^{2} }} + frac{{CL_{i + 1} }}{{sigma_{i + 1}^{2} }}}}{{sigma_{i } + sigma_{i + 1} }} $$


where σi and σi+1 are the uncertainties in coupling length values at i and i + 1. The weighted averaged coupling length values can then be substituted into the previous swing Eq. (15) as

$$ swing = frac{{left| {overline{CL}_{max } – overline{CL}_{min } } right|}}{{CL_{median} }} $$


This updated swing formulation now properly distinguishes between the two trials from the previous example with a higher swing value for trial b than for trial a. Applying this swing formulation (Eq. 18) to BEB-500 S1 yields the plot in Fig. 24.

Fig. 24
figure 24

Normalized swing by trial

In the same fashion as the PRSS parameter, the swing parameter needs to be converted into a fitness parameter for comparison across trials. We use the same approach for swing as we used for PRSS. That is, the lowest swing value among the trials is used as the highest fitness value, and all other trials are compared to it using the deviation significance calculation:

$$ Delta_{swing} = frac{{left| {swing_{i} – swing_{min } } right|}}{{sqrt {sigma_{i}^{2} + sigma_{min }^{2} } }} $$


The computed swing deviation is shown in Fig. 25.

Fig. 25
figure 25

The results of the swing deviation

Solution fitness

Now that we have our two fitness parameters ΔPRSS and Δswing (Eqs. 14 and 19) we combine them by plotting each trial on a scatter plot with the parameters as the values for each axis. To weigh the fitness parameters equally, the values for all of the trials are rescaled to a maximum of 1 (Fig. 26).

Fig. 26
figure 26

The results of the swing deviation

The closer a particular trial is to the origin of this fitness parameter space, the more efficacious a solution the trial is for the trackway. Computing the Euclidean distance for each point creates the final fitness parameter ranking for the trials

$$ fitness = sqrt {Delta_{swing}^{2} + Delta_{PRSS}^{2} } $$


And finally, to more clearly present the quality-of-fitness, the fitness for each trial is subtracted from the poorest fitness to yield the plot in Fig. 11.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit


This article is autogenerated using RSS feeds and has not been created or edited by OA JF.

Click here for Source link (