Skip to content
Publicly Available Published by De Gruyter October 30, 2014

Computational models of upper-limb motion during functional reaching tasks for application in FES-based stroke rehabilitation

Chris Freeman, Tim Exell, Katie Meadmore, Emma Hallewell and Ann-Marie Hughes


Functional electrical stimulation (FES) has been shown to be an effective approach to upper-limb stroke rehabilitation, where it is used to assist arm and shoulder motion. Model-based FES controllers have recently confirmed significant potential to improve accuracy of functional reaching tasks, but they typically require a reference trajectory to track. Few upper-limb FES control schemes embed a computational model of the task; however, this is critical to ensure the controller reinforces the intended movement with high accuracy. This paper derives computational motor control models of functional tasks that can be directly embedded in real-time FES control schemes, removing the need for a predefined reference trajectory. Dynamic models of the electrically stimulated arm are first derived, and constrained optimisation problems are formulated to encapsulate common activities of daily living. These are solved using iterative algorithms, and results are compared with kinematic data from 12 subjects and found to fit closely (mean fitting between 63.2% and 84.0%). The optimisation is performed iteratively using kinematic variables and hence can be transformed into an iterative learning control algorithm by replacing simulation signals with experimental data. The approach is therefore capable of controlling FES in real time to assist tasks in a manner corresponding to unimpaired natural movement. By ensuring that assistance is aligned with voluntary intention, the controller hence maximises the potential effectiveness of future stroke rehabilitation trials.


Stroke is a leading cause of adult disability in the EU, and over 6 million stroke patients require care. Of these, two-thirds have impairment of their affected arm 4 years post-stroke, resulting in an annual cost of 38 billion [33]. Approximately 70% of patients will experience altered arm function after a stroke, and about 40% of survivors will be left with a non-functional arm that is weak and often spastic [24]. This motor deficit limits functional arm use in daily life, but also engagement in the community [50]. Residual impairments and functional deficit in a large percentage of stroke patients indicate that current rehabilitation provision is unsatisfactory [45].

In recent years, there has been growing evidence supporting the effectiveness of rehabilitation robots [30] and functional electrical stimulation (FES) [26] to reduce impairment post-stroke. Both technologies enable a person with limited physical ability to practice tasks, and the resulting sensory feedback is associated with cortical changes that can bring about recovery of functional movement. In particular, there is substantial clinical evidence [8–10, 26] indicating that increased functional recovery is closely related to the accuracy with which FES assists the subject’s own voluntary completion of a task. This finding also has theoretical support from neurophysiology [4, 39] and motor learning research [40].

Model-based FES control is critical to reduce the effect of noise/disturbance, increase accuracy, and enable complex functional tasks to be performed. However, inherent environmental constraints mean that few model-based controllers have transferred to clinical practice [54], where typically open-loop or triggered controllers are employed and resulting performance is limited. This is despite a wide variety of FES upper-limb control techniques having been applied in simulation or laboratory conditions.

Artificial neural networks (ANNs) have been extensively used for upper-limb FES control, especially within neuroprostheses for subjects with spinal cord injury. Here ANNs have been used to represent the inverse dynamics of single-joint, single-muscle systems by creating a mapping from joint kinematics to the required muscle activation. In the case of complex shoulder and arm movements, the inverse dynamics are not unique, and hence sophisticated musculoskeletal models have been developed, such as [5, 11, 47], and combined with optimal cost functions to distribute forces between redundant muscles [3, 18]. ANNs are then trained to approximate the inverse arm dynamics and then used to compute muscle activation levels in real time. This structure is typically combined with an appropriate FES feedback controller to ensure the muscle activation is achieved.

Iterative learning control (ILC) is another model-based approach that has been employed clinically. It has been used in upper-limb stroke rehabilitation [17], to assist movement in the lower limb [1, 2, 32, 41], as well as in the wider biomedical field for drug delivery and heart rate monitoring. When used in FES-based stroke rehabilitation, ILC exploits the repetitive nature of the rehabilitation process where patients attempt the same task multiple times in order to promote re-learning of motor skills. ILC precisely assists the patient’s movement in each attempt of the task and sequentially improves accuracy using data from previous attempts to adjust the FES supplied during the next attempt at the task. By exploiting experience in this way, ILC is able to embed both accuracy and robustness to model uncertainty [17].

ILC update algorithms work by directly using functional mappings that approximate the inverse system dynamics. ILC hence operates in a similar fashion to the aforementioned model-based ANN approaches but implements this mapping directly without employing ANNs. This potentially restricts the available model structure, but has the advantage of not requiring any training procedure, exploiting learning from experience to improve performance, enabling strict control over the amount of FES supplied during each attempt in order to maximise the patient’s volitional effort, and giving rise to transparent analysis of performance, convergence, and robustness to model uncertainty. ILC has yielded a high level of clinical performance in multiple clinical trials [23, 28, 29], most recently assisting functional activities such as closing a drawer, switching on a light switch, stabilising an object, button pressing, and repositioning an object. The setup used in a recent clinical treatment trial is shown in Figure 1.

Figure 1 The components of the GO-SAIL system: (1) control algorithm hardware (including Microsoft Kinect and PrimeSense) and software; (2) touch table virtual and real tasks; (3) FES and multiplexor (4) SaeboMAS® arm support; (5) surface electrode array on forearm, and single electrode pads on anterior deltoid and triceps muscles.

Figure 1

The components of the GO-SAIL system: (1) control algorithm hardware (including Microsoft Kinect and PrimeSense) and software; (2) touch table virtual and real tasks; (3) FES and multiplexor (4) SaeboMAS® arm support; (5) surface electrode array on forearm, and single electrode pads on anterior deltoid and triceps muscles.

