## 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.

## 1 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–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–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–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].

More recently, physics-informed neural networks (PINNs) have been developed [29–32] to approximate the solutions of given physical equations by using feed-forward NNs and minimizing the residuals of the governing partial differential equations (PDEs) and the associated initial and boundary conditions. The PINN method has been successfully applied to problems such as flow and transport in porous media [31], solid mechanics [30], additive manufacturing [33], biomechanics, and biomedical applications [34–37], and inverse problems [29,30,38–40]. In inverse modeling with PINNs, the unknown system characteristics are considered trainable parameters or functions [29,41].

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 muscle-tendon 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 low-frequency 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.

## 2 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.

### 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.

where $A$ is a shape factor [45].

#### 2.1.2 Muscle-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 ($f0M$), the optimal muscle length ($l0M$) corresponding to the maximum isometric force, the maximum contraction velocity ($vmaxM$), the slack length of the tendon ($lsT)$, and the pennation angle ($\varphi )$ to be discussed as follows.

where $fP(l\u0303M)$ is the passive muscle length-dependent force generation function with the specific form given in Appendix B.

#### 2.1.3 Musculoskeletal System Dynamics.

where $q,\u2009q\u02d9,\u2009q\xa8$ 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; $TMT$ is the torque from all muscles in the model calculated by $TMT(a,q,q\u02d9;\kappa )=R(q)FMT(a,q,q\u02d9;\kappa )$, where $R(q)$ are the moment arm's and $FMT(a,q,q\u02d9;\kappa )$ is obtained from Eq. (7). As the muscle lengths $lM$ and velocities $vM$ are functions of the joint motion $q$ and $q\u02d9$, and the total force produced by all the MT complex's in the model is represented as $FMT(a,q,q\u02d9;\kappa )$. Given the muscle activation signals from Eq. (2) and parameters of involved muscle groups, the generalized angular motions $q$ and angular velocities $q\u02d9$ of the joints can be obtained by solving Eq. (8) [50].

### 2.2 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 time-domain training. The enhanced forward dynamics surrogate is then introduced in Sec. 2.2.2, with the feature-encoded training.

where the differential operator $L[\u22c5;\lambda ]$ is parameterized by a set of parameters $\lambda ={\lambda 0,\lambda 1,\u2026,\lambda n}$. The right-hand side $b(t;\omega )$ is parameterized by $\omega =[\omega 1,\omega 2,\u2026,\omega m]$. $B[\u22c5]$ is the operator for initial conditions, and $g$ is the vector of prescribed initial conditions. To simplify notations, the ODE parameters are denoted by $\Gamma ={\lambda ,\omega}$. The solution to the ODE system $q:[0,T]\u2192R$ depends on the choice of parameters $\Gamma $.

#### 2.2.1 Forward Dynamics NN Surrogate With Time-Domain Training.

where $q\u0302$ 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 $kth$ motion trial out of *p* trials of the training data of the MSK forward dynamics with $nk$ data points or time-steps. For the surrogate, the input data of the $kth$ trial is $xk=[xk(1)T,\u2026,\u2009xk(nk)T]T$ with the $ith$ row as $xk(i)T=[tk(i),ek,1(i),\u2026,ek,Na(i)]$ while the output data is $qk=[qk(1),\u2009\u2026,\u2009qk(nk)]T$. Here, $tk(i)$ and $qk(i)$ denote the time and the MSK joint motion at the $ith$ time-step in the $kth$ trial, respectively, and ${ek,j(i)}j=1Na$ denotes the set of the raw sEMG signals of $Na$ muscle groups involved in the MSK joint motion $qk(i)$ at the $ith$ time-step. The entire training dataset has $Ndata=\u2211i=1pni$ data points with the input as $x=x1,\u2009\u2026,\u2009xpT$ and the target output as $q=q1,\u2009\u2026,\u2009qpT$.

where $\beta $ 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, $\beta $ was scaled so that the terms $Jres$ and $\u2009Jdata$ 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 $(\theta ,\Gamma )$ can be obtained efficiently by automatic differentiation [52]. The computational graph is illustrated in Fig. 3.

##### 2.2.1.1 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 time-domain training.

where $WF\u2208Rm\xd7d$ and $bF\u2208Rm$ are weight and bias coefficients, respectively, randomly sampled from a normal distribution $N(0,\sigma 2)$, with $\sigma $ as a tunable hyperparameter and $sin(\u22c5)$ 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 $\theta \u0303$ and the ODE parameters $\Gamma \u0303$ remains the same as described in Eq. (14).

