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 [25]. 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 [913], 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 [2123]. 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 [2932] 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 [3437], and inverse problems [29,30,3840]. 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.12.1.3.

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.
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.
Close modal

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]
u(t)=e(td)
(1)
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
a(t)=exp(Au(t))1exp(A)1
(2)

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 (ϕ) to be discussed as follows.

First, let each muscle group be characterized by a parameter vector
κ=[l0M,vmaxM,f0M,lsT,ϕ]
(3)
The force generated by the muscle group contains an active and a passive component [48,49]. The active force fA component can be expressed as
fA(a,l̃M,ṽM;κ)=afA,L(l̃M;κ)fV(ṽM;κ)l̃M=lM/l0MṽM=vM/vmaxM
(4)
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 fA,L(l̃M) and fV(ṽM) are functions of the length and velocity-dependent force generation properties of the active muscle represented by generic functions of dimensionless quantities as given in Appendix  B [4850], allowing them to be scaled to specific muscles through the parameters described. The total muscle force FM can be expressed as
FM(a,l̃M,ṽM;κ)=f0M(fA(a,l̃M,ṽM;κ)+fP(l̃M;κ))
(5)

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

The muscle force FM obtained in Eq. (5) is transmitted to the joints through the tendon (Fig. 2(b)). The tendon produces force FT only when its length lT is stretched beyond its slack length lsT. In this study, the tendon is assumed to be rigid [6,47] and thus the tendon length lT=lsT is adopted. According to the schematics of the MT complex shown in Fig. 2(b), the total length of MT complex, lMT, is first obtained by
lMT=lsT+lMcosϕ
(6)
where ϕ is the pennation angle. The total force produced by the MT complex, FMT, can be expressed as follows based on the force equilibrium:
FMT(a,l̃M,ṽM,ϕ;κ)=FM(a,l̃M,ṽM;κ)cosϕ
(7)
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).
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).
Close modal

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.

2.1.3 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¨TMT(a,q,q˙;κ)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; TMT is the torque from all muscles in the model calculated by TMT(a,q,q˙;κ)=R(q)FMT(a,q,q˙;κ), where R(q) are the moment arm's and FMT(a,q,q˙;κ) is obtained from Eq. (7). As the muscle lengths lM and velocities vM 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 FMT(a,q,q˙;κ). 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].

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.

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
L[q(t);λ]=b(t;ω)t(0,T],B[q(0)]=g
(9)

