## 2 Materials and Methods
### Gait Tube Analysis  
We employed gait tube analysis to quantify phase-dependent gait stability, a multi-dimensional framework that visualizes COM velocity trajectories in a 3D state space and computes variability via covariance-based ellipsoids [4]. COM velocity trajectories (Vx, Vy, Vz) were extracted and time-normalized to 100 points for each condition. Stride-to-stride variability was assessed by constructing a 3 \(\times\) 3 covariance matrix of Vx, Vy, and Vz at each gait phase, pooling data across participants. Clean strides were defined as those free from marker dropout, tracking artifacts, or timing anomalies. Strides were excluded if their duration or step length deviated more than \(\pm\)2 standard deviations from the participant's mean, or if visual inspection revealed discontinuities in the COM trajectory. Eigenvalues of these matrices were used to derive bivariate ellipsoids, whose volumes quantified local variability, forming a continuous \"gait tube\" that reflects the evolving structure of COM variability throughout the stride. Specifically, for each normalized gait phase (0-100%), we computed a 3 \(\times\) 3 covariance matrix of COM velocity (Vx, Vy, Vz) across all valid strides. The eigenvalues and eigenvectors of this matrix were used to define ellipsoids representing local variability. These ellipsoids were then projected onto the plane perpendicular to the mean trajectory using the Frenet-Serret frame to preserve the directional context of variability [26]. A five-point moving window was applied to smooth temporal fluctuations in ellipsoid volume, improving signal clarity while maintaining phase resolution. This process yielded a temporally continuous map of variability that forms the gait tube structure. This method, implemented via the MATLAB R2024b function (MathWorks, Natick, MA, USA) GaitTubePlot_Supplementary_MATLAB_Code.m (Supplementary Material: MATLAB Code), draws on principles from prior studies of multidimensional gait variability. A detailed mathematical formulation of the gait tube analysis method is provided in the Supplementary Material titled _\"Gait Tube Analysis Mathematical Model\"_, which accompanies this manuscript. Gait tubes were visualized from multiple perspectives: 3D (Vx, Vy, Vz), sagittal (Vx vs. Vy), transverse (Vx vs. Vz), and frontal (Vy vs. Vz), as shown in Figure 1. The mean ellipsoid volume was computed over the entire gait cycle for each condition to enable global comparisons. Direction-specific variability was quantified by calculating the stride-to-stride standard deviation of COM velocity in the ML, AP, and VT directions. A total variability index was defined as the sum of these directional variances. At each normalized time point in the gait cycle (0-100%), all valid strides across participants within a given condition were aligned in time. The COM velocity vectors (Vx, Vy, Vz) at each time point were assembled into a data matrix, and a 3 \(\times\) 3 covariance matrix was computed using standard covariance estimation:  
\[Cov=\frac{1}{n-1}\sum_{i=1}^{n}{(\mathbf{v}_{i}-\overline{\mathbf{v}})(\mathbf{v}_{i}- \overline{\mathbf{v}})^{T}} \tag{1}\]  
where \(\mathbf{v}_{i}\) represents the COM velocity vector for stride \(i\), and \(\overline{\mathbf{v}}\) is the mean vector. This matrix captures the stride-to-stride variability in all three directions. We chose to pool data across participants to generate a generalized group-level representation of the gait stability envelope for each condition. This approach enhances the visibility of consistent, condition-dependent variability patterns while controlling for individual idiosyncrasies. For participant-level comparisons (e.g., ellipsoid volume statistics and correlations), separate analyses were conducted to retain individual-specific features. Z-scores for ellipsoid volume and total variability were computed using the mean and standard deviation of the NoExo condition as a reference, providing a normalized measure of deviation from unassisted walking.