Omnidirectional endpoint force control through functional electrical stimulation

Objective. In recent years, Functional Electrical Stimulation has found many applications both within and outside the medical field. However, most available wearable FES devices are not easily adaptable to different users, and most setups rely on task-specific control schemes. Approach. In this article, we present a peripheral stimulation prototype featuring a compressive jacket which allows to easily modify the electrode arrangement to better fit any body frame. Coupled with a suitable control system, this device can induce the output of arbitrary forces at the end-effector, which is the basis to facilitate universal, task-independent impedance control of the human limbs. Here, the device is validated by having it provide stimulation currents that should induce a desired force output. The forces exerted by the user as a result of stimulation are measured through a 6-axis force-torque sensor, and compared to the desired forces. Furthermore, here we present the offline analysis of a regression algorithm, trained on the data acquired during the aforementioned validation, which is able to reliably predict the force output based on the stimulation currents. Main results. Open-loop control of the output force is possible with correlation coefficients between commanded and measured force output direction up to 0.88. A twitch-based calibration procedure shows significant reduction of the RMS error in the online control. The regression algorithm trained offline is able to predict the force output given the injected stimulation with correlations up to 0.94, and average normalized errors of 0.12 RMS. Significance. A reliable force output control through FES is the first basis towards higher-level FES force controls. This could eventually provide full, general-purpose control of the human neuromuscular system, which would allow to induce any desired movement in the peri-personal space in individuals affected by e.g. spinal cord injury.