where the differential operator L[;λ] is parameterized by a set of parameters λ={λ0,λ1,,λn}. The right-hand side b(t;ω) is parameterized by ω=[ω1,ω2,,ω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 Γ={λ,ω}. The solution to the ODE system q:[0,T]R depends on the choice of parameters Γ.

2.2.1 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, xRnin, to discrete joint motion data outputs, qRnout, by a nonlinear activation function h(·) as follows:
a(k)=h(W(k)a(k1)+b(k)),k=1,,M
(10)
where M is the number of hidden layers; the superscript (k) denotes the layer number. Here, the activation a(k)Rnk is the output of layer k with nk neurons. For the input layer, n0=nin is the input dimension and a(0)=x; W(k)Rnk×nk1 and b(k)Rnk are trainable weight and bias coefficients, respectively; θ={W,b} is used to denote all trainable parameters of the network, where W={W(k)}k=1M+1 and b={b(k)}k=1M+1 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
q̂=W(M+1)a(M)+b(M+1)
(11)

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 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,,xk(nk)T]T with the ith row as xk(i)T=[tk(i),ek,1(i),,ek,Na(i)] while the output data is qk=[qk(1),,qk(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=i=1pni data points with the input as x=x1,,xpT and the target output as q=q1,,qpT.

The NN is trained to learn the mapping from the raw sEMG and time (xk(i)) to the MSK joint motion (qk(i)). The trainable parameters of the NN are obtained by minimizing the following loss function:
Jdata=1Ndatai=1pq̂i(xi;θ)qiL22
(12)
where ·L22 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)
Jres=1Ndatai=1pr(q̂i(xi;θ);Γ)L22
(13)
where r(q̂i(xi;θ);Γ)=L[q̂i(xi;θ);λ]b(t(i);ω) denotes the vector of residuals associated with Eq. (8) for the ith sample; q̂i(xi;θ) is the vector of predicted MSK joint motion from the NN with the trainable parameters θ for the ith sample; Γ={λ,ω} represents the set of ODE parameters relevant to the MSK system. The optimal NN parameters θ̃ and the ODE parameters Γ̃ are obtained by minimizing the composite loss function J as follows:
θ̃,Γ̃=argminθ,Γ(J)=argminθ,Γ(Jdata+βJres)
(14)

where β 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, β was scaled so that the terms Jres and Jdata 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 (θ,Γ) can be obtained efficiently by automatic differentiation [52]. The computational graph is illustrated in Fig. 3.

Fig. 3
A computational graph of time-domain PI-PINNs with Fourier features prefixed to the network
Fig. 3
A computational graph of time-domain PI-PINNs with Fourier features prefixed to the network
Close modal
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.

A Fourier feature xFRm of a given input xRd is defined as [42,43]
xF=sin(WFx+bF)
(15)

where WFRm×d and bFRm are weight and bias coefficients, respectively, randomly sampled from a normal distribution N(0,σ2), with σ 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 θ̃ and the ODE parameters Γ̃ 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.

Due to the periodic nature of the signals, Fourier basis was used along with polynomials up to quadratic terms in the feature-encoded training. The polynomial basis helps capture the signals close to the boundary of the time-domain that have incomplete periodicity. For the jth sEMG signal from the kth trial, ek,j={ek,j(i)}i=1nk, consider the employment of 2nEMG Fourier basis terms with their respective coefficients Ak,j={Ak,j(l)}l=1nEMG{Bk,j(l)}l=1nEMG and three coefficients for the complete quadratic basis terms, i.e., Ck,j(0),Ck,j(1), and Ck,j(2). The set of all input coefficients is denoted by Ik,j, defined as
Ik,j={Ak,j,Bk,j,Ck,j(0),Ck,j(1),Ck,j(2)}
(16)
and the encoded approximation of the jth sEMG signal from the kth trial is given by
ek,jFE(t;Ik,j)=l=1nEMG[Ak,j(l)cos(2πltTk,je)+Bk,j(l)sin(2πltTk,je)]+Ck,j(0)+Ck,j(1)t+Ck,j(2)t2
(17)

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.

Correspondingly, there would be 2nq+3 coefficients for the encoded approximation of the kth motion signal with the Fourier coefficients Fk={Fk(l)}l=1nq,Gk={Gk(l)}l=1nq and polynomial coefficients Hk(0),Hk(1), and Hk(2). The set of coefficients Ok is defined as
Ok={Fk,Gk,Hk(0),Hk(1),Hk(2)}
(18)
qkFE(t;Ok)=l=1nq[Fk(l)cos(2πltTkq)+Gk(l)sin(2πltTkq)]+Hk(0)+Hk(1)t+Hk(2)t2
(19)
where Tkq 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 Ik to the target output coefficients Ok as follows:
Ôk=W(M+1)(h(W(M)(h(W(1)Ik+b(1)))))+b(M+1)
(20)

where W={W(l)}l=1M+1,b={b(k)}l=1M+1 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 Ok, 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 feature-encoded PI-PINN, termed FEPI-PINN, is shown in Fig. 4 
θ̃,Γ̃=argminθ,Γ(JFE)=argminθ,Γ(JdataFE+βJresFE)
(21)
JdataFE=i=1pÔi(Ii;θ)OiL22
(22)
JresFE=i=1pr(qFE(t;Ôi(Ii;θ));Γ)L22
(23)
Fig. 4
The computational graph for the FEPI-PINN framework
Fig. 4
The computational graph for the FEPI-PINN framework
Close modal

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.

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.
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.
Close modal

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.

The model contained two rigid links resembling the upper arm and forearm with lengths lua and lfa, 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-of-freedom 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 κBi and κ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 mfa was attached to one end of the forearm link with a moment arm lfa 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
Iq¨=E(q)+TMT(aBi,aTri,q,q˙;κBi,κTri)
(24)
where
I=mfalfa2
(25)
E(q)=mfaglfasin(q)
(26)
TMT(aBi,aTri,q,q˙;κBi,κTri)=TBiMT(aBi,q,q˙;κBi)TTriMT(aTri,q,q˙;κTri)
(27)
TBiMT(aBi,q,q˙;κBi)=FBiMT(aBi,l̃BiM,ṽBiM,l̃BiT;κBi)l2,Bisin(q)l1,BilBiMT(q)
(28)
TTriMT(aTri,q,q˙;κTri)=FTriMT(aTri,l̃TriM,ṽTriM,l̃TriT;κTri)l2,Trisin(q)l1,TrilTriMT(q)
(29)

with the initial conditions q(0)=π6radians and q˙(0)=0radians/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,atri 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.

Table 1

Parameters involved in the forward dynamics setup of elbow flexion-extension motion

ParameterTypeValueParameterTypeValue
l0,BiMBiceps muscle model0.6 mmfaEquation of motion1.0 kg
vmax,BiMBiceps muscle model6 m/secluaGeometric1.0 m
f0,BiMBiceps muscle model300 NlfaGeometric1.0 m
ls,BiTBiceps muscle model0.55 ml1,BiGeometric0.3 m
ϕBiBiceps muscle model0.0 radiansl2,BiGeometric0.8 m
l0,TriMTriceps muscle model0.4 ml1,TriGeometric0.2 m
vmax,TriMTriceps muscle model4 m/secl2,TriGeometric0.7 m
f0,TriMTriceps muscle model300 NdActivation dynamics0.08 sec
ls,TriTTriceps muscle model0.33 mAActivation dynamics0.2
ϕTriTriceps muscle model0.0 radians
ParameterTypeValueParameterTypeValue
l0,BiMBiceps muscle model0.6 mmfaEquation of motion1.0 kg
vmax,BiMBiceps muscle model6 m/secluaGeometric1.0 m
f0,BiMBiceps muscle model300 NlfaGeometric1.0 m
ls,BiTBiceps muscle model0.55 ml1,BiGeometric0.3 m
ϕBiBiceps muscle model0.0 radiansl2,BiGeometric0.8 m
l0,TriMTriceps muscle model0.4 ml1,TriGeometric0.2 m
vmax,TriMTriceps muscle model4 m/secl2,TriGeometric0.7 m
f0,TriMTriceps muscle model300 NdActivation dynamics0.08 sec
ls,TriTTriceps muscle model0.33 mAActivation dynamics0.2
ϕTriTriceps muscle model0.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]{tk(i),ek,Bi(i),ek,Tri(i),qk(i)}i=1n

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. 6
The dataset with synthetic bicep and tricep sEMG signals with variations in frequency for all trials, and the corresponding motion from the two
Close modal
The set of MSK parameters Γ={Γl}l=14={f0,BiM,vmax,BiM,f0,TriM,vmax,TriM} 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 f0M and vmaxM, which could affect the conditioning of the parameter identification system. To mitigate this issue, normalization [40] was applied to each of the parameters
Γ¯l=ΓlΓl(0)
(30)

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

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 Γ¯ by optimizing Eq. (14), where the residual of the governing equation for the training input xk(i), was expressed as
r(q̂k(xk(i);θq),q̂˙k(xk(i);θq),q̂¨k(xk(i);θq);Γ(Γ;Γ(0)))
=Iq̂̈k(xk(i);θq)E(q̂k(xk(i);θq))TMT(aBi(t(i)),aTri(t(i)),q̂k(xk(i);θq),q̂̇k(xk(i);θq);Γ(Γ¯;Γ(0)))
(31)

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

