Predicting Free Achilles Tendon Strain From Motion Capture Data Using Artificial Intelligence

The Achilles tendon (AT) is sensitive to mechanical loading, with appropriate strain improving tissue mechanical and material properties. Estimating free AT strain is currently possible through personalized neuromusculoskeletal (NMSK) modeling; however, this approach is time-consuming and requires extensive laboratory data. To enable in-field assessments, we developed an artificial intelligence (AI) workflow to predict free AT strain during running from motion capture data. Ten keypoints commonly used in pose estimation algorithms (e.g., OpenPose) were synthesized from motion capture data and noise was added to represent real-world data obtained using video cameras. Two AI workflows were compared: (1) a Long Short-Term Memory (LSTM) neural network that predicted free AT strain directly (called LSTM only workflow); and (2) an LSTM neural network that predicted AT force which was subsequently converted to free AT strain using a personalized force-strain curve (called LSTM+ workflow). AI models were trained and evaluated using estimates of free AT strain obtained from a validated NMSK model with personalized AT force-strain curve. The effect of using different input features (position, velocity, and acceleration of keypoints, height and mass) on free AT strain predictions was also assessed. The LSTM+ workflow significantly improved the predictions of free AT strain compared to the LSTM only workflow (p < 0.001). The best free AT strain predictions were obtained using positions and velocities of keypoints as well as the height and mass of the participants as input, with average time-series root mean square error (RMSE) of 1.72±0.95% strain and r2 of 0.92±0.10, and peak strain RMSE of 2.20% and r2 of 0.54. In conclusion, we showed feasibility of predicting accurate free AT strain during running using low fidelity pose estimation data.


