1 Introduction

Spinal cord injury (SCI) causes paralysis of muscles due to the obstruction of the command signal from the brain to the muscles. Although there is no currently known cure for SCI, electrical stimulation of the paralyzed muscles or a peripheral nerve innervating those muscles can generate muscle contraction for the restoration of movement. One of the advantages of nerve electrodes over muscle electrodes is that a single multiple contact peripheral nerve electrode can activate many muscles innervated by the nerve.

In order to control different muscles with a single multiple contact nerve cuff electrode, selective muscle recruitment capability is desired. A flat interface nerve electrode (FINE) is a type of nerve cuff electrode designed to increase fascicular and subfascicular selectivity by reshaping a nerve trunk and realigning the fascicles inside the nerve (Tyler and Durand 2002). Several studies have shown high spatial selectivity of the FINE in both computational simulation and animal experiments (Leventhal and Durand 2004; Yoo et al. 2004; Tyler and Durand 2002; Schiefer et al. 2008). Recent intraoperative studies with an eight contact FINE on human femoral nerve also demonstrated selective muscle recruitment (Schiefer et al. 2010). However, motion control of neuromuscular systems using a FINE has not been demonstrated so far due to the complexities in neural interface in addition to the inherent complexities in the neuromuscular skeletal systems.

To date, there have been numerous approaches to motion control of FES systems. Open loop control uses predefined stimulation parameters. Since open loop control does not require feedback sensors, it is clinically preferred for its simplicity (Hoshimiya et al. 1989; Kobetic and Marsolais 1994). However, open loop control cannot detect system parameter variations such as muscle fatigue or external disturbance, and thus control performance can degrade in the presence of parameter variations or external disturbances (Zhou et al. 1998). Feedback control has benefits over open loop control because system parameter variations and external disturbances can be compensated for. A simple linear output feedback control, such as proportional integral derivative (PID) control or its modification with co-contraction mapping, showed small output errors in regulation and tracking of step and ramp reference trajectories (Chizeck et al. 1988; Crago et al. 1991; Chizeck et al. 1991; Lan et al. 1991). However, relatively large latency in the musculoskeletal systems limits the bandwidth in feedback controllers. In order to compensate for the relatively large latency in the musculoskeletal system, a feedforward controller (FFC) was adopted (Chang and Kuo 1997). In general, a FFC is an inverse model of the system to be controlled. Artificial neural networks (ANN) have been widely used to approximate an inverse model using measurable system input and output data. However, the redundancy in the musculoskeletal system control with a multiple contact nerve electrode makes it difficult to find an inverse dynamic model for control purposes (Jordan and Rumelhart 1992). The redundancy problem occurs when the number of inputs (number of contacts) is larger than the number of outputs (degree of freedom of the joints). In order to resolve the redundancy problem, optimal control methods can be used (Lan 1997; Heilman et al. 2006; Popovic et al. 1991). In the optimal control approach, the control input is obtained by minimizing a cost function of a model of the system (Popovic et al. 1999; Anderson and Pandy 2001). However, optimal control requires an accurate model of the neuromuscular skeletal system with neural–electrode interface.

In our previous study, we developed a novel non-model based motion control algorithm for musculoskeletal systems and tested it in a computational model of the human ankle joint system (Park and Durand 2008). In the method, the redundancy of the system is reduced by the separation of static and dynamic properties as originally proposed by Kawato (1993). Both the static inverse model and the dynamic inverse model of the simplified system could be obtained using only measurable input and output data with little prior knowledge of the system. However, in our previous study, the computational model did not contain the nervous system model, and it was assumed that each muscle could be activated independently at a level between no activation and full activation.

In this study, our non-model based control algorithm for FES was tested on the motion control of a computational model of the human ankle joint system with a FINE on the sciatic nerve. To achieve this, a realistic multiple level computational model was developed which can serve as a simulation platform to test the performance of various control methods using efferent nerve stimulation for FES. In the model, each individual muscle activation level was hidden to the controller, which increased the difficulty in designing a controller in addition to less than perfect individual muscle selectivity of a multiple contact nerve cuff electrode and increased nonlinearity lying in nerve fiber recruitment property. One of the major differences from the previous study in developing a controller is that the relationship between the stimulus pulse amplitude of each contact of a FINE and the corresponding twitch torque amplitude were obtained in order to find the threshold and maximum pulse amplitude for each contact. The torque simulation model was built with the same musculoskeletal parameters with restraint on the range of motion by a passive joint torque. The rest of the control design procedures are similar to those in the previous study except that the system inputs are the stimulus pulse amplitudes instead of the muscle activation levels, which are hidden variables in this study. The simulation results demonstrate that the non-model based control strategy with the separation of steady state and dynamic properties of the system resulted in small output tracking errors for different reference trajectories under various circumstances.

