In this study, a total of 18 subjects were recruited, including 10 patients after stroke (two women and eight men, age 48.3 ± 12.8 years) and eight healthy subjects (three women and five men, age 29 ± 4.94 years) in the control group. For the paretic and non-paretic groups, the inclusion criteria included the following: (1) chronic stroke survivors (more than 6 months after stroke); (2) the first stroke with unilateral hemiparesis lesions; (3) the ability to walk independently and continuously on a treadmill for at least 5 min; and (4) the ability to follow oral instructions and cooperate with experimental procedures. For the control group, the only inclusion criterion was no history of neuromuscular diseases. All subjects were ambulating independently and were not currently receiving physical therapy. This study was conducted under the approval of the Medical Ethics Committee at the Industrial Injury Rehabilitation Hospital of Guangdong. All subjects had signed a consent form before starting the experiment.

Before the experiment, ten 12-mm reflective markers were affixed to the bilateral lower limbs of each subject and placed at the following anatomical reference positions from bottom to top: the second and third metatarsal space, the lateral malleolus, the midleg, the lateral knee, and the mid-thigh (30). Then, the participants were asked to walk on a treadmill (BH, G6425-F3, Spain) for 5 min wearing specific experimental shoes to familiarize themselves with the experiment and to find the most comfortable speed by adjusting the treadmill speed. During the experiment, subjects were asked to walk on a treadmill for 3 min at a comfortable speed, and this activity was repeated 3 times. There was a 5-min rest period between each walking experiment (31). When the subject was walking on the treadmill, their three-dimensional coordinate trajectory was detected by the 6-cameras motion capture system with a 100 Hz sampling rate (OptiTrack, Natural Point, USA). The coordinate trajectories of all markers were recorded using Tracking Tools software (NaturalPoint, USA) and were processed using MATLAB (MathWorks, Natick, USA) to calculate the joint angle time series for further analysis.

The movement variability in knee and ankle joint angles was quantified by MFDFA. The three groups of the joint angle time series were analyzed as follows: (i) dominant limbs in healthy subjects, (ii) paretic limbs, and (iii) non-paretic limbs of the patients after stroke. The specific steps were as follows (21):

Construct the profile:For the time series x(i), i = 1⋯N, its average value xave=1N∑i=1Nx(i) was subtracted and the profile Y(i) was constructed by the cumulative sum of these differences:

Divide Y(i) into N_s = int(N/s) nonoverlapping segments: Here, s is the length of each equal-length segment. However, the length of the time series N was scarcely divisible by s, and the part that was not divisible might not be calculated. The same calculation process was repeated from the other end. Thus, 2N_s segments were generated in total. Then, the local trend fitting of l−order polynomial was performed on each of the generated segments, and the residual variation between the fluctuations between Y(i) and its l− order fit y_v(i), i = 1, ... , s was calculated as follows:

For the v = 1, ... , N_s segments, we calculated

whereas for the v = N_s+1, ..., 2N_s segments,

By calculating the mean of all segments, the q_th order fluctuation function was obtained:

The exponent variable q was any real number except zero. If q = 0, 1/q was infinity, and a different averaging method had to be used, such as logarithmic averaging. To discover the relationship between the generalized fluctuation functions F_q(s) and the timescale s for different q values, the above steps were repeated for multiple segments. DFA was a special case of MFDFA when q = 2.

Finally, the scaling behavior was determined by plotting F_q vs. s on a logarithmic scale and examining the existence of short and long scales. If the fluctuation function could be described as

it meant that the time series had scale characteristics and that the autocorrelation of the time series had no characteristic time scale. Then, we calculated the slope of the fitted straight line of ln F_q(s) and ln s for different values of q to obtain the scaling exponent h(q). In general, the exponent h(q) was closely related to q, and h(q) was considered to be the generalization of the Hurst exponent.

The generalization of Hurst exponent h(q) had the following relationship with the classical scaling exponent τ(q):

The feature of the mono-fractal series was the existence of only a single Hurst exponent H with a linear relationship between τ(q) and q. The multifractal time series would be characterized by multiple Hurst exponents, and τ(q) depended nonlinearly on q.

The singularity spectrum f(α) had the following relationship with h(q):

where f(α) quantified the dimension of the subset series corresponding to α, which is the singularity strength. A fractal was a repeating pattern that was self-similar across different scales; multifractality referred to patterns that are repeated in multiple ways. Thus, the multifractal spectrum could provide a lot of information about the behavior of fractals in the time series (29). To verify the power law behavior of the fluctuations, the time series were processed by the Fourier spectral analysis first (24, 25, 32). Then, the degree of multifractality was demonstrated by the width of the multifractal spectrum Δα = α_max − α_min. However, for a nonstationary and random walk-like structure time series, values of the degree of correlation h(2) might be above 1. A nonstationary random walk was obtained by integrating stationary noise. Thus, the result h(2) reduced by 1 to obtain the eigenvalues for quantifying the persistent and anti-persistent values of the nonstationary random walk time series (33). The value of h(2) also revealed the characteristics of long-scale correlation in the time series (21). The origin of multifractality could be determined by processing the time series random shuffling. There were two common sources of multifractals: (i) the extensive probability density function and (ii) the existence of many different long-scale correlation fluctuations in the time series (34). By random shuffling, the time series values would be randomly arranged, and the correlations in the original sequence would be destroyed. Therefore, if the multifractality was the result of the long-scale correlations, then the processed series exhibited a non-fractal scaling. When the origin of multifractality was a broad probability density, the reliance of h(q) was unchanged, which was unaffected during the reorganization process. If these two multifractals simultaneously appeared, the processed series would exhibit a decrease in the degree of multifractality.