#### 2.2.2 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.

where $Tk,je$ is the duration of the sEMG signal. Applying the encoding to all sEMG signals of the $kth$ trial gives a set of input coefficients, i.e., $Ik={Ik,j}j=1Na$.

*M*hidden layers is then trained to predict the output coefficients $O\u0302k$ by mapping the input coefficients $Ik$ to the target output coefficients $Ok$ as follows:

where $W={W(l)}l=1M+1$,$\u2009b={b(k)}l=1M+1$ and $h(\u22c5)$ 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.

### 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 non-invasive 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.

## 3 Verification: Elbow Flexion-Extension Motion

To verify the proposed PI-PINN framework, an elbow flexion-extension model was considered and the flowchart of the proposed 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.

with the initial conditions $q(0)=\pi 6\u2009radians$ and $q\u02d9(0)=0\u2009radians/sec$. Given the synthetic sEMG signals ($eBi,eTri$) that were plugged into the sEMG-to-activation dynamics equations (Sec. 2.1.1), leading to muscle activations $aBi,\u2009atri$ 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.

Parameter | Type | Value | Parameter | Type | Value |
---|---|---|---|---|---|

$l0,BiM$ | Biceps muscle model | 0.6 m | $mfa$ | Equation of motion | 1.0 kg |

$vmax,BiM$ | Biceps muscle model | 6 m$/sec$ | $lua$ | Geometric | 1.0 m |

$f0,BiM$ | Biceps muscle model | 300 N | $lfa$ | Geometric | 1.0 m |

$ls,BiT$ | Biceps muscle model | 0.55 m | $l1,Bi$ | Geometric | 0.3 m |

$\varphi Bi$ | Biceps muscle model | 0.0 radians | $l2,Bi$ | Geometric | 0.8 m |

$l0,TriM$ | Triceps muscle model | 0.4 m | $l1,Tri$ | Geometric | 0.2 m |

$vmax,TriM$ | Triceps muscle model | 4 m$/sec$ | $l2,Tri$ | Geometric | 0.7 m |

$f0,TriM$ | Triceps muscle model | 300 N | $d$ | Activation dynamics | 0.08 sec |

$ls,TriT$ | Triceps muscle model | 0.33 m | A | Activation dynamics | 0.2 |

$\varphi Tri$ | Triceps muscle model | 0.0 radians |

Parameter | Type | Value | Parameter | Type | Value |
---|---|---|---|---|---|

$l0,BiM$ | Biceps muscle model | 0.6 m | $mfa$ | Equation of motion | 1.0 kg |

$vmax,BiM$ | Biceps muscle model | 6 m$/sec$ | $lua$ | Geometric | 1.0 m |

$f0,BiM$ | Biceps muscle model | 300 N | $lfa$ | Geometric | 1.0 m |

$ls,BiT$ | Biceps muscle model | 0.55 m | $l1,Bi$ | Geometric | 0.3 m |

$\varphi Bi$ | Biceps muscle model | 0.0 radians | $l2,Bi$ | Geometric | 0.8 m |

$l0,TriM$ | Triceps muscle model | 0.4 m | $l1,Tri$ | Geometric | 0.2 m |

$vmax,TriM$ | Triceps muscle model | 4 m$/sec$ | $l2,Tri$ | Geometric | 0.7 m |

$f0,TriM$ | Triceps muscle model | 300 N | $d$ | Activation dynamics | 0.08 sec |

$ls,TriT$ | Triceps muscle model | 0.33 m | A | Activation dynamics | 0.2 |

$\varphi Tri$ | Triceps muscle model | 0.0 radians |

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 $kth$ sample with $n$ temporal data points contained $[xk,qk]=[tk,ek,Bi,ek,Tri,qk]\u2261{tk(i),ek,Bi(i),ek,Tri(i),qk(i)}i=1n$

where $\Gamma l(0)$ was the initial value of the parameter. Therefore, the set of parameters to be identified became $\Gamma \xaf={\Gamma l}l=14$.

and could be plugged into the residual term $Jres$ in the loss function in Eq. (14).