2 Methods

2.1 Modeling

A computational model of a human ankle joint system with the sciatic nerve was developed for this simulation study (Fig. 1). The ankle joint model contained twelve Hill-type muscles adopted from a lower extremity model (Delp 1990; Heilman et al. 2006), and it is an upgraded version of our previous model, whose details can be found in our previous paper (Park and Durand 2008). For modeling and dynamic simulation, a musculoskeletal modeling and simulation software SIMM (Software for Interactive Musculoskeletal Modeling, MusculoGraphics Inc., Santa Rosa, CA) and SD/FAST (PTC, Needham, MA) were used. The muscles in this model were medial gastrocnemius (MG), lateral gastrocnemius (LG), soleus (Sol), tibialis anterior (TA), tibialis posterior (TP), peroneus brevis (PB), peroneus longus (PL), peroneus tertius (PT), extensor digitorum longus (EDL), extensor hallucis longus (EHL), flexer digitorum longus (FDL), and flexor hallucis longus (FHL). These muscle parameters are summarized in Table 1. The inputs to the musculoskeletal system model were muscle excitation levels between 0 (no excitation) and 1 (full excitation), and the outputs were joint angular displacements; dorsiflexion/plantar flexion of ankle joint and inversion/eversion of subtalar joint. Dorsiflexion and inversion were defined as the positive direction, whereas plantar flexion and eversion were the negative direction. The neutral position was when the foot formed a right angle with the shank; in this position, both the ankle and subtalar joints had angles of 0 degrees. The ranges of motion were restricted by the passive properties of the joint, roughly from \(-40^{\circ }\) to \(30^{\circ }\) for the ankle joint and \(-20^{\circ }\) to \(20^{\circ }\) for the subtalar joint.

Fig. 1
figure 1

Computational model of the human sciatic nerve and the ankle joint system. In the model, the inputs to the system are pulse amplitudes and the outputs are ankle and subtalar joint angles. The output of the nervous system and the input to the musculoskeletal system is the muscle activation level

Full size image
Table 1 List of muscles in the ankle joint system model and their parameters
Full size table

The human sciatic nerve model was based on the histological data (Gustafson et al. 2012). The sciatic nerve divides to the common peroneal (fibular) nerve and the tibial nerve. The muscles innervated by the tibial nerve are MG, LG, Sol, TP, FDL, and FHL. The muscles innervated by the common peroneal nerve are TA, EDL, EHL, PT, PL and PB. The histological cross section of the human sciatic nerve proximal to the branching point to the tibial and the common peroneal nerves is shown in Fig. 2a. The mapping from each fascicle to the target muscle was arbitrarily assigned with the following restrictions. First, the fascicles on medial side were assigned to the muscles innervated by the common peroneal nerve (TA, EDL, EHL, PT, PL and PB), while the fascicles on the lateral side were assigned to the muscles innervated by the tibial nerve (MG, LG, Sol, TP, FDL, and FHL). Second, the total cross sectional area of fascicles assigned to the corresponding muscles was roughly proportional to the maximum muscle strength. In this simulation study, five different configurations in the mapping from the fascicles to the corresponding muscles were used in order to test the robustness of the proposed control method (Fig. 2b–f).

Fig. 2
figure 2

