### Hill-type neuromuscular model of ankle joint

Below, we propose a modified sEMG-US imaging-driven HNM to directly build a relationship between the joint net PF moment and sEMG-US imaging-derived surrogate signals of both LGS and SOL muscles. There are three sub-models involved in the HNM: (a) sEMG-US imaging-derived weighted muscle activation model, (b) muscle-tendon unit geometry model, and (c) muscle contraction dynamic model.

#### Weighted muscle activation model

The neural activation at \(t_{k}\) time instant \(N_{i}(t_{k}),\,i=1,2\), \(k=1,2,3,...\), for SOL and LGS muscles, respectively, considers the electromechanical delay (EMD), \(\tau\), between the onset of an sEMG signal and a muscle contraction and utilizes a second-order recursive filter and is defined as [39]

$$\begin{aligned} N_{i}(t_{k})=\alpha _{i}u_{i}(t_{k-\tau })-\beta _{1i}N_{i}(t_{k-1})-\beta _{2i}N_{i}(t_{k-2}) \end{aligned}$$

(1)

where EMD, \(\tau\), is usually between 30 ms and 120 ms. \(u_{i}(t_{k})\) represents the sEMG’s linear envelope normalized to the peak value of the specific task (with a constant peak value across all walking speeds on each subject in this study). The sEMG’s linear envelope was derived after raw sEMG signals’ band-pass filtering, full-wave rectification, and low-pass filtering. \(\alpha _{i},\) \(\beta _{1i}\), and \(\beta _{2i}\) are coefficients that define the recursive filter’s dynamics of each muscle [39], and the following set of constrains is employed to reach a positive stable solution, i.e.

$$\begin{aligned} \beta _{1i}=\gamma _{1i}+\gamma _{2i},\,\beta _{2i}=\gamma _{1i}\cdot \gamma _{2i},\,\alpha _{i}-\beta _{1i}-\beta _{2i}=1 \end{aligned}$$

(2)

where \(\left| \gamma _{1i}\right| <1\) and \(\left| \gamma _{2i}\right| <1\).

A nonlinear relationship between neural activation \(N_{i}(t_{k})\) and corresponding muscle activation \(a_{1i}(t_{k})\) is given as [39]

$$\begin{aligned} a_{1i}=\frac{e^{A_{i}N_{i}}-1}{e^{A_{i}}-1} \end{aligned}$$

(3)

where \(A_{i}\) represents the nonlinear shape factor for each muscle, which is allowed to vary between − 3 and 0, with \(A_{i}=-3\) being a nonlinear and \(A_{i}=0\) being a linear relationship.

According to the reported results in [40,41,42,43], the MT change was found to correlate with the muscle contraction level or muscle activation through a linear function or piece-wise linear function. Therefore, in this work, the second part of the LGS or SOL muscle activation is calculated from the US imaging-derived MT change. The MT values of LGS and SOL muscles are denoted as \(MT_{i}(t_{k})\). According to the preliminary results in [44], there is a positive relationship between MT change and the targeted muscle contraction level during the walking stance phase. After taking normalization of the MT with respect to the lower bound (at rest state) and upper bound (peak value of the specific task) on each participant, the US imaging-derived muscle activation, \(a_{2i}(t_{k})\), is proposed as

$$\begin{aligned} a_{2i}=\frac{MT_{i}-MT_{i\min }}{MT_{i\max }-MT_{i\min }} \end{aligned}$$

(4)

where \(MT_{i\max }\) and \(MT_{i\min }\) represent the constant subjective MT values when the LGS and SOL muscles are at the walking task-specific maximum voluntary contraction and complete rest condition, respectively. For each individual, \(MT_{i\max }\) and \(MT_{i\min }\) are set as consistent values across different walking speeds. This normalization guarantees that the US imaging-derived muscle activations vary between 0 and 1.

By introducing an allocation gain between the sEMG- and US imaging-derived muscle activations, the synthesized/weighted muscle activation levels for LGS or SOL muscle, \(a_{i}(t_{k})\), can be represented as

$$\begin{aligned} a_{i}=\delta _{i}a_{1i}+(1-\delta _{i})a_{2i} \end{aligned}$$

(5)

where \(\delta _{i}\in [0,\,1]\) represents the muscle activation allocation gain for LGS (\(i=1\)) and SOL (\(i=2\)) muscles.

#### Muscle-tendon unit geometry model