Introduction
Functional electrical stimulation (FES) for artificial generation and support of movements through application of electrical currents represents a promising tool in the rehabilitation of certain neurological patients. Rehabilitation, at its root, has the purpose of forming new neural connections in lieu of damaged ones, typically between the central and the peripheral nervous system, by re-training the patient to perform movements or tasks. FES offers many advantages with respect to rehabilitation facilitated through externally exerted forces, mainly because the patient's muscles are stimulated and thus actively employed for task completion, thus avoiding secondary complications such as muscle atrophy. In the early phase of rehabilitation, FES can be used as an effective tool in a task-specific, restorative therapy program to foster neurological recovery [1]. In the chronic phase after a neurological disease or trauma, FES may still be used as a neuroprosthesis for compensation of completely lost or very weak motor functions. Particularly in individuals with spinal cord injury (SCI) and the associated functional impairments, FES has been successfully employed for assistance in activities of daily living (ADL), both using trans-cutaneous [2] and intramuscular electrodes [3]. Non-invasive FES applied through surface electrodes is also used in applications outside of the medical field, for example VR and AR [4].
Most control schemes focus on restoring functional, task relevant movements, such as reaching and grasping [2,3], and focus on the identification of the dynamics relevant to these [5,6]. Such parameters, however, are not guaranteed to generalize well over different postures. On the one hand, black box approaches typically have to sample the effects of muscle contractions in various postures [7,8]. Musculoskeletal models, on the other hand, can inherently account for at least some effects of posture changes [9,10]. In robotics, impedance-based controls can be used to impose a certain dynamic behaviour between a robot and its environment [11]. Impedance controls are robust in terms of kinematic singularities, and can be well integrated in wider motion-planning algorithms, but the main purpose of impedance control is to facilitate the interaction of a robot with an unpredictable environment with non-linear dynamics [11]. However, the characteristics of impedance control make it robust also with respect to actuators exhibiting these characteristics, as can be human muscles. A suitable impedance could facilitate conversion from a positional error into a desired force output, which is more directly correlated to muscular activity. FESbased force controls have been proposed, among others, in [7,12,13]. These works all present black box models of the endpoint force output as a consequence of FES. In the context of FES, a suitable impedance control could be used to assist movements towards any desired point in the user's peri-personal space, leading to a more general-purpose paradigm, which could be beneficial in rehabilitation, but especially in the case of FES used as a neuroprosthesis. Razavian et al present a demonstration of such a concept [7]. The compliant nature of the human body could allow for the safe inclusion of the positional error's integral over time in the impedance, which would lead to increased robustness with respect to modelling errors. Integrative terms are often excluded from impedance controls in conventional robotics, as they would cause an increase of force output over time, should the robot encounter an obstacle preventing it from reaching the desired pose. In such a scenario, the robot could cause damage to itself or its surroundings if its force output is not limited. This paper should serve as a system description of the FES device, which consists of a wearable surface stimulation device designed to provide proportional force control through FES on the upper limb of a user on up to 10 channels with a resolution of 16 bits, and assess its capability to induce a desired force output in real time. Because of its practicality and versatility, surface FES is widely used in commercially available products, such as the Teslasuit ® platform (VR Electronics Ltd., London, UK). This device has been successfully employed in user studies with able-bodied participants [14], and represents a good commercial benchmark, as it integrates various sensor modalities and gel-less surface stimulation electrodes. However, this system cannot adjust well to different body frames, and the electrode arrangement cannot be modified. The device proposed here, on the other hand, features a Velcro-lined compressive jacket which allows for easy modifications in electrode arrangement to fit any user frame. The jacket also increases repeatability of electrode placement once the ideal arrangement for an individual user has been established.
While this device proposes to be a general-purpose platform designed to test various control algorithms, here the system is driven by a musculoskeletal model presented and validated in [9]. Therein, the musculoskeletal model was tested against a third-party model introduced in [10], which was taken as baseline. The model associates a line of action to each stimulated muscle group, as exemplified in [15]. A core principle and fundamental goal in the design philosophy of this system is the adaptability to different users. In order to achieve this, the musculoskeletal model can be easily modified to better fit each individual without the need for much anatomical expertise. To this end, a calibration procedure able to adjust the model geometry is also introduced here. Canonically, a model based on line of action relies on a line running through the average centroid of the physiological cross section along the whole length of the stimulated muscle groups, as introduced in [16] and, more recently, in [17]. Many studies have demonstrated how inhomogeneities in muscle activation can lead to great effects on joint momenta. In the case of the musculoskeletal model used here, the reconstruction of the line of action is based on the observation of functional effects of muscle contraction on the musculoskeletal system in a given position, as introduced among others in [18]. Here, the model uses a fast and computationally efficient Nearest Neighbour recruitment strategy to calculate the stimulation currents necessary to cause a given endpoint force output.
The system is validated in an online experiment where the FES-induced force output of 3 able-bodied volunteers is compared to a desired force output. Part of these results were published in [19]. In addition, we evaluate the performance of an offline-trained predictor which is able to precisely predict the force output both in task space and in joint space based on the stimulation currents. Such a predictor could be trained based on data from a force-torque sensor prior to the normal operation of the device in a setup similar to the one used during the online system validation.  figure 1(a). The muscle groups stimulated by one electrode pair are considered as a string routed along the line of action's curve, with the muscle force being exerted homogeneously in a tangential direction, which is to say that the force along the line of action has constant module Î  f . Given this, knowing the line of action's routing in 3-D space, we can compute the average force vector¯Î  f 3 through Following the same rationale, we can calculate an average torque vectort Î  3 exerted by the muscle group about the joint's position Î  j Reducing the points on the curve to a finite set of becomes a sum of the form The average forcef and average torquet derived from the line of action characterize a lumped model of the muscle, which behaves like a prismatic joint able to exert a force f determined by the stimulation in the direction¯f f 1 by contraction, with an average effective moment armr, which results in a torquet about the joint at position j with revolute axesē i , as depicted in figure 1(a). Depending on which skeletal segments the muscle group originates from and inserts into, this prismatic actuator approximation of the stimulated muscle group can exert a torque around more than one of the skeleton's rotational axes. In particular, if the i-th degree of freedom is a revolute joint at position j with its axis pointing in a known direction Î  e i 3 , if a muscle group is able to exert a torque about it, we can reduce the expected torque computed as shown in equations (2) and (3) to a scalar torque magnitudet Î  i by computinḡ¯¯( If the muscle group is not able to exert a torque about the i-th degree of freedom, on the other hand, the scalar projection of the expected torque onto joint space is 0. As depicted in figure 1(a), every joint in the musculoskeletal model is modelled as having 3 revolute axes. This is done due to the difficulty of reliably measuring the direction of the anatomically correct revolute axes. This makes every joint in the model as depicted in figure 1(b) defined by a single position j and 3 revolute axes.

Model calibration
The line of action's routing can be initially set based on cursory anatomical expertise and the known electrode placement, but can be further adjusted based on the measured response to stimulation pulses. Here, a simplified calibration procedure is implemented. This is based on the assumption that the average force vectorf should lie in the plane in which the limb moves when a stimulation pulse is applied. The calibration finds the points a i for each line of action which minimize the difference between the expected and the measured torque output in joint space resulting from the stimulation pulse. Although this condition could be satisfied even if the line of action were not to lie entirely on this same plane, this is assumed for simplicity's sake. A stimulation pulse causes an average torquet over the pulse time, which in turn causes the body segments distal to the affected joint at position j to accelerate at an angular acceleration rate¯t a » I distal , where Î I distal 3 3 is the cumulative moment of inertia of the distal body segments, which is assumed to be constant throughout the stimulation, andā Î  3 is the angular acceleration vector due to the twitch. In this study, the moment of inertia was calculated based on the user's mass and the anthropometric tables from [20]. Over the course of this study, the participants were directed to hold their arm in the starting position voluntarily, and therefore the influence of gravity was assumed to be compensated by volitional muscle contraction. Assuming, furthermore, a constant acceleration rate during the stimulation pulse, and an initial rest state of the joint, the distal body segments reach a maximum angular velocity¯( , where t 0 and t end are the times at which the stimulation starts and ends, respectively. Therefore, we have an approximate proportionality between the maximal observed angular velocity and the expected torque of the stimulated muscle group, namelȳ The plane of movement for the limb is defined as passing through the joint's position j and being normal to the angular velocity vector w max of the twitch itself. As shown in figure 2, based on the observed twitch w max , the position of the origin point a 0 and the insertion point a I are set over the q Î  0 and q Î  I coordinates, respectively, in order to minimize the distance of the two points from the twitch plane passing through j and normal to w max . Here, q 0 and q I are the azimuth of the origin and insertion point, respectively, expressed in a cylindrical coordinate frame the height axis of which corresponds to the longitudinal axis of the proximal body segment in the case of the origin point and the distal body segment for the insertion point, as shown in figure 2. The origin and insertion point positions are additionally defined by a radius Î  r 0 and Î  r I , respectively, and a height Î  x 0 and Î  x I . This optimization is usually done based on a series of K observed twitch . A cost function is minimized in order to find the cylindrical coordinate q* I which minimizes the distance between the insertion point a I and the plane in which the body segment moves during the twitch motion. The cost function is represented by the sum of all the observed distances of the line of action routing points from this plane, and it is minimized as follows The cost function can be minimized through a gradient descent. All routing points with freely settable coordinates, such as the origin point a 0 , are adjusted analogously. In addition to routing point placement, the magnitude of the observed twitch angular velocity vector w max leads to an estimation of the proportionality constant g i between the stimulation current s i and the force f . In addition to that, the stimulation is offset by a constant q i so that minimum current is at the edge at which force is exerted as a result of stimulation. This offset is set manually during an initial phase of comfort level setting occurring at the beginning of the experimental session, shown as Comfort level setting in figure 6. By setting these two coefficients, the stimulation currents delivered to the user can be assumed to be within a nearly linear region of the force to stimulation curve. The coefficients were determined by identifying a minimum and maximum threshold for each channel current. The minimum threshold is slightly below the lowest amount of current that the user is able to discern for a given stimulation channel. The maximum is around the highest non-painful stimulation current for which the force output stops rising. These thresholds were maintained throughout the experiment. The expected torque's magnitudet i is computed as in equation (3) by plugging While the present study was performed with the arm locked in a single posture, the calibration procedure had to be performed with the arm free to move, as the twitch resulting from a sharp stimulation signal had to be observed. Since the calibration operates based on several stored twitch vectors, it could place the insertion and origin point in a position that should minimize the distance of the line of action from the twitch plane for all the postures in which the calibration is executed.

Projection matrices
For the purposes of movement control, it is necessary to calculate the projection matrices from the relevant coordinate systems shown in figure 3, which are the muscular Jacobian ≔ Î , computed by differentiating the pose x arm of the arm's end point, which comprises both 3 positional and 3 orientation coordinates, by the same rotational movements q joints . In both cases, due to energy conservation, the Jacobians can project velocities or differentially small shifts in position from the joint space to the Cartesian or muscular space, as well as generalized Cartesian wrenches Î  w hand 6 and muscle forces f mus onto the joint space.
The averaged torque from equation (4) can be used to approximate the muscular Jacobian J mus defined from the muscle space to the joint space, as shown in figure 3, as per where Î  f mus M is the vector of the forces acting on the M actuators, formed by the single scalar force components f i , andt j i , is the scalar projection of expected torque caused by the i-th muscle around the j-th revolute joint, as shown in equation (4).
The arm's Jacobian J arm is similarly calculated and can also be used to project the Cartesian wrench w hand onto the joint space, as per Figure 3. General overview of the coordinate spaces relevant for the control of the human musculoskeletal system. Muscles can be modeled as compliant, prismatic joints. In all spaces, a corresponding set of generalized forces can be defined. These forces are exerted on the musculoskeletal system by both its own actuators and the environment. The momentary positional state of the system and its time derivatives can be linked to a generalized force through an impedance Z , whose inverse is the admittance Y . Positions and forces in either space can be projected onto another space by using forward or inverse kinematics, or employing differential calculus, by the appropriate Jacobians J . These can't always be inverted, and when this is the case, a fitting matrix pseudo-inverse should be put in place. 3 3 is the so-called skew symmetric matrix of the vector a, which satisfies [ ] =á b a b,×being the vector cross product,¯Î  e j 3 represents the direction of the j-th joint's axis in the Cartesian coordinate frame, Î  f hand 3 and t Î  hand 3 are the desired force and the torque to be applied to the endpoint, which are concatenated in a wrench w hand as shown above. These relations are summarized in figure 3. The control system uses a musculoskeletal model shown in figure 1(b), which works based on the principles discussed above. This model is used to compute the stimulation necessary to achieve a certain force output at the endpoint, and this is done by projecting the desired wrench w hand onto the joint space using the transposed arm Jacobian J arm . In order to calculate the necessary muscle forces f mus , and consequently the stimulation currents Î  s mus M , the desired joint torques are approximately projected onto muscle force space using a Nearest Neighbour muscle recruitment strategy, by which only the single muscle group that would cause the torque closest in direction to the desired torque on a given joint are stimulated. Solving for the necessary muscle forces is often nontrivial, as muscle forces are not a conventional vector space, because of the fact that muscles can only actively contract, not expand. The pseudo-code in algorithm 1 explains how this strategy works. Algorithm 1. The muscle force and stimulation solver iterates over all joints shown in figure 1(b). For each joint, only the muscle group which would elicit the torque closest in direction to the desired one is stimulated.

Require:
Î  w hand 6 , which is the desired wrench at the endpoint Require: The algorithm iterates over all joints shown in figure 1(b), which are all modelled as having 3 revolute axes, as shown in figure 1(a). Because of this, both the required and the exerted torques for each joint are three-dimensional. This simplification was put in place because of the difficulty of reliably estimating the direction of certain anatomical musculoskeletal axes. The Nearest Neighbour recruitment strategy can be applied without this simplification, but it does require that every modeled muscle group only exerts a torque about one joint, and that each joint can be defined as having a central position j and up to three revolute axes. The Nearest Neighbour recruitment strategy is efficient in terms of computation and time, in comparison to iterative optimisation algorithms used, for example, in [12]. It also reduces the amount of current delivered to the user, thus improving comfort and minimizing metabolic costs.

Hardware and experimental setup
The FES device's wearable stimulation setup can inject stimulation currents through surface electrodes. The system can provide amplitude-modulated, rectangular current stimulation pulses with 16-bits resolution on up to 10 channels, with a pulse-width of 200 μs, frequency ranging from 0.5 Hz to 100 Hz, and a maximum current amplitude of 70 mA. In addition to the wearable stimulation device, the setup includes a posture tracking sensor system, which in the case of the experiment presented here was the BodyRig [21], an IMU-based posture tracker, visible in figure 4.
In order to generate the stimulation currents, the setup includes 3 FES devices (2 TNS SM2 AKS and 1 TNS SM 2MF, Pierenkemper GmbH, Am Geiersberg 6, 35 630 Ehringshausen, Germany). An intermediate wirelessly controlled modulation box built around the wireless Bluetooth module ESP32 Wroom 32, Espressif systems modulates the generated currents in amplitude from 0 A to the maximum amplitude set on the FES device. The channels are electrically insulated from each other, and are controlled by using analog optocouplers driven by operational amplifier-based driver circuits. The levels of stimulation for each channel are calculated by a remote host running the control model. Exact schematics of the control box are available on request.
The system was validated through an experiment. The setup consisted of the wearable stimulation device, the BodyRig posture tracker [21], as well as a robotic arm [22], which was used here as a measurement device. The robotic arm streamed the measurements of its 6-DoF force-torque sensor, which is visible in figure 5 and served as coupling between the robot's end effector and the cuff holding the user's forearm, to its own host at a rate of 1kHz. The stimulation device's control loop was operating at a rate of 200 Hz on a separate host, as shown in figure 5. The signals to the single stimulation channels, however, were sub-sampled at 20 Hz. This rate for the controller was deemed sufficient for this scenario based on pretests, which showed that more than 99% of the force output signal power lies between 0 Hz and 2 Hz on the frequency spectrum. This is in line with observations published in [23]. The outcome metrics of this study are the measured force outputs compared with the commanded forces in terms of Pearson correlation coefficients. Furthermore, a regression model was trained and cross validated on the acquired data in order to test the feasibility of a machine learning-based controller of the force output via FES.
The surface electrodes were applied in order to stimulate the biceps brachii, the triceps brachii, the deltoid superior, anterior, posterior, the clavicular and the sternocostal head of the pectoralis major, the trapezoid scapular, and the latissimus dorsi.

Experimental protocol
The system is validated here through an experiment employing 3 able-bodied male volunteers (34.3 ± 12.7 years old, 1.76 ± 0.09 m, 77.3 ± 6.67 kg). The experiment was approved by the ethical commission of the university/institution to which the authors are affiliated. This research was conducted in accordance with the principles of the Declaration of Helsinki and in accordance with local statutory requirements. All participants were thoroughly informed about the experimental protocol, and all gave written informed  consent to participate in the study. This research involves no identifiable human subjects, and does not rely on clinical trials. Participants sat in a predetermined position as shown in figure 5 with their right arm coupled to a force-torque sensor attached to a robotic manipulator [22]. Using a robot-mounted force-torque sensor allowed for easy realignment if the robot had to be re-positioned to better fit the size of the user, as the robot's arm made it possible to instantly know the absolute orientation of the forcetorque sensor in space, and therefore to reconstruct the absolute direction of the measured forces and torques in the environment. The session breakdown is shown in figure 6. The participants were first asked to voluntarily exert forces along 6 directions, namely backward, forward, upward, downward, left and right, for 10 repetitions by following visual feedback (C1).
Thereafter, visual feedback was taken away. The device, without any calibration, was then fed desired force output vectors selected randomly among the 6 directions from condition C1 at 2 different magnitudes (namely 50% and 100% of the achievable force output, as determined by the threshold detection during the comfort level setting), over 5 repetitions, for a total of 60 commanded force outputs. All conditions employing FES used such a sequence of commanded forces. The FES device provided stimulation to induce a force output corresponding to the commanded force outputs (C2).
Following this, the calibration procedure described above was performed for all stimulation channels, and the experiment with no visual feedback was then repeated with the calibrated device (C3). Finally, one further condition was tested, where the parameters set through the device's calibration are maintained, except the proportionality constant g i between the stimulation s i and the muscle group's contraction force f i . This proportionality constant is reset to its initial value. This condition has the goal of verifying that, if the calibration procedure improves the performance by correcting the musculoskeletal model or by virtue of increasing, even saturating the current flow through the stimulation channels (C4).

Signal conditioning and data analysis
During the experiment, each sequence consisted of a total of 12 target output forces (6 directions×2 magnitudes), corresponding to a set of stimulation currents delivered by the FES device in all experimental conditions except for the first one, where the participants were directed to exert a given force output through visual feedback. In either case, the contraction lasted 2 seconds, and was always followed by 2 seconds of inactivity before the next contraction. The force output measurement was given by the difference of the wrench measured by the 6-DoF force-torque sensor, minus the wrench measured immediately before the onset of the stimulation currents. The difference was averaged over the whole time of contraction. This was done in order to subtract out of the wrench measurement the effect of the arm's own weight and of the user's volitional actions immediately prior to the onset of stimulation. Because the robotic arm used for measurements and the FES device were using two different coordinate systems, the sets of measured and commanded forces were aligned using the Kabsch algorithm [24] through a rotation around the vertical axis. The reason why the rotation was limited to the vertical axis is that the two coordinate systems can be presumed to be aligned along the vertical direction, as the FES device's coordinate system is based on IMUs, which are able to measure the direction of the gravitational acceleration vector. The measured forces thus calculated were compared to the commanded forces in terms of the Pearson's correlation coefficients of the commanded force direction with the measured force direction.
Besides the analysis of the online system performance, a regression algorithm was trained with the goal of predicting the output wrench Î  w hand 6 based on the stimulation currents. This algorithm was based on ridge regression applied to a so-called Random Fourier Features kernel (RFF) [25]. The predictor was cross-validated over a 10-fold partition of the available data, which consisted of 60 data points for each subject and each of the 3 experimental conditions where the FES device was delivering stimulation currents. The regression's predicted wrench is evaluated in terms of squared Pearson correlation coefficient (R squared) of the prediction w.r.t. the ground truth represented by the measured forces and torques, as well as in terms of root mean square (RMS) error. An analogous predictor was trained to predict the estimated joint torques t Î  joints J as a function of the stimulation currents Î  s M . The joint torques were not directly measured, and had to be estimated by projecting all wrenches acting on the user's body onto the joint space, as shown in equation (8a)-(8d). The estimated joint torques are used as target labels to train the regression algorithm to predict the joint torques based on stimulation currents.

Results
The results are subdivided between the offline analysis of the force output predictor and the evaluation of the force output control in real time, hereafter referred to as online control.
For the offline analysis of the force output predictor, table 1 shows the root mean square error (RMSE) normalized by the maximum range of force or torque output for each subject and experimental condition, and the squared Pearson correlation coefficients (R 2 ) of the regression algorithm trained to predict the force output in Cartesian space and the joint torques based on the stimulation currents, respectively. The joints are as shown in figure 1(b). The predictor'performance achieved a grand average of 0.860 R 2 and 0.124 RMSE for the prediction of the joint torques, and 0.855 R 2 and 0.124 RMSE for the prediction of Cartesian output wrench.
Concerning the results of the online control, figure 7 shows the correlation matrices of the commanded force direction , , 3 after the two have been aligned using the Kabsch algorithm. Each matrix element represents the Pearson's correlation coefficient between one of the components of f m and one of the components of f c . Each row of matrices represents one experimental condition (C1 to C4, as shown in figure 6), and each column of matrices represents one experimental subject (S1 to S3). Within the correlation matrices, the asterisks represent the significance of the correlation coefficient in terms of its p-value, as ascertained through the right-tailed test. The diagonal elements in the correlation matrices were compared across condition 2 and 3 through a paired Wilcoxon signed-rank test, showing a statistically significant improvement in the correlation coefficient with p < 0.05. In condition 2, the diagonal elements of the matrices are 0.615 ± 0.204, and for condition 3 the diagonal coefficients are 0.709 ± 0.128 (mean ± 1SD)). Cohen's d is 0.552. Figure 8 shows the RMSE detected during the online experiment. The error is displayed in its components separately, as well as in its euclidean norm. Next to the RMSE is the achievable range also represented in its three components and its euclidean norm. Seen as the RMSE samples are not normally distributed, and that the individual measurements of the RMSE are assumed to be independent of each other, a non-parametric Mann-Whitney U-test was performed to detect potential effects of the experimental conditions on the euclidean norm of the RMSE across all test participants. Significant effects were detected between conditions C2 and C3 (p < 0.01). Among the directions of force, it seems that the most significant overall effect is on the x-axis (p < 0.001) which, as shown in figure 5, corresponds to the backward/forward direction. However, asignificant effect (p < 0.05) is also detectable on the y-axis between conditions C2 and C3, with a noticeable increase in this component of the RMSE. No significant effect is detectable on the vertical z-axis. Significant effects are also present between conditions C2 and C4 (p < 0.05). On a per-participant basis, the effects vary strongly. Participant S3 shows the strongest effect between conditions C2 and C3 (p < 0.001).
In general, the calibration procedure in both C3 and C4 leads to a consistent decrease in variance between the RMSE components compared to C2. The RMSE for condition C1 can serve as qualitative comparison.

Discussion
The correlation between the commanded and measured force for online control, shown in figure 7, shows a high variation depending on the user. The fourth experimental condition, C4, is associated with significantly worse correlations compared to either C2 or C3. While the purpose of condition 4 is to verify whether any benefit deriving from the FES device's calibration are due to the geometrical adjustments of the musculoskeletal model or due to the adjustment to stimulation intensity, the expected outcome in the latter case would have been a better performance than using the non-calibrated FES device, but worse than the fully calibrated condition. The fact that the correlation is worse than both C2 and C3 would seem to indicate that the most likely cause of this decrease in performance is fatigue and other progressive effects, which are well documented in FES applications [26,27]. In this regard, the study design could have benefited from randomizing the order of conditions C2-C4, or at the very least from allowing more resting time between trials. The calibration procedure leads to a more consistent performance across different users (for this, compare for example the C2 and C3 row of matrices in figure 7) and to a significantly higher correlation of the commanded and measured force output, as well as to a lower overall RMSE. S1 does not show significant improvement between experimental condition 2 and 3. This is likely due to the fact that the non-calibrated state of the system is based on S1ʼs frame. The performance assumed as baseline is the one attained when working in the visual feedback condition C1, where the participants were voluntarily exerting force in the indicated direction. As expected, this condition shows the clearest correlation between commanded and measured force direction.
While the RMSE is relatively high in most cases in relation to the maximum achievable range, this is also true for the control condition C1. Participant 2 shows the highest RMSE relative to the maximum force output across the experimental conditions. This is likely due to the difficulties in stimulating the anterior deltoid group of this specific participant. Calibration seems to sensibly reduce the RMSE along the x direction. However, since the muscle groups available for stimulation cannot adequately cover the needed torque output space, the cumulative RMSE remains high, even though the variance between error components is strongly reduced. Overall, calibration seems to lead to an overall reduction in RMSE in both condition 3 and 4 compared to condition 2. The results of this validation show that an open-loop force control enacted through transcutaneous FES can enable qualitatively good directional control of the endpoint force output. This is significant, as the ability to induce a desired force output is paramount for the implementation of any control, even in position.
The performance of the regression algorithm that was trained and validated to predict the output force, on the other hand, shows better precision and accuracy when compared to the real-time correlation coefficients, as shown in table 1, at least when the arm is static, to the point that the correlation coefficients seem to rival those attained in C1. The prediction of the joint torque output should generalize better over different arm postures, which is an aspect worth investigating in future work. Looking at possible implementations using such a regression machine in order to control a FES setup, the results reported in the works by Schearer et al [12] and Friedrich et al [28] can be helpful, mainly because of their finding that muscle stimulation can be modelled as combining linearly in the force output space. As demonstrated in [27], this fact can be employed to design ridge regression-based FES controllers which can facilitate movement with a time-effective calibration procedure, once the Figure 8. Root mean square error of force output next to the achievable range for each condition (C1-C4) and participant (S1-S3). The lowest row of graphs marked Overall shows the more statistically significant effects over all participants as determined by a nonparametric Mann-Whitney U-test. In this row, the maxima are derived by averaging over the per-subject maxima for each condition and RMSE component. *** p < 0.001, ** p < 0.01, * p < 0.05. electrodes are in place. The work by Schearer et al., in particular, presents an experiment which employs a setup in many ways identical to that adopted in the present study. Crucially, however, both Schearer and Friedrich focus on subjects with implanted electrodes, which are in many ways less challenging from a control point of view, as they allow for a far higher selectivity in muscle stimulation, causing only the firing of the nerve cells in the nerve bundle they are implanted onto, as noted also in [13].
Based on the findings of the offline analysis, it should be possible to use a setup such as the one used in this experiment in order to calibrate the system before normal operation. Such a procedure could entail either the training of a regression algorithm similar to the one used in the offline analysis, or simply using the measured joint torques instead of the maximal angular twitch velocity w max in the calibration used in this study. Either way, the results of the offline analysis clearly indicate that it would be preferable to include torque measurements in the calibration of the system, as this would bypass any error deriving from erroneous modelling of the body's dynamic parameters, which did play a major role in the online calibration used in this experiment.
This study presents some limitations. First and foremost, while many works rely on data from a limited amount of participants [7,8,12,29], future work should focus on including a wider base of participants, including spinal cord injury patients. Furthermore, the presented musculoskeletal model, while being computationally efficient, relies on several simplifying assumptions, and currently does not model biarticular muscle groups. The Nearest Neighbour recruitment strategy can also lead to sudden switches in stimulated muscles with changes in posture or desired wrench. Additionally, as the calibration procedure used in this study does not call for direct measurement of the force output at the endpoint, the calibrated musculoskeletal model only allows for qualitative control of the force output.
Overall, the evidence seems to point to some benefits offered by the proposed calibration procedure. However, this study could not determine definitively whether these improvements in performance stem from the calibration's adjustments in the geometry of the musculoskeletal model or from the adjustments in stimulation intensity, as fatigue is likely to have rendered some effects between condition C4 and others harder to detect. Enlarging the pool of subjects and adjusting for time-dependent effects will hopefully clarify this aspect.

Conclusion
In this study, we described a portable, wearable FES device which, coupled with an appropriate posture tracking system and control architecture, can facilitate qualitatively precise endpoint force output control in three dimensions. We also demonstrated a computationally efficient calibration procedure which can adapt the geometry of the musculoskeletal model and the magnitude of the stimulation currents delivered to the user through simple observation of the twitch, with no need for additional sensors beside the posture tracker.
Future work should focus on investigating the inclusion of force measurement systems for the purpose of calibration as done in [12], and potentially also in order to close the control loop in terms of force during normal operation as done in [7]. In particular, the perspective of using an Exosuit in order to both provide further assistance to the user, as introduced in [30], and to measure the joint torques directly would be worth investigating. The addition of an impedance on top of the force control presented here would allow for general purpose movement control, and should be the focus of future research. In such a case, even if we were to adopt a regression model for the computation of the muscular Jacobian, a model such as the one presented here would still be required to project the desired endpoint force output and the interaction forces measured by the exosuit onto joint torque space.
The precise performance observed during the offline analysis clearly indicates that inclusion of force and torque data in the calibration of the system would improve motion control through FES. The described RFF-based ridge regression algorithm could be able to generalize over non-linearities in the force-stimulation curve. However, including posture in the input vector to the regression could be a challenge, and whether RFF regression will be able to generalize over different body postures remains an open question. As the experiment was conducted in a single posture, the system's adaptability to different limb poses is still theoretical, and not definitively proven by this study.
Future work should also entail user studies involving a larger population of both able-bodied participants and patients suffering from neurological ailments that could benefit from FES-facilitated rehabilitation.

Funding
This work was partially supported by the German Research Society individual research grant HIT-Reha -Human Impedance control for Tailored Rehabilitation, DFG Sachbeihilfe CA1389/3-1, project #505327336.