I. INTRODUCTION
T HE Achilles tendon (AT) transmits force from the triceps surae muscle group to the calcaneus during locomotion [1]. Loading of the AT is a fundamental regulator of tissue health and a main driver of the mechanobiological cascade involved in tissue adaptation [2]. In vivo interventions involving modulation of AT strain, such as variations in strain magnitude, rate, and frequency, have been shown to affect tendon stiffness in humans via changes in morphological and/or material properties [3]. In particular, Arampatzis et al. [4] reported that 14 weeks of AT training (4 times/week, 5 sets of repetitions, 3s loading and 3s relaxation) at high AT strain magnitudes resulted in significantly increased tendon stiffness, whereas training at lower strain magnitudes had no effect on tendon stiffness. Ex vivo bioreactor experiments on animal tendons have also demonstrated that the AT responded to strain magnitudes of 6% with increased expression of collagen type I and type III, altered apoptosis rate, and changes in tissue mechanical properties at intermediate compared to higher (9%) and lower (3%) strain magnitudes [5], [6]. Based on these observations, and with similar findings reported for other musculoskeletal tissues, e.g., bone [7], cartilage [8], it has been proposed that an AT 'sweet spot' (i.e., optimal strain range) may exist for maintaining tissue health and promoting repair [9]. Estimating real-time strain magnitudes could inform AT training and rehabilitation programs to target the AT 'sweet spot' for each individual, possibly maximizing positive tendon adaptation.
The AT is a complex three-dimensional (3D) anatomical structure comprising the free tendon and aponeurosis [10]. The free AT is devoid of muscle fiber attachment along its length and comprises the tendinous tissue between the calcaneus insertion and the soleus muscle-tendon junction. In contrast, the aponeurosis is a proximal extension of the free tendon that attaches to muscle fibers along its length.
Whereas the free AT experiences large longitudinal strains during functional tasks [11], the aponeurosis appears to be subjected to smaller and more complex 3D strains during voluntary muscle contractions [12]. One of the challenges in the field of tendon biomechanics is that AT strain is variously reported in the literature as either free tendon strain or free tendon plus aponeurosis strain (commonly defined as the tendon tissue between the calcaneus and medial gastrocnemius muscle tendon junction), which has given rise to a large disparity in AT strain estimates and makes it difficult to compare AT strain estimates between studies. Generally speaking, the aponeurosis of the AT is about double the length of the free AT and the longitudinal strains experienced by the free AT are at least double those experienced by the free tendon and aponeurosis [13]. Given that the free tendon undergoes large strains and is a common site of tendon injury and degeneration [13], [14], understanding the strains experienced by the free AT is of critical interest to improve tendon function, prevent tendon overstrain injuries during physical performance, and best rehabilitate injured or diseased tendons.
In vivo AT strain estimates during voluntary contraction of the calf muscles are commonly obtained from tendon deformations measured using B-mode ultrasound. This twodimensional (2D) ultrasound modality has provided unique insight into how the triceps surae muscles interact with the AT during dynamic tasks such as walking, running and jumping to enhance power production and movement efficiency of the muscle-tendon unit [15], [16]. However, a limitation of ultrasound is that the measurements can be affected by substantial movement artefact during dynamic movements and cannot account for 'out-of-plane' deformation caused by tissue rotation and/or twisting [17]. The latter problem can be overcome using freehand 3D ultrasound, which integrates 2D B-mode ultrasound with 3D motion capture [18]. However, the freehand 3D ultrasound approach is limited to the assessment of static tendon geometry. An alternative to ultrasound-based approaches is neuromusculoskeletal (NMSK) modeling, which is a method to estimate the internal tissue state from external biomechanical measurements [19]. Electromyogram (EMG)informed NMSK modeling that accounts for the individual's neural patterns has been shown to estimate physiologically plausible joint and tissue loading and muscle forces during a variety of motion tasks and across different populations [20], [21], [22], [23]. The triceps surae force can be estimated and converted to free AT strain by personalized force-strain curve obtained using freehand 3D ultrasound [11], [24]. However, EMG-informed NMSK modeling approach can be highly labor intensive and is typically confined to the laboratory environment.
Computer vision approaches combined with artificial intelligence (AI) have enabled low-cost and rapid tracking of human movement outside of the laboratory. Established tools, such as OpenPose [25], have previously demonstrated high validity and reliability in estimating complex whole-body kinematics of human movement. With these pose estimation techniques, it is possible to use one or more cameras to track human movement in real-world conditions. Various deep learning models have also been used to predict external biomechanical variables. For example, ground reaction forces (GRFs), joint angles and joint moments were accurately predicted using marker data or inertial measurement unit (IMU) signals when compared to traditional physics-based approaches (e.g., inverse kinematic and inverse dynamics) [26], [27], [28]. A neural network was similarly used to predict peak knee adduction moment using synthesized pose estimation data, suggesting the potential to inform the treatment of patients with knee osteoarthritis within clinical settings [29]. An arguably more challenging problem for the field is predicting internal tissue mechanical states such as muscle, tendon and ligament forces and/or strains from pose estimation data. Several recent studies that used the Convolutional Neural Network (CNN) and Long Short-term Memory (LSTM) neural network have suggested that accurate prediction of lower limb joint loadings and muscle forces in a variety of movements is possible [30], [31]. Importantly, the CNN model was real-time capable [30], thereby greatly reducing the computational time compared to NMSK modeling method. At present, however, it is unknown whether an LSTM neural network would be able to accurately predict free AT strain based on low fidelity keypoint data. Free AT strain depends on the mechanical stiffness of the tendon and the force applied to it, which in turn depends on anatomy (e.g., musculotendon units origin and insertion, moment arms), multibody dynamics (e.g., inertial properties of the segment, joint mechanics), and motor control (e.g., muscle coordination and activation). These factors are subject-specific and may not be encoded in pose and anthropometric data. Thus, it is unclear whether an LSTM could accurately predict strain directly from pose estimation data, or whether augmenting an LSTM that predicts tendon force with a personalized constitutive model of the AT would improve strain predictions.
The purpose of this study was to determine if an LSTM model combined with a personalized constitutive model of the free AT force-strain relationship (i.e., LSTM+ workflow) improved the prediction of free AT strain relative to the LSTM only workflow. The effect of different combinations of input features selected from the positions, velocities, and accelerations of keypoints, as well as the anthropometric characteristics of the participants (e.g., height and mass), on the accuracy of free AT strain prediction, was also evaluated. To assess the potential clinical utility of the approach, the accuracy of the best performing prediction model was compared to literature data relating tendon strain to tendon remodeling [4]. We hypothesized that the LSTM+ workflow would result in better free AT strain prediction than the LSTM only model.