As presented in Fig. 1b, consider the ankle joint rotation center in the sagittal plane as *O*, the proximal and distal osteotendinous junction points of each each muscle-tendon unit (MTU) near the knee joint and heel as \(A_{i}\) and \(B_{i}\), and the angle between \(OA_{i}\) and \(OB_{i}\) as \(q(t_{k})\), then each MTU length, \(l_{mt_{i}}(t_{k})\), is represented as the distance between \(A_{i}\) and \(B_{i}\), and is calculated as

$$\begin{aligned} l_{mt_{i}}=\sqrt{l_{OA_{i}}^{2}+l_{OB_{i}}^{2}-2l_{OA_{i}}l_{OB_{i}}\cos (q)} \end{aligned}$$

(6)

where \(l_{OA_{i}}\) and \(l_{OB_{i}}\) represent the distance of \(OA_{i}\) and \(OB_{i}\) obtained from OpenSim (National Institutes of Health for Biomedical Computation, Stanford, USA) [45], respectively. A generic OpenSim model (*gait*2392) was linearly scaled to each participant in OpenSim version 4.1 [45] per [46]. According to the law of sines, the moment arm of each MTU, \(r_{mt_{i}}(t_{k})\), is calculated as

$$\begin{aligned} r_{mt_{i}}=\frac{\partial l_{mt_{i}}(q)}{\partial (q)}=\frac{2l_{OA_{i}}l_{OB_{i}}\sin (q)}{l_{mt_{i}}}. \end{aligned}$$

(7)

The MTU generates contraction force only when it is stretched, which indicates the current tendon length \(l_{t_{i}}(t_{k})\) is equal to or longer than the tendon slack length \(l_{t_{i}}^{sk}\). From the perspective of muscle contraction dynamics that is shown in Fig. 1c, the overall MTU length, \(l_{mt_{i}}(t_{k})\), could also be expressed as

$$\begin{aligned} l_{mt_{i}}=l_{t_{i}}+l_{m_{i}}\cos (\phi _{i}) \end{aligned}$$

(8)

where \(l_{t_{i}}(t_{k})\) and \(l_{m_{i}}(t_{k})\) represent the current tendon length and muscle fascicle length for LGS and SOL muscles, respectively. \(\phi _{i}(t_{k})\) represents the pennation angle that changes with instantaneous \(l_{m_{i}}(t_{k})\). By assuming the individual muscle belly has a relatively small MT change and volume change [47, 48], \(\phi _{i}(t_{k})\) can be approximately calculated as

$$\begin{aligned} \phi _{i}\approx \arcsin \left( \frac{l_{m_{i}}^{0}\sin \phi _{i}^{0}}{l_{m_{i}}}\right) \end{aligned}$$

(9)

where the term \(\phi _{i}^{0}\) denotes the constant pennation angle when the muscle is at optimal fascicle length \(l_{m_{i}}^{0}\). The work in [49, 50] has shown that the optimal fascicle length increases as muscle activation decreases. Due to this coupling, the following relationship between the muscle activation and corresponding optimal fascicle length is given as

$$\begin{aligned} l_{m_{i}}^{0}=l_{mo_{i}}^{0}(\lambda (1-a_{i})+1) \end{aligned}$$

(10)

where \(\lambda\) is the rate of change in the optimal fascicle length, and it is selected as 0.15 [39]. \(l_{mo_{i}}^{0}\) is the optimal fascicle length at the walking task-specific maximum voluntary contraction, and \(l_{m_{i}}^{0}(t_{k})\) is the optimal fascicle length at time \(t_{k}\) and muscle activation \(a_{i}(t_{k})\).

From the preliminary results of plantarflexors’ US imaging in [44], it is very challenging to capture the entire fascicle length \(l_{m_{i}}(t_{k})\) of LGS or SOL muscles by using the current US transducer with a small dimension (width of 38 mm). Therefore, the fascicle length parameters of either LGS or SOL muscles were not measured directly from the US images. Instead, an indirect calculation method was applied as detailed below. According to [51], the nonlinear nominal tendon force-tendon strain relationship is given as