The residual of the equation of motion (Jres) involved in the loss function J (Eq. (14)) was scaled so that the terms Jdata and βJres in Eq. (14) had the same unit. This yielded the scaling parameter βΔt2I, where Δ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×103. In this example, βΔt2I=104. 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 βΔt2I=104 lies within a range leading to optimal motion prediction and parameter identification from the framework.

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. 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
Close modal
Fig. 8
Evolution of the MSK parameters, f0M and vmaxM, of bicep and triceps during the training process
Fig. 8
Evolution of the MSK parameters, f0M and vmaxM, of bicep and triceps during the training process
Close modal

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.

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

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Γ={f0,BiM,vmax,BiM,f0,TriM,vmax,TriM}. The training was performed using the Adam algorithm [56] with an initial learning rate of 5×103.

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 σ=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.

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. 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).
Close modal

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,andITri=[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.

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 feature-encoded 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 training performs better in terms of motion prediction as shown in these figures and Table 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 R2 statistic close to 1 were seen in the case of the feature encoding. Here, the mean squared error (MSE) and R2 statistic are defined as
MSE=1n2qq̂L22
(32)
R2=1i=1n2(q(i)q̂(i))2i=1n2(q(i)q¯)2
(33)

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

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. 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).
Close modal
Table 2

Statistics comparing the testing predictions of the framework with original time-domain training versus the enhanced feature-encoded training

KinematicsFrameworkMSER2 score
Angular motionTime-domain0.50–0.45
Feature-encoded0.030.91
Angular velocityTime-domain1973.75–826.95
Feature-encoded0.370.84
KinematicsFrameworkMSER2 score
Angular motionTime-domain0.50–0.45
Feature-encoded0.030.91
Angular velocityTime-domain1973.75–826.95
Feature-encoded0.370.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 [5759]. Note that since it is difficult to measure the parameter vmaxM in experiments, an estimate of vmaxM, 810l0M/sec, was typically used in hill-type muscle models [50]. Given that l0,BiM=0.1157m and l0,TriM=0.1138m 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.

Fig. 12
Evolution of the identified MSK parameters: (a) f0M and (b) vmaxM 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 Table 3.
Fig. 12
Evolution of the identified MSK parameters: (a) f0M and (b) vmaxM 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 Table 3.
Close modal
Table 3

The identified parameter estimates using feature encoding training, and their values reported in the literature [5759]

ParameterIdentified valuesEstimates from literature
f0,BiM(N)819.56525–849.29
vmax,BiM(m/s)1.241.16–1.32
f0,TriM(N)421.37504–1176
vmax,TriM(m/s)1.1381.14–1.34
ParameterIdentified valuesEstimates from literature
f0,BiM(N)819.56525–849.29
vmax,BiM(m/s)1.241.16–1.32
f0,TriM(N)421.37504–1176
vmax,TriM(m/s)1.1381.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 =

    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̃M =

    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

  •  
  • Ôk =

    predicted output feature-encoded basis coefficients of the kth trial

  •  
  • q =

    generalized angular motion of the MSK system

  •  
  • q̂k =

    predicted angular motion of the kth trial from the PI-PINN framework

  •  
  • TMT =

    torque produced by the muscle-tendon complex

  •  
  • ṽM =

    normalized muscle velocity

  •  
  • vmaxM =

    maximum contraction velocity

  •  
  • u =

    neural excitations

  •  
  • xk =

    input data from the kth motion trial

  •  
  • xF =

    Fourier feature vector

  •  
  • Γi =

    ith parameter characterizing the parameterized ODE system

  •  
  • θ =

    set of weights and biases of the neural network

  •  
  • κ =

    vector of muscle parameters for each muscle group

  •  
  • ϕ =

    pennation angle between muscle and tendon

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