The residual of the equation of motion $(\u2009Jres)$ involved in the loss function $J$ (Eq. (14)) was scaled so that the terms $Jdata$ and $\u2009\beta Jres$ in Eq. (14) had the same unit. This yielded the scaling parameter $\beta \u221d\Delta t2I$, where $\Delta t$ is the time-step size in the time-signal 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\xd710\u22123$. In this example, $\beta \u221d\Delta t2I=10\u22124$. 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, $f0M$ and $vmaxM$, 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 *β* is shown in Appendix C, where the study showed that this estimate of $\beta \u221d\Delta t2I=10\u22124$ lies within a range leading to optimal motion prediction and parameter identification from the framework.

## 4 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 muscle-tendon 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$\u2009\Gamma ={f0,BiM,vmax,BiM,f0,TriM,vmax,TriM}$. The training was performed using the Adam algorithm [56] with an initial learning rate of $5\xd710\u22123$.

### 4.1 Time-Domain PI-PINN Approach.

For the time-domain training, the training data of the $kth$ sample with $nk$ data points contained $[xk,qk]=[tk,ek,Bi,ek,Tri,qk]$. The entire training dataset with a total of $Ndata=n1+n3$ data points could then be defined as follows: the input discrete time and sEMG signals were $t=[t1T,t3T]T$, $eBi=[e1,BiT,e3,BiT]T$, and $eTri=[e1,TriT,e3,TriT]T$. The output was the corresponding angular motion data $q=[q1T,q3T]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 $\sigma =1$.

Figure 10 compares the data of joint kinematics with the predictions from the trained forward dynamics surrogate with time-domain training on both training and testing cases. The noise and oscillations in the raw sEMG signals pose difficulties for the NN surrogate in learning accurate mappings, leading to oscillatory and nonphysiological motion predictions, especially in the testing case. Nevertheless, the motion prediction captured the essential motion pattern, with the mean predicted motion trajectory close to that of the data.

### 4.2 Feature-Encoded PI-PINN (FEPI-PINN) Approach.

For the feature-encoded approach, the set of training data of the $kth$ sample contained ${Ik,Bi,Ik,Tri,Ok}$. The input coefficients of the training data sEMG signals of the biceps and triceps muscles were $IBi=[I1,BiT,I3,BiT]T,\u2009and\u2009ITri=[I1,BiT,I3,TriT]T$. The corresponding output coefficients were $O=[O1T,O3T]T$. The number of Fourier basis terms used for the feature encoding of sEMG signals were $nEMG=1000$ while for the motion, $nq=100$ were used. An NN with two hidden layers, 32 neurons in each layer, was employed.

where $q=[q(1),\u2026,q(n2)]T$ is the motion data of trial 2, $q\u0302=[q\u0302(1),\u2026,q\u0302(n2)]T$ is the prediction from the PI-PINN framework, and $q\xaf$ is the mean of trial 2 s motion data.

Kinematics | Framework | MSE | $R2$ score |

Angular motion | Time-domain | 0.50 | –0.45 |

Feature-encoded | 0.03 | 0.91 | |

Angular velocity | Time-domain | 1973.75 | –826.95 |

Feature-encoded | 0.37 | 0.84 |

Kinematics | Framework | MSE | $R2$ score |

Angular motion | Time-domain | 0.50 | –0.45 |

Feature-encoded | 0.03 | 0.91 | |

Angular velocity | Time-domain | 1973.75 | –826.95 |

Feature-encoded | 0.37 | 0.84 |

### 4.3 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, $f0M$ and $vmaxM$, of both the biceps and the triceps during optimization, with the final converged values of $f0M$ lying within the physiological estimates of these parameters reported in the literature [57–59]. Note that since it is difficult to measure the parameter $vmaxM$ in experiments, an estimate of $vmaxM$, $8\u221210\u2009l0M/sec,$ was typically used in hill-type muscle models [50]. Given that $l0,BiM=0.1157\u2009m$ and $l0,TriM=0.1138\u2009m$ were used in this model, the physiologically reasonable estimated range of $vmax,BiM$ and $vmax,TriM$ were 1.16–1.32 $m/s$ and 1.14–1.34 $m/s$, respectively. The $vmax,BiM$ and $vmax,TriM$ 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 demonstrated the effectiveness of the proposed FEPI-PINN framework and promising potential for real applications.

Parameter | Identified values | Estimates from literature |
---|---|---|

$f0,BiM\u2009(N)$ | 819.56 | 525–849.29 |

$vmax,BiM\u2009(m/s)$ | 1.24 | 1.16–1.32 |

$f0,TriM\u2009(N)$ | 421.37 | 504–1176 |