$$\begin{aligned} F_{t_{i}}(\xi _{i})=\left\{ \begin{array}{ll} 0, &{} \xi \le 0\\ 1480.3F_{i}^{\max }\xi _{i}^{2}, &{} 0<\xi <0.0127\\ (37.5\xi _{i}-0.2375)F_{i}^{\max }, &{} \xi \ge 0.0127 \end{array}\right. \end{aligned}$$

(11)

where \(\xi _{i}=\frac{l_{t_{i}}(t_{k})-l_{t_{i}}^{sk}}{l_{t_{i}}^{sk}}\) represents the tendon strain and \(F_{i}^{\max }\) represents the muscle contraction force at the walking task-specific maximum voluntary contraction. By considering \(F_{t_{i}}(\xi _{i})=F_{mt_{i}}\), where \(F_{mt_{i}}\) is defined in (13), \(l_{t_{i}}(t_{k})\) can be numerically computed based on the Runga–Kutta–Fehlberg algorithm. By substituting (6) and (9) to (8), the muscle fascicle length \(l_{m_{i}}(t_{k})\) will be calculated, as well as the muscle fascicle velocity \(v_{m_{i}}(t_{k})\) by taking the time derivative of \(l_{m_{i}}(t_{k})\).

#### Muscle dynamic contraction model

In the sEMG-US imaging-driven HNM, the individual moment component produced by each muscle, \(M_{i}(t_{k})\), can be represented as

$$\begin{aligned} M_{i}=F_{mt_{i}}r_{mt_{i}} \end{aligned}$$

(12)

where \(r_{mt_{i}}(t_{k})\) is defined in the above geometry model and \(F_{mt_{i}}(t_{k})\) denotes the corresponding contraction force generated on each MTU and is represented as

$$\begin{aligned} F_{mt_{i}}=(F_{ce_{i}}+F_{pe_{i}})\cos (\phi _{i}) \end{aligned}$$

(13)

where \(\phi _{i}(t_{k})\) represents the pennation angle defined above. \(F_{ce_{i}}(t_{k})\) and \(F_{pe_{i}}(t_{k})\) denote the corresponded forces generated by the parallelly located contractile element and the passive element, and can be calculated as

$$\begin{aligned} \begin{array}{lll} F_{ce_{i}} &{} = &{} F_{i}^{\max }f_{l_{i}}(l_{m_{i}})f_{v_{i}}(v_{m_{i}})a_{i}\\ F_{pe_{i}} &{} = &{} F_{i}^{\max }f_{p_{i}}(l_{m_{i}}) \end{array}. \end{aligned}$$

(14)

In (14), \(F_{i}^{\max }\) of the LGS and SOL muscles will be identified based on optimization algorithm in the HNM calibration procedures. \(f_{l_{i}}(l_{m_{i}}(t_{k}))\), \(f_{v_{i}}(v_{m_{i}}(t_{k}))\), and \(f_{p_{i}}(l_{m_{i}}(t_{k}))\) denote the generic muscle contractile force-fascicle length, force-fascicle velocity, and passive elastic force-fascicle velocity curves. These curves were normalized to \(F_{i}^{\max }\), optimal fascicle length \(l_{m_{i}}^{0}\), and maximum fascicle contraction velocity \(v_{m_{i}}^{\max }\). The explicit expressions of \(f_{l_{i}}(l_{m_{i}}(t_{k}))\), \(f_{v_{i}}(v_{m_{i}}(t_{k}))\), and \(f_{p_{i}}(l_{m_{i}}(t_{k}))\) can be found in [25, 32, 52].

### Experimental protocol, data collection and pre-processing

The study was approved by the Institutional Review Board (IRB) at North Carolina State University (IRB number: 20602). Ten young participants (7M/3F, age: 25.4 ± 3.1 years, height: 1.77 ± 0.10 m, mass: 78.0 ± 21.1 kg) without any neuromuscular or orthopedic disorders, were recruited in this study. Every participant got familiar with the experimental procedures and signed an informed consent form before participation.

Figure 2a summarizes the experimental setup for the participants to perform static anatomical poses and dynamic gait trials (walking speeds changing from 0.50 to 1.50 m/s), and Fig. 2b presents the workflow of the data processing and sEMG-US imaging-driven HNM calibration procedures. During all walking trials, three-dimensional coordinates of 39 retro-reflective markers positioned on the participant’s lower extremities following the instructions in [53] were recorded using a 12-camera motion capture system (Vicon Motion Systems Ltd, Los Angeles, CA, USA) at 100 Hz. The GRF signals were collected at 1000 Hz synchronously with makers trajectories using in-ground force plates (AMTI, Watertown, MA, USA) mounted on an instrumented treadmill (Bertec Corp., Columbus, OH, USA) through the commercial real-time data capture software Nexus 2.9. Both GRF signals and markers trajectories were low-pass filtered with a fourth-order Butterworth filter in Visual 3D software (C-Motion, Rockville, MD, USA), and the cut-off frequencies were set as 25 Hz and 6 Hz, respectively. The markers trajectories and GRF signals from static poses and dynamic gait trials were used for joints kinematics and net joint moment calculation by using inverse kinematics and ID algorithms in Visual 3D. During all trials, four wired sEMG sensors (SX230, Biometrics Ltd, Newport, UK) were attached to the shank skin through a double-sided tape to non-invasively record the electrical signals from tibialis anterior (TA), LGS, MGS, and SOL muscles at 1000 Hz synchronously through Nexus 2.9. The locations for sEMG sensors were determined following the instructions in [54]. A US transducer (38-mm of length, 6.4 MHz center frequency, L7.5SC Prodigy Probe, S-Sharp, Taiwan) setup was applied in the walking experiments as described in [34]. The US radio frequency (RF) data were recorded at 1000 frames per second synchronously with markers trajectories, GRF, and sEMG signals, using a pulse sequence trigger signal generated from Nexus 2.9.

sEMG signals were band-pass filtered with the bandwidth between 20 Hz and 450 Hz, and then full-wave rectified and low-pass filtered with a cut-off frequency of 6 Hz. The resulting linear envelopes were normalized to the peak processed sEMG values obtained from all sets of dynamic gait trials with different speeds. The US radio frequency data were converted to B-mode images through beamforming and logarithmic compression, and then a commercial US imaging processing toolbox *UltraTrack* [55] was applied to determine the temporal sequences of the LGS and SOL muscles’ thickness based on the adaptive optical flow tracking algorithms. First, one region of interest (width by height: 336 \(\times\) 400 pixels) that encompassed the LGS and SOL muscles was selected as the area between the superficial and deep aponeuroses. Then, 10 vertical lines were manually defined with evenly distributed distances in the LGS’s region and SOL’s region on the first US imaging frame for each recorded waking trial, as shown in Fig. 2b. Key-frame correction [55] was applied to minimize the time-related drift of MT’s cyclical pattern over multiple gait cycles, where the key frames were selected to be at heel-strike and toe-off time points. After each correction, the new key frames’ vertical lines positions were determined by applying an affine transformation to the key-frame before it. Finally, the temporal sequence of the mean value of those 10 vertical lines’ lengths from LGS and SOL muscles were calculated for LGS’s MT and SOL’s MT, respectively. The MT signals were then low-pass filtered with a cut-off frequency of 30 Hz.

Before data collection, participants practiced dynamic walking steadily with all sensing devices attached to the lower extremities on the treadmill for at least 30 s for each speed. The duration of each walking trial was set as 2 min, and data from the middle 20 s were collected for analysis. In the current study, there were five different walking speeds, 0.50 m/s, 0.75 m/s, 1.00 m/s, 1.25 m/s, and 1.50 m/s, and the order was selected randomly for each participant. The participants were provided with at least 2 min for rest between two successive dynamic walking trials to avoid muscle fatigue.

### Data post-processing

#### HNM model calibration

To determine the values for a set of parameters, including the shape factor \(A_{i}\), the tendon slack length \(l_{t_{i}}^{sk}\), the muscle contraction force at the walking task-specific maximum voluntary contraction \(F_{i}^{\max }\), and the muscle activation allocation gain \(\delta _{i}\), the HNM model calibration was performed with the initial parameters obtained from the literature [56]. The summarized calibration process of the sEMG-US imaging-driven HNM is presented in Fig. 2b. Indeed, the model calibration is a nonlinear optimization problem, where the above parameters set needs to be found for minimizing the designated objective function. During calibration, the parameters were subsequently adjusted within the predefined boundaries, which ensured the muscle-tendon unit always operated within the physiological range [39]. Except for the allocation gains \(\delta _{i}\) were set the range between 0 and 1, other parameters were set with the lower and upper bounds as 50% and 150% of the literature-referred values, while the initial guess of the set of parameters were set as the middle value of each parameter. Other HNM parameters, like muscle-tendon unit length \(l_{mt_{i}}\), optimal muscle fascicle length \(l_{mo_{i}}^{0}\), and optimal pennation angle \(\phi _{i}^{0}\) were assigned using OpenSim and the scaling method as references. The Matlab function ’lsqcurvefit’ with the Levenberg-Marquardt algorithm was applied to solve the nonlinear least-squares optimization problem until the following objective function was minimized during the calibration

$$\begin{aligned} E_{obj}=\frac{1}{N}\sum _{j=1}^{N}\left( \left( \sum _{i=1}^{2}M_{i}(j)\right) -M_{w}(j)\right) ^{2} \end{aligned}$$

(15)

where \(\sum _{i=1}^{2}M_{i}(j)\) represents the net PF moment estimation from both LGS and SOL muscles during stance phase at time instant *j*, \(M_{w}(j)\) represents the net PF moment measurement from ID at time instant *j*, and *N* represents the length of the data used for the calibration.

It should be noted that the HNM calibration is subjective, and for each participant, the calibration procedure was repeated for both single-speed modes and inter-speed mode. In the current study, single-speed modes are defined as the HNM calibration only with data collected from one single speed (five single-speed modes here), while the inter-speed mode is defined as the HNM calibration with all data collected from five different speeds (one inter-speed mode here). In HNM calibration, collected data, including ankle joint’s kinematics, kinetics, sEMG signals, and US imaging from the stance phase of 10 steady gait cycles, were used to minimize the objective function in (15). For the single-speed modes, the 10 steady gait cycles were all from one speed, while for the inter-speed mode, the 10 gait cycles were composed of two steady cycles from each of five speeds. Similarly, by manually setting the muscle activation allocation gain \(\delta _{i}\) to be 1 and 0, we could also conduct the calibration procedures of the sEMG-driven and US imaging-driven HNMs. After the HNM calibration with different modes, new data sets (five gait cycles not involved in the calibration procedures) from walking trials at different speeds were used as input of the calibrated HNMs to predict net PF moment. The prediction results were evaluated through the comparison to the measured net PF moment from ID.

#### Statistical analysis

The validation procedures comprised three tests to assess the sEMG-US imaging-driven, sEMG-driven, and US imaging-driven HNMs’ calibration and prediction abilities of the net ankle joint PF moment during the stance phase. In each type of HNM test, the root mean square error (*RMSE*), *RMSE* normalized to individual peak net PF moment (\(N-RMSE\)), *RMSE* normalized to individual body mass (\(BM-RMSE\)), and coefficient of determination (\(R^{2}\)) values between the calibrated/predicted net PF moment and ID-calculated net PF moment were calculated to evaluate the calibration/prediction performance. The rational of using these three tests is given here. From a mathematical perspective, the proposed sEMG-US imaging-driven HNM is a model-based regression problem. Given that the corresponding data were all collected continuously and synchronously, and the ground truth from the ID calculation and the prediction from different HNMs are all continuous, the most intuitive metric to evaluate the prediction accuracy is to use the *RMSE* between the ground truth and prediction. However, due to the weight, height, and walking pattern’s variations among subjects, the capability (peak value) of net plantarflexion moment generation during treadmill walking at the same speed is different, which would affect the prediction *RMSE* directly. Therefore, the direct *RMSE* values comparison among subjects is likely to introduce high deviations. An effective way to reduce the subjective variation of the prediction *RMSE* is to take the normalization to a subjective characteristic, like the body mass and peak net plantarflexion moment as we selected. The relationship between \(N-RMSE\) and \(BM-RMSE\) is that in the normalization calculation, they have same numerators but different denominators for the same subject. Although \(N-RMSE\) and \(BM-RMSE\) might provide potential redundancy information as can be seen from subsequent calibration and prediction summary results, they show different sensitivities to the walking speeds (\(BM-RMSE\) is more sensitive). Another useful metric to evaluate the prediction accuracy is the coefficient of determination, known as \(R^{2}\), which ranges from 0 to 1 and measures the proportion of variation in the data that is accounted for in the model. In the current study, \(R^{2}\) is used as a parallel evaluation metric as N-RMSE and BM-RMSE, and it is relatively not sensitive to the walking speeds. Although \(R^{2}\) and \(N-RMSE\) might provide potential redundancy information, it is not always true that lower \(N-RMSE\) corresponds to higher \(R^{2}\), so we keep both to provide potential supplementary information.

Shapiro-Wilk parametric hypothesis test was used to determine the normality of the corresponding \(N-RMSE\), \(BM-RMSE\), and \(R^{2}\) values of each calibration step under the single-speed modes and inter-speed mode, and each prediction step. According to the results of the Shapiro-Wilk test, two-way repeated-measure analysis of variance (ANOVA) or Friedman’s tests followed by a Tukey’s honestly significant difference tests (Tukey’s HSD) was applied to evaluate the effect of HNMs’ category and walking speed on different evaluation criteria, including \(N-RMSE\), \(BM-RMSE\), and \(R^{2}\) values, during both calibration and prediction procedures. The significant difference level was chosen as \(p<0.05\) for all statistical tests. Effect sizes were reported as \(\eta _{p}^{2}\) and Cohen’s *d* for main effects from ANOVA or Friedman’s tests and pairwise comparisons from Tukey’s HSD, respectively.