A. Experimental Data Collection
Sixteen trained middle-distance runners (female: 6, age: 25.2±5.0 yr, height: 175.5±7.3 cm, body mass: 64.4±8.4 kg) with no history of AT injuries performed at least 5 trials of overground running at two speeds (3.0±0.3 m/s and 5.0±0.5 m/s) in a biomechanics laboratory. 3D makers trajectories, GRFs, and EMG signals of lower limb muscles were synchronously collected through motion capture system (sampling frequency: 250 Hz; Vicon Vantage Cameras, Vicon Motion Systems Ltd., Oxford, United Kingdom), force plates (sampling frequency: 1000 Hz; Kistler Instrument Corporation, Amherst, NY) and EMG devices (sampling frequency: 1500 Hz; TELEmyo DTS, Noraxon U.S.A. Inc., Scottsdale, AZ). Magnetic resonance imaging (MRI) of the ankle joint and freehand 3D ultrasound of AT during isometric ankle plantarflexion tasks were employed to personalize geometry and mechanical properties of the NMSK model. Biomechanical data from this same cohort were previously published [11]. The study was approved by the institutional human research ethics committee and all participants provided written informed consent prior to their involvement in the study.

B. Generating Reference Free AT Strain From EMG-Informed NMSK Modeling
Free AT strain was obtained through established personalized NMSK modeling pipeline ( Figure 1), with detailed information regarding the processing steps reported in our previous studies [11], [32]. Marker trajectories, GRFs and EMGs were pre-processed prior to analysis using OpenSim [33]. The marker trajectories and GRFs were filtered at 10 Hz cut-off frequency using a zero-lag 4 th order low-pass Butterworth filter. Surface EMGs were bandpass filtered (zero-lag 4 th order bandpass Butterworth filter, 20-450 Hz), full-wave rectified, low-pass filtered and then normalized by maximal voluntary contractions to calculate muscle excitations [34]. A generic OpenSim model (gait2392) [19] was linearly scaled to each participant using anatomical markers acquired during the static trial. To personalize the model, the maximal isometric forces of the muscles in the model were modified based on regression equation for each participant [35]. Then, the moment arms of the medial gastrocnemius, lateral gastrocnemius, and soleus muscles were adjusted to match the experimental free AT moment arm obtained from MRI [36]. The Hill-type muscle model parameters used for NMSK model personalization were optimized via morphological scaling [37]. The force-strain curve for the AT was personalized for each participant by fitting the piecewise function Eq.1 [38] to the experimentally measured AT elongation (i.e., freehand 3D ultrasound) and force data (i.e., ankle joint torque divided by AT moment arm calculated from MRI).
where the a, m, q, and ε 0 ensure C 1 continuity; the f represents force and ε represents strain. The force-strain curve for the other musculotendon units in the model was based on literature data [39]. After completing model preparation, a physics-based modeling pipeline including residual algorithms, inverse kinematics, inverse dynamics, and muscle analysis was run in OpenSim to obtain joint angles, joint moments, musculotendon lengths, and moment arms across all trials [40]. Within the calibrated EMG-informed NMSK modeling toolbox (CEINMS) [41], a calibration procedure was performed to tune model parameters of musculotendon units crossing the ankle, subtalar, and knee joints [41], [42]. The experimental muscle excitations were then mapped to the complete set of muscles in the model [43]. The calibrated model was used to estimate triceps surae muscle forces based on musculotendon kinematics and muscle excitations using the EMG-assisted mode in CEINMS [44]. The AT force was calculated as the sum of medial gastrocnemius, lateral gastrocnemius, and soleus forces. Finally, free AT strain for each participant was obtained using the personalized AT forcestrain curve. AT force and free AT strain estimates from this pipeline were used as target outputs for the following AI workflows.

