The increased risk of cardiovascular disease in people with spinal cord injuries motivates 1 work to identify exercise options that improve health outcomes without causing risk of musculoskele- 2 tal injury. Handcycling is an exercise mode that may be beneficial for wheelchair users but further 3 work is needed to establish appropriate guidelines and requires assessment of the external loads. The 4 goal of this research was to predict the 6-degree-of-freedom external loads during handcycling from 5 data similar to that which can be measured from inertial measurement units (segment accelerations 6 and velocities) using machine learning. Five neural network models and two ensemble models were 7 compared against a statistical model. A temporal convolutional network (TCN) yielded the best pre- 8 dictions. Predictions of forces and moments in-plane with the crank were the most accurate (r = 0.95 - 9 0.97). The TCN model could predict external loads during activities of different intensities, making 10 it viable for different exercise protocols. The ability to predict the loads associated with forward 11 propulsion using wearable-type data enables the development of informed exercise guidelines. 12 1. Introduction 13 Full-time manual wheelchair users, such as persons with spinal cord injuries, are at 14 increased risk of chronic diseases such as cardiovascular disease (CVD), diabetes, and 15 metabolic syndrome [ 1–3]. Although the evidence supporting the benefits of physical 16 activity is clear, many individuals with spinal cord injuries are physically inactive and 17 sedentary leading to worsening symptoms and increased levels of disability [ 4,5]. An 18 optimal exercise for improving cardiopulmonary function requires sufficient taxation of 19 cardio-respiratory fitness levels without risking musculoskeletal injury. Unfortunately, the 20 transition to manual wheelchair use presents challenges as the arms become the source of 21 all activities of daily living. For example, as many as 68% of manual wheelchair users have 22 rotator cuff tendon tears [6]. 23 Handcycling has emerged as a promising exercise mode for wheelchair users [ 7] 24 because the loads on the shoulder are lower than everyday wheelchair propulsion [ 8]. 25 However, exercise guidelines (intensity, frequency, etc.) for handcycling remain to be 26 established and requires assessment of both cardiovascular stimulation and musculoskeletal 27 safety. Musculoskeletal assessments, such as shoulder torque, provide objective measures 28 of shoulder function but requires specialized equipment to measure the external hand 29 loads. 30 Custom-built instrumented crank handles have been used to measure handcycling 31 loads [ 8,9], but are not readily accessible. In contrast, kinematic data collection is commonly 32 used in biomechanics research and can be related to external loads via direct or indirect 33 methods. The ground reaction force during legged locomotion can be solved for directly 34 using Newton’s 2ndlaw if the kinematics of major body segments are measured [ 10,11]. 35 However, tracking all segments can be difficult and may not be feasible for every activity 36 Version July 3, 2024 submitted to Sensors https://www.mdpi.com/journal/sensors Version July 3, 2024 submitted to Sensors 2 of 13 thus motivating the use of statistical or machine learning models to predict external loads 37 from subsets of kinematic data [12,13]. 38 Both inertial measurement units and optical motion capture have been used to estimate 39 ground reaction forces during walking, stair-use, and running [ 10–12,14,15]. Joint-level 40 kinetics have been estimated from kinematic data for everyday wheelchair propulsion [ 13]. 41 However, to our knowledge, no study has reported kinematic-based predictions of the 42 external loads during handcycling and is therefore the objective of this study. To mimic data 43 that could be measured via an inertial measurement unit we based our analyses on data 44 from the radius segment. We evaluated two classes of machine learning models (neural 45 network and ensemble) and compared their prediction performance against a statistical 46 regression model. 47 2. Materials and Methods 48 2.1. Dataset and Model inputs 49 Biomechanics data of full-time wheelchair users (n=20) completing a handcycling 50 exercise protocol as part of a separate study [ 9] were used. For each participant, six 51 propulsion cycles were extracted during moderate-intensity (45% peak power) handcycling 52 and six cycles during high-intensity (90% peak power). A custom instrumented handle 53 was used to measure the three forces and three moments applied at the right handcycle 54 handle (sampling frequency = 2 kHz). Kinematic data of the right arm was recorded using 55 a 10-camera motion capture system (Vicon, sampling frequency = 100 Hz). A rigid-body 56 model [ 16] was scaled to each participant and used to determine the radius segment center 57 of mass position and orientation using inverse kinematics (Figure 1). Data from each 58 propulsion cycle was interpolated into 360 data points corresponding to the 360 degrees of 59 crank rotation. 60 Figure 1. Motion-capture markers (pink spheres) used with subject-scaled rigid-body musculoskeletal model of upper arm. Radius segment (red circle) center of mass data derived from the radius position ⃗pand orientation ⃗ewith respect to the global reference frame. The segment positions were numerically differentiated twice using the central differ- 61 ence formula to obtain linear accelerations (Eq. 1):. 62 ⃗a(t) =⃗p(t+τ)−2⃗p(t) +⃗p(t−τ) τ2(1) where ⃗pis the radius center of mass position, tis the current time, and τis the interval 63 between data points. Orientations (Euler angles) were differentiated once to obtain angular 64 velocities (Eq. 2): 65 ⃗ω(t) =⃗e(t+τ)−⃗e(t−τ) 2τ(2) where ⃗eis the vector of Euler angles denoting the radius center of mass in the global 66 coordinate system. 67 The training data for all models was based on 90% of participants (n=18, randomly 68 selected) while data from the remaining two participants were used as the test set. Data for 69 model validation was based on a further split of the training dataset using n=2 participants’ 70 data. Participant weight (lbs), wingspan (in), and exercise intensity were included with 71 Version July 3, 2024 submitted to Sensors 3 of 13 the kinematic data as model input features. Exercise intensity was defined as as 0.5 or 72 0.9 to denote the varying protocol intensities. Prior to model training, each feature was 73 standardized by subtracting the mean and dividing by the standard deviation. This process 74 was used to standardize the training, test, and validation datasets. 75 2.2. Model Descriptions 76 Linearregression: We first assessed the predictive performance of an elastic net linear 77 regression model (scikit-learn [ 17]) herein referred to as "Linear". The Linear model was 78 trained with least absolute shrinkage and selection operator (L1) and Tikhonov (ridge) 79 (L2) regularization of coefficients. The Linear model parameters were α=0.003 ,ρ=0.5, 80 and maximum training iterations = 2000. The model was fit by minimizing the objective 81 function (Eq. 3): 82 min ⃗β1 2nsamples X⃗β−⃗y 2 2+αρ ⃗β +α(1−ρ) 2 ⃗β 2 (3) where ⃗βis the vector of coefficients, Xis the design matrix of predictors, ⃗yis the vector 83 of outputs, αis a constant that multiplies the penalty terms, and ρis the ratio of L1 84 regularization to L2 regularization. 85 Five neural network and two ensemble models (described below) were evaluated. Our 86 choice of models was motivated by the work of others in predicting external loads [12,14] 87 in addition to testing models of varying complexity. Each of the neural network models 88 was compiled using mean squared error as the loss function and the Adam algorithm as the 89 optimizer (learning rate = 0.0001). Models were trained for a maximum of 1000 epochs with 90 an early stop callback, which stopped training early if the loss function did not decrease by 91 0.01 in 50 epochs. A batch size of 16 propulsion cycles was used during training. All neural 92 network models were created using Tensorflow[ 18], the random forest model was created 93 with scikit-learn [ 17], and the gradient boosting model was created with XGBoost [ 19]. All 94 models were implemented using Python (v 3.9) . 95 Dense: A dense neural network (Figure 2A) model was created with one hidden 96 dense layer containing 200 nodes and L1/L2 regularization of 0.00001 applied to the kernel 97 weights matrix. Batch normalization and dropout (45%) were applied after the hidden 98 layer. The output layer was a dense layer with six nodes denoting the six output features. 99 Bidirectional long short -term mem ory(BiLSTM): two layer BiLSTM model (Figure 2B). 100 Each BiLSTM layer had 200 nodes with L1/L2 regularization (0.001) on the kernel weights 101 matrix, recurrent kernel weights matrix, and on the bias vector. After each BiLSTM layer, 102 batch normalization and dropout (45%) were applied. The output layer was a dense layer 103 with six nodes. 104 Version July 3, 2024 submitted to Sensors 4 of 13 Figure 2. Neural network model architectures: (A) Dense (B) Bi-directonal long short-term memory (BiLSTM) Convolutional neuralnetwork (CNN): Two 1D convolutional layers (Figure 3A) used 105 to mimic the structure of the BiLSTM model. Each convolutional layer had 200 filters, 106 a kernel size of 64, causal padding, rectified linear unit (ReLU) activation function, and 107 L1/L2 kernel regularization (0.00001). After each convolutional layer, batch normalization 108 and dropout (50%) were applied. The output layer was a time-distributed dense layer with 109 six nodes. 110 Figure 3. Convolution neural network architectures (A) Convolutional neural network (CNN) (B) Temporal convolutional network (C) Residual neural network (ResNet) Tem poralconvolutional network (TCN): Six hidden layers with 2n−1dilation rate, 111 where nis the layer number (Figure 3B). Each layer had 64 filters with a kernel size of 32 112 and L1/L2 regularization (0.01). After each layer was an ReLU activation layer followed by 113 a dropout layer (40%). The output layer was a time-distributed dense layer with 6 nodes 114 and L1/L2 kernel regularization (0.01). 115 Resid ualneuralnetwork (ResNet): 40 layers of residual blocks (Figure 3C) . Each 116 residual block was made of a 1D convolutional layer starting with 32 filters and a kernel 117 size of 32. The kernel size in the convolutional layer doubled every five layers. The 118 Version July 3, 2024 submitted to Sensors 5 of 13 convolutional layers used ’same’ padding and had L2 regularization (0.001). After each 119 convolutional layer, batch normalization and an ReLU activation layer were used. The 120 output layer was a time-distributed dense layer with 6 nodes and L2 kernel regularization 121 (0.001). 122 Random Forest(RF) and GradientBoost ingModel (GBM): Finally, two ensemble 123 models were evaluated. A random forest regressor [ 17] model was created with 400 124 decision trees, with each tree having a maximum depth of 10 levels. The second ensemble 125 model was a Gradient Boosting Machine (GBM) implemented using XGBoost [ 19]. This 126 model sequentially built 2000 decision trees, each with a maximum depth of 10, using a 127 learning rate of 0.01. It employed both row and column subsampling at 70% and included 128 L1 and L2 regularization (both set to 0.2). The GBM was trained for 1000 boosting rounds. 129 2.3. Statistical Analysis 130 For each of the 24 propulsion cycles in the test set, the median absolute percentage 131 error (MdAPE) was calculated for each kinetic degree of freedom predicted. This analysis 132 was repeated for the eight models tested. The Linear model was used as the baseline for 133 comparison as it was the simplest model tested. 134 The distribution of model errors was non-normal for some degrees of freedom based 135 on Shapiro-Wilkes tests for normality. Therefore, for each degree of freedom, a Kruskal- 136 Wallis test was used to determine whether the errors differed between the models tested. 137 If the differences were significant, a Dunn’s test for multiple comparisons was used to 138 determine which models differed in performance compared to the linear model. 139 The correspondence of the model prediction to the ground truth was also compared 140 qualitatively by evaluating the shape of the predictions and quantitatively using Pearson’s 141 correlation coefficients. Kruskal-Wallis tests were used to determine whether the distribu- 142 tions of correlation coefficients for each degree of freedom differed between models and a 143 post-hoc Dunn’s test was used to determine differences compared to the linear model. 144 Wilcoxon Rank Sum tests were used to compare predictions made by the best per- 145 forming model between the moderate and high-intensity propulsion cycles for each of 146 the output features. Statistical analyses were performed using the SciPy package [ 20] in 147 Python. 148 3. Results 149 Overall, errors in predictions of kinetics in the x and y directions, which contribute 150 to the tangential and radial loads during handcycling, were lower (range: 32 to 56%) than 151 kinetics in the out-of-plane z-direction (range: 45 to 101%). The Dense and RF models had 152 significantly worse prediction performance (errors of 101 and 100%, respectively) than the 153 Linear model for Fz(66% error), and the RF model predictions were also significantly worse 154 forTz. Most of the neural network models performed significantly better (error range: 38 - 155 42 %) than the Linear model (error = 54%) for predicting Fx(Table 1). The ResNet and TCN 156 models performed better than the Linear model for Fy(errors of 38 and 36%, respectively) 157 andTx(errors of 32%). The GBM model predicted better than the linear model for Tyonly 158 (error = 39%). There were no other significant differences in prediction error compared to 159 the Linear model. 160 Version July 3, 2024 submitted to Sensors 6 of 13 Table 1. Comparison of model performance (median absolute percentage error (MdAPE)) across output features. Significance levels shown for models with lower error than the linear regression model. Model TypeModel NameKinetic Prediction Error (MdAPE [%]) Fx Fy Fz Tx Ty Tz Statistical Linear 53.97 (8.02) 51.65 (4.77) 66.36 (9.99) 44.84 (4.42) 51.05 (11.31) 44.86 (12.52) Neural NetworkDense 41.72 (7.16)* 42.44 (24.87) 101.45 (22.24) 40.6 (10.31) 55.74 (9.2) 55.74 (20.42) BiLSTM 45.02 (6.22) 41.97 (12.07) 84.66 (14.89) 35.14 (10.55) 45.11 (7.28) 61.67 (19.5) CNN 39.83 (5.56)*** 43.18 (5.16) 85.64 (12.83) 40.27 (13.55) 45.32 (7.72) 51.53 (9.41) TCN 37.53 (6.26)*** 35.88 (5.17)*** 80.79 (6.66) 31.64 (4.31)** 40.44 (6.9) 46.79 (16.77) ResNet 41.58 (6.1)* 37.83 (7.58)** 88.82 (18.71) 31.96 (8.68)* 41.57 (6.88) 49.23 (8.85) EnsembleRF 47.35 (7.26) 43.35 (18.97) 99.53 (16.05) 41.65 (6.81) 49.45 (8.42) 65.4 (17.05) GBM 42.51 (6.36) 47.71 (6.73) 85.15 (5.72) 35.6 (8.12) 38.9 (7.92)* 53.45 (9.61) ***p<0.001,**p<0.01,*p<0.05. Data are expressed as: Median (Median Absolute Deviation). Of the neural network models, the TCN model had the lowest overall error and the 161 GBM model had the lowest error of the ensemble models. The range of error values for the 162 best statistical, neural network (TCN), and ensemble (GBM)) models were smaller in the x 163 and y-directions compared to the z-direction (Figure 4). 164 Figure 4. Comparison of model prediction median absolute percentage error (MdAPE) for each output feature across the 24 propulsion cycles of the test set. *** p<0.001, ** p<0.01, * p<0.05 The shape of the average prediction from the TCN model matched well with the 165 ground truth for the in-plane kinetic degrees of freedom (Figure 5, A1) with correlation 166 coefficients ranging from 0.95 - 0.97 (Table 2). However, the model tended to underestimate 167 the minima and maxima of each feature. For the x and y degrees of freedom, the errors 168 tended to be highest at the beginning and end of the propulsion cycle. Force and torque 169 profiles in the z-direction were poorly predicted. 170 Version July 3, 2024 submitted to Sensors 7 of 13 Figure 5. Average ground truth and average prediction of a propulsion cycle from the TCN model across each output feature. Standard deviation ribbons are shown. Version July 3, 2024 submitted to Sensors 8 of 13 Table 2. Comparison of model performance (Pearson’s Correlation Coefficient (r)) across output features. Significance levels shown are compared to the linear regression model. Model TypeModel NameOutput Feature Prediction Pearson’s Correlation Coefficient (r) Fx Fy Fz Tx Ty Tz Statistical Linear 0.81 (0.05) 0.86 (0.02) 0.78 (0.06) 0.89 (0.01) 0.84 (0.05) 0.3 (0.09) Neural NetworkDense 0.92 (0.02)*** 0.95 (0.01)*** 0.47 (0.24) 0.94 (0.02)** 0.88 (0.04) 0.41 (0.16) BiLSTM 0.91 (0.03)*** 0.94 (0.02)*** 0.4 (0.14) 0.94 (0.02)*** 0.89 (0.04) 0.3 (0.23) CNN 0.94 (0.02)*** 0.97 (0.01)*** 0.49 (0.21) 0.94 (0.03)*** 0.93 (0.03)*** 0.43 (0.34) TCN 0.97 (0.01)*** 0.95 (0.02)*** 0.56 (0.09) 0.95 (0.01)*** 0.96 (0.02)*** 0.51 (0.2)* ResNet 0.91 (0.04)*** 0.96 (0.01)*** 0.41 (0.41) 0.97 (0.01)*** 0.93 (0.01)*** 0.55 (0.16) EnsembleRF 0.83 (0.06) 0.86 (0.04) -0.01 (0.24) 0.86 (0.03) 0.86 (0.04) 0.14 (0.5) GBM 0.92 (0.03)*** 0.93 (0.02)*** 0.63 (0.1) 0.96 (0.01) 0.93 (0.01)*** 0.54 (0.18) ***p<0.001,**p<0.01,*p<0.05. Data are expressed as: Median (Median Absolute Deviation). Finally, we evaluated whether model performance was sensitive to exercise intensity 171 level. Errors in forces in the y-direction were moderately higher during high-intensity 172 exercise (mean error = 41%) compared to moderate intensity (mean error = 32%) (Figure 6). 173 There were no significant differences in prediction performance between the two exercise 174 intensities for the remaining degrees of freedom. 175 Figure 6. Comparison of TCN model prediction median absolute percentage error (MdAPE) between moderate intensity and high intensity propulsion cycles for each output feature. **: p=0.00323) 4. Discussion 176 We have shown that machine learning models can improve predictions of external 177 loads during handcycling compared to a simpler elastic net linear regression model. Neural 178 network models tended to perform better than ensemble models. Neural network pre- 179 dictions can leverage a longer time range of values, whereas tree-based ensemble models 180 make instantaneous predictions at each timestep. Among the neural network models, 181 convolutional neural networks outperformed recurrent (BiLSTM model) and dense neural 182 network models. 183 Of all models tested, the predictive performance of the TCN was the most accurate 184 for the kinetic parameters Fxand Fy. During handcycling, most force is applied in the 185 tangential and radial directions with respect to the crank center, and is a combination of 186 theFxand Fyforces. In contrast, forces applied in the lateral direction ( Fz) are lower in 187 magnitude. Thus, while the lateral forces and torques had higher percent errors, the error 188 Version July 3, 2024 submitted to Sensors 9 of 13 magnitude was relatively small. The poor ability to predict lateral forces and torques may 189 be due to the variability in propulsion technique (shaded blue regions in Figure 5) [ 9] 190 despite the fact that handcycling is a closed chain task and relatively constrained in terms 191 of motion. Individual technique has been classified for manual wheelchair propulsion but 192 has not been used in handcycling. The inclusion of propulsion style, either as a categorical 193 variable or via the inclusion of additional segment kinematics as inputs, in future machine 194 learning models may prove useful for improving kinetic predictions. Another option 195 for improving model performance is to tune the hyperparameters of the neural network 196 via grid searching though requires significant computational resources. Finally, although 197 there were some differences in prediction performance between moderate intensity and 198 high intensity propulsion cycles, a single model is likely sufficient for multiple exercise 199 intensities provided the intensity data is used as an input to the model. 200 The prediction performance for handcycling loads was similar to other studies that 201 predicted external loads from kinematic data during locomotion. Johnson et al. used a 202 combination of optical motion capture kinematics and inertial measurement unit data to 203 predict the ground reaction forces and moments during various athletic maneuvers and 204 reported Pearson’s correlation coefficients from r=0.55 - 0.90 [ 12], similar to our results 205 (r-range: 0.51 and 0.97) using the TCN model (Table 2). Liu et al. used joint angles from 206 rigid body models to predict normalized ground reaction forces during stair use and 207 reported Pearson’s correlation coefficients between r=0.908 and r=0.991. Normalizing the 208 external loads reduces individual variability, which may increase model predictions, but is 209 challenging to apply to data from persons with spinal cord injury. 210 There are a few limitations of this study. The size of the dataset used was relatively 211 small with 18 independent datasets used for training models. Remarkably, despite this 212 small dataset, our best predictions were similar to others that used data from 80 subjects 213 [14] for legged locomotion. Although the kinematic data used in this study was created to 214 mimic that of a smartwatch-based inertial measurement unit, the data was derived from 215 optical motion capture. Validation of the model predictions using IMU-derived kinematic 216 data is needed. We also restricted our input data to the radius segment of the arm. It is 217 possible that with additional training data from the adjacent segments (upper arm or hand), 218 model performance could be improved. However, our goal was to evaluate the predictive 219 performance of different models using a minimal dataset. 220 The results of this study indicate that a TCN model could be used in place of an instru- 221 mented crank handle in future studies of handcycling propulsion biomechanics. However, 222 it should be noted that the model is best suited for predicting loads that contribute to 223 forward propulsion with more work needed to improve predictions in the lateral direction. 224 Nonetheless, the ability to predict propulsive loads during handcycling now enables future 225 research aimed at providing evidence-based exercise recommendations for people with 226 spinal cord injuries by quantifying the external loads of different protocols. Importantly, 227 the use of wearable-type data combined with machine learning predictions of hand loads 228 allows for the assessment of handcycling outside of the lab environment. 229 Author Contributions: Conceptualization, Griffin C. Sipes, Ian Rice, Mariana E. Kersh; Methodology, 230 Griffin C. Sipes, Mariana E. Kersh; Software, Griffin C. Sipes, Matthew Lee; Validation, Griffin C. 231 Sipes, Matthew Lee; Formal analysis, Griffin C. Sipes; Investigation, Griffin C. Sipes, Matthew Lee; 232 Resources, Griffin C. Sipes, Kellie Halloran; Data curation, Griffin C. Sipes, Matthew Lee; Writing 233 - original draft preparation, Griffin C. Sipes, Matthew Lee; Writing - review and editing, Griffin 234 C. Sipes, Matthew Lee, Kellie Halloran, Ian Rice, Mariana E. Kersh; Visualization, Griffin C. Sipes, 235 Mariana E. Kersh; Supervision, Ian Rice, Mariana E. Kersh; Project administration, Mariana E. Kersh; 236 Funding acquisition, Griffin C. Sipes, Mariana E. Kersh. All authors have read and agreed to the 237 published version of the manuscript. 238 Funding: This research was funded by the National Science Foundation Graduate Research Fellow- 239 ship Program and the De Luca Foundation Seed for Science Initiative. 240 Version July 3, 2024 submitted to Sensors 10 of 13 Data Availability Statement: The raw data supporting the conclusions of this article will be made 241 available by the authors on request. 242 Acknowledgments: We thank Prof. Sanmi Koyejo for helpful guidance regarding best practices for 243 machine learning during the early stages of this study. We are also grateful for the assistance of Dr. 244 Joey Peters and Michael Focht in data collection. 245 Conflicts of Interest: The authors declare no conflicts of interest. 246 Abbreviations 247 The following abbreviations are used in this manuscript: 248 249 CVD Cardiovascular disease PPO Peak power output IMU Inertial measurement unit BiLSTM Bi-directional long short-term memory CNN Convolutional neural network ReLU Rectified linear unit ResNet Residual neural network TCN Temporal convolutional network RF Random forest regressor GBM Gradient boosting machine MdAPE Median absolute percentage error250 Version July 3, 2024 submitted to Sensors 11 of 13 Appendix A 251 Figure A1. Representative (median error) propulsion cycle showing ground truth and TCN model prediction.