A drawback to previous use of ILC is that it requires joint reference trajectories to follow, with these provided from unimpaired subject data. However, since optimisation-based ILC works by iteratively minimising a cost function using experimental data, there is clear scope for it to automatically generate such trajectories if a suitable model of unimpaired movement is embedded within it. This expansion in the scope of ILC has been facilitated by recent extensions to the ILC framework, such as [14, 15, 38], in which more general classes of constrained optimisation problem are addressed. There is hence potential for ILC to precisely support natural motion if models of human movement can be expressed as suitable optimisation problems.

Many studies have been reported characterising human motion during upper-limb reaching tasks. The majority extract relationships between key variables (e.g., timing and amplitude of kinematic, kinetic or electromyographic data) in order to examine the effect of task conditions and/or participant groups on task execution. These include effect of age [55], task conditions [48], and compensatory strategies post stroke [52]. For example [7, 42, 44], find significant inter-patient differences in the performance of reaching tasks following stroke compared with unimpaired participants. Human sensorimotor control has also been expressed computationally to more fully capture the underlying dynamics of movement [22, 51]. Approaches can be divided into those that attempt to simulate the internal feedback/feedforward mechanisms present in the central nervous system, and those that try only to model the resulting kinematic motion at the task level. The latter have traditionally posed reaching tasks as optimisation problems, involving, for example, the minimisation of jerk [13], torque change [46], variance [19], interaction torques, or combinations of these [34]. The form of an optimisation holds the possibility of application within FES-based control strategies such as ILC, however, the focus has so far overwhelmingly been on planar point to point tasks. The case of functional tasks necessary to complete activities of daily living has not yet been addressed. The first step in modelling functional tasks is the selection of an appropriate underlying musculoskeletal model.

Musculoskeletal models linking FES and resultant motion have been employed in the design of FES control systems by several authors [12, 18, 21, 25]. However, the majority are purely employed for offline analysis of control performance and behaviour. An exception is the simulation environment of Chadwick et al. [5], which is designed for non-invasive real-time simulation and user feedback. One musculoskeletal model that has been embedded within a real-time upper-limb FES control scheme is that of Freeman et al. [16], which was used within clinical trials with stroke participants. This form has been shown to accurately capture the dynamics of movement in stroke while also yielding computationally tractable solutions.

This paper extends the upper-limb model of Freeman et al. [16] to incorporate 7-degree-of-freedom (dof) shoulder and elbow motion. It then develops a framework in which functional tasks are represented as optimal objective functions, together with additional constraints acting upon variables within the model. The resulting optimisation can be solved using iterative optimisation approaches to yield computational models of functional movement, together with corresponding FES signals. Moreover, it can also be solved using experimental data within the ILC update structure, thereby enabling high-precision functional task completion to maximise effectiveness of stroke rehabilitation.

In addition to enabling more accurate assistance, the development of computational models of natural movement also provides measures of stroke patients’ movement that correspond closely to functional ability. This paper thereby addresses the significant gap, highlighted in [42], that exists between performance measures based on kinematic and kinetic data (generated by, e.g., robotic therapy systems), and clinical measures based on functional ability (e.g., Fugl-Meyer, action research arm test).



Following University of Southampton, Faculty of Health Sciences ethical approval, 12 unimpaired volunteers were recruited to the study. Inclusion criteria were that the participants had to be 30–80 years old, able to comply with study protocol, able to communicate effectively, and able to provide written informed consent. Exclusion criteria were the requirement of an interpreter, uncorrected visual impairment, a skin disease or allergy to sticky tape, severe pain in the arm, shoulder, or hand, and, for the control participants, a neurological condition that affects movement in the arm. Participants were recruited and tested over a 2-month period.

The unimpaired participants (six men and six women) were aged between 49 and 77 years (mean 64, SD 10). All unimpaired participants, except for one, were right-handed. The side tested for each participant was randomised; six participants completed the tasks using their right hand, and the other six using their left hand. Five participants were tested using their dominant hand. The testing side for unimpaired participants was randomised so as to be more representative of the stroke population, which are almost equally distributed between left- and right-sided incidence [20].

Experimental protocol

All participants attended one testing session lasting 2–3 h in which the kinematic movements of the upper limb and hand during three functional reaching tasks were recorded. The tasks were closing a drawer, turning off a light switch, and picking up a can to drink from. These tasks were selected as they offered ecologically valid reach to grasp tasks of varying demands and were similar to tasks used in many studies investigating stroke patients’ movements [6, 43]. Participants performed tasks both freely and when supported by an exoskeleton unweighing system (ArmeoSpring, Hocoma, Zurich, Switzerland); however, only data collection and results for the former case are presented in this paper.

Position data were recorded using a Vicon MX T-Series motion capture system (Vicon, Oxford, UK) using 12 cameras (6× T40 and 6× T160) sampling at 100 Hz. Reflective anatomical markers were positioned on key landmarks of the torso, shoulder complex, upper limb, wrist, and hand, as shown in Figure 2. Marker clusters were attached to the sternum and acromion on the side that was being tested, with additional markers placed on the radial and ulnar styloids and the second and fifth carpometacarpal (CMC) and metacarpophalangeal (MCP) joints of the hand. A marker wand was used to locate specific anatomical landmarks with respect to the relevant marker cluster (see [49]). These additional bony landmarks consisted of the sternal notch, xyphoid process, C7 vertebra, T8 vertebra, sternoclavicular joint, acromioclavicular joint, acromion angle of the scapula, medial spine of the scapula, inferior angle of the scapula, and medial and lateral epicondyles of the elbow. The elbow joint centre was estimated as the midpoint between the medial and lateral epicondyles. The glenohumeral joint centre was defined according to the regression method of Meskers et al. [31]. Additional markers were positioned on the task objects to aid movement identification. Four markers were positioned in a square on the front edge of the drawer and light switch. Markers were located so that they did not inhibit participants’ movements during the tasks.

Figure 2 (A) Positioning of Vicon marker clusters used during trials and (B) additional marker wand locations used for calibration.

Figure 2