C. Predicting Free AT Strain From AI Workflows
Two AI workflows (i.e., LSTM only and LSTM+) were developed and evaluated in their ability to predict free AT strain (Figure 1). The LSTM only workflow was investigated firstly to predict the free AT strain from synthesized pose estimation data. The LSTM+ workflow was developed integrating the LSTM model with a constitutive model of AT (i.e., forcestrain curve) previously described in Section B. Specifically, the LSTM+ workflow included an LSTM model to predict the AT force output from NMSK modeling pipeline, which was subsequently used as input to the personalized AT force-strain curve to obtain free AT strain. For both workflows, ten frames of time window data were used to predict free AT strain value in the frame following the selected time window.
1) Data Preparation: Maker trajectory data were spacenormalized to remove positional offset in the anterior-posterior and left-right directions, therefore translated by subtracting the initial values from the landmark of 7 th cervical (C7) in each trial. Ten keypoints selected from pose estimation model were synthesized from relevant maker trajectories [25]. Specifically, the bilateral knee, ankle joint centers calculated as the middle point of corresponding anatomical markers on each joint. The hip keypoints were created from bilateral anterior and posterior superior iliac spines via linear regression [45]. The bilateral toe keypoints were directly extracted from corresponding markers. Neck and pelvis keypoints were calculated as the middle of anatomical markers on acromia and hip keypoints, respectively. The dataset was augmented 10 times by adding noise to the x, y, and z coordinates of each keypoint. Specifically, the noise was applied using axis-specific and keypoint-specific gaussian distributions representing the errors of OpenPose compared to gold standard optical motion capture [46]. As the errors of neck, pelvis, and toe keypoints were not reported, the noise distribution of adjacent keypoints were used. Noisy position data was low-pass filtered using a 2 nd order Butterworth low-pass filter with a cut-off frequency of 10 Hz, which is appropriate for real-time applications [47]. Velocities and accelerations were calculated using backward numerical differentiation of filtered keypoint positions. Through this approach, a total of 1510 trials containing positions, velocities, and accelerations of synthetized keypoints were generated for subsequent model training.
Three combinations of input features were investigated. The first combination of input features comprised of 3D positions, velocities, and accelerations of the keypoints combined with height and mass of each participant (i.e., Pos-Vel-Acc), for a total of 92 features. Pos-Vel and Pos combinations stepwise removed accelerations and both accelerations and velocities of the keypoints from Pos-Vel-Acc, resulting in 62 and 32 total features respectively. For each trial, the input features were cropped by a sliding window of length of ten frames to generate input matrixes (Pos-Vel-Acc: 10 frames × 92 features; Pos-Vel: 10 frames × 62 features; Pos: 10 frames × 32 features) for model training. All the input features were mean-removed and scaled to unit variance to facilitate convergence of gradient descent [48]. Target values (i.e., free AT strain for LSTM only workflow or AT force for LSTM+ workflow) were extracted from the NMSK workflow outputs to pair each input time window.
2) AI Workflow Development: Both workflows implemented hyperparameter tuning and Leave-One-Out Cross Validation (LOOCV) for model development and evaluation. The LSTM model with two bi-directional LSTM layers was selected from preliminary model training. Each bi-directional LSTM layer followed with one BatchNormalisation layer and one DropOut layer was introduced to minimize overfitting and improving model prediction performance [49], [50]. Hyperparameter tuning was implemented using hyperband algorithm [51] to tune the number of neurons of each layer, dropout rate for DropOut layers, learning rate and batch size in reasonable ranges. The reduction factor for the hyperband algorithm was set to 3, and the maximal epochs of the hyperband tuning was 500. The weights of the LSTM model optimized by Adam [52] and the early stopping with 10 epochs of patience was used to monitor the validation loss during tuning process. The hyperparameter tuning was run using TensorFlow (v2.6.0) in Machine Learning eResearch Platform (MLeRP) (32GB RAM and Tesla A100 GPU) on cloud. The augmented dataset was split into 7:1 by randomly choosing 14 participants for training and 2 participants for validation. The hyperparameter tuning was run to search the best hyperparameters for free AT strain of LSTM only prediction and AT force of LSTM+ prediction among three different combinations of input features (Pos-Vel-Acc, Pos-Vel and Pos). The hyperparameters used in the final model are summarized in Table I. Finally, six models with three input combinations across two types of targeted outputs were generated and evaluated. The hyperparameters that achieved the best model performance were reloaded to each model for running LOOCV. In this step, early stopping with 20 epochs of patience was used for monitoring validation loss. A callback function was applied to save the best model parameters and was reloaded to make predictions for validation. For the LSTM+ workflow, the predicted AT forces were converted to free AT strain using personalized force-strain curve for further evaluation.

