Quantifying changes in individual-specific template-based representations of center-of-mass dynamics during walking with ankle exoskeletons using Hybrid-SINDy

Ankle exoskeletons alter whole-body walking mechanics, energetics, and stability by altering center-of-mass (CoM) motion. Controlling the dynamics governing CoM motion is, therefore, critical for maintaining efficient and stable gait. However, how CoM dynamics change with ankle exoskeletons is unknown, and how to optimally model individual-specific CoM dynamics, especially in individuals with neurological injuries, remains a challenge. Here, we evaluated individual-specific changes in CoM dynamics in unimpaired adults and one individual with post-stroke hemiparesis while walking in shoes-only and with zero-stiffness and high-stiffness passive ankle exoskeletons. To identify optimal sets of physically interpretable mechanisms describing CoM dynamics, termed template signatures, we leveraged hybrid sparse identification of nonlinear dynamics (Hybrid-SINDy), an equation-free data-driven method for inferring sparse hybrid dynamics from a library of candidate functional forms. In unimpaired adults, Hybrid-SINDy automatically identified spring-loaded inverted pendulum-like template signatures, which did not change with exoskeletons (p > 0.16), except for small changes in leg resting length (p < 0.001). Conversely, post-stroke paretic-leg rotary stiffness mechanisms increased by 37–50% with zero-stiffness exoskeletons. While unimpaired CoM dynamics appear robust to passive ankle exoskeletons, how neurological injuries alter exoskeleton impacts on CoM dynamics merits further investigation. Our findings support Hybrid-SINDy’s potential to discover mechanisms describing individual-specific CoM dynamics with assistive devices.


SLIP dynamics
SLIP parameters are described in Figure S1.Variables  and  denote sagittal and frontal-plane leg angles with respect to vertical, respectively.The center-of-mass (COM) has mass M. Left: Before normalization, model parameters are denoted as: Leg radial stiffness (kL) and damping (cL), leg resting length (L0), rotary stiffness in the sagittal (ks) and frontal (kf) planes, and rotary damping in the sagittal (cs) and frontal (cf) planes.Non-normalized parameters were estimated for each participant.Right: After normalization, parameters are denoted: Leg radial stiffness (  ) and damping (  ), leg resting length ( ̃0), rotary stiffness in the sagittal (  ) and frontal (  ) planes, and rotary damping in the sagittal (  ) and frontal (  ) planes [1][2][3][4] .Normalization was applied to compare participants after model fitting.
We solve the 3D SLIP dynamics with respect to linear accelerations as follows.Let SLIP model mass be denoted by M, and 3D positions of the center-of-mass (COM) relative to the feet be denoted x, y, & z.We solve for linear accelerations of the COM for each foot on the ground as: ̈=  ̈ −  ̇2 + 2 ̇̇  +  ̈ ̈= − ̈ −  ̇2 − 2 ̇̇  +  ̈ ̈=  ̈ −  ̇2 + 2 ̇̇  +  ̈ 1 where single dots denote first derivatives and double dots denote second derivatives.Other variables are defined as in Figure S1.Note that we assume that anterior-posterior and vertical dynamics are decoupled from frontal-plane dynamics.This simplification is supported by prior studies of COM dynamics during gait 5,6 .We solve radial ( ̈) and rotary accelerations ( ̈,  ̈) as a function of spring and damper forcing terms and gravity, g:  where variables correspond to leg radial stiffness (kL) and damping (cL), leg resting length (L0), rotary stiffness in the sagittal (ks) and frontal (kf) planes, and rotary damping in the sagittal (cs) and frontal (cf) planes.Note that we assume that the resting angle of the rotary stiffness mechanisms is zero 2 .Substituting forcing terms into equation (1) yields: Omitting Coriolis terms and adopting the notation  = [, , ]  , with vectors and matrices in boldface, we can write the 3D bipedal SLIP dynamics as: where the summation corresponds to the right and left legs.The left-most brackets contain mechanisms that impart forces radially along the leg:   is the leg stiffness,  is the instantaneous leg length,  0 is the leg resting length,   is the leg damping,  ̇ is the instantaneous leg velocity.The middle bracket contains mechanisms that impart forces transverse to the leg axis in the sagittal plane:   and   are the sagittalplane rotary stiffness and damping, respectively.  denotes the sagittal-plane leg projection.Analogously in the right-most brackets,   and   represent the frontal-plane rotary stiffness and damping, respectively, and   denotes the frontal-plane leg projection.
Note that we add gravitational acceleration vector, g, to both sides of equation ( 4) to remove it from the right-hand side of the formulation.We can rewrite the dynamics as a linear map from the nonlinear transformations of the SLIP states to linear accelerations: ̈ = (, ̇)

