The bionic neural network for external simulation of human locomotor
  system
Muscle forces and joint kinematics estimated with musculoskeletal (MSK)
modeling techniques offer useful metrics describing movement quality.
Model-based computational MSK models can interpret the dynamic interaction
between the neural drive to muscles, muscle dynamics, body and joint
kinematics, and kinetics. Still, such a set of solutions suffers from high
computational time and muscle recruitment problems, especially in complex
modeling. In recent years, data-driven methods have emerged as a promising
alternative due to the benefits of flexibility and adaptability. However, a
large amount of labeled training data is not easy to be acquired. This paper
proposes a physics-informed deep learning method based on MSK modeling to
predict joint motion and muscle forces. The MSK model is embedded into the
neural network as an ordinary differential equation (ODE) loss function with
physiological parameters of muscle activation dynamics and muscle contraction
dynamics to be identified. These parameters are automatically estimated during
the training process which guides the prediction of muscle forces combined with
the MSK forward dynamics model. Experimental validations on two groups of data,
including one benchmark dataset and one self-collected dataset from six healthy
subjects, are performed. The results demonstrate that the proposed deep
learning method can effectively identify subject-specific MSK physiological
parameters and the trained physics-informed forward-dynamics surrogate yields
accurate motion and muscle forces predictions.
Yue Shi, Shuhao Ma, Yihui Zhao
# The bionic neural network for external simulation of human locomotor system

###### Abstract

Muscle forces and joint kinematics estimated with musculoskeletal (MSK) modeling techniques offer useful metrics describing movement quality. Model-based computational MSK models can interpret the dynamic interaction between the neural drive to muscles, muscle dynamics, body and joint kinematics, and kinetics. Still, such a set of solutions suffers from high computational time and muscle recruitment problems, especially in complex modeling. In recent years, data-driven methods have emerged as a promising alternative due to the benefits of flexibility and adaptability. However, a large amount of labeled training data is not easy to be acquired. This paper proposes a physics-informed deep learning method based on MSK modeling to predict joint motion and muscle forces. The MSK model is embedded into the neural network as an ordinary differential equation (ODE) loss function with physiological parameters of muscle activation dynamics and muscle contraction dynamics to be identified. These parameters are automatically estimated during the training process which guides the prediction of muscle forces combined with the MSK forward dynamics model. Experimental validations on two groups of data, including one benchmark dataset and one self-collected dataset from six healthy subjects, are performed. The results demonstrate that the proposed deep learning method can effectively identify subject-specific MSK physiological parameters and the trained physics-informed forward-dynamics surrogate yields accurate motion and muscle forces predictions.

Electromyography (EMG) exoskeleton online learning adversarial learning edge-computing parallel computing

## 1 Introduction

Accurate estimation of muscle forces holds pivotal significance in diverse domains. As a powerful computational simulation tool, the musculoskeletal (MSK) model can be applied for detailed biomechanical analysis to understand muscle-generated forces thoroughly, which would be beneficial to various applications ranging from designing efficacious rehabilitation protocols Dong et al. (2006), optimizing motion control algorithms Ding et al. (2010), Luo et al.
[2023], enhancing clinical decision-making Azimi et al. [2020], Karthick et al. [2018], Sekiya et al. [2019] and the performance of athletes Zaman et al. [2022], McErlain-Naylor et al. [2021]. Thus far, the majority of the MSK models are based on physics-based modeling techniques to interpret transformation among neural excitation, muscle dynamics, joint kinematics, and kinetics.

It is challenging to offer a biologically consistent rationale for the choice of any objective function Kim et al. [2009], Modenese et al. [2011], given our lack of knowledge about the method used by the central nervous system (CNS) Rane et al. [2019]. Moreover, estimation of muscle forces and joint motion during inverse dynamics analysis suffers from high computational time and muscle recruitment problems, especially in complex modeling Trinler et al. [2019], Lee et al. [2009], Lund et al. [2015]. The surface electromyography (sEMG) signal is a non-invasive technique that measures the electrical activity of muscles. The sEMG signal is able to reflect the motor intention of the human 20 ms to 200 ms prior to the actual initiation of the joint motion and muscle activation Campanini et al. [2020], Teramae et al. [2018]. MSK forward dynamics model is an effective method for mapping sEMG to muscle force and joint motion. The sEMG-based forward dynamics approaches calculate the muscle forces based on the sEMG signals and muscle properties which can directly reflect the force and velocity of muscle contraction without assuming optimization criteria or constraints. For instance, Zhao et al. proposed an MSK model driven by sEMG for estimating the continuous motion of the wrist joint and a genetic algorithm for parameters optimization Zhao et al. [2020a]. Thomas et al. presented an MSK model to predict joint moment and related muscle forces of knee and elbowBuchanan et al. [2004]. Although MSK forward dynamics approaches are effective, access to individualised physiological parameters is challenging. Optimisation algorithms are used commonly for parameter identification which needs the muscle forces or joint moments calculated from inverse dynamics as the targets. However, the process of obtaining the targets is time-consuming Simonetti et al. [2021], Silvestros et al. [2019]. To address the time-consuming and uncertainty issues of model-based methods, data-driven methods have also been explored to establish relationships between movement variables and neuromuscular status, i.e., from sEMG to joint kinematics and muscle forces Su et al. [2021], Wu et al. [2021]. Burton et al. Burton et al. [2021] implemented four machine/deep learning methods, including recurrent neural network (RNN), convolutional neural network (CNN), fully-connected neural network (FNN), and principal component regression, to predict trend and magnitude of the estimated joint contact and muscle forces. Zhang et al. proposed a physics-informed deep learning framework for the prediction of the muscle forces and joint Kinematics where a new loss function from motion equation was designed as soft constraints to regularize the data-driven model Zhang et al. [2023]. Shi et al. seamlessly integrated Lagrange's equation of motion and inverse dynamic muscle model into the generative adversarial network (GAN) framework for sEMG-based estimation of muscle force and joint kinematics Shi et al. [2023].