D. AI Workflow Evaluation and Statistical Analysis
For each trial, coefficient of determination (r 2 ), root mean square error (RMSE) and normalized RMSE (nRMSE) [53] were calculated for the time-series between the predicted free AT strain and estimates of free AT strain obtained from a NMSK model with personalized AT force-strain curve. As the RMSE values were not normally distributed, they were logarithmically transformed and then tested through a two-way ANOVA on two workflows by three combinations of input features at significance level 0.05 using SPSS 22.0, followed by the post hoc testing using Bonferroni correction. Also, the best predictions with the lowest RMSE for each workflow were subsequently selected to run 1D statistical parametric mapping (1D-SPM) to compare the predicted free AT strain with the NMSK workflow outputs in subject level. In details, as the raw dataset was augmented 10 times with noise, the average RMSE of each sub-dataset was calculated for different conditions. The predictions with lowest average RMSE for each workflow were selected to run two-tailed paired sample t-test using 1D-SPM and the differences were considered statistically significant for p-values < 0.01. To evaluate the scalar value prediction of AI workflows, the free AT peak strain values from the best LSTM+ and LSTM only workflow were extracted for each trial. Linear regression models were fitted between the best predicted and the ground truth peak free AT strain values, and the corresponding r 2 and RMSE for each trial were computed.
To address whether the best workflow was sufficiently accurate to be useful as a tool for studying tendon remodeling, the prediction error of peak strain magnitude was evaluated in the context of existing literature. Arampatzis et al. [4] reported that a strain magnitude during training of 4.55% was associated with positive tendon adaptation. As the strain threshold of 4.55% was for the free tendon plus aponeurosis (i.e., calcaneus to medial gastrocnemius muscle tendon junction), to facilitate comparison, the free AT tendon strain from the present study was converted to a corresponding free tendon plus aponeurosis strain. This conversion was achieved by dividing free tendon strain estimates from the present study by 2.2, which is the ratio of free AT strain to free AT plus aponeurosis strain reported by Magnussen et al. [13]. The error of the prediction from the best workflow in the present study was computed from the difference between the predicted free AT peak strain and ground truth from NMSK modeling and converted to corresponding values for free AT plus aponeurosis strain, expressed as the upper and lower limits of the 95% confidence interval. The error distribution was then applied to the mean of the free AT strain prediction and compared to the 4.55% AT strain threshold from Arampatzis et al. [4].