Akaike Information Criterion (AIC)
In the following description of the AIC, we consider a nonlinear system with dynamics described by equation 5, where the estimated system dynamics are described by  ̂(), as the formulation below is not specific to models identified by Hybrid-SINDy.
The AIC is widely used to compare candidate models according to their number of free parameters (i.e., model complexity), , and log-likelihoods (i.e., representativeness of the observations), (, ).The AIC has the form:  = 2 − 2 ln((, )).6 Mangan and colleagues (2019) assumed that model errors are independently, identically, and normally distributed, enabling the AIC to be written in terms of the number of samples, , and the sum of squared residuals 3 : where in the outer summation, P denotes the number of samples and in the inner summation, I denotes the number of output states.A correction may be applied to the AIC formulation to account for bias due to finite sample sizes (AICc): The AICc score provides a level of support for each candidate model, with lower AICc scores indicating stronger support for a model.To compare competing models, a relative AICc score (Δ) may be calculated as the difference between the AICc score of the minimum-AICc model and that of competing models: where the subscript j, denotes the j th model.Burnham and

Comparison to alternative modeling frameworks
Here we present a brief analysis comparing Hybrid-SINDy to alternative modeling dimensionality reduction frameworks: Stepwise regression 8 and principal components analysis (PCA).

Stepwise regression
Stepwise regression is similar to Hybrid-SINDy in that it identifies sparse system models, often by eliminating or adding terms until terms no longer improve accuracy or are statistically significant.However, unlike Hybrid-SINDy, stepwise regression only selects a single model, even if multiple models are equally plausible representations of the system 9 .Further, stepwise regression is known to select different models depending on the data, especially when predictor variables are colinear 9,10 .However, it is unclear if stepwise regression could produce similar models as Hybrid-SINDy.
Approach: To test this, we replaced the SINDy algorithm with a similar approach using stepwise regression.Our stepwise regression approach is similar to that of Shamaei and Colleagues (2013), who estimated ankle quasi-stiffness from mechanics-based functional forms 8 .
First, to provide an example of how stepwise regression produces different and less-sparse results than SINDy, we compared stepwise regression, as described in 8 , to SINDy for the synthetic SLIP without noise and with moderate noise ( = 1 mm).We report results for a single cluster in single-limb support.
Next, we applied a similar framework to human walking data and produced template signatures for walking in all three exoskeleton conditions, following the methodology described in the main manuscript.
Results: For synthetic SLIP data, Even in the noise-free condition, stepwise regression selected artificially high-dimensional dynamics (Table S1.1).The stepwise regression model had rank = 13, compared to the SINDy-based model, which identified 4 models ranging from 5-to 0-dimensional.Note that the SLIP model has at most 6 non-zero terms in each gait phase (see Supplemental -S2, Table S2.1 for synthetic SLIP parameters).With noise (Table S1.2), SINDy identified higher-dimensional solutions but also identified additional low-dimensional solutions that would likely be more plausible according to the AIC.In human data, stepwise regression also identified rank = 13 dynamics, while SINDy identified lowerrank dynamics (Table S1.3).Consistent with single-cluster results in the synthetic SLIP, stepwise regression-based template signatures were not sparse (Figure S1.2; compare to Figure 3 in the main manuscript).Figure S1.3 shows that the post-stroke template signatures are also not sparse, making it unclear which mechanisms are most important for template representations of CoM dynamics (compare to Figure 6 in the main manuscript).

Interpretation:
These results showcase that stepwise regression tends to select less parsimonious models than Hybrid-SINDy, hindering interpretation of template signatures.While regression hyperparameters may be tuned to reduce this dimensionality, Hybrid-SINDy's use of clustering and model selection within clusters likely enables more sparse template representations of CoM dynamics to be identified.Therefore, Hybrid-SINDy was appropriate for our work.