The total length of the MT system lMT is given by
lMT=lMcosϕ+lT
(A1)

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

Fig. 13
Simplified model to explain geometrical relations between lM and lMT
Fig. 13
Simplified model to explain geometrical relations between lM and lMT
Close modal
The total length and velocity of the MT complex can then be written as
lMT=l12+l22+2l1l2cosqvMT=l̇MT=q̇(l1l2sinq)/lMT
(A2)
where l1,l2 are the projected lengths of the MT complex on the forearm and upper-arm. By using the fixed-height assumption to maintain a constant muscle volume
lMsinϕ=l0Msinϕ0
(A3)
the muscle length can then be found as
lM=(lMTlST)2+(l0Msinϕ0)2
(A4)
using which the pennation angle ϕ can also be found. The velocity of the MT system can then be found as
vM=vMTcosϕ
(A5)

Appendix B: Hill-Type Muscle Models

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

The active muscle force dependent on variation in length is given as
fA,L(l̃M)=9(l̃M0.4)2,l̃M0.614(1l̃M)2,0.6l̃M1.49(l̃M1.6)2,l̃M>1.4
(B1)
fP(l̃M)=0,l̃M1γ1(exp(γ2(l̃M1))1),1l̃M1.4(γ1γ2exp(0.4γ2))l̃M+γ1((11.4γ2)exp(0.4γ2)1),l̃M>1.4
(B2)

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

Appendix C: Sensitivity Study of Weight Parameter β in Composite Loss Function

The scaling parameter β in the loss function in Eq. (14) was estimated as βΔ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 β for stable and accurate solution of optimization problem in Eq. (14). Six different parameters β=100,101,102,103,104,105 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 β. The estimated βΔt2I=104 falls in the range of β=[105,103] where good optimization solutions are obtained.

Fig. 14
Results of the sensitivity study showing the range of β within which the framework shows stable results
Fig. 14
Results of the sensitivity study showing the range of β within which the framework shows stable results
Close modal

References

