Functional electrical stimulation (FES) has received significant attention in recent years for its positive impact on motor recovery in post-stroke patients or those with spinal cord injury. By applying electrodes to the skin of the forearm, FES can elicit hand movements that otherwise might be impossible. However, achieving precise and physiological hand movements is still challenging (Westerveld et al., 2012; Malesevic et al., 2012; Imatz-Ojanguren et al., 2016); common applications are therefore limited to simple movements such as hand opening and closing. In order to stimulate complex and functional hand movements, different wrist and finger muscles must be selectively activated in a specific pattern (Imatz-Ojanguren et al., 2016). This is particularly difficult at the forearm due to the anatomy—numerous small muscles arranged in several layers (Westerveld et al., 2012; Bao et al., 2018; Baker et al., 1993). A muscle is activated selectively when it is activated while minimizing the activation of neighboring muscles. In addition to the stimulation parameters (current amplitude, pulse width and shape, and frequency) and chosen electrodes (size and geometry), the selective activation of a muscle depends particularly on the stimulation site (Westerveld et al., 2012; Bao et al., 2018; Baker et al., 1993). The best position to stimulate a single muscle is generally referred to as the “motor point” (MP) (Baker et al., 1993; Imatz-Ojanguren et al., 2016; RaviChandran et al., 2020b).

Current research approaches in the context of selective muscle activation often attempt to enable selective stimulation at different MPs by increasing the number of electrodes. Arrays with many small electrodes are placed on the forearm in order to determine the best stimulation positions and settings using manual or semi-automated algorithms (Koutsou et al., 2016; Salchow et al., 2016a; Crema et al., 2018). The evaluation of the stimulated movement is performed either manually by an expert (Salchow et al., 2016a) or automatically by measuring movements (Malesevic et al., 2012; Bao et al., 2018; Salchow et al., 2016b) or forces (Keller et al., 2007). All previous studies use non-specific arrays (electrodes arranged in regular grids) and are based on a limited number of experiments because these are very time-consuming and are further limited by premature muscle fatigue during FES. The search algorithms are therefore not optimized—they require many repetitions and either take a long time or only cover a small range of the possible settings. Despite extensive research, these approaches have not yet succeeded in developing a system suitable for everyday use that reliably and without great additional effort enables the stimulation of selective muscle activation in the forearm as a basis for complex hand movements. It is thus evident that other or complementary methods are needed. Particularly in the field of biomedical engineering, there is great potential in coupling human-related models with digital twins using biomechanical or bioelectronic simulations (Neumann et al., 2020).

In the context of FES, simulation studies have been successfully used to enhance understanding and support the development of new FES applications (Imatz-Ojanguren et al., 2016). Existing simulation studies have primarily focused on two influencing factors: electrode configuration and fat-layer thickness. Evaluation criteria in simulation-based analyses often cover either the activation of single modeled nerve fibers or a general volume area in which fiber activation occurs (activation volume). For instance, Gomez-Tames et al. (2011) analyzed the influence of fat-layer thickness on activation volume, and five other studies have investigated the interdependency of fat thickness and electrode size (Doheny et al., 2010; Gomez-Tames et al., 2012; Gomez-Tames et al., 2013b; Kuhn and Keller, 2005; Kuhn et al., 2010), inter-electrode distance (Doheny et al., 2010; Gomez-Tames et al., 2012; 2013b), and electrode shape (Gomez-Tames et al., 2012). With respect to the electrode configuration, Kuhn et al. (2009b) analyzed the influence of electrode material properties (hydrogel resistivity) and composed electrodes (as often used with array electrodes) on the required stimulation current to achieve a defined activation (depth of activation volume). Cooper et al. (2011) also addressed electrode material properties and simulated resistivity changes during prolonged use. Two other studies have addressed just the electrode shape with respect to the nerve fiber orientation: Goffredo et al. (2014) analyzed different patterns of an electrode array; Loitz et al. (2015) rotated a rectangular single electrode. In contrast, RaviChandran et al. (2020a) analyzed the influence of electrode shape in terms of different edge length but same area. However, none of the studies analyzed the influence of stimulation settings and person-specific characteristics on the selective activation of individual muscles, taking into account the surrounding muscles. Such studies are critical to support experimental improvement and enable the transfer of these approaches to daily practice.