Principal components analysis (PCA)
Principal component analysis (PCA) could identify interpretable basis vectors describing CoM dynamics.
For example, we could apply (e.g.) PCA or NMF to the nonlinear template states and use information criteria to select an optimal number of modes in the model 11 .
However, PCA would not satisfy our goal of interpreting individual differences in template-based representations of CoM dynamics.PCA-based models would not be sparse in the space of template parameters: If the optimal (e.g., according to the AIC) PCA or NMF-based model was unidimensional, the first PC would contain multiple template states.The higher-dimensional modes of these models would likely also contain multiple non-zero template states.Therefore, PCA would not provide the same mechanically intuitive interpretability as the sparse regression of template states used by Hybrid-SINDy.
Approach: To exemplify the limitation of PCA (with potential generalization to other dimensionality reduction approaches), we show that a sparse PCA-based model is still dense in the template mechanism space.Using data from our synthetic SLIP, we applied PCA to the template states, then used SINDy to identify the linear map,   , from PCs to CoM accelerations.We then projected the map back into the original template basis (  in Table S1.4).
Results: Table S1. 4 shows an example of a 1-D linear map identified by SINDy.While the linear map is sparse in the PC space, it becomes dense when projected back into the original template mechanism space (Table S1.4).Validating Hybrid-SINDy using a synthetic SLIP walking model We evaluated Hybrid-SINDy's ability to accurately identify and select walking dynamics, using synthetic data generated from simulations of a bipedal SLIP walker with known dynamics (Table S2.5 and Figure S2.4) 1 .The bipedal SLIP had dynamics matching those used in human data (Supplemental -S1), with the addition of foot mass to enable simulation of a full gait cycle.The synthetic SLIP had asymmetric leg stiffness and damping, and masses representing the pelvis (M = 56 kg) and feet (Mf = 7 kg) (Table S2.5; Figure S2.4).We omitted rotary mechanisms from the SLIP dynamics but used the full function library described in Supplemental -S1, equation 4. To replicate the human walking datasets, we simulated 120 gait cycles each from randomly perturbed initial conditions with an average initial velocity of 1.20 m/s.
The SLIP walked at approximately 1 stride/s.We conducted two analyses using the simulations of a synthetic SLIP model: First, we evaluated the effects of measurement noise on Hybrid-SINDy's ability to estimate synthetic SLIP parameters in the presence of sensor noise.Second, we evaluated Hybrid-SINDy's ability to estimate synthetic SLIP parameters when a subset of the terms constituting the true system dynamics is missing from the function library.