$vmax,TriM\u2009(m/s)$ | 1.138 | 1.14–1.34 |

Parameter | Identified values | Estimates from literature |
---|---|---|

$f0,BiM\u2009(N)$ | 819.56 | 525–849.29 |

$vmax,BiM\u2009(m/s)$ | 1.24 | 1.16–1.32 |

$f0,TriM\u2009(N)$ | 421.37 | 504–1176 |

$vmax,TriM\u2009(m/s)$ | 1.138 | 1.14–1.34 |

## 5 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, 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 feature-encoded 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 flexion-extension 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 velocity-dependent 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.

## Acknowledgment

All at the University of California San Diego, are very much appreciated.

## Funding Data

National Institute of Health (Grant No. 5R01AG056999-04; Funder ID: 10.13039/100000002).

U.S. Office of Naval Research (Grant No. N00014-20-1-2329; Funder ID: 10.13039/100000006).

## Nomenclature

- $a$ =
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

- $fA$ =
active force component of hill-type muscle model

- $fA,L$ =
length-dependent active force generation component

- $fP\u2009$ =
passive force component of hill-type muscle model

- $fV$ =
velocity-dependent active force generation component

- $f0M$ =
maximum isometric force in the muscle

- $FM$ =
total muscle force

- $FMT$ =
total force produced by muscle-tendon complex

- $Ik,j$ =
set of all input feature-encoded basis coefficients of the $jth$ sEMG signal from the $kth$ trial

- $J$ =
composite loss function of the PI-PINN minimizing data and ODE residual

- $l\u0303M$ =
normalized muscle length

- $lMT$ =
total length of muscle-tendon complex

- $l0M$ =
optimal muscle length corresponding to the maximum isometric force

- $lsT$ =
slack length of the tendon

- $lfa$ =
forearm length

- $lua$ =
upper arm length

- $mfa$ =
mass of forearm

- $Ok$ =
set of all output feature-encoded basis coefficients of the $kth$ trial

- $O\u0302k$ =
predicted output feature-encoded basis coefficients of the $kth$ trial

- $q$ =
generalized angular motion of the MSK system

- $q\u0302k$ =
predicted angular motion of the $kth$ trial from the PI-PINN framework

- $TMT$ =
torque produced by the muscle-tendon complex

- $v\u0303M$ =
normalized muscle velocity

- $vmaxM$ =
maximum contraction velocity

- $u$ =
neural excitations

- $xk$ =
input data from the $kth$ motion trial

- $xF$ =
Fourier feature vector

- $\Gamma i$ =
$ith$ parameter characterizing the parameterized ODE system

- $\theta $ =
set of weights and biases of the neural network

- $\kappa $ =
vector of muscle parameters for each muscle group

- $\varphi $ =
pennation angle between muscle and tendon

### Appendix A: Calculating Muscle Length and Velocity Using a Rigid Tendon Model

Given the current joint angle $q$ and the angular velocity $q\u0307$, the current length, $lMT$ of the MT system can be calculated using trigonometric relations. As an example, consider the following simplified model (Fig. 13).

### Appendix B: Hill-Type Muscle Models

For the length-dependent muscle force relations, this work uses the equations given in Ref. [48,49].

where $\gamma 1=0.075$ and $\gamma 2=6.6$ correspond to parameters in the passive muscle force model related to an adult human. The muscle force velocity relationship $fV(v\u0303M)$ is used directly from Ref. [50].

### Appendix C: Sensitivity Study of Weight Parameter $\beta $ in Composite Loss Function

The scaling parameter $\beta $ in the loss function in Eq. (14) was estimated as $\beta \u221d\Delta t2I$, as described in Secs. 2.2.1 and 3. A sensitivity study performed in this appendix uses the verification problem in Sec. 3 to demonstrate the proper range of $\beta $ for stable and accurate solution of optimization problem in Eq. (14). Six different parameters $\beta =100,10\u22121,10\u22122,10\u22123,10\u22124,10\u22125$ were chosen as shown in Fig. 14 and the results show the motion prediction of one representative training trial and the optimal MSK parameters identified for each $\beta $. The estimated $\beta \u221d\Delta t2I=10\u22124$ falls in the range of $\beta =[10\u22125,10\u22123]$ where good optimization solutions are obtained.

## References

**399**, p.

**449**, p.

**33**, pp.

**8**(2), pp.

**11590**, pp.