Existing simulation models only represent the “outer” anatomy (i.e., the different homogeneous tissue layers) and are not suitable for analyses for selective muscle activation due to missing or randomly arranged innervation areas. All models cover a 3D volume conductor with one or various embedded nerve fibers which are coupled to a nerve model. This two-step approach was first introduced by McNeal (1976). The volume conductor is used to simulate the extracellular potential distribution at the nerve fibers. On this basis, the activation of the nerve fibers is then calculated with the nerve model. All existing 3D volume conductors consist of three to five homogeneous tissue layers, mostly arranged in symmetrical cylinders. In contrast, past models differed in the arrangement of the nerve fibers: single nerve fibers at specified locations (Doheny et al., 2010; Cooper et al., 2011; Gomez-Tames et al., 2013b; RaviChandran et al., 2020a), nerve fiber bundles at specified locations (Goffredo et al., 2014; Kuhn et al., 2009b; 2010), many nerve fibers homogeneously distributed in the whole muscle layer (Kuhn and Keller, 2007), and many nerve fibers located in a specified volume representing one single muscle (Gomez-Tames et al., 2013c; a; Loitz et al., 2015). However, none of the existing models provides muscle-specific activation resolution, which is required to use the simulation to develop new applications and products that enable the stimulation of complex hand movements.

This paper introduces the design and evaluation of a new, person-specific simulation model that enables the analysis and evaluation of simulation results at the level of selective muscle activations. In contrast to previous models, our volume conductor is anatomy-based, which means that an individual anatomy given by MRI data is approximated by regular geometries. Furthermore, the nerve fibers in our model are arranged in various muscle-specific regions. This allows conclusions to be drawn about the activation of each individual muscle as the result of a stimulation pulse at a specific electrode position. The developed model is evaluated with experimental measurements and shows good alignment for strength–duration and force-intensity curves. Thus, the new model enables the simulation of single muscle activation, which forms the basis for complex movements and further simulation-based analysis.

The two main objectives of this paper are as follows:

• The introduction of a new modeling method to develop a digital twin of the forearm that accurately and reliably predicts real stimulated muscle activations for different electrode configurations and stimulation settings.

The new individual model basically consists of three parts: an anatomically based 3D volume conductor with four homogeneous tissue layers and two stimulation electrodes placed on the skin surface (Figure 1A), a muscle-specific nerve fiber arrangement in various regions of interest (ROIs) (Figure 1A, B), and a standard nerve model connected to the single nerve fibers (Figure 1C). According to the two-step approach mentioned above, the volume conductor was used to simulate the extracellular potential distribution caused by a single stimulation pulse, and the nerve model was used to simulate the activation state of the single fibers based on this potential distribution. In addition, individual muscle activations are calculated based on the simulated fiber activations, which is enabled by the muscle-specific arrangement of nerve fibers. This model and simulation setup is then evaluated by comparing simulated strength–duration (SD) diagrams with experimental measurements.

Previous models have been limited in simulating the activation of several single muscles, which is necessary for using them as a digital twin and supporting the development of individualized applications and new products. To demonstrate that our simulation model maps real muscle activations correctly and reliably, we generated strength–duration (SD) diagrams for various muscles, both for simulations and the experimental measurements. Typically, in experimental studies, the lowest muscle response is chosen for the creation of the SD diagram. However, in simulations, even at low intensities, individual nerve fibers are activated; such activations may not be measurable in experiments and could lead to an overestimation of activation in the simulation compared to the experiments. To ensure reproducible and comparable thresholds from both experimental measurements and simulations, we opted for the onset of the characteristic sigmoidal progression of the intensity-force curve instead of the usual first muscle twitch. For the example shown in Figure 4A, this means that the activation thresholds used to create the SD diagrams were shifted by force increases from 3–4 mA to 6–7 mA. This applied to all four repetitions of the experiment. Figure 4B shows, for the same muscle, the fitted mean curve that was used for the further comparisons to simulated SD diagrams. Additionally, the maximal variance for the four repetitions is shown. Although the three electrode positions differed (<1 cm) when repeating the marking procedure, the resulting SD plots show only small differences for all muscles with relative standard deviations (RSDs) between 0.04 and 0.09, except for the ECRL, which is slightly higher at 0.14.