Effects of measurement noise on algorithm performance
Approach: As SINDy's performance is sensitive to sensor noise, we identified template signatures using simulation measurements with Gaussian noise ranging logarithmically from 10 nm to 10 cm added to the position measurements 2 .After adding noise, we differentiated each gait cycle separately before concatenating the strides to avoid errors in the estimated time derivatives due to changes in the initial conditions of each simulation.We quantified the Hybrid-SINDy's ability to reconstruct model dynamics using the number of correctly selected mechanisms and the accuracy of estimated coefficients.We also evaluated the ability of each model to predict the ground truth (noise-free) CoM accelerations.
Results: In the noise-free condition, the Hybrid-SINDy algorithm correctly identified the salient features of the simulated SLIP's template signature during single-limb support (Figure S2.5).Specifically, Hybrid-SINDy identified non-zero leg stiffness (  ), resting length ( ̃), and damping (  ) terms onlyconsistent with the SLIP's dynamics.In the noise-free condition, 86% of terms were correctly selected and only one term, leg damping in first double-limb support, was incorrectly selected (Figure S2  ).B) Effects of noise on model coefficient percent error (left) and the percent of correctly selected mechanisms (right).All horizontal axes represent noise magnitudes on a log scale.For coefficient error, percent errors were computed for stiffness (solid line), leg length (dashed line), and damping (dotted line) separately.Green lines denote the  = 1 mm noise level.
When measurement noise with a standard deviation of 1 mm-similar to that of our motion capture system-was added to the system measurements, Hybrid-SINDy correctly selected 83% of mechanisms, with incorrect mechanisms selection occurring primarily in double-limb support.Template signature coefficients were estimated with an average error of less than 29% in single-limb support and swing (Figure S2.5B).In single-limb support, Hybrid-SINDy was still able to estimate large coefficients (i.e., leg stiffness and resting length) with less than 1% error.Even at this noise level, the estimated template signature could still predict the noise-free data with similar accuracy to the noise-free data (r 2 = 0.83).The accuracy of coefficient estimates and model predictions rapidly declined at noise levels with standard deviations larger than 1 mm.
Interpretation: Hybrid-SINDy's ability to identify template signatures from forward simulations of a synthetic SLIP model with known dynamics at measurement noise levels comparable to marker-based motion capture suggests that individual differences in template signatures of human gait data likely reflect some differences in the underlying CoM dynamics.Hybrid-SINDy performed as expected in the noisefree condition, accurately identifying single-limb support and swing dynamics, except for small damping terms.Terms with small contributions to system behavior, such as damping in this study, may be omitted by the AICc in favor of parsimony 3,4 .Therefore, Hybrid-SINDy will not always select the true system dynamics but will select terms critical to describing the salient aspects of CoM dynamics.In our human walking data, Hybrid-SINDy may, therefore, omit mechanisms describing subtle idiosyncrasies in CoM dynamics that could be important for gait stability or efficiency.Future experiments using different constraints or perturbations may identify additional mechanisms needed to describe CoM dynamics.
The reduced accuracy of synthetic SLIP template signatures in double-limb support and with increasing noise magnitude is consistent with the results of Hybrid-SINDy applied to a spring-mass hopper 5 .However, accurate estimates of single-limb support CoM dynamics and longer double-limb support times in human walking data support our use of 800 samples per cluster for human template signatures.Longer time-series datasets would help to ensure that clusters are small enough to span only a single gait phase but large enough to remain robust to sensor noise.
Interpretation: Small differences in kinematics in humans compared to the synthetic SLIP simulations suggest that large differences in template signature coefficients are unlikely due entirely to differences in kinematics.However, at least 5-10% differences in template signature coefficients can be attributed to changes in kinematics rather than CoM dynamics.Larger differences in coefficient estimates may be needed for human data, depending on how well the model captures the underlying dynamics.Approach: To investigate the effect of covariation between state variables for the human data, we correlated SLIP states within each gait phase separately (Figure S2.7; average coefficient of determination across participants).
Results: Leg stiffness and resting length were most strongly correlated, as expected (r 2 = 0.99).These states were both selected by Hybrid-SINDy.Most correlations were less than r 2 = 0.70 (Figure S2.8 highlights strong correlations: r 2 > 0.70).While frontal-plane rotary stiffness (  ) and damping (  ) in double-limb support were strongly correlated, only frontal-plane rotary stiffness was selected for a subset of participants (see Figure 3, main manuscript).Interpretation: Hybrid-SINDy's performance appears robust to moderate correlations in the state variables.Covariation in the state variables did not result in unexpected variables being included in the template signatures.However, different combinations of parameters could have similar statistical plausibility.In this case, our use of the Akaike Information Criterion and multi-model inference acknowledges the existence of multiple plausible template signatures 3,4 .While we selected mechanisms with distinct relationships to CoM accelerations, studies using larger or non-mechanistic function libraries may be more strongly impacted by covariation in the state variables.Such studies should evaluate covariation among state variables before model fitting.

Figure S1 :
Figure S1: A 2D schematic of the 3D SLIP with model parameters (left) and normalized parameters (right).Variables  and  denote sagittal and frontal-plane leg angles with respect to vertical, respectively.The center-of-mass (COM) has mass M. Left: Before normalization, model parameters are denoted as: Leg radial stiffness (kL) and damping (cL), leg resting length (L0), rotary stiffness in the sagittal (ks) and frontal (kf) planes, and rotary damping in the sagittal (cs) and frontal (cf) planes.Non-normalized parameters were estimated for each participant.Right: After normalization, parameters are denoted: Leg radial stiffness (  ) and damping (  ), leg resting length ( ̃0), rotary stiffness in the sagittal (  ) and frontal (  ) planes, and rotary damping in the sagittal (  ) and frontal (  ) planes[1][2][3][4] .Normalization was applied to compare participants after model fitting.