(A) Positioning of Vicon marker clusters used during trials and (B) additional marker wand locations used for calibration.

Participants were seated at a table that was adjusted so that the underside was 10 cm above their knee. For each task, the participant was asked to start with their hand (palm down) on their knee. On a “start” command, the participant completed the task and then placed their hand back on their knee. Tasks were performed at both self-selected and maximal speeds. Five successful trials were collected for each task and each speed; trials were repeated if any of the reflective markers were occluded for more than 25 epochs during the trial. Participants were given a 15- to 30-s break between each trial. All participants completed the drawer closing task first. Maximum reach for each participant was measured from the anterior edge of the acromion (AA) to the end of the index finger.

Drawer closing task

A custom-made cabinet that had a drawer with a large, round, central knob was placed on the table in front of the participant so that the drawer knob was directly in line with the participant’s shoulder for the side being tested. The cabinet was placed at a distance corresponding to 100% of the participant’s maximum reach, with the drawer knob directly in line with the participant’s shoulder for the side being tested. When opened, the drawer knob was at a distance of 75% of the participant’s maximum reach. Participants were asked to move their hand from their knee to push the drawer closed using the knob and to return their hand back to their knee.

Light switch task

On the opposite side of the cabinet, a standard light switch had been mounted. The cabinet was positioned so that the light switch was directly in line with the participant’s shoulder, at 75% of their maximum reach. Participants were asked to move their hand from their knee to turn off the light switch and to return their hand back to their knee. Participants were always required to push in the top of the light switch, which required more control and stability of the arm than pressing the bottom of the switch.

Data analysis

Kinematic data, collected using the Vicon system, were reconstructed using Vicon Nexus (Version 1.8) software to provide three-dimensional position data for all the markers. Marker position data were filtered using a fourth-order low-pass Butterworth filter, with a typical cutoff frequency of 10 Hz. Anatomical landmarks, recreated from their known positions with respect to the marker clusters, were used to define local coordinate systems for the thorax, scapular, humorous, and ulna following ISB guidelines [53]. The hand was defined with an origin at the midpoint between the second and fifth CMC markers, the Y axis was defined as being parallel with the line formed by the midpoint of the CMC markers to the midpoint of the MCP markers pointing dorsally, the Z axis parallel to the line from the second CMC marker to the fifth CMC marker pointing laterally (with the hand supinated), with the X axis orthogonal to the Y and Z axes.

Positional data were averaged at each time point across the five repetitions, for each participant, task, and speed. Key timings were then extracted, comprising the start and end of the movement, defined by the initial hand movement from the participant’s knee and the hand returning to the knee, respectively. For the light task, an additional timing was when the light switch was pressed. For the drawer task, two additional timings were at the start and end of the drawer movement.

Computational model development

A simplified biomechanical model of the electrically stimulated human arm is shown in Figure 3 and is an extension of the planar model [16] validated with both unimpaired and stroke subjects. Parameters li and ai denote link and centre of mass lengths, respectively, for i∈{1, 2, 3}, and mi and Ii are mass and inertias, respectively. Joint angles θi, i∈{1, 2, 3}, relate to a 3-dof shoulder joint, θ4 to a 1-dof elbow joint, θ5 to forearm rotation, and θi, i∈{6, 7}, corresponds to wrist flexion/extension, abduction/adduction. This rigid body representation has been shown to adequately capture the dynamic properties of the upper limb, while giving rise to optimisations that are computationally tractable (a key issue when subsequently used in real-time control) [16]. Note that the experimental results reported in this paper omit joint angles θi, i∈{5, 6, 7}; however, they are included here for completeness.

Figure 3 (A) Unimpaired participant performing task (with exoskeletal support) and (B) 7 dof model of the upper limb.

Figure 3

(A) Unimpaired participant performing task (with exoskeletal support) and (B) 7 dof model of the upper limb.

The dynamics are represented by

(1)B(θ(t))θ¨(t)+C(θ(t),θ˙(t))θ˙(t)+F(θ(t),θ˙(t))+G(θ(t))=τ(u(t),θ(t),θ˙(t))-JΤ(θ(t))h(t), (1)

where the joint angle vector θ(t)=[θ1(t), θ2(t), θ3(t), θ4(t), θ5(t), θ6(t), θ7(t)]T and B(·) and C(·) are the 7×7 inertial and Corelis matrices, respectively. In addition, J(·) is the system Jacobian, h(t) is the externally applied force/torque, and F(·) and G(·) are friction and gravitational vectors. The general form Fi(θ(t),θ˙(t))=Fi(θi(t),θ˙i(t)) is used and identification procedures for the terms in Eq. 1 can be found in [16]. It is assumed that n muscles are involved in the movement, with elements of vector u(t) specifying the innervation input applied to each (either through voluntary action or application of electrical stimulation). The ith element of the muscle torque vector τ(·) is the sum of moments generated by each of these muscles, which each may impart a moment about the ith joint. From [27], it is assumed that each can be represented as a Hill-type model of the form

(2)τi,j(uj(t),θi(t),θ˙i(t))=hi,j(uj(t),t)×Fm,i,j(θi(t),θ˙i(t)),j=1,,n,i=1,,7, (2)

where hi,j(uj(t), t) is a Hammerstein structure incorporating a static non-linearity, hIRC,i,j(uj(t)), representing the isometric recruitment curve, cascaded with linear activation dynamics, hLAD,i,j(t). Here hLAD,i,j(t)=L-1{hLAD,i,j}(s), where L denotes the Laplace transform. The term Fm,i,j(θi(t),θ˙i(t)) models the multiplicative effect of the joint angle and joint angular velocity on the active torque developed by the muscle. It is assumed for simplicity that hIRC,i,j(·)=hIRC,i(·) and hLAD,i,j(·)=hLAD,j(·) ∀i, so that hi,j(·)=hj(·) ∀i. This leads to the structure