1.
Pandy
,
M. G.
,
2001
, “
Computer Modeling and Simulation of Human Movement
,”
Annu. Rev. Biomed. Eng.
,
3
(
1
), pp.
245
273
.10.1146/annurev.bioeng.3.1.245
2.
Buchanan
,
T. S.
,
Lloyd
,
D. G.
,
Manal
,
K.
, and
Besier
,
T. F.
,
2004
, “
Neuromusculoskeletal Modeling: Estimation of Muscle Forces and Joint Moments and Movements From Measurements of Neural Command
,”
ASME J. Appl. Biomech.
,
20
(
4
), pp.
367
395
.10.1123/jab.20.4.367
3.
Lloyd
,
D. G.
, and
Besier
,
T. F.
,
2003
, “
An EMG-Driven Musculoskeletal Model to Estimate Muscle Forces and Knee Joint Moments In Vivo
,”
J. Biomech.
,
36
(
6
), pp.
765
776
.10.1016/S0021-9290(03)00010-1
4.
Pau
,
J. W. L.
,
Xie
,
S. S. Q.
, and
Pullan
,
A. J.
,
2012
, “
Neuromuscular Interfacing: Establishing an EMG-Driven Model for the Human Elbow Joint
,”
IEEE Trans. Biomed. Eng.
,
59
(
9
), pp.
2586
2593
.10.1109/TBME.2012.2206389
5.
Buffi
,
J. H.
,
Werner
,
K.
,
Kepple
,
T.
, and
Murray
,
W. M.
,
2015
, “
Computing Muscle, Ligament, and Osseous Contributions to the Elbow Varus Moment During Baseball Pitching
,”
Ann. Biomed. Eng.
,
43
(
2
), pp.
404
415
.10.1007/s10439-014-1144-z
6.
Zhao
,
Y.
,
Zhang
,
Z.
,
Li
,
Z.
,
Yang
,
Z.
,
Dehghani-Sanij
,
A. A.
, and
Xie
,
S.
,
2020
, “
An EMG-Driven Musculoskeletal Model for Estimating Continuous Wrist Motion
,”
IEEE Trans. Neural Syst. Rehabil. Eng.
,
28
(
12
), pp.
3113
3120
.10.1109/TNSRE.2020.3038051
7.
Kian
,
A.
,
Pizzolato
,
C.
,
Halaki
,
M.
,
Ginn
,
K.
,
Lloyd
,
D.
,
Reed
,
D.
, and
Ackland
,
D.
,
2019
, “
Static Optimization Underestimates Antagonist Muscle Activity at the Glenohumeral Joint: A Musculoskeletal Modeling Study
,”
J. Biomech.
,
97
, p.
109348
.10.1016/j.jbiomech.2019.109348
8.
Goodfellow
,
I.
,
Bengio
,
Y.
, and
Courville
,
A.
,
2016
,
Deep Learning
,
MIT Press
,
Cambridge, MA
.
9.
Kaneko
,
S.
,
Wei
,
H.
,
He
,
Q.
,
Chen
,
J.-S.
, and
Yoshimura
,
S.
,
2021
, “
A Hyper-Reduction Computational Method for Accelerated Modeling of Thermal Cycling-Induced Plastic Deformations
,”
J. Mech. Phys. Solids
,
151
, p.
104385
.10.1016/j.jmps.2021.104385
10.
Kim
,
B.
,
Azevedo
,
V. C.
,
Thuerey
,
N.
,
Kim
,
T.
,
Gross
,
M.
, and
Solenthaler
,
B.
,
2019
, “
Deep Fluids: A Generative Network for Parameterized Fluid Simulations
,”
Comput. Graph. Forum
,
38
(
2
), pp.
59
70
.10.1111/cgf.13619
11.
Xie
,
X.
,
Zhang
,
G.
, and
Webster
,
C.
,
2019
, “
Non-Intrusive Inference Reduced Order Model for Fluids Using Deep Multistep Neural Network
,”
Mathematics
,
7
(
8
), p.
757
.10.3390/math7080757
12.
Fries
,
W.
,
He
,
X.
, and
Choi
,
Y.
,
2022
, “LaSDI: Parametric Latent Space Dynamics Identification,”
Comput. Methods Appl. Mech. Eng.
, 399, p.
115436
.10.1016/j.cma.2022.115436
13.
He
,
X.
,
Choi
,
Y.
,
Fries
,
W. D.
,
Belof
,
J.
, and
Chen
,
J.-S.
,
2022
, “GLaSDI: Parametric Physics-Informed Greedy Latent Space Dynamics Identification,”
arXiv:2204.12005
[physics].10.48550/arXiv.2204.12005
14.
Ghaboussi
,
J.
,
Garrett
,
J. H.
, and
Wu
,
X.
,
1991
, “
Knowledge‐Based Modeling of Material Behavior With Neural Networks
,”
J. Eng. Mech. ASCE
,
117
(
1
), pp.
132
153
.10.1061/(ASCE)0733-9399(1991)117:1(132)
15.
Lefik
,
M.
,
Boso
,
D. P.
, and
Schrefler
,
B. A.
,
2009
, “
Artificial Neural Networks in Numerical Modelling of Composites
,”
Comput. Methods Appl. Mech. Eng.
,
198
(
21–26
), pp.
1785
1804
.10.1016/j.cma.2008.12.036
16.
Au
,
A. T. C.
, and
Kirsch
,
R. F.
,
2000
, “
EMG-Based Prediction of Shoulder and Elbow Kinematics in Able-Bodied and Spinal Cord Injured Individuals
,”
IEEE Trans. Rehabil. Eng.
,
8
(
4
), pp.
471
480
.10.1109/86.895950
17.
Wu
,
W.
,
Saul
,
K. R.
, and
Huang
,
H.
,
2021
, “
Using Reinforcement Learning to Estimate Human Joint Moments From Electromyography or Joint Kinematics: An Alternative Solution to Musculoskeletal-Based Biomechanics
,”
ASME J. Biomech. Eng.
,
143
(
4
), p.
044502
.10.1115/1.4049333
18.
Leserri
,
D.
,
Grimmelsmann
,
N.
,
Mechtenberg
,
M.
,
Meyer
,
H. G.
, and
Schneider
,
A.
,
2022
, “
Evaluation of SEMG Signal Features and Segmentation Parameters for Limb Movement Prediction Using a Feedforward Neural Network
,”
Mathematics
,
10
(
6
), p.
932
.10.3390/math10060932
19.
Ma
,
C.
,
Lin
,
C.
,
Samuel
,
O. W.
,
Xu
,
L.
, and
Li
,
G.
,
2020
, “
Continuous Estimation of Upper Limb Joint Angle From SEMG Signals Based on SCA-LSTM Deep Learning Approach
,”
Biomed. Signal Process. Control
,
61
, p.
102024
.10.1016/j.bspc.2020.102024
20.
Ren
,
J.-L.
,
Chien
,
Y.-H.
,
Chia
,
E.-Y.
,
Fu
,
L.-C.
, and
Lai
,
J.-S.
,
2019
, “
Deep Learning Based Motion Prediction for Exoskeleton Robot Control in Upper Limb Rehabilitation
,”
International Conference on Robotics and Automation (ICRA)
, Montreal, QC, Canada, May 20–24, pp.
5076
5082
.10.1109/ICRA.2019.8794187
21.
Kirchdoerfer
,
T.
, and
Ortiz
,
M.
,
2017
, “
Data Driven Computing With Noisy Material Data Sets
,”
Comput. Methods Appl. Mech. Eng.
,
326
, pp.
622
641
.10.1016/j.cma.2017.07.039
22.
Eggersmann
,
R.
,
Kirchdoerfer
,
T.
,
Reese
,
S.
,
Stainier
,
L.
, and
Ortiz
,
M.
,
2019
, “
Model-Free Data-Driven Inelasticity
,”
Comput. Methods Appl. Mech. Eng.
,
350
, pp.
81
99
.10.1016/j.cma.2019.02.016
23.
He
,
Q.
, and
Chen
,
J.-S.
,
2020
, “
A Physics-Constrained Data-Driven Approach Based on Locally Convex Reconstruction for Noisy Database
,”
Comput. Methods Appl. Mech. Eng.
,
363
, p.
112791
.10.1016/j.cma.2019.112791
24.
He
,
Q.
,
Laurence
,
D. W.
,
Lee
,
C.-H.
, and
Chen
,
J.-S.
,
2021
, “
Manifold Learning Based Data-Driven Modeling for Soft Biological Tissues
,”
J. Biomech.
,
117
, p.
110124
.10.1016/j.jbiomech.2020.110124
25.
He
,
X.
,
He
,
Q.
,
Chen
,
J.-S.
,
Sinha
,
U.
, and
Sinha
,
S.
,
2020
, “
Physics-Constrained Local Convexity Data-Driven Modeling of Anisotropic Nonlinear Elastic Solids
,”
DCE
,
1
, p.
e19
.10.1017/dce.2020.20
26.
Carrara
,
P.
,
De Lorenzis
,
L.
,
Stainier
,
L.
, and
Ortiz
,
M.
,
2020
, “
Data-Driven Fracture Mechanics
,”
Comput. Methods Appl. Mech. Eng.
,
372
, p.
113390
.10.1016/j.cma.2020.113390
27.
He
,
X.
,
He
,
Q.
, and
Chen
,
J.-S.
,
2021
, “
Deep Autoencoders for Physics-Constrained Data-Driven Nonlinear Materials Modeling
,”
Comput. Methods Appl. Mech. Eng.
,
385
, p.
114034
.10.1016/j.cma.2021.114034
28.
Bahmani
,
B.
, and
Sun
,
W.
,
2022
, “
Manifold Embedding Data-Driven Mechanics
,”
J. Mech. Phys. Solids
,
166
, p.
104927
.10.1016/j.jmps.2022.104927
29.
Raissi
,
M.
,
Perdikaris
,
P.
, and
Karniadakis
,
G. E.
,
2019
, “
Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations
,”
J. Comput. Phys.
,
378
, pp.
686
707
.10.1016/j.jcp.2018.10.045
30.
Haghighat
,
E.
,
Raissi
,
M.
,
Moure
,
A.
,
Gomez
,
H.
, and
Juanes
,
R.
,
2021
, “
A Physics-Informed Deep Learning Framework for Inversion and Surrogate Modeling in Solid Mechanics
,”
Comput. Methods Appl. Mech. Eng.
,
379
, p.
113741
.10.1016/j.cma.2021.113741
31.
He
,
Q.
, and
Tartakovsky
,
A. M.
,
2021
, “
Physics‐Informed Neural Network Method for Forward and Backward Advection‐Dispersion Equations
,”
Water Resour. Res.
,
57
(
7
), p. e2020WR029479.10.1029/2020WR029479
32.
Karniadakis
,
G. E.
,
Kevrekidis
,
I. G.
,
Lu
,
L.
,
Perdikaris
,
P.
,
Wang
,
S.
, and
Yang
,
L.
,
2021
, “
Physics-Informed Machine Learning
,”
Nat. Rev. Phys.
,
3
(
6
), pp.
422
440
.10.1038/s42254-021-00314-5
33.
Zhu
,
Q.
,
Liu
,
Z.
, and
Yan
,
J.
,
2021
, “
Machine Learning for Metal Additive Manufacturing: Predicting Temperature and Melt Pool Fluid Dynamics Using Physics-Informed Neural Networks
,”
Comput. Mech.
,
67
(
2
), pp.
619
635
.10.1007/s00466-020-01952-9
34.
Kissas
,
G.
,
Yang
,
Y.
,
Hwuang
,
E.
,
Witschey
,
W. R.
,
Detre
,
J. A.
, and
Perdikaris
,
P.
,
2020
, “
Machine Learning in Cardiovascular Flows Modeling: Predicting Arterial Blood Pressure From Non-Invasive 4D Flow MRI Data Using Physics-Informed Neural Networks
,”
Comput. Methods Appl. Mech. Eng.
,
358
, p.
112623
.10.1016/j.cma.2019.112623
35.
Sahli Costabal
,
F.
,
Yang
,
Y.
,
Perdikaris
,
P.
,
Hurtado
,
D. E.
, and
Kuhl
,
E.
,
2020
, “
Physics-Informed Neural Networks for Cardiac Activation Mapping
,”
Front. Phys.
,
8
, p. 42.10.3389/fphy.2020.00042
36.
Alber
,
M.
,
Buganza Tepole
,
A.
,
Cannon
,
W. R.
,
De
,
S.
,
Dura-Bernal
,
S.
,
Garikipati
,
K.
,
Karniadakis
,
G.
,
Lytton
,
W. W.
,
Perdikaris
,
P.
,
Petzold
,
L.
, and
Kuhl
,
E.
,
2019
, “
Integrating Machine Learning and Multiscale Modeling—Perspectives, Challenges, and Opportunities in the Biological, Biomedical, and Behavioral Sciences
,”
NPJ Digit. Med.
,
2
(
1
), pp.
1
11
.10.1038/s41746-019-0193-y
37.
Yazdani
,
A.
,
Lu
,
L.
,
Raissi
,
M.
, and
Karniadakis
,
G. E.
,
2020
, “
Systems Biology Informed Deep Learning for Inferring Parameters and Hidden Dynamics
,”
PLoS Comput. Biol.
,
16
(
11
), p.
e1007575
.10.1371/journal.pcbi.1007575
38.
Xu
,
K.
,
Tartakovsky
,
A. M.
,
Burghardt
,
J.
, and
Darve
,
E.
,
2021
, “
Learning Viscoelasticity Models From Indirect Data Using Deep Neural Networks
,”
Comput. Methods Appl. Mech. Eng.
,
387
, p.
114124
.10.1016/j.cma.2021.114124
39.
Tartakovsky
,
A. M.
,
Marrero
,
C. O.
,
Perdikaris
,
P.
,
Tartakovsky
,
G. D.
, and
Barajas-Solano
,
D.
,
2020
, “
Physics-Informed Deep Neural Networks for Learning Parameters and Constitutive Relationships in Subsurface Flow Problems
,”
Water Resour. Res.
,
56
(
5
), p. e2019WR026731.10.1029/2019WR026731
40.
He
,
Q.
,
Stinis
,
P.
, and
Tartakovsky
,
A. M.
,
2022
, “
Physics-Constrained Deep Neural Network Method for Estimating Parameters in a Redox Flow Battery
,”
J. Power Sources
,
528
, p.
231147
.10.1016/j.jpowsour.2022.231147
41.
He
,
Q.
,
Barajas-Solano
,
D.
,
Tartakovsky
,
G.
, and
Tartakovsky
,
A. M.
,
2020
, “
Physics-Informed Neural Networks for Multiphysics Data Assimilation With Application to Subsurface Transport
,”
Adv. Water Resour.
,
141
, p.
103610
.10.1016/j.advwatres.2020.103610
42.
Wang
,
S.
,
Wang
,
H.
, and
Perdikaris
,
P.
,
2021
, “
On the Eigenvector Bias of Fourier Feature Networks: From Regression to Solving Multi-Scale PDEs With Physics-Informed Neural Networks
,”
Comput. Methods Appl. Mech. Eng.
,
384
, p.
113938
.10.1016/j.cma.2021.113938
43.
Wong
,
J. C.
,
Ooi
,
C.
,
Gupta
,
A.
, and
Ong
,
Y.-S.
,
2022
, “Learning in Sinusoidal Spaces with Physics- Informed Neural Networks,”
IEEE Trans. Artif. Intell.
, pp. 1–15.10.1109/TAI.2022.3192362
44.
Sartori
,
M.
,
Reggiani
,
M.
,
Farina
,
D.
, and
Lloyd
,
D. G.
,
2012
, “
EMG-Driven Forward-Dynamic Estimation of Muscle Force and Joint Moment About Multiple Degrees of Freedom in the Human Lower Extremity
,”
PLoS One
,
7
(
12
), p.
e52618
.10.1371/journal.pone.0052618
45.
Buongiorno
,
D.
,
Barsotti
,
M.
,
Barone
,
F.
,
Bevilacqua
,
V.
, and
Frisoli
,
A.
,
2018
, “
A Linear Approach to Optimize an EMG-Driven Neuromusculoskeletal Model for Movement Intention Detection in Myo-Control: A Case Study on Shoulder and Elbow Joints
,”
Front. Neurorobot.
,
12
, p.
74
.10.3389/fnbot.2018.00074
46.
Winters
,
J. M.
,
1990
, “
Hill-Based Muscle Models: A Systems Engineering Perspective
,”
Multiple Muscle Systems
,
Springer
,
New York,
pp.
69
93
.
47.
Millard
,
M.
,
Uchida
,
T.
,
Seth
,
A.
, and
Delp
,
S. L.
,
2013
, “
Flexing Computational Muscle: Modeling and Simulation of Musculotendon Dynamics
,”
ASME J. Biomech. Eng.
,
135
(
2
), p.
021005
.10.1115/1.4023390
48.
Zhang
,
Y.
,
Chen
,
J.
,
He
,
Q.
,
He
,
X.
,
Basava
,
R. R.
,
Hodgson
,
J.
,
Sinha
,
U.
, and
Sinha
,
S.
,
2020
, “
Microstructural Analysis of Skeletal Muscle Force Generation During Aging
,”
Int. J. Numer. Methods Biomed. Eng.
,
36
(
1
), p. e3295.10.1002/cnm.3295
49.
Chen
,
J.-S.
,
Basava
,
R. R.
,
Zhang
,
Y.
,
Csapo
,
R.
,
Malis
,
V.
,
Sinha
,
U.
,
Hodgson
,
J.
, and
Sinha
,
S.
,
2016
, “
Pixel-Based Meshfree Modelling of Skeletal Muscles
,”
Comput. Methods Biomech. Biomed. Eng. Imaging Vis.
,
4
(
2
), pp.
73
85
.10.1080/21681163.2015.1049712
50.
Thelen
,
D. G.
,
2003
, “
Adjustment of Muscle Mechanics Model Parameters to Simulate Dynamic Contractions in Older Adults
,”
ASME J. Biomech. Eng.
,
125
(
1
), pp.
70
77
.10.1115/1.1531112
51.
Wang
,
S.
,
Yu
,
X.
, and
Perdikaris
,
P.
,
2022
, “When and Why PINNs Fail to Train: A Neural Tangent Kernel Perspective,”
J. Comp. Phys.
, 449, p.
110768
.10.1016/j.jcp.2021.110768
52.
Baydin
,
A. G.
,
Pearlmutter
,
B. A.
,
Radul
,
A. A.
, and
Siskind
,
J. M.
,
2018
, “
Automatic Differentiation in Machine Learning: A Survey
,”
J. Mach. Learn. Res.
,
18
, p.
43
.
53.
Tancik
,
M.
,
Srinivasan
,
P. P.
,
Mildenhall
,
B.
,
Fridovich-Keil
,
S.
,
Raghavan
,
N.
,
Singhal
,
U.
,
Ramamoorthi
,
R.
,
Barron
,
J. T.
, and
Ng
,
R.
,
2020
, “Fourier Features Let Networks Learn High Frequency Functions in Low Dimensional Domains,”
Adv. Neural. Inf. Process. Sys.
, 33, pp.
7537
7547
.https://proceedings.neurips.cc/paper/2020/file/55053683268957697aa39fba6f231c68-Paper.pdf
54.
Hermens
,
H. J.
,
Freriks
,
B.
,
Merletti
,
R.
,
Stegeman
,
D.
,
Blok
,
J.
,
Rau
,
G.
,
Disselhorst-Klug
,
C.
, and
Hägg
,
G.
,
1999
, “European Recommendations for Surface ElectroMyoGraphy,” Roessingh Res. Develop., 8(2), pp.
13
54
.
55.
Lin
,
Y.-A.
,
Zhao
,
X.
,
Silder
,
A.
,
Sessoms
,
P.
,
Fraser
,
J. J.
, and
Loh
,
K. J.
,
2021
, “
Graphene Kinesiology Tape for Monitoring Distributed Human Movements of the Ankle-Foot Complex
,”
Nano-, Bio-, Info-Tech Sens. Wear. Sys.
, 11590, pp.
57
69
.10.1117/12.2582873
56.
Kingma
,
D. P.
, and
Ba
,
J.
,
2017
, “Adam: A Method for Stochastic Optimization,”
arXiv:1412.6980
[cs].https://arxiv.org/abs/1412.6980
57.
Holzbaur
,
K. R. S.
,
Murray
,
W. M.
, and
Delp
,
S. L.
,
2005
, “
A Model of the Upper Extremity for Simulating Musculoskeletal Surgery and Analyzing Neuromuscular Control
,”
Ann. Biomed. Eng.
,
33
(
6
), pp.
829
840
.10.1007/s10439-005-3320-7
58.
Garner
,
B. A.
, and
Pandy
,
M. G.
,
2003
, “
Estimation of Musculotendon Properties in the Human Upper Limb
,”
Ann. Biomed. Eng.
,
31
(
2
), pp.
207
220
.10.1114/1.1540105
59.
An
,
K. N.
,
Hui
,
F. C.
,
Morrey
,
B. F.
,
Linscheid
,
R. L.
, and
Chao
,
E. Y.
,
1981
, “
Muscles Across the Elbow Joint: A Biomechanical Analysis
,”
J. Biomech.
,
14
(
10
), pp.
659
669
.10.1016/0021-9290(81)90048-8
60.
Hajian
,
G.
, and
Morin
,
E.
,
2022
, “
Deep Multi-Scale Fusion of Convolutional Neural Networks for EMG-Based Movement Estimation
,”
IEEE Trans. Neural Sys. Rehabil. Eng.
,
30
, pp.
486
495
.10.1109/TNSRE.2022.3153252
61.
Xiong
,
D.
,
Zhang
,
D.
,
Zhao
,
X.
, and
Zhao
,
Y.
,
2021
, “
Deep Learning for EMG-Based Human-Machine Interaction: A Review
,”
IEEE/CAA J. Autom. Sin.
,
8
(
3
), pp.
512
533
.10.1109/JAS.2021.1003865
62.
He
,
X.
,
Taneja
,
K.
,
Chen
,
J.-S.
,
Lee
,
C.-H.
,
Hodgson
,
J.
,
Malis
,
V.
,
Sinha
,
U.
, and
Sinha
,
S.
,
2022
, “
Multiscale Modeling of Passive Material Influences on Deformation and Force Output of Skeletal Muscles
,”
Int. J. Numer. Methods Biomed. Eng.
,
38
(
4
), p.
e3571
.10.1002/cnm.3571
63.
Weiss
,
K.
,
Khoshgoftaar
,
T. M.
, and
Wang
,
D.
,
2016
, “
A Survey of Transfer Learning
,”
J. Big Data
,
3
(
1
), p.
9
.10.1186/s40537-016-0043-6