Human sciatic nerve modeling. a The cross section of human sciatic nerve. Medial side has fascicles connected to the common peroneal nerve and lateral side has fascicles connected to the tibial nerve (Gustafson et al. 2012). b-f Five different configurations of mapping from each fascicle to muscles (mapping #1– #5, respectively). Each fascicle is modeled as a regular dodecagon. Fascicles branching tibial nerve (MG, LG, Sol, TP, FDL, and FHL) are located lateral (left in the figure) to fascicles branching to common peroneal nerve (TA, EDL, EHL, PT, PL and PB). The numbers in each fascicle are the fascicle indexes. g Color codes for the muscles. h Realignment of fascicles with a FINE after applying realignment algorithm in Appendix A (height: 2.5 mm). The black rectangles are platinum electrode contacts with contact number from 1 to 20

Full size image

A FINE reshapes the nerve trunk and realigns the fascicles inside the nerve according to the dimension of the FINE. In order to predict the fascicular geometry of the nerve after the FINE placement, a nerve-reshaping algorithm was developed based on geometrical information. In this algorithm, the fascicular shape remains unchanged and the total area inside the nerve trunk is preserved. In addition, the vertical position of a fascicle is shifted down in proportion to the decrease in the nerve trunk height and the horizontal position of the fascicle is shifted right proportional to the increase in the nerve trunk width. In case that there is an overlap of fascicles, the involved fascicles are moved in the horizontally opposite direction in proportion to the fascicular size.

In the simulation model of the FINE, the opening height was 2.5 mm and there were 20 platinum contacts with 10 contacts on each side. Each contact size was 1 mm \(\times \) 1 mm, and the distance between two adjacent contacts was 1.6 mm. After the fascicular realignment was determined by the reshaping algorithm (Fig. 2h), a finite element method was used to find extracellular potential inside the fascicles using an electromagnetic field simulation software (Maxwell, ANSYS Inc, Canonsburg, PA). The conductivity of endoneurium, epineurium, and perineurium was adopted from a human femoral nerve model (Schiefer et al. 2008), and the perineurium thickness was 3 % of the fascicle diameter (Grinberg et al. 2008). For simplicity, the cross sectional area of each fascicle was modeled as a regular polygon with 12 sides and all materials in this model were assumed to be purely resistive. Since the DC conduction model is linear, the extracellular potential obtained using a unit current (1 mA) can be linearly scaled to find the extracellular potential for arbitrary current amplitudes.

The axons inside the fascicles were evenly distributed with \(50\, \upmu \hbox {m}\) distance between the adjacent axons in the transversal plane. The positions of the nodes of Ranvier in relation to the electrode in the longitudinal direction were randomly selected with uniform random distribution. The axonal diameters followed a Gaussian distribution between \(10\,\upmu \hbox {m}\) and \(20\,\upmu \hbox {m}\), and a single cable model was adopted for the axonal dynamics model (Leventhal and Durand 2004). The axon model has 15 nodes, and NEURON software (Carnevale and Hines 2006) was used to run simulation with cathodic pulses with pulse width of \(50\,\upmu \hbox {s}\). Each contact was assumed to generate stimulus pulse in sequence, and the activated axons were assumed to be within a refractory period when the neighboring contacts generate stimulus pulses. It was also assumed that the time-shifted summation of subthreshold extracellular potential generated by stimulation of neighboring contacts was not able to generate suprathreshold extracellular potential. Therefore, axonal activation can be simply calculated by comparing the input stimulus amplitude of each contact to the threshold amplitude for the axon. The threshold pulse amplitude at each contact for each axon were found using a binary search algorithm with a resolution of \(0.1\,\upmu \mathrm{A}\).

The outputs of the nervous system model were the muscle excitation levels, which were the inputs to the musculoskeletal model. The muscle excitation level of each muscle was assumed to be the weighted sum of excited axons over the total axons as in the following equation.

$$\begin{aligned} e_i =\frac{\sum _{j\in {\text {Activated}}} {A_{i,j} } }{\sum _{j=1}^{N_i } {A_{i,j} } } \end{aligned}$$
(1)

In this equation, \(e_{i}\) is the excitation level of \(i\)th muscle, and \(A_{i,j}\) is the cross sectional area of \(j\)th axon innervating \(i\)th muscle. \(N_{i}\) is the total number of axons innervating \(i\)th muscle. As a result, the nervous system model was simplified as the mapping from the stimulation amplitudes to the muscle excitation levels.

For the muscle activation–deactivation dynamic model between the muscle excitation and the muscle activation, a first order dynamics with default activation time constant (13 ms) and deactivation time constant (69 ms) of SIMM was used (Thelen 2003). The stimulation frequency was set to 25 Hz, and thus the muscle excitation level changed every 40 ms. The stimulation frequency was selected because it is within typical stimulation frequency for FES system (Sheffler and Chae 2007), and is high enough to reduce ripples in muscle force response for pulse train stimulation. In summary, in the current ankle joint system model with a FINE on the sciatic nerve, the inputs to the complete model were the stimulation amplitudes and the outputs of the model were ankle and subtalar joint angles. During the motion control simulation, the internal variables such as the muscle excitation levels were hidden to the controller, and only the input stimulation amplitudes and the output joint angles were accessible by the controller.

Individual muscle selectivity using a FINE on the sciatic nerve was also tested in order to characterize the system. Individual muscle selectivity was defined as the muscle excitation level of each muscle either without activating other muscles (0 % cost) or allowing maximum 10 % activation of the other muscles (10 % cost).

2.2 Controller design

As described in our previous work (Park and Durand 2008), the controller is composed of an inverse steady state controller (ISSC), a FFC and a feedback controller (FBC) (Fig. 3a). The system here is the ankle joint system model with a FINE on the sciatic nerve, where the inputs are the stimulus pulse amplitudes at each contact and the outputs are ankle/subtalar joint angles. ISSC is an inverse model of the system at steady states, where the inputs are the desired ankle/subtalar joint angles, and the outputs are the stimulus pulse amplitudes. The FFC is a dynamic inverse model of the combination of the system and ISSC in series so that the combination of FFC and ISSC becomes a dynamic inverse model of the system. A PID controller was used for the feedback controller to compensate for external disturbances and system parameter variations. Details about the controller design can be found in our previous paper (Park and Durand 2008). In this study, for simplification, only type I controller of our previous work is discussed, where the FBC output is fed to the system through ISSC. Similar results are expected for type II controller with proper selection of FBC, where the FBC output is directly fed to the system in parallel to ISSC output.

Fig. 3
figure 3

Controller structure. a The controller is composed of Inverse steady state controller (ISSC), artificial neural network (ANN) feedforward controller (FFC), and PID feedback controller (FBC). \(\mathbf{r}[k]\) is the desired ankle and subtalar angles, \(\mathbf{y}[k]\) is the actual ankle and subtalar joint angles, \(\mathbf{v}[k]\) is the input to ISSC, \(\mathbf{v}_\mathrm{ff}[k]\) is the feedforward controller output, \(\mathbf{v}_\mathrm{fb}[k]\) is the feedback controller output, and \(\mathbf{u}[k]\) is the input to the neuromuscular skeletal system at time step \(k\). \(\mathbf{r}[k], \mathbf{y}[k], \mathbf{v}[k], \mathbf{v}_\mathrm{ff}[k]\) and \(\mathbf{v}_\mathrm{fb}[k]\) are vectors with 2 elements (ankle and subtalar joint angle). The number of the elements of \(\mathbf{u}[k]\) is the same as the number of contacts. Since, the future desired trajectory \(\mathbf{r}[k+1], {\ldots }, \mathbf{r}[k+p]\) were used for the feedforward controller, the desired trajectory should be determined p steps in prior. b Training of the feedforward controller is done offline. During the training, the time series of the system output becomes the input to the feedforward controller, and the input to ISSC becomes the desired output of the feedforward controller. The feedforward controller is trained to minimized the error

Full size image

2.2.1 Inverse steady state controller (ISSC)

ISSC is a static inverse model of the system composed of a lookup table whose entries are evenly spaced ankle and subtalar joint angles, which are the inputs to ISSC, and their corresponding stimulus pulse amplitudes, which are the outputs of ISSC. The steady state values were chosen as the pulse amplitudes and their corresponding output angles at 2 s after the initiation of stimulation (25 Hz). The convergence of the joint angles for constant inputs is described in our previous study (Park and Durand 2008). Compared to our previous study where the input to the system was the muscle excitation levels between 0 and 1, the incorporation of the nervous system model required determining lower bound and upper bound of the stimulation amplitudes of each contact. In order to find these lower and upper bounds, instead of measuring the motion response to the stimulation, twitch torque responses at a fixed position were used. Measuring the twitch torque response has several benefits over measuring the motion responses. First, it requires less time since single pulse stimulations are used. Second, the balancing maximum torque for each direction can be attainable, which is comparable to maximum muscle activation level in our previous study (Park and Durand 2008). While the threshold amplitudes were used for the lower bound, the upper bound could be selected more flexibly based on monotonicity and balance of the torque responses. In addition, the electrode contacts were categorized into four functional groups according to their corresponding primary torque responses: plantar flexion (PF), dorsiflexion (DF), inversion (INV) and eversion (EV). In order to generate motion covering the entire range, at least one contact should be selected from each functional group. Otherwise, the range of motion can be limited. Once the contacts, threshold amplitudes and maximum amplitudes were determined, steady state input-output data to build an ISSC were obtained by repetitively applying the interpolated or extrapolated stimulus parameters to the neuromuscular skeletal system and measuring the ankle and subtalar joint output angles. For the interpolation and extrapolation of the stimulus amplitudes, the generalized barycentric method (Floater et al. 2006) was used in the output joint coordinate as described in the following equations.

$$\begin{aligned}&\lambda _i (\mathbf{v})\ge 0,\quad i=1,2,\ldots ,n \nonumber \\&\sum _{i=1}^n \lambda _i (\mathbf{v}) =1 \nonumber \\&\sum _{i=1}^n \lambda _i (\mathbf{v}) \mathbf{v}_\mathbf{i} = \mathbf{v} \nonumber \\&\mathbf{u}=\sum _{i=1}^n \lambda _i (\mathbf{v}) \mathbf{u}_\mathbf{i}, \end{aligned}$$
(2)

where v is the joint angles at which the stimulation amplitudes are to be obtained, and \(\mathbf{u}_{\mathbf{i}}\) and \(\mathbf{v}_{\mathbf{i}}\) are the stimulus pulse amplitudes and their corresponding joint angles of the selected steady state data for interpolation or extrapolation. \(n\) is the number of data used for the interpolation or extrapolation. For interpolation, nearest four data points were used, and for extrapolation, the convex hull of the existing data were used. When enough number of steady state data points covering the entire output joint coordinate space of interest were obtained, the ISSC table was built by obtaining the corresponding stimulus amplitudes for the target joint angles in the entries of the ISSC table using interpolation or extrapolation. During control stage, bilinear interpolation was used to calculate stimulation amplitudes when the desired output angles did not exactly match the entries in the lookup table.

2.2.2 Feedforward controller (FFC) and Feedback controller (FBC)

The FFC was implemented with a time-delayed multilayer perceptron (MLP) network with one hidden layer of 30 nodes to model the inverse dynamics of the combination of the system and ISSC in series. A significant reduction in redundancy due to the ISSC made it possible to apply the direct inverse training method to train the FFC. In the direct inverse training method, the system output trajectory becomes the input trajectory to the FFC and the input trajectory to the ISSC becomes the desired output trajectory of FFC (Fig. 3b). For the training of the ANN, the Levenberg-Marquardt algorithm in MATLAB was used. A PID controller was adopted for the feedback controller. The PID gain matrices were roughly tuned by trial and error beginning from small gain.

3 Results

3.1 Individual muscle selectivity

The ability of the FINE to recruit individual muscles was first tested. The muscle recruitment selectivity of each of 12 muscles in 5 different configurations is shown in Fig. 4. The average individual muscle selectivity was \(66\pm 42\,\%\,(\hbox {mean}\pm \hbox {s.d.}\)) for 0 % cost (no activation of untargeted muscles), and \(73\pm 38\,\%\) at 10 % cost (allowing maximum 10 % activation of untargeted muscles). The results indicate large variation in the individual muscle selectivity depending on the fascicular configuration inside the nerve and electrode contact geometry.

Fig. 4
figure 4

Individual muscle selectivity of FINE. a Individual muscle selectivity in five different mappings for 0 % cost (gray) and for 10 % cost (gray + black). 0 % cost means no other muscle activation except the target muscle, and 10 % cost means allowance of untargeted muscle activation up to 10 %. The individual muscle selectivity depends highly on the mapping configuration. b The individual muscle selectivity averaged over the five different mapping configurations. The averaged muscle selectivities for 0 % cost and 10 % cost are \(66\pm 24\,\%\) and \(73\pm 21\,\%\) (\(\hbox {mean}\pm \hbox {s.d.}\)), respectively. These results indicate large variations in the individual muscle selectivity

Full size image

3.2 Steady state error of the inverse steady state controller (ISSC)

First, the twitch torque responses at neutral position were measured by applying single pulse stimulation with the amplitude from \(0\, \upmu \hbox {A}\) to \(200\,\upmu \hbox {A}\) with the resolution of \(5\, \upmu \hbox {A}\) for each of 20 contacts (Fig. 5a). The twitch torque response indicates that the contacts can be categorized into four functional groups (DF, PF, INV and EV). Then the lower bound and upper bound of the amplitude for each contact were determined (Fig. 5b). The thresholds, maximum amplitudes and primary torque direction are summarized in Table 2. For each of five mapping, three different subsets of electrode contacts were tested in building the ISSCs to demonstrate the robustness of the controller with respect to the contact selection. The first and second combinations used only one contact from each functional contact group, and the third combination used two contacts. The simulation results hereafter include all 15 configurations (5 fascicle-muscle mappings \(\times \) 3 contact sets).

Table 2 Thresholds, maximum amplitudes (in parentheses), in \(\mu A\) and primary torque direction for each contact in 5 different mappings. The four main directions are dorsiflexion (DF), plantar flexion (PF), inversion (INV), and eversion (EV)
Full size table
Fig. 5
figure 5

Twitch torque response for single pulse stimulation. a The twitch torque amplitude is plotted for the pulse amplitudes from \(0\, \upmu \hbox {A}\) up to \(200\, \upmu \hbox {A}\) with \(5\,\upmu \hbox {A}\) resolution for each electrode contact from 1 to 20 (Mapping #1). b An example of the twitch torque response for contact 3, where the lower bound and upper bound were chosen to be \(45\,\upmu \hbox {A}\) and \(70\,\upmu \hbox {A}\). c The number of contacts in each primary torque direction: dorsiflexion (DF), plantar flexion (PF), inversion (INV) and eversion (EV). d The average lower and upper bounds for each primary torque direction

Full size image

An example of the steady state input–output data and the corresponding ISSC are shown in Fig. 6. The average steady state errors for 15 ISSCs at 651 lookup table entry points (31 for ankle joint \(\times \) 21 for subtalar joint) were \(1.2^{\circ }\pm 1.0^{\circ }\) and \(1.4^{\circ } \pm 1.2^{\circ }\) for ankle and subtalar joints, respectively. The maximum average steady state errors occurred for mapping #2 and contact set #2, and they were \(1.9^{\circ }\pm 2.4^{\circ }\) and \(1.8^{\circ }\pm 1.9^{\circ }\) for the ankle and subtalar joint, respectively. These results indicate the steady state error by ISSC is small regardless of the selection of contacts and fascicle configurations.

Fig. 6
figure 6

Example of steady state data and ISSC. FFC An example of steady state data displayed in the output joint coordinate (Mapping #1 and Contact set #1). Each unfilled circle represents the stimulus pulse amplitudes (not shown in the figure), and their corresponding ankle and subtalar joint angles. The steady state input–output data were obtained by iterative interpolation process. ISSC table is obtained by interpolating the steady state data to find the corresponding stimulus pulse amplitudes for the entries of the ISSC table. b ISSC table obtained from the steady state data in a. The inputs to the ISSC are the desired ankle and subtalar joint angles, and the outputs of the ISSC are the corresponding stimulus pulse amplitudes for each contact to obtain the desired joint angles. In this ISSC example, the amplitudes at the other contacts except contact 2, 8, 15 and 19 are zero

Full size image

3.3 Controller performance and the effects of FFC and FBC

The performance of the proposed controller was tested with 20 sinusoidal reference trajectories with different frequencies, amplitudes and phases. These reference trajectories had frequency ranging between 0.5 and 1.0 Hz, and their amplitudes and phases were randomly selected. In order to determine the effect of each of the components of the controller on the overall performance, the three different control structures were compared using sinusoidal signals: ISSC alone, ISSC with FFC, and ISSC with feedforward and feedback controller. In the simulation, random sensor noise between \(-0.5^{\circ }\) and \(0.5^{\circ }\) and random input noise between -\(0.5\,\upmu \hbox {A}\) and \(0.5\,\upmu \hbox {A}\) were also added to approximate real environment. The average ankle and subtalar joint root mean square (RMS) tracking errors for 15 different configurations and for 20 sinusoidal reference trajectories were \(12.6^{\circ }\pm 6.3^{\circ }\) and \(9.3^{\circ }\pm 4.4^{\circ }\) for ISSC alone, \(3.6^{\circ }\pm 2.6^{\circ }\) and \(3.4^{\circ }\pm 2.6^{\circ }\) for ISSC with FFC, \(1.7^{\circ }\pm 0.9^{\circ }\) and \(1.9^{\circ }\pm 1.3^{\circ }\) for the proposed controller (ISSC+FFC+FBC). Statistical analysis demonstrated that the proposed controller showed better trajectory tracking performance than ISSC alone or ISSC with FFC both in the ankle and the subtalar joint angles (\(\hbox {p}<0.001\); one-way ANOVA with Bonferroni correction). In addition, ISSC with FFC demonstrated smaller output errors than ISSC alone (Fig. 7). An example of the proposed controller performance results is shown in Fig. 8. The output error was the smallest when both the feedforward and feedback controller were used. The small output errors by the proposed controller indicate that the controller could find an inverse model of the system to be controlled by separating static and dynamic properties.

Fig. 7
figure 7

Results for sinusoidal reference trajectories. Performance of the different controller configurations is shown for 20 sinusoidal reference trajectories between 0.5 and 1.0 Hz. The ISSC+FFC+FFB demonstrated the least output errors

Full size image
Fig. 8
figure 8

Example of tracking sinusoidal reference trajectories. Performance of the combined controller for 1.0 Hz sinusoidal reference trajectory is demonstrated (Mapping #1 and Contact set #1). The thin blue line is the command trajectory and the thick red line is the actual output trajectory. The RMS errors for ankle and subtalar joint angles are \(3.4^{\circ }\) and \(2.3^{\circ }\), respectively. The pulse amplitudes at the contacts are shown in the bottom figure

Full size image

3.4 Effects of the system parameter variations

The robustness of the ISSC+FBC+FFC method was tested for the system parameter variations such as muscle fatigue. In this simulation, the maximum muscle strength of each muscle was decreased maximally to 50 % from its original value and 20 different combinations of maximum muscle strength change were tested using 20 s pseudo-random reference trajectory. The average RMS output errors increased slightly, compared with the performance on unmodified muscle strengths: from \(1.0^{\circ }\pm 0.2^{\circ }\) and \(1.1^{\circ }\pm 0.3^{\circ }\) to \(1.3^{\circ }\pm 0.4^{\circ }\) and \(1.5^{\circ }\pm 0.9^{\circ }\), respectively. The performance of the proposed controller did not degrade significantly with the muscle parameter changes because feedback control compensated for the system parameter variation.

3.5 Effects of external disturbances

In order to test the robustness of the ISSC+FBC+FFC method in the presence of external disturbance, sinusoidal external force was applied at the point 20 cm from the subtalar joint along the foot. Two different maximum magnitudes were used to test the controller performance in the presence of external disturbance: 10 N and 25 N. A filtered pseudo-random trajectory in the previous section was used as a command signal.

The output RMS errors were \(1.1^{\circ }\pm 0.2^{\circ }\) and \(1.1^{\circ }\pm 0.4^{\circ }\) for the ankle and subtalar joint, respectively (\(\hbox {mean}\pm \hbox {s.d.}\)) without external disturbance. The output RMS errors for the sinusoidal disturbance with maximum magnitude of 10 N were \(1.9^{\circ }\pm 0.5^{\circ }\) and \(1.2^{\circ }\pm 0.4^{\circ }\) for the ankle and subtalar joint, respectively. The RMS errors for a 25 N disturbance were \(3.7^{\circ }\pm 1.0^{\circ }\) and \(1.7^{\circ }\pm 0.7^{\circ }\) for the ankle and the subtalar joint, respectively. An example of the effect of external disturbance on the combined controller performance is shown in Fig. 9. Although the error increased due to the external disturbance, the average RMS output error was less than 10 % of the range of the motion. The results indicate that the feedback controller could effectively compensate the external disturbance to minimize the output error.

Fig. 9
figure 9

Example of tracking filtered random trajectories with external disturbance. Performance of the combined controller for filtered random trajectories in the presence of sinusoidal external disturbance is demonstrated (Mapping #1 and Contact set #3). The RMS errors for the ankle and the subtalar joint angles are \(2.6^{\circ }\) and \(0.8^{\circ }\) for 25 N disturbance, respectively. The bottom graph shows sinusoidal external disturbances with the frequency of 0.5 Hz applied at the point 20 cm from the subtalar joint along the foot

Full size image

4 Discussion

The computational neuromuscular skeletal system model described above is the first implementation (to the best of our knowledge) used for trajectory tracking motion control simulation study with multiple scale levels of axons, nerve cuff electrode, and musculoskeletal implementation of the human ankle joint system and the sciatic nerve. Due to the difficulties in obtaining an accurate model for control using multiple contact nerve electrodes, the stimulation parameters for control were usually obtained by trial and error (Mushahwar et al. 2007).

An ideal FFC is a dynamic inverse model of the system to be controlled. In the absence of an accurate system model, finding an inverse model using only measurable input and output data is very difficult, especially when the system has redundancy. One method to reduce the system redundancy is to separate static inverse model and dynamic inverse model (Kawato et al. 1993). In order to find an inverse steady state model of a system with redundancy, a novel method was presented utilizing iterative interpolation. Since our method only requires measurable input and output data, it is suitable for neuromuscular skeletal system control with a multiple contact nerve electrode.

In this study, a total of 15 different configurations were used to demonstrate the robustness of the control method: 5 mappings \(\times \) 3 contact sets. The individual muscle selectivity test showed that the average individual muscle selectivity was \(66\pm 41\,\%\) without activating any other muscles and increased up to \(75\pm 37\,\%\) when the other muscles were allowed to be activated up to 10 %. The individual muscle selectivity level is close to the results in the human femoral nerve simulation study (Schiefer et al. 2008). The large variation in the selectivity is due to the fact that each fascicular selectivity depends on the distance between the electrode contact and the fascicle.

The simulation results for sinusoidal reference signals (\(0.5\, \hbox {Hz}\sim 1.0\,\hbox {Hz}\)) showed good trajectory tracking performance. The results indicate that the controller could appropriately approximate an inverse dynamic model of the system to be controlled without resorting to an analytical model of the system even with less than perfect individual muscle selectivity of the nerve cuff electrode. In the worst case of all 15 configurations, the average ankle joint angle error was \(3.2^{\circ }\pm 3.1^{\circ }\), which is a little above 5 % of the range of motion of the ankle joint. These results demonstrate the robustness of the algorithm in finding ISSCs and the FFC.

The controller performance was also tested in the presence of the system parameter variations and external disturbance. Although the output error increased in these conditions, the FBC compensated for the system changes and disturbances. However, since the FBC output was applied to the system through the ISSC, the error compensation was limited by the ISSC. Direct output feedback to the system can be used, but in that case, the large number of PID gains should be obtained: A difficult task with a system that is nonlinear and strongly coupled (Park and Durand 2008).

In developing the neuromuscular skeletal model, several assumptions were made. Since there was no available information about the mapping between the fascicles inside the sciatic nerve to the muscles, five different hypothetical mappings were used. In addition, a geometrical algorithm was used to reshape the nerve and relocate the fascicles inside the nerve trunk when a FINE was hypothetically placed on the sciatic nerve. There were also assumptions regarding the axonal distribution inside each fascicle and the muscle activation level proportional to the cross sectional area of each axon. Since the position of the nodes and diameters are unknown, the simulations were carried out with random assignment of these variables.

In the ankle joint system model presented above, the foot was not supported, and therefore, ground reaction force was not considered. For the application to walking, the ground reaction force should also be considered and the system characteristics including the joint stiffness will be drastically changed due to the body weight. In addition, the coupling between the joints may affect the controller performance. In order to accommodate different situations, multiple ISSCs may need to be built to account for variations in the system and environment. In addition, the incorporation of torque control in addition to motion control using joint stiffness information may enhance the accuracy of the gait control. Although a dynamic inverse model of a redundant multiple- input multiple-output (MIMO) system with unknown dynamics could be obtained using the proposed method, as the system output dimension increases, the number of trials in building an ISSC increases exponentially. In such a system with large output dimension, accurate control may not be possible without an accurate model of the system. Another limitation in the proposed control method is that either ISSC or FFC does not account for the system change over time. Although in this simulation study, the FBC could compensate for the maximum muscle strength change very well, the actual neuromuscular skeletal system change may be different from the simple maximum muscle strength change; for example, muscle recruitment property changes. Therefore, adaptive ISSC or adaptive FFC may improve the controller performance in physiological systems (Kim et al. 2009). The main limitation of this simulation study is the lack of experimental validation to support the findings. Although rabbit ankle joint motion control experiments with one degree of freedom demonstrated satisfactory reference trajectory tracking performance (Park and Durand 2009), higher degree of freedom experiments on animal models or humans will validate the control algorithm for the practical implementation for FES systems in the future.

In conclusion, a novel non-model based control algorithm was presented for the human ankle joint system control with a multiple contact nerve cuff electrode on the sciatic nerve. For the simulation, a complete neuromuscular skeletal system model was developed. The controller could track different reference trajectories even with system parameter variations and external disturbances. This control technology could be useful for FES system using nerve stimulation.