Human Neuromuscular System Identification Using Functional Electrical Stimulation for the Development of a Digital Twin of the Locomotor System

Introduction Our research group has been developing and applying a human digital twin of the locomotor system and has proposed a simple method for estimating the dynamics of the neuromusculoskeletal system using functional electrical stimulation based on the equilibrium point hypothesis, which focuses on coordination between extensor and flexor muscles. This method defines two parameters: the electrical agonist-antagonist ratio (<i>r</i> _E) and sum (<i>s</i> _E), representing the ratio and sum of the stimulation intensities applied to the extensor and flexor muscles, respectively. Our previous study showed that the relationship between <i>r</i> _E and the evoked force, i.e., the neuromuscular system (NMS), can be approximated by a second-order system with dead time under isometric conditions, and that the NMS parameters vary with <i>s</i> _E. However, this variation has not yet been modeled. This study investigates how <i>s</i> _E influences the parameters of isometric elbow joint motion with one degree of freedom. Methods Under isometric contraction, we conducted experiments to estimate the parameters of a second-order system, such as proportional gain (<i>K</i> _p), natural frequency (<i>ω</i> _n), and damping ratio (<i>ζ</i>), at 15 different <i>s</i> _E levels. Data were collected from 10 participants (nine males and one female; mean age: 22.7 ± 0.8 years; all right-handed). For group-averaged and individual data, we fitted models describing the relationship between <i>s</i> _E and each parameter. Model performance was evaluated using the corrected Akaike information criterion (AICc) across linear, quadratic, and exponential models. Results For <i>K</i> _p, the quadratic model with a concave shape best fit the group mean data as indicated by the AICc values (linear: -2.78, quadratic: -15.8, and exponential: -14.9). For <i>ω</i> _n and <i>ζ</i>, the convex quadratic models best described the group mean (for <i>ω</i> _n, linear: 3.30, quadratic: -2.74, and exponential: 3.46; for <i>ζ</i>, linear: -49.8, quadratic: -53.5, and exponential: -49.8). However, at the individual level, some participants exhibited monotonic trends. Discussion For <i>K</i> _p, although the quadratic model provided the best fit for the group mean, the exponential model showed comparable AICc values. Moreover, when summing AICc values across individuals, the exponential model yielded the lowest AICc sum, suggesting that the relationship between <i>K</i> _p and <i>s</i> _E can be reasonably approximated by an exponential function. For <i>ω</i> _n and <i>ζ</i>, the overall trend with <i>s</i> _E was best described by a convex quadratic function. However, due to interindividual differences in muscle properties, some participants did not exhibit a turning point within the tested <i>s</i> _E range, resulting in monotonic trends. These convex patterns may be explained by the influence of the refractory period of skeletal muscle fibers. Conclusions The clinical significance of the model obtained in this study lies in its potential to contribute to the development of the human digital twin of the locomotor system. By incorporating dynamics in which each parameter changes in real time with <i>s</i> _E​, it may become possible to estimate human movement from electromyographic (EMG) signals. However, because the stimulation frequency used in this study was higher than the EMG frequency, the influence of the refractory period may have been amplified. Future studies should investigate whether similar parameter trends are observed at stimulation frequencies closer to those of EMG signals.