Figure 5A shows the simulated SD curves for one exemplary muscle (extensor carpi ulnaris, or ECU) using the linear cable model (LC) as well as the activating function (AF) compared to the experimental SD curve. The simulated SD curves vary depending on the chosen activation threshold (ath), which is defined as the percentage activation of the corresponding region of interest (ROI). To determine which activation thresholds are appropriate, we conducted the following analysis. Figure 5 includes uniformly distributed ath between 0.1 and 0.6 to illustrate the influence of this factor and the range of the resulting SD curves. In general, the graphs generated with the LC model follow the same trajectory as the experimental curve, whereas the graphs generated using the AF especially fail to map the characteristic curve for low pulse widths (PW) up to 200 µs. Therefore, we excluded the AF in further analysis. To analyze the deviations of the simulated and experimental SD curves, the relative differences were plotted for each PW in 1 µs steps. As illustrated in Figure 5B, the SD graphs for the LC model fit for either low PW or those higher. The curves of the ath with comparable overall good fit intersect at PW 200 ± 50 μs. We therefore determined muscle-specific activation thresholds for two PW ranges, 20–200 μs and 201–500 μs, based on the mean of the relative differences between the SD curves at each PW in 1 μs steps. Table 6 summarizes these values for the different muscles and activation thresholds, with the best fitting cases highlighted in bold. For the muscles that could not be selectively/experimentally validated, we chose the activation thresholds based on the average of the validated muscles for the respective pulse width ranges: 0.20 for PW = 20–200 μs and 0.12 for PW = 201–500 μs.

Figure 6 shows the distribution of the relative differences between simulated PW to the experimental mean values (A) compared to the distribution of the inter-experimental differences (single experiments to experimental mean). The comparison shows that for both pulse width (PW) intervals, the remaining differences in the simulated values are in the same range as the inter-experimental differences for all muscles. It can be assumed that for a higher number of experiments, the experimental mean converges against the simulated curve.

Figure 7 shows the comparison of the simulated intensity-force curves with (corrected) and without (default) inclusion of the muscle-specific activation threshold (ath) with the experimentally measured intensity-force curves for different PWs for one muscle (ECRB) as an example. The experiment in this case was performed additionally and was not part of the determination of the activation thresholds. As can be seen, the corrected simulated curve fits well with the experimental curve for all PWs. Both curves show the characteristic onset of the sigmoidal progression of intensity-force curves. The characteristic late bending cannot be seen as the amplitude was not increased until the maximum muscle force was reached. The default simulated curve instead differs clearly from the experimental curve and does not show the characteristic sigmoidal progression.

Despite the established use of FES for many years and extensive research, current applications at the forearm to stimulate hand movements are still limited to simple movements due to the difficulties of achieving selective muscle activations. Most research in this area has focused on experimentally determining motor points (MPs) using multichannel systems with electrode arrays and manual or semi-automatic search algorithms (Salchow et al., 2016a; Malesevic et al., 2012; Keller et al., 2007; Koutsou et al., 2016). However, the results and their impact are limited, mainly due to the small number of experiments (time-consuming, premature muscle fatigue). Simulations have been used successfully in the past to complement experimental research, but existing simulation models are not suitable for analyses of selective muscle activation.

This study used a new modelling approach to develop an individual forearm model that allows the analysis of FES at the level of individual muscle activation through simulations. Like existing models, our simulation model represents the “outer” anatomy by homogeneous tissue layers (bone, muscle, fat, and skin), as in, for example, Kuhn and Keller (2005), Loitz et al. (2015), RaviChandran et al. (2020b), and Gomez-Tames et al. (2013b). In contrast to these models, the shape of our tissue layers is neither fully simplified (as e.g. in Kuhn and Keller, 2005) nor fully anatomic (as e.g. in Gomez-Tames et al., 2013b). The shape of our model is anatomically based, meaning that an individual anatomy given by MRI data has been approximated with regular geometries—mostly ellipses. The advantage of this choice is that the resulting volumes and surfaces are still mathematically defined and, therefore, the model can be easily extended or modified. This is important:

• to allow the modeling of electrodes of any size and geometry (including electrode arrays in future applications) at any location on the skin;

As all changes to the electrode configuration are automatically updated via scripts and macros, the model can be easily integrated into optimization loops. Optimization loops can be used, for example, to find the best individual stimulation positions, which could then be integrated into an individual forearm sleeve.

In addition to the “outer” anatomy, our model also includes aspects of the “inner” anatomy; more precisely, it includes muscle-specific regions of interest (ROIs) that represent the muscle areas where the nerve entry point (NEP) and the innervation zone (IZ) of individual muscles are located. The nerve fibers were arranged in these various muscle-specific ROIs, which differ from existing models, which cover three main types of arrangement: single fiber or one group of fibers placed at a specific depth (e.g., Kuhn et al., 2010; Goffredo et al., 2014; RaviChandran et al., 2020b), uniform distribution of various fibers covering the whole muscle layer (e.g., Kuhn and Keller, 2007), and fibers arranged in one muscle-specific region (e.g., Gomez-Tames et al., 2013a; Loitz et al., 2015). In the latter case, this region was placed in the muscle layer based only on general anatomical information; in our model, the muscle-specific regions were placed based on the transfer of literature data regarding the position of NEPs to the underlying anatomy (MRI data). Due to the muscle-specific nerve fiber arrangement, our model enables the simulation of single muscle activations and, consequently, simulation-based analyses of selective muscle activation.