However, despite their adaptability, and flexibility, deep learning methods still suffer from potential limitations. It is a 'black box' method, the intermediate processes of data-driven models do not consider the physical significance underlying the modeling process. More importantly, the above deep learning methods used labeled input and output data for training and it is difficult for us to get the value of the muscle forces corresponding to the sEMG. Given the considerations of the limitations concerning both model-based and data-driven models, in this work, we propose a novel PINN framework combining the information from MSK forward dynamics for the estimation of joint angle and related muscle forces, simultaneously identifying the physiological parameters. We find that the entire MSK forward dynamics system can be regarded as an ordinary differential equation (ODE) system with unknown physiological parameters, where the independent variable is the joint angle. Specifically, a fully connected neural network (FNN) estimates the MSK forward dynamics model. The physiological parameters are identified during back-propagation. In addition, muscle forces are predicted through the guide of muscle contraction dynamics.

The main contributions of this paper include: The PINN framework is proposed for the estimation of joint motion and related muscle forces, simultaneously identifying the physiological parameters. In the absence of true values, muscle forces are estimated based on embedded muscle contraction dynamics through sEMG signals and joint motion signals.

## 2 Methods

In this section, an introduction to each sub-process of the MSK forward dynamics system is first presented in Part. A. Then, we demonstrate a novel PINN framework in Part. B.
### Musculoskeletalelat forward dynamics

#### 2.1.1 Muscle activation dynamics

Muscle activation dynamics refer to the process of transforming sEMG signals \(e_{n}\) into muscle activation signals \(a_{n}\) which can be estimated through Eq. 1.

\[a_{n}=\frac{e^{Ae_{n}}-1}{e^{A}-1} \tag{1}\]

where preprocessed sEMG signals \(e_{n}\) are utilized as neural activation signals \(u_{n}\) for estimating muscle activation signals Pau et al. (2012) to simplify the model. \(A\) is a nonlinear factor which has the range of highly non-linearity (-3) to linearity (0.01) Buchanan et al..

#### 2.1.2 Muscle-Tendon Model

Hill's modelling technique is used to compute the muscle-tendon force \(F^{mt}\), which consists of an elastic tendon in series with a muscle fibre. The muscle fibre includes a contractile element (CE) in parallel with a passive elastic element (PE), as shown in Fig. 1. \(l_{n}^{mt},l_{n}^{m},l_{n}^{t}\) is the muscle-tendon length, muscle fibre length and tendon length of the nth muscle respectively. The model is parameterized by the physiological parameters maximum isometric muscle force \(F_{0,n}^{m}\), the optimal muscle fibre length \(l_{0,n}^{m}\), the maximum contraction velocity \(v_{0,n}\), the slack length of the tendon \(l_{s,n}^{t}\) and the pennation angle \(\varphi_{0,n}\), which are difficult to measure in vivo and varies between the age, gender. We employ the following vector \(\chi_{n}\) to represent the personalized physiological parameters of the nth muscle that require identification.

\[\chi_{n}=[l_{0,n}^{m},v_{0,n},F_{0,n}^{m},l_{s,n}^{t},\varphi_{0,n}]\]

Pennation angle \(\varphi_{n}\) is the angle between the orientation of the muscle fibre and tendon, and the pennation angle at the current muscle fibre length \(l_{n}^{m}\) is calculated by

\[\varphi_{n}(l_{n}^{m};\chi_{n})=\sin^{-1}(\frac{l_{0,n}^{m}\sin\varphi_{0,n}} {l_{n}^{m}}) \tag{2}\]

In this study, the tendon is assumed to be rigid Millard et al. (2013) and thus the tendon length \(l_{s,n}^{t}=l_{n}^{t}\) is adopted. The consequent fibre length can be calculated as follows:

\[l_{n}^{m}=(l_{n}^{mt}-l_{n}^{t}){\cos^{-1}\varphi_{n}} \tag{3}\]

The \(F_{CE,n}\), is the active force generated by CE, which can be written as

\[F_{CE,n}(a_{n},l_{n}^{m},v_{n};\chi_{n})=a_{n}f_{v}(\overline{v}_{n})f_{a}( \overline{l}_{n,a}^{m})F_{0,n}^{m} \tag{4}\]

The function \(f_{a}(\cdot)\) represents the active force-length relationship at different muscle fibre lengths and muscle activations, which are written as