Figure S1. 2 :
Figure S1.2:The percentage of unimpaired legs whose stepwise regression-based template signatures contained each mechanism in each gait phase.For 24 legs, the percentage of legs for which each template signature mechanism was selected by the Hybrid-SINDy algorithm.Colors denote each gait phase.Mechanisms selected in a larger percentage of legs suggest common representations of CoM dynamics, while less frequently selected mechanisms reflect individual-specific template features describing CoM dynamics.

Figure S1. 3 :
Figure S1.3:Inter-leg differences in template signatures in one stroke survivor identified using stepwise regression.Non-paretic (colored bars) and paretic (white bars) template signatures for one individual with post-stroke hemiparesis.Each plot represents a different gait phase.

Figure S2. 4 :
Figure S2.4:Synthetic SLIP model depiction and simulated CoM states.A) Two-dimensional depictions of the simulated SLIP.The simulated SLIP used leg springs and dampers, as well as foot masses to simulate full strides.B) Time-series measurements of simulated CoM position, velocity, and acceleration for the simulated SLIP (left; average ± 1SD).Template states (right) show trajectories of the template variables used in the function library to estimate template signature coefficients (average ± 1SD).
.5B).Both double-limb support phases were very short in the simulated SLIP: 1.3% and 4.6% of samples in first and second double-limb supportsmaller than the 7.4% of the training data in each cluster.Consequently, all clusters with centroids in the double-limb support phases contained samples spanning double-and single-limb support.Hybrid-SINDy incorrectly omitted only damping terms with small coefficients (FigureS2.5).In single-limb support and swing, coefficients were within 2% of the ground truth values.While prediction accuracy (averaged across output dimensions) was lower than expected (r 2 = 0.88), errors occurred primarily at the transitions between gait phases and in the double-limb support phases (FigureS2.5A,left: deviations in the predictions (colors) relative to ground truth simulation results (gray); transitions occur where colors change).

Figure S2. 5 :
Figure S2.5:Parameter estimates, prediction accuracy of simulated CoM kinematics, and the effects of noise on model identification accuracy.A) Model coefficients (top) can be compared between walking conditions.Bar plots show ground truth synthetic SLIP parameters (gray) and the corresponding estimated template signature coefficients (colors).The range of percent errors in coefficients is shown above each plot.The average (±1SD) ground truth CoM accelerations of the synthetic SLIP (bottom; gray) and model predictions in each gait phase (colors).Model accuracies are quantified using the

Figure S2. 7 :
Figure S2.7:Correlations between SLIP states in human walking data for each gait phase.Colormaps of R-squared (R 2 ) values denoting correlations between each pair of state variables in each gait phase (one panel for each phase, denoted by titles above each panel).Black boxes denote R 2 = 0 and white boxes denote R 2 = 1.The lower triangle is omitted for clarity (the correlation matrix is symmetric).Along the axes, the second subscript on each variable denotes the right (R) or left (L) leg.Only axial stiffness (  ) and leg resting length ( 0 ) states were, as expected, very strongly correlated.

Figure S2. 8 :
Figure S2.8:Strong correlations between SLIP states in human walking data for each gait phase.Colormaps of R-squared (R 2 ) values highlight strong correlations (R 2 > 0.70) between each pair of state variables in each gait phase (one panel for each phase, denoted by titles above each panel).Black boxes denote R 2 < 0 and white boxes denote R 2 > 0.70.The lower triangle is omitted for clarity (the correlation matrix is symmetric).Along the axes, the second subscript on each variable denotes the right (R) or left (L) leg.While frontal-plane rotary stiffness (  ) and damping (  ) in double-limb support were strongly correlated, only frontal-plane rotary stiffness was selected for some participants.However, Hybrid-SINDy consistently selected frontal-plane rotary stiffness, suggesting that stronger correlations are needed to select different mechanisms for different participants. ,

Table S1
.1: Noise-free synthetic SLIP template signatures in one cluster in single-limb support.

Table S1 .
3: Example human template signatures in one cluster in single-limb support.

Table S1
12: 1-dimensional PCA-based model of CoM dynamics.describingCoMdynamics using template models.This limitation would likely generalize to other techniques like non-negative matrix factorization12.
Covariation of template signature state variablesCovariation among the state variables corresponding to mechanisms in the Hybrid-SINDy function library can hinder the interpretation of template signatures.Specifically, covariation in state variables can alter coefficient estimates if both mechanisms corresponding to two co-varying states are selected, or if redundant mechanisms are selected at random for different participants.