The experimental evaluation shows an overall good coincidence of the simulated and experimentally measured muscle activations. The small RSD of the different measurements indicate that the chosen reference system provides a robust method for transferring electrode positions between the digital twin and real forearm (Figure 4). The comparison of the simulated strength–duration (SD) curves with the measured curves shows that the simulations with the LC model are generally suitable for reproducing the characteristic course of the SD curves, unlike the simulations with the AF. The best fit is achieved when splitting the considered pulse width (PW) range in two intervals: 20–200 μs and 201–500 μs. For each muscle, different activation thresholds showed the best results. The definition of a global activation threshold would always lead to high differences in some muscles and is only used for the muscles not yet considered experimentally. Considering these limitations for low PWs up to 200 μs, the maximal mean relative difference between simulated and experimental curve is 0.045 and, for PWs higher than 200 μs, 0.013.

The main reason for the need to split the considered PW range in two intervals is that the measured amplitudes still decrease for high PWs. With increasing PW, the amplitude normally converges toward a minimum amplitude required to stimulate a muscle response—the rheobase (Baker et al., 1993). According to Baker et al. (1993), PW higher than 300 μs does not usually affect the stimulation amplitude. Therefore, we are currently working on evaluating a second existing person-specific model to understand whether this effect is person-specific or occurs more often and that perhaps the reporting studies have a bias, such as due to less precise measurement systems.

In addition to this experimental anomaly, there exist different simulation-based reasons that could explain the lack of the simulated curves to represent the measured curves over the full PW range.

• The selected nerve fiber orientation: all nerve fibers are modeled parallel to the ROI, which means parallel to the muscle course. This may not represent the actual fiber orientations with sufficient accuracy, as the terminal nerve branch splits after the NEP, and the individual fibers pull in all directions to reach all muscle fibers. It has been shown that the relative orientation of the nerve fibers to the electrode (if not rotationally symmetrical) influences activation (Loitz et al., 2015).

The first two of these simulation-based reasons may affect different muscles in different ways, and therefore might also affect the “global activation threshold” aspect mentioned above. However, an important question here is, “Is there one global activation threshold for all muscles?”. Measured SD curves for different muscles show differences (Baker et al., 1993), which indicates it is not. This is consistent with anatomical studies showing that fiber type distribution differs between muscles (Hwang et al., 2007). Based on this hypothesis, the definition of mean values as realized for the muscles not measured will lead to inappropriate muscle activations for these muscles. Therefore, we are currently working on individual measurements of these models to determine muscle-specific activation thresholds.

Considering the muscle-specific activation thresholds for the two PW intervals leads to promising results in predicting muscle activation independently of the stimulation intensity, which is determined by PW and amplitude. The simulated strength-intensity curves fit additional experimental data not used in prior determination of the activation thresholds and shows the characteristic force increase at the same amplitude. As muscle activation patterns are composed of single muscle activations, it can be assumed that the model is also suitable for simulating these and enabling further conclusions on stimulated movements or forces when several muscles are activated simultaneously, as is usual.

Overall, the approach presented holds significant promise of being suitable for designing digital twins of the forearm that can enable accurate and reliable prediction of real stimulated muscle activations. The next steps cover the two aspects mentioned above: experimental evaluation of a second individual model, which is already designed, and experimental evaluation of the remaining ROIs for the model presented. Furthermore, we aim to use the activation patterns simulated with the current model to predict the resulting movement or force. In a first step, each measured force could be related to a weighted sum of the activations of all muscles affecting this force location and direction (e.g., all extensors analyzed in these studies contribute to an extension in the wrist). In a second step, the model presented could be coupled with a biomechanical musculo-skeletal model to predict and visualize real movements. The current model is already suitable for integration into product development processes due to the automated model and simulation workflow modification: it can be integrated in optimization loops to find the best stimulation positions and design a patient-specific forearm sleeve for easy stimulation setup, or it can be used for developing optimized search algorithms for common electrode arrays with regular grids. In order to increase the benefits of this digital twin for product development, a future parameterization is planned. Such a parameterized digital twin will enable the prediction of muscle-specific activations by FES on the forearm for different forearm anatomies, which will particularly further improve the development of universal solutions using optimization or AI methods.