A Feature-Encoded Physics-Informed Parameter Identification Neural Network for Musculoskeletal Systems

Abstract Identification of muscle-tendon force generation properties and muscle activities from physiological measurements, e.g., motion data and raw surface electromyography (sEMG), offers opportunities to construct a subject-specific musculoskeletal (MSK) digital twin system for health condition assessment and motion prediction. While machine learning approaches with capabilities in extracting complex features and patterns from a large amount of data have been applied to motion prediction given sEMG signals, the learned data-driven mapping is black-box and may not satisfy the underlying physics and has reduced generality. In this work, we propose a feature-encoded physics-informed parameter identification neural network (FEPI-PINN) for simultaneous prediction of motion and parameter identification of human MSK systems. In this approach, features of high-dimensional noisy sEMG signals are projected onto a low-dimensional noise-filtered embedding space for the enhancement of forwarding dynamics prediction. This FEPI-PINN model can be trained to relate sEMG signals to joint motion and simultaneously identify key MSK parameters. The numerical examples demonstrate that the proposed framework can effectively identify subject-specific muscle parameters and the trained physics-informed forward-dynamics surrogate yields accurate motion predictions of elbow flexion-extension motion that are in good agreement with the measured joint motion data.


Introduction
Human movement resulting from the interaction of various subsystems within the human body is controlled by the excitation signals from the central nervous system that lead to joint motion in the musculoskeletal (MSK) system.These subsystems, governed by parameterized nonlinear differential equations [1,2], form the forward dynamics problem [2].Given information on muscle activations, the joint motion of a subject-specific MSK system can be obtained by solving a forward dynamics problem.The muscle activations can be estimated by the surface electromyography (sEMG) signals through a noninvasive procedure [2,3].Surface electromyography driven forward dynamics have been widely used for predictions of joint kinetics or kinematics [2][3][4][5].For rehabilitation applications in assessing muscular pathologies or weakened muscle groups, Pau et al. [4] used a simplified geometric model of the MSK system for elbow flexion to predict the motion of the elbow joint, given an sEMG signal.Zhao et al. [6] utilized sEMG signals to simulate wrist kinematics for various flexion/extension trials.Standard optimization techniques such as genetic algorithms [4,6], simulated annealing [7], and nonlinear least squares [2,3] have been used for parameter identification.
In recent years, machine learning (ML) or deep-learning-based approaches have become a viable alternative due to their flexibility and capability in extracting complex features and patterns from data [8] and have been successfully applied to various problems in engineering applications such as reduced-order modeling [9][10][11][12][13], and materials modeling [14,15], among others.For motion prediction, data-driven approaches have been introduced to directly map the input sEMG signal to joint kinetics/kinematics, bypassing the forward dynamics equations and the need for parameter estimation.For example, Au et al. [16] used only sEMG data as input to a time-delayed neural network (NN) to predict shoulder motion.Wu et al. [17] applied reinforcement learning for the sEMG and joint moment mapping.Leserri et al. [18] used features of the raw sEMG signals to map them to the elbow joint motion.Ma et al. and Ren et al. [19,20] used enhanced recurrent neural networks and convolutional neural networks, respectively, to map the raw sEMG signal to the motion of the upper limb.While these ML approaches do not require calibrating physiological parameters, the resulting ML-based surrogate models lack interpretability and may not satisfy the underlying physics.
To integrate data and physical models by using ML approaches, data-driven computing that enforces constraints of conservation laws in the learning algorithms of a material database has been developed in the field of computational mechanics [21][22][23].This paradigm has been applied to other engineering problems, such as nonlinear material modeling [22,24,25] and, fracture mechanics [26], among others.Furthermore, deep manifold embedding techniques have been introduced in data-driven computing for extracting low-dimensional feature space [27,28].
In this study, we propose a physics-informed parameter identification neural network (PI-PINN) for the simultaneous prediction of motion and parameter identification with application to MSK systems.Using the raw transient sEMG signals obtained from the sensors and the corresponding joint motion data, the PI-PINN learns to predict the motion and identifies the parameters of the hill-type muscle models representing the contractile muscletendon complex.
The standard PINNs with a fully connected architecture present difficulties in learning the high-frequency components of the solution, which is known as spectral bias [32,42,43].Wong et al. [43] pointed out that PINNs tend to have a strong bias toward lowfrequency patterns which are considered trivial solutions to the governing equations.To mitigate the issue of spectral bias, in this work a Fourier feature transformation is applied to modulate the signals input to embedding space [42,43].A similar approach was proposed in Ref. [37] where a feature layer has been introduced for encoding certain dynamic patterns, such as periodicity.Inspired by these studies, we further propose a feature-encoded approach for the PI-PINN framework, termed feature-encoded physics-informed parameter identification neural network (FEPI-PINN), in dealing with the highly oscillatory sEMG signals.Here, the sEMG signals and joint motion data are first projected onto a low-dimensional space consisting of Fourier and polynomial bases.The NN is then used to map the associated basis coefficients between the input signal and the target output signal.Subsequently, the mapped coefficients are used to reconstruct the joint motion using the bases.
This paper is organized as follows.Section 2 introduces the subsystems and mathematical formulations of MSK forward dynamics, including the EMG-to-activation dynamics, the muscle-tendon contraction dynamics, and the MSK system dynamics, followed by an introduction of the proposed FEPI-PINN framework for simultaneous motion prediction and system parameter identification.Section 3 discusses the verification of the proposed algorithms using synthetic motion data.In Sec. 4, the proposed frameworks are validated by modeling the elbow flexion-extension movement using subject-specific sEMG signals and recorded motion data.Concluding remarks and future work are summarized in Sec. 5.