(3)τi(u(t),θi(t),θ˙i(t))=j=1nτi,j(uj(t),θi(t),θ˙i(t))=j=1n{hi,j(uj(t),t)×Fm,i,j(θi(t),θ˙i(t))}=j=1n{hj(uj(t),t)×Fm,i,j(θi(t),θ˙i(t))}. (3)

Now let hLAD,j(t), have continuous time state-space model matrices {Am,j, Bm,j, Cm,j} (state, input, and output, respectively), and states xm,j(t), so that

(4)τi(u(t),θi(t),θ˙i(t))=j=1n{Cm,jxm,jhj(uj)Fm,i,j(θi(t),θ˙i(t))}. (4)

System 1 can then be expressed in state-space form as the first-order ordinary differential equation:

(5)x˙(t)=[θ˙(t)B(θ(t))-1X(θ(t),θ˙(t))Am,1xm,1Am,nxm,n]f(x(t))+[00Bm,1hIRC,1(u1(t))Bm,nhIRC,n(un(t))]g(u(t)),θ(t)=[I00]x(t)h(x(t)),θ(0)=θ0,t[0,T], (5)

where x(t)=[θ(t)Τ,θ˙(t)Τ,xm,1(t)Τxm,n(t)Τ]Τ, and X(θ(t),θ˙(t)) has elements

(6)Xi(θ(t),θ˙(t))=j=1n{Cm,jxm,jFm,i,j(θi(t),θ˙i(t))}-Ci(θ(t),θ˙(t))θ˙(t)-Fi(θ(t),θ˙(t))-Gi(θ(t))+(JΤ(θ(t)))ih(t). (6)

A variety of models of human movement have been proposed in the literature, and here the minimum input energy is considered [22, 51]. The problem of coming to rest at the light switch at time t=T is hence expressed as

