Edited by: Mitsuhiro Hayashibe, University of Montpellier, France
Reviewed by: Eric Jon Perreault, Northwestern University, USA; Dingguo Zhang, Shanghai Jiao Tong University, China
*Correspondence: Ian Williams, Center for Bio-Inspired Technology, Institute of Biomedical Engineering, Imperial College London, B422 Bessemer Building, South Kensington Campus, London, SW7 2AZ, UK e-mail:
This article was submitted to Neuroprosthetics, a section of the journal Frontiers in Neuroscience.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Accurate models of proprioceptive neural patterns could 1 day play an important role in the creation of an intuitive proprioceptive neural prosthesis for amputees. This paper looks at combining efficient implementations of biomechanical and proprioceptor models in order to generate signals that mimic human muscular proprioceptive patterns for future experimental work in prosthesis feedback. A neuro-musculoskeletal model of the upper limb with 7 degrees of freedom and 17 muscles is presented and generates real time estimates of muscle spindle and Golgi Tendon Organ neural firing patterns. Unlike previous neuro-musculoskeletal models, muscle activation and excitation levels are unknowns in this application and an inverse dynamics tool (static optimization) is integrated to estimate these variables. A proprioceptive prosthesis will need to be portable and this is incompatible with the computationally demanding nature of standard biomechanical and proprioceptor modeling. This paper uses and proposes a number of approximations and optimizations to make real time operation on portable hardware feasible. Finally technical obstacles to mimicking natural feedback for an intuitive proprioceptive prosthesis, as well as issues and limitations with existing models, are identified and discussed.
A device capable of giving an amputee a sense of feeling back from their prosthetic limb could help millions of people live happier, more productive lives (Blank et al.,
Direct neural feedback in the form of a neural prosthesis has the potential to provide high quality and intuitive feedback. The incredible capability of neural prostheses to transform lives has already been vividly demonstrated in recent years by the rise of cochlear implants for the deaf and the tantalizing progress in retinal implants for the blind. A proprioceptive prosthesis on the other hand could in theory provide a user with feedback of their limb's position, motion and the forces it is exerting, as well as potentially providing therapeutic benefit for phantom limb issues (Dhillon and Horch,
The ideal for a sensory neural prosthesis would be to mimic naturally occurring neural patterns and stimulate the appropriate neurons with those patterns—providing the user with comprehensive feedback that is as intuitive as possible. However, major obstacles remain to be overcome (see section 4.1), and it seems likely that neural prostheses will rely on the brain's ability to interpret limited and abnormal feedback for some time yet.
Mimicking the function and signals of specialized neurons is an active area of focus for cochlear and retinal prostheses in order to enhance the user's ability to interpret the feedback. The brain's ability to adapt and learn is impressive, but fitting in with its pre-existing neural processing may offer better performance. However, this progression, (from simple graded stimulation to systems of modulation that mimic natural patterns) has not yet been addressed for proprioception, despite tantalizing indications that limited but appropriate neural stimulation can generate limb state representations in the brain (Weber et al.,
The aim of this paper is to create a real time model of proprioceptive signals from specific receptors to demonstrate its feasibility and to support future work investigating the possible benefits of mimicking natural signals. Our concept for developing a proprioceptive prosthesis for a transhumeral amputee is shown in Figure
In the human body the generation of proprioceptive neural signals is implicitly linked with musculoskeletal biomechanics as well as muscle and proprioceptor dynamics (Proske and Gandevia,
The integration of sensory feedback models with representations of musculoskeletal components is still a relatively new field and most publications have focused on the lower limb and locomotion. Upper limb models considering only 1 degree of freedom have previously been proposed (Lan et al.,
The addition of static optimization (an inverse dynamics tool) is proposed to estimate muscle forces and activations. However, standard implementations are computationally demanding—unsuited to the real-time, portable, and low-power nature of a proprioceptive prosthesis—and as such approximations are proposed to address this.
Numerous models of muscle spindles have been proposed in the literature [see Prochazka and Gorassini (
There are relatively fewer GTO models in the literature [see Mileusnic et al. (
This paper proposes a system for modeling ensemble average proprioceptor signals (see section 4.1 for discussion of this approach) for a simplified representation of the upper limb with 7 degrees of freedom and 17 muscles. Approximations to existing models and tools are proposed with the aim of creating a real time system capable of running on portable hardware.
The system described here is shown in Figure
The biomechanical modeling is underpinned by data from a 3D musculoskeletal model of the upper limb in OpenSim (Delp et al.,
The model covers 7 degrees of freedom in the upper limb: 3 at the shoulder (describing elevation angle, shoulder elevation and shoulder rotation), 2 at the elbow (covering elbow flexion and forearm pronation), and 2 at the wrist (covering flexion and deviation). Joint kinematics and ranges of motion were unchanged from the original model.
A standard three component dimensionless Hill type muscle model (0° pennation angle) was used and scaled to fit individual muscles as proposed by Zajac (
The torque (
In OpenSim the musculotendon length and muscle moment arms are calculated based on the muscle's origin and insertion points as well as anatomical wrapping points and constraints. However, running a 3D model is computationally intensive. A more efficient (although less accurate) approach based on fitting a polynomial surface to the length-joint angles relationship and another for the moment arm-joint angles relationship was described in van den Bogert et al. (
In order to determine the polynomial coefficients for this relationship, the OpenSim musculoskeletal model was swept through the full range of motion of the various joints and at each pose the lengths and moment arms of all the muscles were recorded. This data was then processed in Matlab with the polyfitn function to generate polynomial surfaces fitted to this data. The polyfitn function outputs the polynomial surface coefficients (
The length of the muscle (
It has been widely noted that there is redundancy in the human musculoskeletal system and hence there is typically not a unique combination of muscle forces to generate any particular motion. The situation is further complicated by the fact that muscles are often multi-articular and produce moments around each of the joints they span. Methodologies such as static or dynamic optimization [which rely on minimizing or maximizing some optimization criteria (Erdemir et al.,
Probably the simplest proposed optimization criteria is to try to minimize the total amount of muscle activation (
Minimize β subject to:
A weakness of this optimization criteria is that once a minimum β value has been calculated, the optimization process does not try to reduce muscle activations below this value, e.g., if muscle
The open source simplex package lp_solve was used to solve this linear programme around the four joints in the elbow and wrist. Given the limited subset of shoulder spanning muscles being modeled (and the target application being a transhumeral amputee with extant shoulder musculature), it was not possible to resolve the torques at the three shoulder joints (
Muscle activation dynamics were accounted for by using a method similar to that used by Thelen (
In our work we do not have the muscle excitation (
The feasible activation range was calculated for each muscle and included as constraints in the linear programme solver.
The muscle spindle outputs were simulated using a model based on that proposed by Mileusnic et al. (
The parameter “G” in the Mileusnic model is a scaling term—mapping ideal normalized spindle firing rates to feline data in the paper—and was estimated based on changes in spindle firing rates of up to 150 pulses per second (pps), that occur due to fusimotor stimulation in a feline muscle. There is limited data about the fusimotor sensitivity of human muscle spindle, but the maximum observed change in spindle output due to fusimotor signals has been observed to be <30 pps (Prochazka and Hulliger,
The GTO model used here is based on the model described by Lin and Crago (
Firstly a non-linearity:
Secondly, the output of the non-linearity is then fed into a linear dynamics transfer function:
For efficient implementation this transfer function was transformed into the z-domain in Matlab using a bilinear approximation with a sample frequency of 1 kHz and warped to fit at 6 Hz giving a z-domain transfer function of:
The focus of the work presented here is on choosing and modifying existing validated models to create a real time system. As such the approximations will be validated against the original models and the computational efficiency compared. To generate a dataset for realistic comparison, 30 s (at 120 Hz) of motion capture data from the mocapdata.com website (product_id = 15,044 showing an actor swinging his arms while walking across a room, then brushing his teeth before walking back to the original spot) was scaled, fed into OpenSim and the resulting joint angles for the upper limb were used for all simulations. This data set was chosen because it has a range of fast and slow upper limb motions and because tooth brushing represents an example of where a prosthesis user would not be able to visually monitor their limb and so feedback could provide significant benefit.
The polynomial approximation for estimating length showed close conformance with the values generated by OpenSim's 3D model throughout the dataset; giving a coefficient of determination (
Figure
To test how well the equilibrium spindle approximation corresponded to the original Mileusnic et al model, a number of the validation runs from the original paper were repeated and the results are shown in Figures
There is little data available to assess the models performance for human spindles, however, the limited validation possible did highlight a potential limitation of the Mileusnic model which is discussed below. Figure
Spindle lengths vary significantly less from muscle to muscle than muscle fasicle lengths do, this discrepancy is possible because the spindle attachment to muscle endpoints or perimysium is varied to provide consistent proprioceptive acuity across a range of joints and muscles (Proske et al.,
The system was run to generate primary afferent and GTO neural signals to enable 3rd party validation of this work, and the neural firing patterns are shown in Figure
For the purposes of this analysis the code was broken down into two parts—a deterministic part and a non-deterministic part. The deterministic part consisting of the newly written code (calculating muscle lengths, moment arms, force-velocity/length relationships, muscle spindle output, etc.), while the non-deterministic part of the code (formulating and solving the optimization linear programme) used an open source library.
The 30 s 120 Hz dataset (consisting of 3600 samples) was processed by a 2.1 GHz laptop in under 1.09 s (i.e., 27.5 times faster than real time). Profiling showed that 15% of that time was spent in the deterministic part of the code and 85% in the non-deterministic optimization code. A manual estimate of the number of instructions required to process each time sample (based on the source code) was conducted for the deterministic part of the code—yielding a value of approximately 77,000 instructions. A conservative estimate of the total number of instructions (deterministic and non-deterministic) necessary for processing each sample was made based on scaling the estimated number of instructions by the relative processing duration of the deterministic and non-deterministic parts and then multiplying by a factor of 2—this provided an estimate of 1.03 million instructions per sample.
Comparisons between the proposed spindle model and a C-code implementation of the standard Mileusnic model (using a standard euler method solver) showed that the majority of the improvement in processing speed was due to differences in the time steps that can be used (rather than reduction in the number of calculations per time step). The standard Mileusnic model can become unstable if too large a time step is chosen (experimentation showed that a maximum timestep in the order of 0.1–1 ms is required, depending on the dataset), whereas the proposed solution is stable regardless of the timestep. As such it was necessary to upsample the 120 Hz dataset used here, by a factor of 8, to obtain results with the standard model but not for the proposed model. However, the maximum timestep for the proposed model will be upper bounded by the limb position update frequency and the maximum firing frequency of the spindle—meaning that the efficiency improvement is data dependent—but in this situation was in the region of an eightfold improvement.
The fitting of cubic polynomials to calculate length and moment arm make these elements of the processing almost negligible and appears to represent a reduction in required calculations by multiple orders of magnitude compared to 3D modeling. Observed execution time for all the biomechanical modeling (length, moment arm and also the static optimization) presented here was around five orders of magnitude faster than in OpenSim, however, this is a deeply unfair comparison—pitting optimized C-code against the performance of the general purpose OpenSim package.
Providing intuitive and comprehensive feedback which is familiar or trivially easy to interpret is the ultimate goal for any neural prosthesis. However, it is unknown what the most effective and achievable format for providing proprioceptive neural feedback to prosthesis users is in the near term and there are numerous challenges in implementing, comparing and optimizing competing methods.
Stimulating a small number of neurons with a pattern that is linearly related to prosthetic limb parameters (e.g., joint angle or end effector force) is possibly the simplest approach and has been demonstrated to provide benefit in a laboratory control task (Clippinger et al.,
An alternative feedback modulation method is one that aims to approximate all the naturally occurring neural feedback patterns in the human limb. This approach is in its infancy and there is a need for physiological experimentation to verify the suitability of this approach and investigate how close to this ideal the feedback needs to be in order to demonstrate benefit compared to simpler forms of modulation, however, it seems a logical, albeit distant, ideal to aim for. Major obstacles remain such as limitations in our understanding of proprioceptors and the modulating signals from the brain, as well as our limited ability to interface with and selectively modulate large numbers of neurons. Our concept for implementing this approach breaks the process down into three steps:
Prosthetic limb properties can differ substantially from those of a human limb. For the purposes of this mapping the differences between the human and prosthesis can be grouped into three main categories: (a) physical properties—weight, moments of inertia, size and shape (including for instance the number of digits); (b) actuation properties—strength, speed, joint coupling and actuator non-linearities; and (c) kinematic properties—degrees of freedom, range of motion, axis of rotation and joint structure (including for instance joint complexes). These differences were most evident in the days when cable and hook prostheses dominated the market. However, under the twin pressures of prosthesis users' desire for cosmesis and functionality (in a world of tools and equipment designed to be operated by the human hand), there has been a strong trend toward anthropomorphic convergence. Prototype upper limbs such as the DEKA arm or commercially available prostheses such as the i-Limb, clearly demonstrate the progress that has been made toward approximating the human upper limb. Even the rise in underactuation for finger joints—which was largely driven by actuator weight considerations—moves prosthetics closer to the human form and with anthropomorphic design as a guiding principle, it is a trend that looks set to continue. This has important implications because the closer the match between prosthetic limbs and human limbs, the easier the process of mapping the state of one to the other becomes. Assuming a modeled human limb can match the states and motion of the prosthesis, then modeling the neural signals can be consdered a problem of biomechanics and proprioceptor modeling (assuming feedback is applied in the Peripheral Nervous System). Receptors in the muscles, joints and skin all sense tissue deformation and provide proprioceptive feedback to the CNS. Ideally all these receptors would be modeled for a proprioceptive prosthesis so that appropriate stimulation could be applied to any receptors interfaced with. However, in a system constrained by power, portability and (as a result) complexity, it is necessary to prioritize. Discriminating criteria include the importance of each receptor type to motor control, our ability to model the underlying tissue deformation and our understanding of the receptor firing patterns. Here we elected to focus on muscle spindles and GTOs, both of which stand out on the grounds of the quality of information they provide to the CNS and the quality of the models available for musculo-tendon and proprioception modeling. However, even for these relatively well understood receptors modeling limitations are evident such as a lack of parameters for human receptors, difficulties fitting parameterized models to different muscles and uncertainty regarding fusimotor input (see section 4.3 for further discussion). When electrodes are implanted in or around a nerve, it is unknown which neurons will be stimulated and what sensation or motor effect they can elicit. The main factors in determining which neurons are stimulated are the electrode-neuron distance, the stimulus strength and the neuronal diameter (with larger neurons recruited at lower thresholds). There is typically a trade-off between the number of neurons stimulated and the selectivity achieved; with stimulation of non-target neurons with a synchronous barrage leading to unusual or potentially noxious sensations (Smith and Leslie, The stimulation pattern to apply potentially depends on how selectively individual neurons can be targeted. Our approach is to target fascicle level selectivity because higher selectivity electrodes typically need to penetrate the perineurium which introduces a break in the blood-nerve barrier and has been observed to cause endoneurial accumulation, fibrous build up due to tissue rejection and neural damage caused by relative motion (boring at the electrode tip) (Biran et al.,
The system presented here is focused on real time prediction of neural stimulation patterns and firing rates that would be suitable for human nerve stimulation and which are based on data that could be available from a prosthetic limb. A number of models and approximations were used to achieve this aim. The cubic polynomial approximations for muscle length curves fitted the OpenSim data closely, however, the equivalent moment arm approximations showed substantially worse correlation. It was observed that the quality of a muscle's moment arm approximation decreased as the number of joints the muscle spanned increased and that the surface fit for moment arms about the shoulder joints were particularly poor (although that did not matter for this application). Considering these differences in input data and the difference in optimization criteria employed, the output of the static optimization stage showed a reasonably good match with the standard OpenSim tool. However, substantial differences were visually evident in the distribution of load between the three tricep muscle branches and this is reflected in the very low R2 values for these branches. This was likely due to the linear nature of the optimization, which becomes increasingly poor at load sharing as the number of joints and muscles increases. Examination of the results showed that simply averaging the activations of the three tricep branches would have closely fitted the OpenSim results (R2 values of 0.912, 0.970, and 0.906 for the lateral, long and medial heads respectively). In more complex systems (with greater numbers of joints and muscles) it may be necessary to compartmentalize the optimization process and run multiple iterations or implement some alternative method of sharing load. Results for bicep long and short head activation were also well below average and may be a result of the poorer moment arm fit observed for these muscle branches.
The simplified version of the spindle model closely matched the outputs of the original model for the validations proposed in the original paper as well as for some real movement data. Discrepancies were visible during rapid movements, however, given the duration of these transient differences and peak firing rates in humans of approximately 100 Hz, these discrepancies represent only a low number of missed action potentials. As mentioned in the results, the efficiency improvement provided by the proposed model appears to be largely data dependent and related to the sample frequency of the system, the maximum spindle output frequency and the maximum step size for stable solving of the differential equations in the standard model. It should be noted that the analysis here assumed a standard euler method for solving these equations, but that many alternative numerical methods exist and could improve or guarantee stability.
The proposed parameter change and adjustment to rest and threshold lengths allowed the model to estimate human spindle recordings to within a standard deviation, but without a significantly greater quantity of human data it is unclear how widely applicable these adjustments are.
The proposed approach of mimicking naturally occurring neural signals seems logical, however, physiological experimental to demonstrate and quantify benefit remains an important area of future work.
The inaccuracy of the polynomial fitting of moment arms for shoulder and highly multi-articular muscles indicates that alternative fitting models may need to be investigated to achieve good fits if the system is to be extended to include other joints.
The scaling of animal spindle and GTO models to fit human data is an area requiring further investigation. Here we have proposed simple modifications to better align the models with human firing rates, however, further data, validation and modification is required. The shape of the modeled human firing patterns was qualitatively different from the recorded data—most notably in that it lacked an initial burst or peak. Further analysis of recorded human spindle data from Cordo et al. (
In addition the Mileusnic spindle model is formulated on a limited subset of muscles and may need modification for wider applicability. In this work we noted a need for further research on how well the model copes with muscles whose optimal length is above or below the range of lengths the muscle is physiologically limited to. In our work looking at the recordings in Edin and Vallbo (
The work presented here assumed zero fusimotor activity, producing signals that are analogous to those experienced during passive motion of the limb. However, if the user is able to control their prosthesis naturally (i.e., the prosthesis responds to physiologically appropriate neural commands—e.g., following Targeted Muscle Reinnervation), then there will be a descending fusimotor signal that will act in the CNS on the pathways carrying the stimulated proprioceptive feedback, potentially interfering with it and reducing effectiveness. Ultimately it would be preferential to integrate fusimotor behavior and the Mileusnic spindle model was in part chosen to enable future implementations to easily introduce this. However, current proposed models of fusimotor action are in their infancy and largely based on recordings performed on decerebrate cats or surrogate outcomes like obtaining a linear relationship between joint angle and spindle firing (Taylor et al.,
The integration of motion capture, a musculoskeletal model, static optimization and proprioceptor models enables some proprioceptive signals to be non-invasively modeled for real movements. Further integration of inverse dynamic modeling to estimate joint torques, as well as a suitable musculoskeletal model of the feline hind limb (a preferred experimentation model), would enhance this system, allowing novel experimentation and extensive validation.
Realistic models that link human motion to proprioceptor signals could 1 day form the basis for a proprioceptive neural prosthesis in much the same way retinal and cochlear implants aim to mimic auditory and retinal cells. In contrast to previous neuro-musculoskeletal models, this work has proposed: the integration of static optimization; modifications to approximate human proprioceptors; and a variety of approximations and optimizations to reduce computational complexity without substantial degradation of the output. A key uncertainty in aiming to provide natural feeling proprioceptive feedback to a prosthesis user is how close to normal it needs to be in order to provide benefit over simpler forms of feedback modulation. This work aims to build capability to explore this question.
The model presented here is able to simulate muscle lengths, moment arms and activations as well as the corresponding muscle spindle and GTO neural signals in real time on low power hardware. This system potentially enables physiological experimentation into intuitive proprioceptive feedback as well as novel forms of proprioceptive and motor control and maybe 1 day could form part of a system capable of giving amputees feeling in their prosthetic limbs.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.