Methods
In this section, the mathematical formulations of subsystems for the forward dynamics of the human MSK system are described, followed by an introduction of the proposed PI-PINN framework designed for simultaneous forward dynamics learning and parameter identification.
Fig. 1 A flowchart depicting the interaction of the different subsystems related to the motion of the MSK system.Excitation from neurons is transmitted to muscle fibers (activation dynamics) that contract to produce force (muscle-tendon contraction dynamics).These forces generate torques at the joints (structural level MSK dynamics), leading to joint motion.The boxed entities are the nonlinear differential equations that relate to the different states of the system.
2.1 Musculoskeletal Forward Dynamics.The hierarchical interaction of various subsystems of MSK forward dynamics is illustrated in Fig. 1, where the activation dynamics transform neural excitation (which can be measured by sEMG signals) to muscle activation that drives muscle fibers to produce force through the muscle-tendon (MT) contraction dynamics, leading to joint motion (translation and rotation) of MSK systems via the MSK system dynamics.The mathematical formulations of these subsystems are introduced in Secs.2.1.1-2.1.3.

2.
1.1 EMG-to-Activation Dynamics.The raw sEMG signals eðtÞ are measures of the neural excitation u t ð Þ from the central nervous system.The EMG-to-Activation dynamics describes a nonlinear transformation from sEMG signals to muscle activations a t ð Þ [4,6,44], which activate muscle fibers in a muscle group to produce active force.The neural excitation signal uðtÞ can be transformed from the sEMG signals eðtÞ by [2,3,45] where d is the electromechanical delay between the neural excitation originating and reaching the muscle group [45].The muscle activation signal aðtÞ can then be obtained from the neural excitation u t ð Þ in Eq. ( 1) by where A is a shape factor [45].
2.1.2Muscle-Tendon Contraction Dynamics.The MT contraction dynamics is introduced to relate the muscle contraction to force production via a hill-type muscle model [46,47] parameterized by the maximum isometric force in the muscle (f M 0 ), the optimal muscle length (l M 0 ) corresponding to the maximum isometric force, the maximum contraction velocity (v M max ), the slack length of the tendon (l T s Þ, and the pennation angle (/Þ to be discussed as follows.
First, let each muscle group be characterized by a parameter vector The force generated by the muscle group contains an active and a passive component [48,49].The active force f A component can be expressed as where a is the activation function in Eq. ( 2), l M is the normalized muscle length, ṽM is the normalized velocity of the muscle, and f A;L l M ð Þ and f V ṽM ð Þ are functions of the length and velocitydependent force generation properties of the active muscle represented by generic functions of dimensionless quantities as given in Appendix B [48][49][50], allowing them to be scaled to specific muscles through the parameters described.The total muscle force F M can be expressed as where f P lM ð Þ is the passive muscle length-dependent force generation function with the specific form given in Appendix B.
The muscle force F M obtained in Eq. ( 5) is transmitted to the joints through the tendon (Fig. 2(b)).The tendon produces force F T only when its length l T is stretched beyond its slack length l T s .In this study, the tendon is assumed to be rigid [6,47] and thus the tendon length l T ¼ l T s is adopted.According to the schematics of the MT complex shown in Fig. 2(b), the total length of MT complex, l MT , is first obtained by where / is the pennation angle.The total force produced by the MT complex, F MT , can be expressed as follows based on the force equilibrium: The rigid-tendon model simplifies the MT contraction dynamics [46,47,50] which accounts for the interaction of the activation, force-length, and force-velocity properties of the MT complex.More details on how muscle length and velocity are calculated can be found in Appendix A.

Musculoskeletal System
Dynamics.The motion of the MSK system is modeled as the rotational movement of rigid links (bones) with force producing actuators (MT complex's).The force produced by these actuators is converted to torques at the joints of the body which ultimately leads to movement.This torque equilibrium can be expressed as I q ð Þ € q À T MT ða; q; _ q; jÞ À EðqÞ ¼ 0 (8) where q; _ q; € q are the vectors of generalized angular motions, angular velocities, and angular accelerations, respectively; EðqÞ is the torque from the external forces acting on the MSK system, e.g., ground reactions, gravitational loads, etc.; I q ð Þ is the inertial matrix; T MT is the torque from all muscles in the model calculated by T MT ða; q; _ q; jÞ ¼ R q ð Þ F MT ða; q; _ q; jÞ, where R q ð Þ are the moment arm's and F MT ða; q; _ q; jÞ is obtained from Eq. ( 7).As the muscle lengths l M and velocities v M are functions of the joint motion q and _ q, and the total force produced by all the MT complex's in the model is represented as F MT ða; q; _ q; jÞ.Given the muscle activation signals from Eq. ( 2) and parameters of involved muscle groups, the generalized angular motions q and angular velocities _ q of the joints can be obtained by solving Eq. ( 8) [50].

Simultaneous Forward Dynamics
Learning and Parameter Identification.The proposed PI-PINN framework for simultaneous prediction of motion and parameter identification of human MSK systems is first introduced in Sec.2.2.1, where the MSK forward dynamics are learned by an NN surrogate that predicts motion for identification of MSK properties using timedomain training.The enhanced forward dynamics surrogate is then introduced in Sec.2.2.2, with the feature-encoded training.
With the governing equations for the MSK forward dynamics introduced in Sec.2.1, the following parameterized ordinary differential equation (ODE) system is defined as where the differential operator L½Á; k is parameterized by a set of parameters k ¼ fk 0 ; k 1 ; …; k n g.The right-hand side b t; x ð Þ is parameterized by x ¼ ½x 1 ; x 2 ; …; x m .B½Á is the operator for initial conditions, and g is the vector of prescribed initial conditions.To simplify notations, the ODE parameters are denoted by C ¼ fk; xg.The solution to the ODE system q : ½0; T !R depends on the choice of parameters C.

Forward Dynamics NN Surrogate With Time-Domain
Training.Here, a multilayer NN is used to relate data inputs containing discrete sEMG signals and discrete time, x 2 R nin , to discrete joint motion data outputs, q 2 R nout , by a nonlinear activation function h Á ð Þ as follows: where M is the number of hidden layers; the superscript ðkÞ denotes the layer number.Here, the activation R nk are trainable weight and bias coefficients, respectively; h ¼ fW; bg is used to denote all trainable parameters of the network, where are the set of all weights and biases, respectively.A hyperbolic tangent activation function h Á ð Þ is used for the hidden layers, where the activation is applied to all components of its input vector.As the application of interest is a regression task, a linear activation is used for the output layer, which means the output of the last hidden layer is mapped to the prediction of the output layer where q is the approximation of the output training data q.
An NN approximates the MSK forward dynamics, which predicts MSK motion given the time history of raw sEMG signals of muscle groups.The raw sEMG signals inform the motion prediction with all transient information of sensor measured sEMG signals without any processing.Since we use the signals in their time-domain, we denote this as time-domain training.
Let us consider the k th motion trial out of p trials of the training data of the MSK forward dynamics with n k data points or timesteps.For the surrogate, the input data of the k th trial is T and the target output as q ¼ q 1 ; …; q p b c T .The NN is trained to learn the mapping from the raw sEMG and time The trainable parameters of the NN are obtained by minimizing the following loss function: where k Á k 2 L2 denotes the L2 norm.In addition to training an MSK forward dynamics surrogate, the proposed framework aims to simultaneously identify important MSK parameters from the training data by minimizing residual of the governing equation of MSK system dynamics in Eq. ( 8) denotes the vector of residuals associated with Eq. ( 8) for the i th sample; qi x i ; h ð Þis the vector of predicted MSK joint motion from the NN with the trainable parameters h for the i th sample; C ¼ fk; xg represents the set of ODE parameters relevant to the MSK system.The optimal NN parameters h and the ODE parameters C are obtained by minimizing the composite loss function J as follows: where b is the parameter to regularize the loss contribution from the ODE residual term in the loss function and is estimated analytically during the training [51].For this study, b was scaled so that the terms J res and J data in the loss function given in Eq. ( 14) were with the same unit.The details are described in Sec. 3. The gradients of the NN outputs with respect to the NN MSK parameters h; C ð Þ can be obtained efficiently by automatic differentiation [52].The computational graph is illustrated in Fig. 3.

Fourier features.
It has been pointed out in the literature that training PINNs can be challenging as the solution could be trapped in a local minima, leading to large training data errors [42,43,51].Wang et al. [42,51] showed that these networks tend to learn low-frequency components of the solution rather than the high-frequency counterparts, which is described as the spectral bias of PINNs.Wong et al. [43] further showed that PINNs are biased toward flat outputs which are considered as trivial solutions to the governing equations.It has been shown that employing Fourier features enables PINNs to learn the higher frequency components and avoid local minima, leading to improved learning performance [42].PINNs with Fourier features have also been applied to image regression tasks in biomedical applications with high-dimensional data sources as magnetic resonance imaging [53].Due to the oscillatory nature of the sEMG signals in the time-domain, the mapping to the MSK motion could contain features of various frequencies.Here, PINNs with Fourier features are employed in the forward dynamics surrogate with timedomain training.
A Fourier feature x F 2 R m of a given input x 2 R d is defined as [42,43] x where W F 2 R mÂd and b F 2 R m are weight and bias coefficients, respectively, randomly sampled from a normal distribution N 0; r 2 À Á , with r as a tunable hyperparameter and sinðÁÞ is the sinusoidal activation function.They can be easily added to the computational graph of PINN, as shown in Fig. 3.The equations for forward propagation remain the same with an additional Fourier feature transformation before the input layer.The optimization problem to obtain the optimal trainable network parameters h and the ODE parameters C remains the same as described in Eq. ( 14).

Forward Dynamics NN Surrogate With Feature-Encoded Training.
Mapping oscillations in the input signal of the network to a smooth less-oscillatory output signal is susceptible to spectral bias.Given the high degree of oscillations in the transient raw sEMG signals, it was observed in the tests with time-domain training (Sec.4) that the motion predictions were affected by artifacts of those oscillations.To this end, the framework was enhanced by encoding the features of the raw sEMG signals and motion signals to a low-dimensional space with Fourier and polynomial basis functions.This dimensionality reduction enables noise filtering and provides a smooth low-dimensional representation of the high-dimensional noisy data.A forward dynamics NN surrogate is then trained in the low-dimensional feature embedding space to perform mapping and simultaneously identify muscle parameters in the MSK system.
Due to the periodic nature of the signals, Fourier basis was used along with polynomials up to quadratic terms in the featureencoded training.The polynomial basis helps capture the signals close to the boundary of the time-domain that have incomplete periodicity.For the j th sEMG signal from the k th trial, , consider the employment of 2n EMG Fourier basis terms with their respective coefficients and three coefficients for the complete quadratic basis terms, i.e., The set of all input coefficients is denoted by I k;j , defined as and the encoded approximation of the j th sEMG signal from the k th trial is given by e FE k;j t; where T e k;j is the duration of the sEMG signal.Applying the encoding to all sEMG signals of the k th trial gives a set of input coefficients, i.e., I k ¼ I k;j f g Na j¼1 .Correspondingly, there would be 2n q þ 3 coefficients for the encoded approximation of the k th motion signal with the Fourier where T q k is the duration of the motion signal.The encoded coefficients of the approximations can be obtained by a least-squares minimization between the approximate and the original signals.
An NN with M hidden layers is then trained to predict the output coefficients Ôk by mapping the input coefficients I k to the target output coefficients O k as follows: where and hðÁÞ are the set of weights, set of biases, and activation functions (hyperbolic tangent) of the NN, respectively.Note that rather than learning the mapping related to transformed input features [37,42,43] in the time-domain, this feature-encoding approach is designed to learn the mapping of the basis coefficients, which can reduce mapping complexity and improve prediction accuracy as shown in our validation example.
Given the predicted output coefficients O k , the predictions of angular motion, velocity, and acceleration are obtained from Eq. ( 19) and its time derivatives.The loss function in Eq. ( 14) is modified as follows.The computational graph for the featureencoded PI-PINN, termed FEPI-PINN, is shown in Fig. 4 h; C ¼ argmin Journal of Biomechanical Engineering DECEMBER 2022, Vol.144 / 121006-5 2.3 Sensor Data Acquisition.The subject performs three trials of varying elbow flexion and extension motions as shown in Fig. 5, with each trial lasting 10 s.Three retroreflective markers with a diameter of 14 mm were attached to the acromion (shoulder), the humeral lateral epicondyle (elbow), and the radial styloid (wrist) of the right arm of the subject.Two delsys trigno surface sEMG sensors were affixed onto the biceps and triceps muscle groups, following surface electromyography for the noninvasive assessment of muscles project recommendations [54].To minimize signal artifacts, the subject's arm was wiped with isopropyl alcohol before sEMG sensors were attached.The marker positions during the duration of the motion were captured by a Vicon motion capture system (Vicon Motion Systems Ltd., UK) at 150 Hz [55], and the sEMG signals were recorded at 2250 Hz.
The kinematic data and the sEMG data were time-synchronized using the Vicon system via a trigger module.To obtain the muscle activations needed to calculate the MSK ODE residual, the raw sEMG signals were first centered with respect to their medians to center them around zero.The centered sEMG signals were processed using a Hampel filter to remove outliers resulting from skin artifacts, before going through a full-wave rectification.The rectified signals then go through a second-order Butterworth filter with a cutoff frequency of 3 Hz.Finally, the maximum voluntary contraction for each muscle was then taken as the maximum voltage from all trials for each muscle.The centered, rectified, and filtered sEMG signals were then normalized by the maximum voluntary contraction voltage to complete processing.They were then passed to the EMG-to-activation dynamics model for each muscle group to get the muscle activations as described in Sec.2.1.The marker trajectories were directly used as is without any processing to estimate the elbow angle.

Verification: Elbow Flexion-Extension Motion
To verify the proposed PI-PINN framework, an elbow flexionextension model was considered and the flowchart of the proposed Fig. 5 Overview of the application of this framework to the recorded motion data.The location of the motion capture markers is circled in red and the EMG sensors in blue.The simplified rigid body model is used in the forward dynamics equations within the framework with appropriately scaled anthropometric properties (for geometry) and physiological parameters (for muscle tendon material models).The raw sEMG signals are mapped to the target angular motion of the elbow and used to simultaneously characterize the MSK system using the proposed frameworks for PI-PINN.Transactions of the ASME PI-PINN computational framework for simultaneous forward dynamics prediction and parameter identification of MSK parameters is shown in Fig. 5. Synthetic sEMG signals and associated motion responses were considered and the time-domain PI-PINN approach was examined.The FEPI-PINN approach was employed in dealing with the recorded oscillatory sEMG data in the validation problem in Sec. 4. The model contained two rigid links resembling the upper arm and forearm with lengths l ua and l fa , respectively.Both rigid links were connected at a hinge resembling the elbow joint "A," and the upper arm link was fixed at the top joint "B," while the lines connecting the two links represented the biceps and triceps muscle-tendon complex, as shown in Fig. 2(b).The degree-offreedom of the model was the elbow flexion angle q.The biceps and triceps muscle-tendon complex assemblies were modeled by two hill-type models with parameters j Bi and j Tri , as discussed in Sec.2.1.It was assumed that the mass in the forehand is concentrated at the wrist location, hence, a mass m fa was attached to one end of the forearm link with a moment arm l fa from the elbow joint.We assumed that the tendons were rigid as discussed in Ref. [47] for ease of computation.The equation of motion for this rigid body system was where T MT a Bi ; a Tri ; q; _ q; j Bi ; T MT Bi a Bi ; q; _ q; j Bi ð Þ¼ with the initial conditions q 0 ð Þ ¼ p 6 radians and _ q 0 ð Þ ¼ 0 radians=sec.Given the synthetic sEMG signals (e Bi ; e Tri ) that were plugged into the sEMG-to-activation dynamics equations (Sec.2.1.1),leading to muscle activations a Bi ; a tri and using the parameters summarized in Table 1, the motion of the elbow joint, q, was obtained by solving the MSK forward dynamics problem using a synthetic solver.
Five synthetic samples of the elbow flexion-extension motion were generated by solving Eq. ( 24) given synthetic muscle sEMG signals, as shown in Fig. 6.The training dataset contained the data of trials 1, 2, 4, and 5, while the data of trial 3 was used for testing, each trial with n ¼500 data points.The data of the k th sample with n temporal data points contained ½x k ; The set of MSK parameters o were considered to be identified from the training data using the proposed framework.Due to different units and physiological nature of the parameters, they could be very different in terms of scale.For example, there could be a difference of more than two orders of magnitude between f M 0 and v M max , which could affect the conditioning of the parameter identification system.To mitigate this issue, normalization [40] was applied to each of the parameters where was the initial value of the parameter.Therefore, the set of parameters to be identified became The proposed PI-PINN framework, as described in Sec.2.2, was applied to simultaneously learn the MSK forward dynamics surrogate and identify the MSK parameters C by optimizing Eq. ( 14), where the residual of the governing equation for the training input x i ð Þ k , was expressed as and could be plugged into the residual term J res in the loss function in Eq. ( 14).
The residual of the equation of motion ð J res Þ involved in the loss function J (Eq. ( 14)) was scaled so that the terms J data and bJ res in Eq. ( 14) had the same unit.This yielded the scaling parameter b / Dt 2 I , where Dt is the time-step size in the timesignal and I is the moment of inertia of the mass around the elbow.
A three-layer feedforward NN with 32 neurons in each layer was used.The training was performed using the Adam algorithm [56] with an initial learning rate of 5Â10 À3 .In this example, b / Dt 2 I ¼ 10 À4 .The trained forward dynamics surrogate could accurately predict the joint kinematics given an input signal from the testing data, as shown in Fig. 7(b).Meanwhile, the MSK parameters, f M 0 and v M max , of both the biceps and the triceps were accurately identified from the motion data, with an error of less than 1%, as shown in Fig. 8.A sensitivity study of the results of this verification problem with respect to b is shown in Appendix C, where the study showed that this estimate of b / Dt 2 I ¼ 10 À4 lies within a range leading to optimal motion prediction and parameter identification from the framework.

Validation: Elbow Flexion-Extension Motion
The recorded motion data and sEMG signals were collected as per the data acquisition protocols mentioned in Sec.2.3.The subject performed three trials of the elbow flexion-extension motion, where the sEMG sensors were fixed on two muscle groups, namely, the biceps and triceps.The sEMG signals were processed as described in Sec.2.1 to obtain muscle activation signals for each muscle group that was used to calculate the MSK ODE residual.We used the same simplified rigid body model as in Sec. 3 and appropriately scaled the anthropometric properties (for the geometry of the model) and physiological parameters (for muscletendon material models used for the muscle groups) based on the generic upper body model defined in Refs.[50] and [57].Figure 9 shows the measured data of the three trials, including the transient raw sEMG signals and the corresponding angular motion of the elbow flexion-extension of the subject.
In this example, the raw sEMG signals were used as input.The time-domain and the feature-encoded PI-PINN frameworks were examined and compared.The data of trials 1 and 3 were used for training, while trial 2 was used for validation, where each signal contained 1000 temporal data points.The normalization described in Eqs. ( 30) and (31) were adopted.The set of muscle parameters to be identified by the framework include the maximum isometric force and the maximum contraction velocity from both muscle groups, which are denoted as The training was performed using the Adam algorithm [56] with an initial learning rate of 5Â10 À3 .
4.1 Time-Domain PI-PINN Approach.For the time-domain training, the training data of the k th sample with n k data points contained ½x k ; q k ¼ ½t k ; e k;Bi ; e k;Tri ; q k .The entire training dataset with a total of N data ¼ n 1 þ n 3 data points could then be defined as follows: the input discrete time and sEMG signals were , and e Tri ¼ e T 1;Tri ; e T 3;Tri h i T .The output was the corresponding angular motion data q ¼ ½q T 1 ; q T 3 T .An NN with three hidden layers, with 64 neurons in each layer, was used.The Fourier features used in the time-domain training were controlled by parameter r ¼ 1.
Figure 10 compares the data of joint kinematics with the predictions from the trained forward dynamics surrogate with time- Figure 11 shows that a good agreement between the predictions and motion data was achieved for both training and testing cases using the feature encoding method.It demonstrated that the trained forward dynamics surrogate could learn the MSK forward dynamics and effectively predict the joint kinematics given a sEMG signal of the same type of motion used for training.Compared with the results obtained from the time-domain training, as shown in Fig. 10, the feature-encoded training showed a significant improvement in terms of prediction accuracy, which was attributed to the dimensionality reduction that filters out unnecessary noise and oscillations of original high-dimensional data.Figures 10(c) and 10(d) and Figs.11(c), and 11(d) compare the data with the predictions of both methods in terms of angular motion q and velocity _ q.Using exact derivatives in the featureencoded training also led to smoother _ q values as compared to the ones obtained through the NN in time-domain training.It is noted that the motion prediction on the testing case was with a slight error, which could be attributed to the small size of the training dataset (i.e., using two trials).Overall, the feature-encoded  2 where the error statistics are reported.The testing predictions from both training methods were compared to the true data and lower mean squared error and R 2 statistic close to 1 were seen in the case of the feature encoding.Here, the mean squared error (MSE) and R 2 statistic are defined as where is the prediction from the PI-PINN framework, and q is the mean of trial 2 s motion data.

Parameter Identification.
The identified MSK parameters from the FEPI-PINN that yielded higher accuracy of motion predictions are summarized in Table 3. Figure 12 shows the evolution of the MSK parameters, f M 0 and v M max , of both the biceps and the triceps during optimization, with the final converged values of f M 0 lying within the physiological estimates of these parameters reported in the literature [57][58][59].Note that since it is difficult to measure the parameter v M max in experiments, an estimate of v M max , 8 À 10 l M 0 =sec; was typically used in hill-type muscle models [50].
Given that l M 0;Bi ¼ 0:1157 m and l M 0;Tri ¼ 0:1138 m were used in this model, the physiologically reasonable estimated range of v M max;Bi and v M max;Tri were 1.16-1.32m=s and 1.14-1.34m=s, respectively.The v M max;Bi and v M max;Tri identified by the proposed framework were 1.24 m=s and 1.138 m=s, respectively, as summarized in Table 3, which shows good agreement.The results

Discussion and Conclusions
In this work, a physics-informed parameter identification neural network (PI-PINN) framework, as well as its feature-encoded version (FEPI-PINN), was proposed for motion prediction and parameter identification of MSK systems.In this approach, the recorded marker data for positional information and raw signals from the sEMG sensors (for estimation of the activity of the two muscle groups involved in the motion) were utilized for training.The framework was built using a neural network that learns the mapping between the raw sEMG signals and joint kinematics,   3.
minimizing a loss function that consists of the error in the training data and the residual of the MSK forward dynamics equilibrium.Under this PI-PINN framework, time-domain training was first investigated.In this approach, sEMG and motion signals were mapped for their entire time duration, and Fourier features were used to learn the relationship between them as a forward dynamics surrogate.Next, an alternative FEPI-PINN approach was proposed to enhance the model performance by introducing a feature transformation encoder to project the noisy and oscillatory raw data onto a low-dimensional feature embedding space.In this featureencoded PI-PINN (i.e., FEPI-PINN) framework, the forward dynamics surrogate relates basis coefficients of inputs to those of target outputs encoded in noise-filtered embedding space, yielding an enhanced prediction accuracy.
The time-domain method was first verified with synthetic data where PI-PINN was trained to learn the motion of elbow flexionextension given synthetic sEMG signals.Simultaneously, the parameters related to the maximum isometric force and maximum contraction velocity in the hill-type models of the biceps and triceps muscle groups were also identified.In the validation problem, the elbow flexion-extension model was scaled according to the subject's anthropometric geometry.The recorded raw sEMG signal was employed as input to predict the angular motion.It was found that the time-domain PI-PINN training led to nonphysiological motion predictions due to the noise and oscillations in the raw sEMG signal that increases mapping complexities.FEPI-PINN, on the other hand, yielded smoother and physiologically accurate mappings of the angular motion and angular velocity.It was noticeable that, despite the limited size of the dataset, the FEPI-PINN approach outperformed the time-domain PI-PINN approach.A reasonable agreement was found between the identified parameters and the range of their physiological parameter estimates.
The applications of this technique are aimed at assistive devices or exoskeletons [20] that rely on guiding the motion of the subject based on control signals such as sEMG that inform muscle activity.Accurate prediction of the intended motion of the subject, based on their muscle activity is crucial for the successful implementation of these technologies.While in the earlier studies, methods have been proposed to map the sEMG signals to the joint motion by using deep learning techniques such as convolutional neural networks, recurrent neural networks, and auto-encoders [19,60,61], these black-box ML-based models do not satisfy the underlying physics and often lack interpretability.The proposed FEPI-PINN framework, on the other hand, provides a systematic approach to enforce the underlying physics into the mapping, improving the interpretability of ML-based motion.This study also suggested the potential extension of the FEPI-PINN framework for MSK digital twin applications.Under this framework, upon completion of model training with sufficient data from the subject(s), the characterized MSK and NN parameters provide a real-time model for joint motion prediction using sEMG signals without solving the nonlinear ODEs, as long as the sEMG signals are within the range of the training data.The characterized MSK physiology parameters also offer the potential diagnosis for muscle injury and disease development.
While a simplified MSK model was used to model the elbow flexion-extension movement, a more physiologically accurate representation of muscle tissues with fats, connective tissues, and muscle fibers could be used to inform the length and velocitydependent force generation capacity of the muscles [62] in future work.To incorporate different types of motions and account for subject variability, transfer learning [63] could be applied to carry over the features from the previously learned motions for model enhancement.and ODE residual l M ¼ normalized muscle length l MT ¼ total length of muscle-tendon complex l M 0 ¼ optimal muscle length corresponding to the maximum isometric force l T s ¼ slack length of the tendon l fa ¼ forearm length l ua ¼ upper arm length m fa ¼ mass of forearm O k ¼ set of all output feature-encoded basis coefficients of the k th trial Ôk ¼ predicted output feature-encoded basis coefficients of the k th trial q ¼ generalized angular motion of the MSK system qk ¼ predicted angular motion of the k th trial from the PI-PINN framework T MT ¼ torque produced by the muscle-tendon complex ṽM ¼ normalized muscle velocity v M max ¼ maximum contraction velocity u ¼ neural excitations x k ¼ input data from the k th motion trial x F ¼ Fourier feature vector C i ¼ i th parameter characterizing the parameterized ODE system h ¼ set of weights and biases of the neural network j ¼ vector of muscle parameters for each muscle group / ¼ pennation angle between muscle and tendon

Fig. 2 A
Fig. 2 A representation of the muscle-tendon complex in the arm modeled by homogenized hill-type models.Each dark grey line in (a) is a homogenized MT complex described by the model shown in (b).
denote the time and the MSK joint motion at the i th time-step in the k th trial, respectively, and e i ð Þ k;j n o Na j¼1 denotes the set of the raw sEMG signals of N a muscle groups involved in the MSK joint motion q i ð Þ k at the i th time-step.The entire training dataset has N data ¼ P p i¼1 n i data points with the input as x ¼ x 1 ; …; x p b c

Fig. 3 A
Fig. 3 A computational graph of time-domain PI-PINNs with Fourier features prefixed to the network

Fig. 4
Fig. 4 The computational graph for the FEPI-PINN framework

Table 1 Fig. 6
Fig. 6 The dataset with synthetic bicep and tricep sEMG signals with variations in frequency for all trials, and the corresponding motion from the two

Fig. 7
Fig. 7 Comparison between the data of joint kinematics and the predictions from the trained forward dynamics surrogate: (a) the training cases and (b) the testing case

Fig. 9 Fig. 8
Fig. 9 The measured transient raw sEMG signals and the corresponding angular motion of the elbow flexion-extension of the subject

Fig. 10
Fig. 10 Comparison of the angular motion data with the predictions from the time-domain PI-PINN: (a) training case: trial 1; (b) training case: trial 3; (c) testing case: trial 2. The comparison of angular velocity of the testing case is shown in (d).

Fig. 11
Fig. 11 Comparison of the angular motion data with the predictions from the FEPI-PINN: (a) training case: trial 1; (b) training case: trial 3; (c) testing case: trial 2. The comparison of angular velocity of the testing case is shown in (d).

Fig. 12
Fig. 12 Evolution of the identified MSK parameters: (a) f M 0 and (b) v M max of bicep and triceps during the training process.It shows the identified parameters converge to estimates that are in the physiological range of those parameters, as summarized in Table3.

Nomenclaturea
¼ muscle activations A ¼ shape factor used to relate neural excitation to muscle activation d ¼ electromechanical delay between origin of neural excitation from central nervous system and reaching the muscle group e ¼ raw sEMG signals captured by sensors f A ¼ active force component of hill-type muscle model f A;L ¼ length-dependent active force generation component f P ¼ passive force component of hill-type muscle model f V ¼ velocity-dependent active force generation component f M 0 ¼ maximum isometric force in the muscle F M ¼ total muscle force F MT ¼ total force produced by muscle-tendon complex I k;j ¼ set of all input feature-encoded basis coefficients of the j th sEMG signal from the k th trial J ¼ composite loss function of the PI-PINN minimizing data

Table 2
Statistics comparing the testing predictions of the framework with original time-domain training versus the enhanced feature-encoded training 121006-10 / Vol.144, DECEMBER 2022 Transactions of the ASME demonstrated the effectiveness of the proposed FEPI-PINN framework and promising potential for real applications.