(7)minu||u||2,suchthat{y(T)=k(θ(T))=py˙(T)=(ddt(k(θ(t))))t=T=J(θ(T))θ˙(T)=0 (7)

subject to the dynamics (5) over t∈[0, T]. Here y=k(·) is the direct kinematics equation of the system, with y the hand position in Cartesian space and p is the light switch position in Cartesian space. Problem 7 can be solved through iterative optimisation methods, which are based on linear approximations of the system dynamics [35]. Given an operating point (u¯,x¯), define

(8)A(t)=(ddxf(x(t)))x¯(t),B(t)=(ddug(u(t)))u¯(t), (8)
(9)C(t)=(ddxk(h(x)))x¯(t)=[J(θ¯(t))00],t[0,T], (9)

then the relationship between position, y(t), and stimulation, u(t), signals over t∈[0,T] can be approximated by the linear time-varying system

(10)y=g1(x¯,u¯)u:z˙(t)=A(t)z(t)+B(t)u(t)y(t)=C(t)z(t),t[0,T]. (10)

Likewise, the relationship between velocity, y˙(t), and stimulation, u(t), can be approximated by

(11)y˙=g2(x¯,u¯)u:z˙(t)=A(t)z(t)+B(t)u(t)y˙(t)=C(t)A(t)z(t),t[0,T]. (11)

The linear operators g1(x¯,u¯):L2n[0,T]3:uy(T) and g2(x¯,u¯):L2n[0,T]3:uy˙(T) in Eqs. 10 and 11 are defined as

(12)g1(x¯,u¯)u:=0TC(T)Φ(T,τ)B(τ)u(τ)dτ,g2(x¯,u¯)u:=0TC(T)A(t)Φ(T,τ)B(τ)u(τ)dτ, (12)

where Φ(t, τ) is the state transition matrix for systems 10 and 11. Around operating point (u¯,x¯), the optimisation problem (7) can therefore be replaced by the linear counterpart

(13)minu||u||2, such that {g1(x¯,u¯)u=pg2(x¯,u¯)u=0 (13)

The associated Lagrangian of Eq. 13 is hence

(14)L(u)=||u||2+2<λ1,g1(x¯,u¯)u-p>+2<λ2,g2(x¯,u¯)u>, (14)

with multipliers λ1, λ2. Here the inner product is defined by

(15)<u1,u2>=0Tu1Τ(t)u2(t)dt, with <u1,u1>=||u1||2. (15)

It follows that the relevant conditions for a stationary point (u, λ1, λ2) are

(16)g1(x¯,u¯)u-p=0,g2(x¯,u¯)u=0,u+g1(x¯,u¯)λ1+g2(x¯,u¯)λ2=0, (16)

where the adjoint operators g1(x¯,u¯):3L2n[0,T] and g2(x¯,u¯):3L2n[0,T] are computed as

(17)w=g1(x¯,u¯)v:z˙(t)=-AΤ(t)z(t)w(t)=BΤ(t)z(t),z(T)=CΤ(T)v,t[0,T] (17)


(18)w=g2(x¯,u¯)v:z˙(t)=-AΤ(t)z(t)w(t)=BΤ(t)z(t),z(T)=AΤ(t)CΤ(T)v,t[0,T], (18)

respectively. From Eq. 16, the solution to Eq. 13 is given by

(19)u=g12(x¯,u¯)(g12(x¯,u¯)g12(x¯,u¯))1[p0], (19)

where the compound operators are defined by

(20)g12(x¯,u¯)u:=[g1(x¯,u¯)ug2(x¯,u¯)u],g12(x¯,u¯)(v1v2):=g1(x¯,u¯)v1+g2(x¯,u¯)v2 (20)

To solve the original problem (7), the solution (19) is employed in an iterative update procedure that corresponds to an implementation of the Newton method approach. This necessaritates modifying Eq. 19 to take the form

(21)uk+1=uk+αg12(xk,uk)(g12(xk,uk)g12(xk,uk))-1([p0]-[yk(T)y˙k(T)]), (21)

where α∈(0,1] is a scalar gain chosen to affect a compromise between robustness and convergence speed. The final update procedure is given in Table 1. In this paper, the model (5) is used in simulation within steps (c) and (d) to model voluntary task completion. However, the procedure in Table 1 can equally be employed experimentally to control FES applied to the real human arm by redefining uk+1 to be the FES input vector on trial k+1, with element uk+1,i(t) the stimulation signal applied to the ith muscle at time t. This requires replacing the associated components of the muscle model (3) by identified parameters that capture the dynamics of applied FES [27] and replacing the simulation outputs yk, y˙k by their experimentally recorded alternatives within step (d). Note that the procedure of Table 1 is similar to ILC approaches in [14, 15, 36–38] and inherits many of their favourable robustness and convergence properties. The problem definition for the drawer closing task is given in Appendix A.1, and extensions to include the return component of both tasks are given in Appendix A.2 and A.3.

Table 1

Update procedure to solve light switch task.

a.Determine task parameters T, p, and parameters within model (5).
b.Set k=0 and choose an initial input u0=0.
c.Apply input uk to the system.
d.Record the outputs yk(T), y˙k(T), and states xk.
e.Stop if, for sufficiently small positive scalars δ, ε, ||uk-uk1||<δ, and [p0]-[y(T)y˙(T)]<ε.
f.Linearise the system about (u¯,x¯)=(uk,xk) using Eqs. 9–11, 17, and 18.
g.Calculate update uk+1 using Eq. 21.
h.Increment k and go to step (c).


To examine how accurately the motor control models developed represent the experimentally collected unimpaired movement, it is necessary to perform the modelling procedure set out in Table 1. The first step is to extract parameters appearing within the arm dynamic model (5) and task descriptions (7) and (22). The former uses estimated arm lengths, masses, and inertias, with non-conservative forces taking the simple decoupled form Fi(θi(t),θ˙i(t))=biθ˙i(t)+ki(θi(t)-θ˜i), where bi, ki, and set-point θ˜i are scalars for the joint angles i=1, …, 4 employed in the model.

Task parameters are defined by the placement of the participant and manipulated objects. The initial joint variables are taken from the experimentally recorded data set, which is averaged at each time point over the five repetitions of each task to yield θ^. In the case of the light switch task, the time taken to press the switch is calculated and used in optimisation (7) to model the reaching component. Likewise, p is set as the light switch position.

Performing the procedure of Table 1 in simulation with initial condition θ0=θ^(0) produces the simulated joint angle signal θ. Results expressing the difference in experimental and simulated joint angles appear in Table 2 for the self-selected speed where the normalised error, expressed as a percentage, is ||θ^-θ||||θ^-θ0||×100. For the reach component of the light switch task at self-selected speed, the mean fitting across the 12 participants is 71.437%, confirming high accuracy. The procedure has been repeated for the maximal speed light switch task, with fitting results shown in Table 2. The fitting is now 75.892%. Figure 4A and B shows the signals θ and θ^ (solid and dotted lines, respectively) at both speeds for a single participant (P1). Figure 5A shows the corresponding paths in Cartesian space for the reach component of the maximal speed light switch task. These figures confirm accurate fitting of the experimentally recorded signals.

Table 2

Light switch pressing task fitting error.

Maximal speedSelf-selected speed
Reach onlyReachReturnReach and returnReach onlyReachReturnReach and return
Figure 4 Simulated and experimental joint angles for light switch task (A) maximal speed reach component, (B) self-selected speed reach component, (C) maximal speed reach and return components, (D) self-selected speed reach and return components.

Figure 4

Simulated and experimental joint angles for light switch task (A) maximal speed reach component, (B) self-selected speed reach component, (C) maximal speed reach and return components, (D) self-selected speed reach and return components.

Figure 5 Simulated and experimental paths in Cartesian space for maximal speed reach component of the (A) light switch reach task and (B) drawer closing task.

Figure 5

Simulated and experimental paths in Cartesian space for maximal speed reach component of the (A) light switch reach task and (B) drawer closing task.

Fitting results for the full reach and return light switch task are shown in Table 2 at both speeds, with values also given for the reach and return sub-components of the task. The mean fitting for the reach and return light switch task across the 12 participants is now 69.166% and 64.626% for maximal and self-selected speeds, respectively. Figure 4C and D shows the signals θ and θ^ (solid and dotted lines, respectively) at both speeds for participant P1.

For the drawer closing task, the time when the participant first makes contact with the drawer knob and the time when the drawer is fully closed are directly extracted from the experimental data. In addition, p is set equal to the position of the drawer knob when fully closed, and D^ and d define the constraint imposed by the drawer’s runners. These are used in Eq. 22 to provide a model of the reach component (which includes the drawer closing action). Fitting results are shown in Table 3 and confirm a mean accuracy of 77.194% and 83.961% for the self-selected and maximal speeds, respectively. Figure 6A and B shows the experimental and simulated joint angles at both speeds for participant P1. Figure 5B shows the corresponding paths in Cartesian space for the reach component of the maximal speed drawer closing switch task. These results again confirm that the model can accurately fit the experimental data.

Table 3

Drawer closing task fitting error.

Maximal speedSelf-selected speed
Reach onlyReachReturnReach and returnReach onlyReachReturnReach and return
Figure 6 Simulated and experimental joint angles for drawer closing task: (A) maximal speed reach component, (B) self-selected speed reach component, (C) maximal speed reach and return components, (D) self-selected speed reach and return components.

Figure 6

Simulated and experimental joint angles for drawer closing task: (A) maximal speed reach component, (B) self-selected speed reach component, (C) maximal speed reach and return components, (D) self-selected speed reach and return components.

Mean fitting results for the reach and return drawer closing task are shown in Table 2 and confirm a total accuracy of 69.417% and 73.698% for the self-selected and maximal speeds, respectively. Figure 6C and D shows the signals θ and θ^ (solid and dotted lines, respectively) at both speeds for participant P1.


Fitting results confirm that unimpaired human motion can be accurately predicted using a simple underlying model of the upper limb, together with a suitable optimisation procedure. This means that functional motion can be posed as optimisation problems and embedded in FES control schemes to enable clinically relevant tasks to be assisted. Furthermore, since a reference trajectory is no longer required, this can be achieved without the need to collect data for each new movement (as would be necessary to provide joint reference trajectories). Moreover, the patient’s own estimated dynamic parameters in the model mean the task is tailored to each patient. The close link with volitional control this facilitates is anticipated to lead to more effective FES-based rehabilitation.

Corresponding author: Chris Freeman, Electronics and Computer Science, University of Southampton, Southampton SO17 1BJ, UK, Phone: +442380593486, Fax: +442380593709, E-mail:


This work is supported by the Engineering and Physical Sciences Research Council grant number EP/I01909X/1.

Appendix A. Additional operator derivations

A.1 Drawer closing task

In the case of the drawer closing task, Eq. 7 is replaced by

(22)minu||u||2,such that {Dy=dy(T)=k(θ(T))=py˙(T)=(ddt(k(θ(t))))t=T=J(θ(T))θ˙(T)=0 (22)

The operator D:L23[0,T]L21[T1,T2] is defined by Dy(t)=(Dy)(t)=D^y(t)=D^k(θ(t)),t∈[T1, T2] with D^ a 1×3 matrix with full row rank. Matrix D^ embeds the tracking requirement D^y(t)=d on [T1, T2], which is used to ensure that the hand follows a straight line in Cartesian space (representing the draw knob trajectory). To solve Eq. 22, introduce the linear operator

(23)Dg3(x¯,u¯)u:=D^0tC(t)Φ(t,τ)B(τ)u(τ)dτ,t[T1,T2] (23)

so that Eq. 13 is replaced by

(24)minu||u||2, such that {g1(x¯,u¯)u=pg2(x¯,u¯)u=0Dg3(x¯,u¯)u=d (24)

and the input solution (19) is replaced by

(25)u=g123(x¯,u¯)(g123(x¯,u¯)g123(x¯,u¯))-1[p0d] (25)

with the operators

(26)g123(x¯,u¯)u:=[g12(x¯,u¯)uDg3(x¯,u¯)u],g123(x¯,u¯)(v1v2):=g12(x¯,u¯)v1+g3(x¯,u¯)Dv2. (26)

The iterative update (21) used in step (g) of Table 1 is replaced by

(27)uk+1=uk+αg123(xk,uk)(g123(xk,uk)g123(xk,uk))-1([p0d]-[yk(T)y˙k(T)D^yk]) (27)

Again, this can either be applied to the arm model (5), or experimentally, to control FES applied to a patient. In the latter case, yk and y˙k are replaced by the corresponding experimentally recorded signals.

A.2 Light switch task return component

To model the return component of the light switch task, Eq. 7 must be extended by redefining the linear operators as g1(x¯,u¯):u(y(T)y(T1)), and g2(x¯,u¯):u(y˙(T)y˙(T1)) using

(28)g1(x¯,u¯)u:=[0TC(T)Φ(T,τ)B(τ)u(τ)dτ0T1C(T1)Φ(T1,τ)B(τ)u(τ)dτ], (28)
(29)g2(x¯,u¯)u:=[0TC(T)A(T)Φ(T,τ)B(τ)u(τ)dτ0T1C(T1)A(T1)Φ(T1,τ)B(τ)u(τ)dτ], (29)

and replacing p by (pk(θ^(T1))). Here T is the time taken to press the light switch, and T1 is the time taken for the hand to return to the starting position.

A.3 Drawer closing task return component

To model the return component of the drawer closing task, Eq. 22 must again be extended by redefining the linear operators as g1(x¯,u¯):u(y(T)y(T3)), and g2(x¯,u¯):u(y˙(T)y˙(T3)) using Eq. 29.


[1] Aubin PM, Cowley MS, Ledoux WR. Gait simulation via a 6-DOF parallel robot with iterative learning control. IEEE Trans Biomed Eng 2008; 55: 1237–1240.10.1109/TBME.2007.908072Search in Google Scholar

[2] Bae J, Tomizuka M. A gait rehabilitation strategy inspired by an iterative learning algorithm. Mechatronics 2012; 22: 213–221.10.1016/j.mechatronics.2012.01.009Search in Google Scholar

[3] Blana D, Kirsch RF, Chadwick EK. Combined feedforward and feedback control of a redundant, nonlinear, dynamic musculoskeletal system. Med Biol Eng Comput 2009; 47: 533–542.10.1007/s11517-009-0479-3Search in Google Scholar

[4] Burridge JH, Ladouceur M. Clinical and therapeutic applications of neuromuscular stimulation: a review of current use and speculation into future developments. Neuromodulation 2001; 4: 147–154.10.1046/j.1525-1403.2001.00147.xSearch in Google Scholar

[5] Chadwick EK, Blana D, van den Bogert AJ, Kirsch RF. A real-time, 3-D musculoskeletal model for dynamic simulation of arm movements. IEEE Trans Biomed Eng 2009; 56: 941–948.10.1109/TBME.2008.2005946Search in Google Scholar

[6] Cirstea M, Ptito A, Levin M. Arm reaching improvements with short-term practice depend on the severity of the motor deficit in stroke. Exp Brain Res 2003; 152: 476–488.10.1007/s00221-003-1568-4Search in Google Scholar

[7] Cirstea MC, Levin MF. Compensatory strategies for reaching in stroke. Brain 2000; 123: 940–953.10.1093/brain/123.5.940Search in Google Scholar

[8] de Kroon JR, IJzerman MJ. Electrical stimulation of the upper extremity in stroke: cyclic versus EMG-triggered stimulation. Clin Rehabil 2008; 22: 690–697.10.1177/0269215508088984Search in Google Scholar

[9] de Kroon JR, IJzerman MJ, Chae JJ, Lankhorst GJ, Zilvold G. Relation between stimulation characteristics and clinical outcome in studies using electrical stimulation to improve motor control of the upper extremity in stroke. J Rehabil Med 2005; 37: 65–74.10.1080/16501970410024190Search in Google Scholar

[10] de Kroon JR, van der Lee JH, IJzerman MJ, Lankhorst GJ. Therapeutic electrical stimulation to improve motor control and functional abilities of the upper extremity after stroke: a systematic review. Clin Rehabil 2002; 16: 350–360.10.1191/0269215502cr504oaSearch in Google Scholar

[11] Delp SL, Loan JP. A graphics-based software system to develop and analyze models of musculoskeletal structures. Comput Biol Med 1995; 25: 21–34.10.1016/0010-4825(95)98882-ESearch in Google Scholar

[12] Ferrarin M, Palazzo F, Riener R, Quintern J. Model-based control of FES-induced single joint movements. IEEE Trans Neural Syst Rehabil Eng 2001; 9: 245–257.10.1109/7333.948452Search in Google Scholar

[13] Flash T, Hogan N. The coordination of arm movements: an experimentally confirmed mathematical model. J Neurosci 1985; 5: 1688–1703.10.1523/JNEUROSCI.05-07-01688.1985Search in Google Scholar

[14] Freeman CT. Constrained point-to-point iterative learning control with experimental verification. Control Eng Pract 2012; 20: 489–498.10.1016/j.conengprac.2012.01.003Search in Google Scholar

[15] Freeman CT, Tan Y. Iterative learning control with mixed constraints for point-to-point tracking. IEEE Trans Control Syst Technol 2012; 21: 604–616.Search in Google Scholar

[16] Freeman CT, Hughes AM, Burridge JH, Chappell PH, Lewin PL, Rogers E. A model of the upper extremity using FES for stroke rehabilitation. ASME J Biomech Eng 2009; 131: 031006-1–031006-10.10.1115/1.3005332Search in Google Scholar

[17] Freeman CT, Rogers E, Hughes AM, Burridge JH, Meadmore KL. Iterative learning control in healthcare: electrical stimulation and robotic-assisted upper limb stroke rehabilitation. IEEE Control Syst Mag 2012; 32: 18–43.10.1109/MCS.2011.2173261Search in Google Scholar

[18] Giuffrida JP, Crago PE. Functional restoration of elbow extension after spinal-cord injury using a neural network-based synergistic FES controller. IEEE Trans Neural Syst Rehabil Eng 2005; 13: 147–152.10.1109/TNSRE.2005.847375Search in Google Scholar

[19] Harris CM, Wolpert DM. Signal-dependent noise determines motor planning. Nature 1998; 20: 780–784.Search in Google Scholar

[20] Hedna VS, Bodhit AN, Ansari S, et al. Hemispheric dfferences in ischemic stroke: is left-hemisphere stroke more common? J Clin Neurol 2013; 9: 97–102.10.3988/jcn.2013.9.2.97Search in Google Scholar

[21] Heilman BP, Kirsch RF. Model-based development of a user control algorithm for postural control via a FES-based standing neuroprosthesis. In: Conf Proc IEEE Eng Med Biol Soc 2004; 4614–4617.Search in Google Scholar

[22] Huang VS, Krakauer JW. Robotic neurorehabilitation: a computational motor learning perspective. J Neuroeng Rehabil 2009; 6.10.1186/1743-0003-6-5Search in Google Scholar

[23] Hughes AM, Freeman CT, Burridge JH, Chappell PH, Lewin P, Rogers P. Feasibility of iterative learning control mediated by functional electrical stimulation for reaching after stroke. J Neurorehab Neural Repair 2009; 23: 559–568.10.1177/1545968308328718Search in Google Scholar

[24] Intercollegiate Working Party for Stroke. National clinical guidelines for stroke 2008. Technical report. London: Royal College of Physicians, 2008.Search in Google Scholar

[25] Kirsch R, Acosta A, van der Helm F, Rotteveel R, Cash L. Model-based development of neuroprostheses for restoring proximal arm function. J Rehabil Res Dev 2001; 38: 619–626.10.1109/IEMBS.2001.1019748Search in Google Scholar

[26] Langhorne P, Coupar F, Pollock A. Motor recovery after stroke: a systematic review. Lancet Neurol 2009; 8: 741–754.10.1016/S1474-4422(09)70150-4Search in Google Scholar

[27] Le F, Markovsky I, Freeman CT, Rogers E. Identification of electrically stimulated muscle models of stroke patients. Control Eng Pract 2010; 18: 396–407.10.1016/j.conengprac.2009.12.007Search in Google Scholar

[28] Meadmore KL, Exell T, Hallewell E, et al. The application of precisely controlled functional electrical stimulation to the shoulder, elbow and wrist for upper limb stroke rehabilitation: a feasibility study. J Neuroeng Rehabil 2014; 11: 105.10.1186/1743-0003-11-105Search in Google Scholar

[29] Meadmore KL, Hughes A-M, Freeman CT, et al. Functional electrical stimulation mediated by iterative learning control and 3d robotics reduces motor impairment in chronic stroke. J Neuroeng Rehabilil 2012; 32: 1–11.10.1186/1743-0003-9-32Search in Google Scholar

[30] Mehrholz J, Platz T, Kugler J, Pohl M. Electromechanical and robot-assisted arm training for improving arm function and activities of daily living after stroke. Cochrane Database Syst Rev 2008; 4: CD006876.10.1002/14651858.CD006876.pub2Search in Google Scholar

[31] Meskers C, van der Helm F, Rozendaal L, Rozing P. In vivo estimation of the glenohumeral joint rotation center from scapular bony landmarks by linear regression. J Biomech 1998; 31: 93–96.10.1016/S0021-9290(97)00101-2Search in Google Scholar

[32] Nahrstaedt H, Schauer T, Hesse S, Raisch J. Iterative learning control of a gait neuroprosthesis. Automatisierungs 2008; 56: 494–501.10.1524/auto.2008.0728Search in Google Scholar

[33] Nichols M, Townsend N, Luengo-Fernandez R, et al. European cardiovascular disease statistics 2012. Technical report, European Heart Network, Brussels. Sophia Antipolis: European Society of Cardiology, 2012.Search in Google Scholar

[34] Ohta K, Svinin MM, Luo Z, Hosoe S, Laboissiere RL. Optimal trajectory formation of constrained human arm reaching movements. Biol Cybern 2004; 91: 23–36.10.1007/s00422-004-0491-5Search in Google Scholar

[35] Ortega JM, Rheinboldt WC. Iterative solution of nonlinear equations in several variables. London: Academic Press, 1970.Search in Google Scholar

[36] Owens DH, Freeman CT, Chu B. Multivariable norm optimal iterative learning control with auxiliary optimization. Int J Control 2013; 86: 1026–1045.10.1080/00207179.2013.771822Search in Google Scholar

[37] Owens DH, Freeman CT, Chu B. An inverse model approach to multivariable norm optimal iterative learning control with auxiliary optimization. Int J Control 2014; 87: 1646–1671.10.1080/00207179.2014.880951Search in Google Scholar

[38] Owens DH, Freeman CT, Van Dinh T. Norm-optimal iterative learning control with intermediate point weighting: theory, algorithms, and experimental evaluation. IEEE Trans Control Syst Technol 2013; 21: 999–1007.10.1109/TCST.2012.2196281Search in Google Scholar

[39] Rushton DN. Functional electrical stimulation and rehabilitation – an hypothesis. Med Eng Phys 2003; 25: 75–78.10.1016/S1350-4533(02)00040-1Search in Google Scholar

[40] Schmidt RA, Lee TD. Motor control and learning: a behavioural emphasis. 4th ed. Leeds, UK: Human Kinetics Europe, 2005.Search in Google Scholar

[41] Seel T, Schauer T, Raisch J. Iterative learning control for variable pass length systems. In: IFAC World Congress 2011: 4880–4885.10.3182/20110828-6-IT-1002.02180Search in Google Scholar

[42] Thies SB, Tresadern PA, Kenney LP, et al. Movement variability in stroke patients and controls performing two upper limb functional tasks: a new assessment methodology. J Neuroeng Rehabil 2009; 6: 2.10.1186/1743-0003-6-2Search in Google Scholar

[43] Thrasher TA, Zivanovic V, McIlroy W, Popovic MR. Rehabilitation of reaching and grasping function in severe hemiplegic patients using functional electrical stimulation therapy. Neurorehabil Neural Repair 2008; 22: 706–714.10.1177/1545968308317436Search in Google Scholar

[44] Tippett W, Alexander L, Rizkalla M, Sergio L, Black S. True functional ability of chronic stroke patients. J Neuroeng Rehabil 2013; 10: 1–13.10.1186/1743-0003-10-20Search in Google Scholar

[45] Turton A, Pomeroy V. When should upper limb function be trained after stroke? Evidence for and against early intervention. NeuroRehabilitation 2002; 17: 215–224.10.3233/NRE-2002-17307Search in Google Scholar

[46] Uno Y, Kawato M, Suzuki R. Formation and control of optimal trajectory in human multijoint arm movement. Minimum torque-change model. Biol Cybern 1989; 61: 89–101.10.1007/BF00204593Search in Google Scholar

[47] van der Helm FC. A finite element musculoskeletal model of the shoulder mechanism. J Biomech 1994; 27: 551–569.10.1016/0021-9290(94)90065-5Search in Google Scholar

[48] Viau A, Feldman AG, McFadyen BJ, Levin MF. Reaching in reality and virtual reality: a comparison of movement kinematics in healthy subjects and in adults with hemiparesis. J Neuroeng Rehabil 2004; 1: 1–7.10.1186/1743-0003-1-11Search in Google Scholar

[49] Warner M, Chappell P, Stokes M. Measuring scapular kinematics during arm lowering using the acromion marker cluster. Hum Movement Sci 2012; 31: 386–396.10.1016/j.humov.2011.07.004Search in Google Scholar

[50] Wolfe CD. The impact of stroke. Br Med Bull 2000; 56: 275–286.10.1258/0007142001903120Search in Google Scholar

[51] Wolpert DM, Ghahramani Z, Flanagan JR. Perspectives and problems in motor learning. Trends Cogn Sci 2001; 5: 487–494.10.1016/S1364-6613(00)01773-3Search in Google Scholar

[52] Wu C, Trombly CA, Lin K, Tickle-Degnen L. A kinematic study of contextual effects on reaching performance in persons with and without stroke: influences of object availability. Arch Phys Med Rehabil 2000; 81: 95–101.10.1016/S0003-9993(00)90228-4Search in Google Scholar

[53] Wu G, van der Helm F, Veeger H, et al. ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion – part ii: shoulder, elbow, wrist and hand. J Biomech 2005; 38: 981–992.10.1016/j.jbiomech.2004.05.042Search in Google Scholar

[54] Zhang D, Guan TH, Widjaja F, Ang WT. Functional electrical stimulation in rehabilitation engineering: a survey. In: 1st International Convention on Rehabilitation Engineering and Assistive Technology, Singapore, 2007.Search in Google Scholar

[55] Zoia S, Pezzetta E, Blason L, et al. A comparison of the reach-to-grasp movement between children and adults: a kinematic study. Dev Neuropsychol 2006; 30: 719–738.10.1207/s15326942dn3002_4Search in Google Scholar

Received: 2014-2-4
Accepted: 2014-10-6
Published Online: 2014-10-30
Published in Print: 2015-6-1

©2015 by De Gruyter