IV. DISCUSSION
Consistent with our hypothesis, the performance of LSTM+ workflow (i.e., LSTM model combined with a personalized constitutive model of the free AT force-strain relationship) yielded better predictions of free AT strain than the LSTM only workflow. Furthermore, the best predictions of free AT strain were obtained using input features that consisted of positions and velocities of keypoints as well as the height and mass of the participant. Prior studies reported good prediction of external biomechanics (i.e., joint angles and moments) from motion capture [27], [54]. Our study extended these findings by demonstrating that predictions of internal tissue mechanics (i.e., the free AT strain) may also be achieved from synthetized keypoints using an AI workflow. Given that free AT strain is a critical mechanical stimulus for tendon adaptation [9], [24], these findings provide some preliminary evidence in support of using strain prediction from motion capture data through AI for the purpose of guiding AT training and rehabilitation. LSTM+ better predicted free AT strain for both timeseries and peak value compared with the LSTM workflow, suggesting that the personalized force-strain curves used in the LSTM+ model contained unique information that could not be extracted from keypoint kinematics, height, and mass. Importantly, AT mechanical properties, such as Young's modulus and stiffness, vary with age [55], level of physical activity [32], and disease [56]. Previous modeling studies have highlighted the sensitivity of free AT strain prediction to subject-specific variation in AT geometry [57], suggesting that future implementations of our proposed LSTM workflow would likely require input describing personalized geometry (e.g., cross-sectional area) and/or mechanical properties (e.g., stiffness) of the AT. The best overall model predictions (i.e., lowest RMSE) were obtained using the Pos-Vel combination of input features in the LSTM+ workflow, indicating that employing derivatives of keypoint positional data as input features improved overall model performance. To reduce the risk of overfitting, we selected a shallow network architecture with only two bi-directional LSTM layers. It is therefore plausible that our neural network was unable to automatically extract velocity and acceleration features from the keypoint positional data, thus resulting in accuracy improvements when velocity and acceleration were explicitly provided as input. Similar feature engineering approaches are not uncommon, and have been applied to prediction of GRF from inertial measurement units, where mean, standard deviation, and range values from sensor data were explicitly provided in input to the neural network to increase prediction accuracy [58], [59].
The current study also evaluated whether the accuracy of the best performing prediction model (LSTM+ workflow) was sufficient to be useful in a practical context for studying tendon remodeling. The error in free tendon plus aponeurosis strain associated with the LSTM+ prediction of ±1.84% for a mean peak strain of 5.84% resulted in 90% of running trials experiencing strains that exceed the strain threshold of 4.55%, which has been reported to result in positive tendon remodeling [4]. While it is acknowledged that the exact modes of strain stimuli (i.e., magnitude, rate, frequency) that would cause adaptation is still a matter of ongoing investigation, the results from the present study suggest that our AT strain predictions are sufficiently accurate to distinguish the strain during running from those that have been reported to be too low to induce positive tendon remodeling. Caution however is warranted in generalizing these findings to other tasks, especially those associated with smaller tendon strains than running, as it might be expected that the error associated with our approach would overlap the 4.55% strain threshold to a greater extent than reported here for running. We therefore conclude that further refinement and evaluation would be necessary for the approach to have widespread clinical application.
An important strength of current study is that axis-and keypoint-specific gaussian noise were added to the keypoint data prior to input into the LSTM. Introducing noise to input data can improve the training of AI models and their overall robustness against noise [60]. This is particularly relevant in real-world scenarios where noise is always present when estimating human pose. Several limitations of our study should be acknowledged. Analyses were limited to running, thereby reducing generalizability of findings to other motor tasks. To enhance generalizability, a multi-task dataset relevant to AT training and rehabilitation (e.g., walking, heel rise and heel drop) is recommended for future investigations. Keypoints were synthesized from laboratory-based motion capture data; consequently, evaluation of performance when using keypoint data obtained in the field via video or depth cameras is yet to be performed. The best prediction of free AT strain was from the LSTM+ workflow, which required estimation of AT mechanical properties. However, performing freehand 3D ultrasound measurement to obtain AT mechanical properties relies on laboratory-based equipment and extensive data post processing. As such, future investigations should explore approaches for rapidly estimating AT properties and to evaluate model performance in real-time. These steps are critical to improving the efficiency and clinical utility of the proposed approach.

V. CONCLUSION
This study demonstrated the feasibility of using an LSTM model to predict time-series free AT strain through synthesized pose estimation data containing noise. Including personalized information of the AT further increased accuracy of free AT strain prediction compared to a ground truth estimate obtained using a validated NMSK modeling framework. The error associated with the tendon strain predictions was sufficient to distinguish the strains in running from those that are likely too low to lead to tendon remodeling, but the accuracy would likely need to be improved for the approach to have broad clinical utility. It is foreseeable that AT training could be soon performed in ecologically valid settings, enabling training approaches based on mechanobiology and removing the guesswork in exercise therapy.