Statistical differences in the values of Δα and h(2) for each scaling region between subjects with paretic limbs and healthy subjects, and between subjects with non-paretic limbs and healthy subjects were assessed using an independent t-test. Paired t-tests were applied between the results of Δα and h(2) in those with paretic and those with non-paretic limbs to investigate whether there were statistical differences between the two conditions. For all statistical analyses, the p-value of 0.05 significance level was set. The SPSS 24.0 software package was used for all statistical calculations (SPSS Inc., Chicago, USA).

The representative examples of the variation of the fluctuation function F with scale s are shown in Figure 1. First, the third-order polynomials (MFDFA1-3) were used in the detrending procedure. F decreased with the increase of the fitting order l because the increase of the order would result in a smaller residual. Then, values of h(q) for −10 ≤ q ≤ 10 were displayed to reveal the effect of different magnitudes of joint angles anomalies on a behavior scale. In Figure 2, the fluctuation function F_q of the knee angles' time series of a healthy subject and its behavior with the change of s were similar to F₂(s) in Figure 1. By performing linear fitting in two independent regions, the characteristics of multifractal scales on the short and long scales were characterized. The vertical red and blue dashed lines represent the boundaries of the short-scale and long-scale regions, respectively. Specifically, the short-scale region started after the DFA transient artifact with the smallest s values corresponded to 20 samples and this ended with 100 samples. The long-scale region started with 170 samples and usually ended with N/8, and the specific figure corresponded to 1,000 samples. In some joint angles, the phenomenon of crossover was less obvious in the scaling description, and the slopes h(q) were also estimated over the all-scale region. Furthermore, the all-scale region started the same way as the short-scale region, commencing with 20 samples and usually ending with 1,000 samples. This included the change in slope at the crossover, where the change in slope increased as q decreased. By performing linear fitting in the three-scale regions, the characteristics of multifractal scales on the short, long, and all scales were characterized. The vertical red and blue dashed lines are the boundaries of the short-scale and long-scale. The different scale ranges for each participant were the same in this study. The representative fluctuation behavior started with a linear segment with a slope >1, and the slope range was usually between 1.36 and 1.75.

The Fourier power spectrum P(v) of the knee angle time series is shown in Figure 3A (v is the frequency, axes are shown in logarithm scale). Statistical self-similarity is manifested as a power-law scaling of the frequency distribution. The power law coefficient (slope) is changed throughout the entire v range (Figure 3A), suggesting multifractality in the knee angle time series. Multifractals reflected the different scale behavior of time series anomalies. An increase in the moment q could realize the transition from small anomalies to large anomalies. Thus, the representative changes of h(q) with q in the three cases are presented in Figure 3B. The values of h(q) decreased continuously as the moment q increased in Figure 3B, which indicated a multifractal behavior. The change of h(q) with q and the non-linear reliance of τ(q) on q reflected the multifractal characteristics of human gait in these three cases (Figure 3C). The degree of multifractality was quantitatively determined by the distribution range of fractal dimension f(α), which was characterized by Holder exponent α. The typical inverted parabola form could be observed in the singularity spectrum f(α) vs. α (Figure 3D).

Results were reported separately for the degree of correlation and the degree of multifractality. The box plot of the multifractal analysis results showed that, when the order was l = 2, the time series of knee and ankle joint angles exhibited correlations and multifractality under the three fitting regions (Figures 4, 5). For the knee angles, the statistical results of h(2) indicated that there was persistence in all-scale and short-scale regions and anti-persistence in the long-scale region. Especially, h(2) corresponding to the knee angle was significantly higher in the paretic limb and the non-paretic limb of patients than in healthy subjects (p < 0.05). For the ankle angles, the statistical results of h(2) indicated that there was persistence in all-scale and anti-persistence in long-scale regions. In the short-scale region, healthy subjects and the non-paretic limb of subjects were persistence, while the paretic limb of the patients wase anti-persistence. The value of h(2) corresponding to the ankle angle was also significantly higher in the paretic limb and the non-paretic limb of patients with stroke than in the healthy subjects in all-scale and short-scale regions, but in the long-scale, only the non-paretic limbs of the patients were significantly higher than the healthy subjects (p < 0.05), although there was no significant difference in h(2) between the paretic and non-paretic limbs of patients after stroke. As presented in Figure 4, compared with the ankle joint angles, the knee joint angles were characterized by larger h(2) in the short-scale region, and the difference between the knee–ankle angles was greater in the patients. The values of Δα in knee angle were significantly higher in the paretic limb and the non-paretic limb of patients than in healthy subjects on all-scale, and the non-paretic limb of patients was higher than the healthy subjects on the short-scale region (p < 0.05). Although a significant difference occurred in the values of Δα in the ankle angles between the paretic limb and the non-paretic limb of patients with stroke compared with healthy subjects only in the long-scale region (p < 0.05), there was no significant difference on all-scale and short-scale regions.

Moreover, the corresponding random shuffle series were analyzed for all six data sets (knee and ankle angles of three groups) to determine the cause of the multifractality. As is evident in Tables 1, 2, the statistical differences of Δα and h(2) of each group of shuffled series were significantly smaller. The results showed that the shuffled data of the three groups all had scale exponents of about 0.5 and a multifractal strength of about 0.2. The results suggested that the multifractality in lower limb joint angles was predominantly the result of the long-range correlations.