\[f_{a}(\overline{l}_{n,a}^{m})=e^{-(\overline{l}_{n,a}^{m}-1)^{2}k^{-1}} \tag{5}\]

where the \(\overline{l}_{n,a}^{m}=l_{n}^{m}/(l_{0,n}^{m}(\lambda(1-a_{n})+1)\) is the normalized muscle fibre length concerning the corresponding activation levels and \(\lambda\) is a constant, which is set to 0.15 Lloyd and Besier (2003). The \(k\) is a constant to approximate the force-length relationship, which is set to 0.45 Thelen (2003). The function \(f_{v}(\overline{v}_{n})\) represents the force-velocity relationship between the \(F_{n}^{m}\) and the normalized contraction velocity \(\overline{v}_{n}\) Zhao et al. (2020)

\[f_{v}(\overline{v}_{n})=\left\{\begin{array}{cc}\frac{0.3(\overline{v}_{n} +1)}{-\overline{v}_{n}^{t}+0.3}&\overline{v}_{n}\leq 0\\ \frac{2.3\overline{v}_{n}+0.039}{1.3\overline{v}_{n}+0.039}&\overline{v}_{n}> 0\end{array}\right. \tag{6}\]

Figure 1: Wrist model with the primary muscle units associated with wrist flexion and extension
where \(\overline{v}_{n}=v_{n}/v_{0,n}\). \(v_{0,n}\) is typically set as 10 \(l_{0,n}^{m}\)/sec Zajac (1989). The \(v_{n}\) is the derivative of the muscle fibre length with respect to time \(t\).

Note that the passive force \(F_{PE,n}\) is the force produced by PE which can be calculated as

\[F_{PE,n}(l_{n}^{m};\chi_{n})=\left\{\begin{array}{cc}0&l_{n}^{m}\leq l_{0,n} ^{m}\\ f_{P}(\overline{l}_{n}^{m})F_{0,n}^{m}&l_{n}^{m}>l_{0,n}^{m}\end{array}\right. \tag{7}\]

where \(\overline{l}_{n}^{m}=l_{n}^{m}/l_{0,n}^{m}\) the normalized muscle fibre length. The \(f_{p}(\cdot)\) is

\[f_{P}(\overline{l}_{n}^{m})=\frac{e^{10\overline{l}_{n}^{m}-1)}}{e^{5}} \tag{8}\]

The \(F_{n}^{mt}\) is the summation of the active force \(F_{CE,n}\) and the passive force \(F_{PE,n}\), which can be written as

\[F_{n}^{mt}(a_{n},l_{n}^{m},v_{n};\chi_{n})=(F_{CE,n}+F_{PE,n})\cos\varphi_{n} \tag{9}\]

#### 2.1.3 Joint Kinematic Modelling Technique

The single joint configuration is presented to estimate the wrist's continuous joint motion. The muscle-tendon length \(l_{n}^{mt}\) and moment arm \(r_{n}\) against wrist joint angle \(q\) are calculated using the polynomial equation and the scale coefficient Ramsay et al. (2009). The total joint torque is calculated as

\[\tau(\mathbf{a},\mathbf{r},\mathbf{l}^{m},\mathbf{v};\tilde{\chi})=\sum_{n=1} ^{N}F_{n}^{mt}(a_{n},l_{n}^{m},v_{n};\chi_{n})r_{n} \tag{10}\]

Where N represents the number of muscles included. \(\mathbf{a},\mathbf{r},\mathbf{l}^{m},\mathbf{v},\tilde{\chi}\) represent the sets of muscle activation, muscle arm length, muscle fibre length, muscle contraction velocity and the physiological parameters of all included muscles, respectively. Since the muscle activation level is not directly related to the joint motion, it is necessary to compute joint acceleration using forward dynamics. In this research, we only consider the simple single hinge joint's flexion and extension and take the model of the wrist joint's flexion as an example. The wrist joint is assumed to be a single hinge joint, and the palm and fingers are assumed to be a rigid segment rotating around the wrist joint in the sagittal plane. Thus, we can have the following relationship based on the Lagrange equation

\[\ddot{q}=\frac{\tau-mgL\sin q-c\dot{q}}{I} \tag{11}\]

where \(q,\dot{q},\ddot{q}\) are the joint motion, joint angular velocity and joint angular acceleration. \(I\) is the moment of inertia of hand, which is equal to \(mL^{2}+I_{p}\). \(I_{p}\) is the moment of inertia at the principal axis which is parallel to the flexion/extension axis. \(m\) and \(L\) are the mass of the hand and the length between the rotation centre to the hand's centre of mass, which are measured from subjects. \(C\) is the damping coefficient representing the elastic and viscous effects from tendons and ligaments. \(\tau\) calculated from Eq. 10. The three parts mentioned above constitute the main components of the MSK forward dynamics model, which establishes the process from sEMG to joint angle motion from a model-based method.

The reason why the model-based approach is introduced here is because it is the foundation on which the PINN framework is built later. We embedded the entire MSK forward dynamics approach into the neural network

The reason why the model-based method is introduced here is because this is the basis on which the PINN framework is built. We embed the entire MSK forward dynamics method into the neural network

### PINN framework

The proposed PINN framework for the estimation of joint angle and muscle forces together with physiological parameters identification is introduced. The computational graph is illustrated in Fig. 2. To be specific, in the data-driven component, a Fully connected Neural Network (FNN) is utilised to automatically extract the high-level features and build the relationship between sEMG signals and the joint angle/muscle forces, while the physics-informed component entails the underlying physical relationship between the joint angle and muscle forces. In this manner, in the data-driven component, the recorded sEMG signals and discrete time steps \(t\) are first fed into FNN. With the features extracted by FNN and integration of the MSK forward dynamics model, the predicted muscle forces and joint angles could be achieved.
#### 2.2.1 Architecture and Training of FNN

To demonstrate the effectiveness of the proposed physics-informed deep learning framework, we chose an architecture of FNN which was composed of four fully-connected blocks and two regression blocks. Specifically, each the fully-connected block has one Linear layer, one ReLU layer and one dropout layer. Similar to the fully-connected blocks, there are one ReLU layer and one dropout layer in each regression block. Between these two regression blocks, one is for the prediction of joint angle, the other is for the prediction of muscle forces. The training was performed using the Adam algorithm with an initial learning rate of \(1*10^{-4}\). In the model training phase, the batch size is set as 1, and FNN is trained by stochastic gradient descent with momentum. Additionally, the maximum iteration is 1000 and the dropout rate is 0.3.

#### 2.2.2 Design of Loss Functions

Unlike state-of-the-art methods, the loss function of the proposed framework consists of the MSE loss and the physics-based losses. The MSE loss is to minimise the mean square error of the ground truth and prediction, while the physics-informed losses are for the parameters identification and guidance on the muscle force prediction. The total loss shows below:

\[L_{total}=L_{q}+L_{r,1}+L_{r,2} \tag{12}\]

\[L_{q}=MSE(q) \tag{13}\]

\[L_{r,1}=R_{1}(q) \tag{14}\]

\[L_{r,2}=R_{2}(q) \tag{15}\]

where \(L_{q}\) denotes the MSE loss function of the joint angle since it has the ground truth from measurement, while \(L_{r,1}\) represents the first residual loss function imposed by the MSK forward dynamics introduced in Sec. 2.1, which helps identify physiological parameters simultaneously. \(L_{r,2}\) denotes the second residual loss function which is utilized for the prediction of muscle forces. \(R_{1}(q),R_{2}(q)\) both indicate the function of the predicted joint angle. The mathematical formulas for each part of the loss will be specified next.

MSE LossMSE loss is calculated by

\[MSE(q)=\frac{1}{T}\sum_{t=1}^{T}(\hat{q}_{t}-q_{t})^{2} \tag{16}\]

where \(q_{t}\) is denoted as the ground truth of the joint angle at time step t, and \(\hat{q}_{t}\) is the predicted joint angle from the network. Additionally, \(T\) denotes the total sample number.

Physics-informed LossPhysics-informed governing laws, reflecting underlying relationships among the muscle force and kinematics in human motion, are converted to constraints during the Neural Network's training phase. The first residual loss function designed from the Lagrange equation which expressed the motion state of the MSK system is

Figure 2: A computational graph of the proposed PINN.
shown below

\[\begin{split} R_{1}(q)&=\frac{1}{T}\sum_{t=1}^{T}(I \dot{q}_{t}+C\dot{q}_{t}+mgL\sin\hat{q}_{t}\\ &-\tau_{t}(\mathbf{a},\mathbf{r},\mathbf{l}^{m},\mathbf{v}^{m}; \tilde{\chi}))\end{split} \tag{17}\]

\(\tau_{t}\) represents the joint torque which is calculated by Eq. 10. The moment arm \(\mathbf{r}\) against the joint angle is calculated using the polynomial equation and the scale coefficient Ramsay et al. (2009), which can be represented as \(\mathbf{r}(q)\). We use \(\mathbf{l}^{mt}(q)\) to represent the relationship between joint angle and the muscle-tendon length which is exported from Opensim. Therefore, muscle fibre length \(\mathbf{l}^{m}(q)\) can also be expressed in terms of the joint angle \(q\) from Eq. 3. The \(\mathbf{v}(q)\) is the derivative of the muscle fibre length \(\mathbf{l}^{m}(q)\), which can be represented in terms of \(q,\tilde{q}\). As mentioned in Eq. 11, this equation expresses the dynamic equilibrium equation that the dynamical system should satisfy at time \(t\). Therefore, for the MSK forward dynamics surrogate model from the neural network, this equation is still satisfied. When the predicted angle of the neural network is \(\hat{q}_{t}\) at time \(t\), the joint torque at this time can be expressed as \(\tau_{t}(\mathbf{a},\mathbf{r}(\hat{q}_{t}),\mathbf{l}^{mt}(\hat{q}_{t}), \mathbf{v}^{m}(\hat{q}_{t});\tilde{\chi})\) = \(\tau_{t}(\mathbf{a},\hat{q}_{t},\hat{q}_{t};\tilde{\chi})\), which can be regarded as an ordinary differential equation (ODE) with unknown physiological parameters \(\tilde{\chi}\) in respect of the predicted joint angle \(\hat{q}_{t}\). During the training phase, \(\hat{q}_{t},\hat{q}_{t}\) can be automatically derived through backpropagation. These physiological parameters \(\tilde{\chi}\) are defined as the weights of the surrogate neural network which are then automatically updated as the neural network minimises the residual loss function \(L_{r,1}\) during backpropagation. The estimation of \(\hat{\chi}\) can be written as

\[\hat{\chi}=\underset{\tilde{\chi}}{\arg\min}L_{r,1} \tag{18}\]

Due to the time-varying nature and difficulty in measurement of muscle forces, it is challenging for us to get the measured muscle forces as the ground truth for training. The second residual loss function \(L_{r,2}\) is designed for the guidance of the muscle forces prediction which is based on \(L_{r,1}\). During each iteration, we identify physiological parameters through \(L_{r,1}\). These potential physiological parameters are used to calculate the muscle forces through the MSK forward dynamics model embedded as the learning target of the neural network. The loss function \(L_{r,2}\) is given by

\[R_{2}(q)=\frac{1}{T}\sum_{t=1}^{T}\sum_{n=1}^{N}(\hat{F}_{t,n}^{mt}-F_{t,n}^{ mt}(a_{t,n},l_{t,n}^{m},v_{t,n}^{m};\hat{\chi}_{n})) \tag{19}\]

where \(\hat{F}_{t,n}^{mt}\) is the estimated muscle force from the network at the time step \(t\). \(F_{t,n}^{mt}\) represents the muscle force calculated by Eq. 9. \(l_{t,n}^{m},v_{t,n}^{m}\)indicates the muscle fibre length and muscle contraction speed of the nth muscle at time \(t\). They all can be expressed in terms of the \(q,\hat{q}\). \(\hat{\chi}_{n}\) stands for the identified physiological parameters for the nth muscle by Eq. 18. Specifically, muscle force of the nth muscle at time \(t\) can be expressed as follows \(F_{t,n}^{mt}(a_{t,n},\hat{q}_{t},\hat{q}_{t};\hat{\chi}_{n})=F_{t,n}^{mt}(a_{t,n},l_{t,n}^{m},v_{t,n}^{m};\hat{\chi}_{n})\) which can be regarded as another ODE with respect to the joint angle \(\hat{q}_{t}\) as well, guiding the muscle forces prediction in the absence of the ground truth.

## 3 Dataset and Experimental Methods

A self-collected dataset of wrist motions is utilised to demonstrate the feasibility of the proposed framework.

#### 3.0.1 Selection and initialisation of physiological parameters to be identified

Among the physiological parameters of the muscle-tendon units involved, we only chose the maximum isometric muscle force \(F_{0,n}^{m}\) and the optimal muscle fibre length \(l_{0,n}^{m}\) for the identification in order to increase the generality of the model. The nonlinear shape factors \(A\) in the activation dynamics also need to be identified. Physiological parameters other than these were obtained by linear scaling based on the initial values of the generic model from Opensim. The selection and initialisation of all the physiological parameters are summarized table. 1. Since there may be differences in terms of magnitude and scale between each parameter because of their different physiological nature, it is necessary to deflate them in about the same interval before training the network.

#### 3.0.2 Self-Collected Dataset

Approved by the MaPS and Engineering Joint Faculty Research Ethics Committee of the University of Leeds (MEEC18-002), 6 subjects participated in this experiment. The consent forms are signed by all subjects. We take the subject's weight data and the length of their hand in order to calculate the moment of inertia of the hand.
In the experiment, subjects were informed to maintain a fully straight torso with the \(90\) abducted shoulder and the \(90\) flexed elbow joint. The continuous wrist flexion/extension motion was recorded using the VICON motion capture system. The joint angle was computed at 250 Hz through the upper limb model using 16 reflective markers. Meanwhile, sEMG signals were recorded by Avanti Sensors at 2000 Hz from the main wrist muscles, including FCR, FCU, ECRL, ECRB and ECU. Moreover, the sEMG signals and motion data were synchronised and re-sampled at 1000 Hz. 5 repetitive trials were performed for each subject, and a three-minute break was given between trials to prevent muscle fatigue.

The measured sEMG signals were band-pass filtered (20 Hz and 450 Hz), fully rectified, and low-pass filtered (6 Hz). Then, they were normalised concerning the maximum voluntary contraction recorded before the experiment, resulting in the enveloped sEMG signals. Each wrist motion trial, consisting of time steps \(t\), pre-processed sEMG signals and wrist joint angles, is formed into at by 7 matrix.

#### 3.0.3 Evaluation Criteria

To quantify the estimation performance of the proposed framework, root mean square error (RMSE) is first used as the metric. RMSE is the response variable. Specifically, RMSE indicates the discrepancies in the amplitude and between the estimated variables and the ground truth, which can be calculated by

\[RMSE=\sqrt{\frac{1}{T}\sum_{t=1}^{T}(\hat{y_{t}}-y_{t})^{2}} \tag{20}\]

where \(y_{t}\) and \(\hat{y_{t}}\) indicate the ground truth and the corresponding predicted value, respectively.

The coefficient of determination which denoted \(R^{2}\) is used as another metric, which could be calculated by

\[R^{2}=1-\frac{\sum_{i=1}^{T}(y_{t}-\hat{y}_{t})^{2}}{\sum_{t=1}^{T}(y_{t}-\hat {y}_{t})^{2}} \tag{21}\]

Where \(T\) is the number of samples, \(y_{t}\) is the measured ground truth at the time \(t\), \(\hat{y}_{t}\) is the corresponding estimation of the model, and \(\bar{y}\) is the mean value of all the samples.

#### 3.0.4 Baseline Methods and Parameters Setting

To verify the effectiveness of the proposed physics-informed deep learning framework in the prediction of the muscle forces and joint angle, several state-of-the-art methods, including FNN, SVR Zhang et al. (2020), are considered as the baseline methods for the comparison. Specifically, FNN has the same neural network architecture as our proposed model which has four fully-connected blocks and two regression blocks but without the physics-informed component. Stochastic gradient descent with Adam optimiser is employed for FNN training, the batch size is set as 1, the maximum iteration is set as 1000, and the initial learning rate is 0.001. The radial basis function (RBF) is selected as the kernel function of the SVR method and the parameter \(C\) which controls the tolerance of the training samples is set as 100 and the kernel function parameters \(\gamma\) which controls the range of the kernel function influence is set to 1.

## 4 Results

In this section, we evaluate the performance of the proposed framework on a self-collected dataset which is detailed in Sec. 3 by comparing it with selected baseline methods. Specifically, we first present the process of parameter

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline Parameters & FCR & FCU & ECRL & ECRB & ECU \\ \hline \(F_{m}^{m}\)(N) & 407 & 479 & 337 & 252 & 192 \\ \(I_{0,n}^{o}\)(m) & 0.062 & 0.051 & 0.081 & 0.058 & 0.062 \\ \(v_{0,n}\)(m/s) & 0.62 & 0.51 & 0.81 & 0.58 & 0.62 \\ \(l_{s,n}^{t}\)(m) & 0.24 & 0.26 & 0.24 & 0.22 & 0.2285 \\ \(\varphi_{0,n}\)(rad) & 0.05 & 0.2 & 0 & 0.16 & 0.06 \\ \hline A & & & 0.01 & & \\ \hline \hline \end{tabular}
\end{table}
Table 1: Physiological parameters involved in the forward dynamics setup of wrist flexion-extension motion for subject 1
identification by our model and verify the reliability of the identification results. Then, the overall comparisons depict the outcomes of both the proposed framework and the baseline methods. This includes representative results showcasing the included muscles' force and joint angles, alongside comprehensive predictions for six healthy subjects. In addition, we investigate the robustness and generalization performance of the proposed framework in the intrasession scenario. Finally, the effects of the network architectures and hyperparameters on our model are evaluated separately. The proposed framework and baseline methods are trained using PyTorch on a laptop with GeForce RTX 3070 Ti graphic cards and 32G RAM.

### Parameters Identification

The subject-specific physiological parameters are identified during the training. Table. 2 presents the estimation and the variation of the parameters of subject one. The variation pertains to the extent of disparity between the initial guess and the determined value derived from our model. The initial guess of the parameters is given in Table. 1 in Sec. 2. Physiological boundaries of the parameters are chosen according to Saul et al. (2015). The boundaries of maximum

Figure 4: Evolution of the optimal muscle fiber length \(l_{0}^{m}\) identified during the training process of subject 1.

Figure 3: Evolution of the maximum isometric muscle force \(F_{0}^{m}\) identified during the training process of subject 1.
isometric force \(F_{0}^{m}\) are set to \(50\%\) of the initial guess, while the boundaries of optimal muscle fibre length \(l_{0}^{m}\) are set to \(\pm 0.001\) of the initial guess. As shown in the Table. 2, the parameters identified through our proposed framework are all within the physiological range and possess physiological consistency. Furthermore, the deviations of the optimal fibre length exhibit minimal disparity from the initial value, while the maximum isometric forces diverge significantly from the initial value. As well as the two main muscle-tendon physiological parameters, the identified non-linear shape factor A muscle activation dynamics parameter A is -2.29 which is physiologically acceptable in the range of -3 to 0.01.

Fig. 4 shows the evolution of the physiological parameters identified during the training process. In conjunction with the comparison of results obtained through the PINN framework within physiological ranges, we establish the credibility of the results from another perspective. We set the estimation of the physiological parameters from the MSK forward dynamics model optimized by Genetic Algorithm (GA) as a reference for comparison Zhao et al. (2020). As depicted in Fig. 4, the red dashed lines represent the reference values estimated by the GA optimization method, the black dashed lines indicate the estimates from our proposed framework, and the blue solid lines illustrate the variation process of the parameters. The comparison revealed small differences between the values gained from two different methods.

### Overall Comparisons

The overall comparisons between the proposed framework and baseline methods are then performed. Fig. 5 depicts the representative results of the proposed framework for the wrist joint prediction, including the wrist angle, muscle force

\begin{table}
\begin{tabular}{l c c c} \hline \hline  & \multicolumn{3}{c}{Parameter index} \\ \cline{2-4}  & \multicolumn{2}{c}{\(l_{0}^{lm}\)(m)} & \multicolumn{2}{c}{\(F_{0}^{m}\)(N)} \\ \cline{2-4} Muscle index & Estimation\(\&\)Ariation & Estimatio\(\&\)Ariation \\ \hline FCR & 0.056 & 90.78\% & 475.2 & 116.71\% \\ FCU & 0.061 & 120.62\% & 644.1 & 129.37\% \\ ECRB & 0.050 & 86.15\% & 166.7 & 66.05\% \\ ECRL & 0.082 & 101.85\% & 475.2 & 140.91\% \\ ECU & 0.052 & 84.51\% & 286.6 & 148.58\% \\ \hline \hline \end{tabular}
\end{table}
Table 2: Identified physiological parameters of Subject one

Figure 5: Representative results of the wrist joint through the proposed data-driven model. The predicted outputs include the wrist angle, FCR muscle force, FCU muscle force, ECRL muscle force, ECRB muscle force, and ECU muscle force.
\begin{table}
\begin{tabular}{c c c c c c c c c} \hline \hline Subject & Methods & FCR(N) & FCU(N) & ECRL(N) & ECRB(N) & ECU(N) & Angle(rad) \\ \hline \multirow{3}{*}{S1} & Ours & 4.99/0.98 & 4.62/0.97 & 4.35/0.98 & 5.33/0.96 & 3.12/0.96 & 0.09/0.98 \\  & FNN & 3.27/0.99 & 5.37/0.96 & 4.20/0.98 & 4.29/0.98 & 3.03/0.96 & 0.09/0.98 \\  & SVR & 6.02/0.97 & 5.59/0.96 & 6.05/0.96 & 4.50/0.98 & 5.21/0.94 & 0.11/0.97 \\ \hline \multirow{3}{*}{S2} & Ours & 6.01/0.96 & 4.47/0.97 & 3.71/0.98 & 4.95/0.97 & 2.56/0.97 & 0.08/0.98 \\  & FNN & 5.57/0.97 & 4.69/0.97 & 2.79/0.99 & 5.76/0.96 & 2.04/0.98 & 0.06/0.99 \\  & SVR & 7.35/0.95 & 5.58/0.96 & 5.93/0.96 & 6.01/0.96 & 4.31/0.96 & 0.09/0.98 \\ \hline \multirow{3}{*}{S3} & Ours & 5.43/0.98 & 3.74/0.98 & 3.49/0.97 & 3.91/0.98 & 3.20/0.96 & 0.13/0.96 \\  & FNN & 6.25/0.97 & 4.27/0.98 & 5.29/0.97 & 5.59/0.96 & 2.97/0.96 & 0.13/0.96 \\  & SVR & 7.63/0.95 & 6.91/0.96 & 7.97/0.95 & 6.25/0.96 & 3.21/0.94 & 0.14/0.96 \\ \hline \multirow{3}{*}{S4} & Ours & 5.54/0.97 & 5.21/0.97 & 5.31/0.97 & 3.58/0.97 & 3.92/0.97 & 0.04/0.99 \\  & FNN & 5.38/0.97 & 5.82/0.96 & 5.57/0.97 & 4.77/0.96 & 4.28/0.96 & 0.08/0.98 \\  & SVR & 6.27/0.96 & 6.35/0.96 & 6.45/0.96 & 5.42/0.96 & 5.39/0.95 & 0.10/0.97 \\ \hline \multirow{3}{*}{S5} & Ours & 5.41/0.96 & 5.23/0.97 & 3.98/0.98 & 3.87/0.97 & 4.14/0.95 & 0.11/0.97 \\  & FNN & 4.57/0.97 & 4.61/0.97 & 4.53/0.97 & 4.41/0.97 & 3.82/0.96 & 0.14/0.96 \\  & SVR & 8.65/0.95 & 5.34/0.97 & 6.16/0.96 & 5.79/0.96 & 5.75/0.94 & 0.18/0.94 \\ \hline \multirow{3}{*}{S6} & Ours & 3.21/0.99 & 3.27/0.98 & 5.61/0.97 & 6.28/0.94 & 3.97/0.95 & 0.06/0.99 \\  & FNN & 3.97/0.98 & 5.11/0.97 & 5.50/0.97 & 6.24/0.94 & 4.32/0.94 & 0.05/0.99 \\ \cline{1-1}  & SVR & 7.09/0.94 & 6.19/0.96 & 7.15/0.95 & 5.72/0.95 & 4.97/0.93 & 0.07/0.98 \\ \hline \hline \end{tabular}
\end{table}
Table 3: RMSE AND R\({}^{2}\) of the proposed framework and baseline methods of wrist joint

Figure 6: Average RMSEs across all the subjects in (a) included muscles’ force and (b) joint angle, respectively.
of FCR, muscle force of FCU, muscle force of ECRL, muscle force of ECRB, muscle force of ECRB, and muscle force of ECU, respectively. For the prediction of the joint angle, we use the measured value as the ground truth. However, for the prediction of muscle forces, since it is challenging to measure directly, we use the muscle forces calculated from the sEMG-driven MSK forward dynamics model as the ground truthZhao et al. (2020). As presented in Fig. 5, the predicted values of muscle forces and joint angle can fit the ground truths well, indicating the great dynamic tracking capability of the proposed framework.

To quantitatively evaluate the performance of our model, detailed comparisons of all the subjects between the PINN and the baseline methods are presented in Table. 3. We use the data with the same flexion speed to train and test the proposed framework and baseline methods in the wrist joint, according to Table. 3, deep learning-based methods, including the proposed framework, FNN, achieve better-predicted performance than the machine learning-based method, the SVR, as evidenced by smaller RMSEs and higher \(R^{2}\) in most cases. Because these deep learning-based methods could automatically extract high-level features from the collected data.

Fig. 6 (a) illustrates the average RMSEs of the muscle forces of the PINN framework and baseline methods across all the subjects. Among deep learning-based methods, our proposed framework achieves an overall performance similar to that of the FNN without direct reliance on actual muscle force labels. Fig. 6 (b) illustrates the average RMSEs of the joint angle. Our model achieves satisfactory performance with lower standard deviations, and its predicted results are with smaller fluctuations compared with the pure FNN method and the machine learning method. In the training process, the labelled joint angles are not only trained by the MSE loss but also enhanced by the physics-informed losses. The embedded physics laws provide the potential relationships between the output variables as learnable features for the training process.

### Evaluation of Intression Scenario

The intraression scenario is also considered to validate the robustness of the proposed framework. For each subject, we train the model with data of one flexion speed and then test the performance using data of a different flexion speed from the training data which can be seen as a dataset that our trained model is unfamiliar with. Fig. 7 depicts the corresponding experimental results, Our model demonstrates exceptional performance in datasets with different distributions, but the predicted results of some baseline methods are degraded. In particular, concerning the prediction results of the muscle forces of the ECRB and ECU, the prediction results yielded by our proposed framework demonstrate a notably enhanced congruence with the underlying ground truth. The SVR method also demonstrates its efficacy in accurately tracking specific predictions, such as the prediction of the muscle forces of the FCU and ECRL. Nevertheless, errors remain in other prediction results. Hence, the application of this method is constrained by various factors, including the

Figure 7: Evolution of the physiological parameters identified during the training process.
nature of the prediction target, leading to both instability and inherent limitations. The SVR method also exhibits the proficiency of dynamical tracking in part of the prediction results, such as the prediction of the muscle forces of the FCU and ECRL. However, the error remains in other prediction results which reflects the limitation of the stability. Although the FNN method shows favourable performance on the training data, when it comes to the intrasession scenario, the performance got much worse than that on the training data. Our proposed framework manifests a discernible capability to predict muscle forces and joint angles on data characterized by the diverse distribution in the absence of the real value of the muscle forces.

### Effects of Hyperparameters

To investigate the effects of hyperparameters, including learning rate, types of activation functions, and batch size, on the performance of the proposed framework, the detailed results are shown in Table. 4, Table. 5, and Table. 6. Specifically, we consider three different learning rates, i.e., 0.01, 0.001 and 0.0001, and the maximum iteration is set as 1000. According to Table. 4, we can find that the proposed framework achieves better performance with a smaller learning rate. As observed in Table. 5, the proposed framework with ReLU and Tanh could achieve comparable performance, but better than the one using Sigmoid as the activation function. In Table. 6, the performance of the proposed framework achieves better performance with the decrease of the batch size and achieves its best performance when the batch size is 1. Even though setting the batch size to 1 for training the deep learning framework can still retain the drawbacks such as over-fitting, poor generalisation, unstable training, etc., in previous experiments, they have been proven not to be issues.
series sequences, such as Recurrent Neural Network (RNN) and Long Short-Term Memory (LsTM) is another possible solution for the problem.

For our proposed framework, the number of physiological parameters that can be efficiently identified is limited considering the network's performance and the time cost of model training. However, a more complex but physiologically consistent model with more parameters will directly benefit the accuracy of the muscle force prediction in our proposed model. Therefore, in future work, we will devoted to enhancing the model's performance to effectively identify a broader range of physiological parameters which will enable a deeper integration of the MSK forward dynamics model into our framework.

In order to improve the experimental feasibility and generalisation of the model, we partially simplified the MSK forward dynamics model by reducing the number of individualised physiological parameters and by choosing a rigid rather than an elastic tendon. A more physiologically accurate representation of muscle tissues with connective tissues and muscle fibres could inform the muscles' length and velocity-dependent force generation capacity. In addition, We will achieve precise physiological parameter identification and muscle force prediction on more accurate, complex MSK models in the future.

## 6 Conclusion

This paper presents a novel physics-informed deep-learning framework for joint angle and muscle force estimation. Specifically, our proposed model uses the ODE based on the MSK forward dynamics as the residual loss for the identification of personalized physiological parameters and another residual constraint based on the muscle contraction dynamics for the estimation of muscle forces in the absence of the ground truth. Comprehensive experiments on the wrist joint for muscle forces and joint angles' prediction indicate the feasibility of the proposed framework. However, it is worth noting that the entire MSK forward dynamics model embedded in the framework needs to be adjusted when the proposed approach is extended to other application cases.
