## Abstract

This paper presents a novel multiscale finite element-based framework for modelling electromyographic (EMG) signals. The framework combines (i) a biophysical description of the excitation–contraction coupling at the half-sarcomere level, (ii) a model of the action potential (AP) propagation along muscle fibres, (iii) a continuum-mechanical formulation of force generation and deformation of the muscle, and (iv) a model for predicting the intramuscular and surface EMG. Owing to the biophysical description of the half-sarcomere, the model inherently accounts for physiological properties of skeletal muscle. To demonstrate this, the influence of membrane fatigue on the EMG signal during sustained contractions is investigated. During a stimulation period of 500 ms at 100 Hz, the predicted EMG amplitude decreases by 40% and the AP propagation velocity decreases by 15%. Further, the model can take into account contraction-induced deformations of the muscle. This is demonstrated by simulating fixed-length contractions of an idealized geometry and a model of the human tibialis anterior muscle (TA). The model of the TA furthermore demonstrates that the proposed finite element model is capable of simulating realistic geometries, complex fibre architectures, and can include different types of heterogeneities. In addition, the TA model accounts for a distributed innervation zone, different fibre types and appeals to motor unit discharge times that are based on a biophysical description of the *α* motor neurons.

## 1. Introduction

One of the main advantages of electromyographic (EMG) signals is that they can be measured *in vivo*, providing information on the physiology of the motor system. Thus, it is not surprising that EMG signals are widely used as a source of data in scientific research such as movement analysis, rehabilitation medicine, neurology or development of biofeedback techniques. The EMG signal reflects the electrical activity of skeletal muscles. Muscles consist of electro-active biological cells called muscle fibres. Upon activation through their corresponding motor neuron, action potentials (APs, short-term depolarizations of the membrane potential) are generated in the muscle fibres at the neuromuscular junction. The APs then propagate along the fibres and induce the contraction in the underlying sarcomeres. The EMG is the superposition of the APs of all fibres of all motor units (MUs, the muscle fibres innervated by a single motor neuron and the motor neuron itself) in the three-dimensional muscle tissue.

One of the major drawbacks of EMG signals is that they are hard to interpret and analyse [1]. A problem of surface EMG (sEMG) signals is that they are recorded on the skin surface separated from the signal source through subcutaneous and skin tissue. Besides the difficulty of identifying individual MU discharges from the sEMG signal [2], the distance between the signal source and the detection electrode might change due to the contraction of the muscle. Furthermore, the electrode might also pick up signals from neighbouring muscles, for example, due to the movement of the muscles relative to the detection electrode, and it is often difficult to separate the signals in the EMG [3,4].

In contrast to sEMG, intramuscular EMG measurements provide information only about a very small region, picking up signals from very few MUs [5]. Moreover, intramuscular detections are very sensitive to movements of the muscle, loosing easily their recording position. While these problems also exist in isometric contractions, they are especially pronounced in non-isometric contractions, which explains the fact that EMG signals are often recorded under isometric conditions [2]. Furthermore, based on the limited information provided by the EMG, it is often impossible to explain changes observed in recorded signals, for example during sustained contractions. To accurately and reliably interpret the EMG, further information on the processes taking place in the muscle (e.g. tissue movements, changes in the AP shape and amplitude) is required. Therefore, computational models predicting the EMG have great potential to improve comprehension of recorded signals.

A large number of EMG models, both analytical [6,7] and numerical [8–10], have been developed in the past decades, see [11] for a recent review. While (semi-)analytical models have the advantage of short computation times, they are restricted to simplified geometries (cylindrical and plane layer volume conductors) with idealized, layered tissue conductivities [12,13].

Real muscles, however, are often not well described by such simplified geometries, and the geometry has a considerable effect on the surface-detected AP shape [14]. To simulate realistic muscle geometries with complex fibre architectures and different heterogeneities such as locally varying material properties, numerical models based on the finite element method are required. Previous finite element models described the conductivity tensor using an analytical fibre path function [15], which, however, still does not allow for arbitrarily complex geometries. Furthermore, existing numerical models of the EMG [9,14,15] are single-field models that cannot account for the electro-mechanical feedback, i.e. the contraction-induced deformation of the muscle that is triggered by the propagating APs.

A further shortcoming of existing EMG models is that neither the numerical nor the analytical models are capable of simulating changes in the AP. This is due to the fact that these models prescribe phenomenologically the AP shape based on an algebraic description (impulse response or Rosenfalck approximation) and a constant propagation velocity instead of simulating the underlying biophysical processes. Experimental data, however, suggest that the EMG amplitude and the AP propagation velocity change during sustained contractions [16,17], which increases the variability in the EMG. Neglecting these changes of the AP might lead to inaccurate model predictions.

The aim of this work is to overcome the above-mentioned limitations by establishing a multiscale model that is capable of simulating different physiological and (bio-)physical factors that influence the EMG. Owing to the fact that the EMG is the result of cascades of complex biophysical processes at length scales ranging from the subcellular level to the whole organ level, a multiscale modelling framework is required. In this sense, the multiscale chemo–electro–mechanical skeletal muscle model proposed in [18] serves as the basis for the proposed EMG modelling work. In this model, the generation and propagation of APs along single muscle fibres is described by a biophysical Hodgkin–Huxley-type membrane model in conjunction with a transient diffusion equation. The advantage of the multiscale formulation over phenomenological descriptions is that it intrinsically accounts for changes in the AP shape and propagation velocity, which are observed in the EMG, for example during sustained contractions. Furthermore, the model considers the electro-mechanical feedback by simulating the excitation–contraction coupling in the sarcomeres [19] and linking the sarcomere-based contraction-induced stresses to a continuum-mechanical formulation of skeletal muscle behaviour. That way, the chemical and electrophysiological processes responsible for the AP generation are coupled to force generation and muscle deformation. This allows the EMG to be simulated during non-isometric contractions and provides a framework to investigate the influence of (sub-)cellular properties on the EMG signal. Furthermore, the proposed multiscale model also provides a framework to investigate the effect of positional changes of the source, innervation zone and tendon relative to the detection electrode, which might considerably influence the sEMG. Thus, the multiscale model eventually yields a better understanding of recorded EMG signals during isometric and non-isometric contractions.

## 2. Material and methods

The propagation of electrical signals through biological tissue can be described using the bidomain model [20].

### 2.1. Bidomain model

The bidomain model, see figure 1, is a continuum model describing the electrophysiology of excitable biological tissues. It relates currents crossing the cell membrane to intracellular and extracellular potentials [20]. The bidomain model distinguishes between electro-active materials existing in the muscle domain, *Ω*^{M} and electrically inactive tissues existing in the body domain, *Ω*^{B}. Neglecting the conductivity of the surrounding air, there is no current crossing the outer surface of the body region, *Γ*^{B}. The muscle and body domains can interact through the current *I*_{body} across the muscle–body interface, *Γ*^{M}: = ∂*Ω*^{M}. The muscle domain is modelled as two interpenetrating domains—the intracellular domain and the extracellular domain—that coexist at each spatial point * x* at all times

*t*. The intracellular domain and extracellular domain are separated by a selective-permeable cell membrane and can interact through a current crossing the cell membrane,

*I*

_{m}. Each domain has its own conductivity tensor

*and potential*

*σ**ϕ*.

Assuming quasi-static conditions, the bidomain equations are given by 2.1 2.2

The bidomain equations are a system of two coupled partial differential equations (PDEs) linked to a set of ordinary differential equations (ODEs). The PDEs solve for the extracellular potential *ϕ*_{e}(** x**,

*t*) and the membrane potential

*V*

_{m}(

**,**

*x**t*), which is defined as the difference of the intracellular and the extracellular potential, i.e.

*V*

_{m}(

*,*

**x***t*) =

*ϕ*

_{i}(

*,*

**x***t*) −

*ϕ*

_{e}(

*,*

**x***t*). Further,

*σ*_{i}(

*) and*

**x**

*σ*_{e}(

*) are the intracellular and extracellular conductivity tensors, respectively,*

**x***A*

_{m}(

*) denotes the surface-area-to-volume ratio of the membrane and*

**x***C*

_{m}(

*) is the capacitance of the cell membrane per unit area. The first bidomain equation (2.1) relates changes in the potential to the currents crossing the cell membrane, which consist of time-dependent capacitive currents and ionic currents, i.e. . In general, the ionic currents depend on the membrane potential, further state variables, and the time, i.e.*

**x***I*

_{ion}=

*I*

_{ion}(

*V*

_{m},

*,*

**y***t*). The state variables require to solve an ODE system, which is coupled to the PDEs via

*I*

_{ion}. Note that the second bidomain equation (2.2) corresponds to the Poisson equation used in other EMG models [8,9,12,21], in which the right-hand side of equation (2.2) is given by a signal source term.

For electrically inactive tissue, the first bidomain equation (2.1) vanishes, and the extracellular bidomain equation (2.2) reduces to the generalized Laplace equation [20]
2.3
where *σ*_{o}(* x*) denotes the conductivity tensor, and

*ϕ*

_{o}(

*,*

**x***t*) is the potential of the electrically inactive tissue.

The bidomain equations provide a strong coupling between the membrane potential and the extracellular potential. However, under certain conditions, the bidomain model can be simplified to the monodomain model, which only depends on the membrane potential (§2.3.1). Once the membrane potential distribution has been determined, the extracellular potential can be computed, as described in §2.2.

### 2.2. Propagation of electrical signals through a volume conductor

For a given distribution of the membrane potential, the extracellular potential in the muscle region can be computed using the extracellular bidomain equation (2.2). Further, the potential in the body region is obtained by solving the generalized Laplace equation (2.3).

The electrically inactive body region can be subdivided into several regions with different tissue behaviour, e.g. fat and skin. Using the finite element method, there is no theoretical limitation on the geometry of the regions or on the forms and numbers of heterogeneities. The practical limitation is given through computational limits due to the resolution of the finite element mesh. Furthermore, the approach allows for spatially varying fibre directions throughout the muscle by introducing a local fibre coordinate system [20]. This way, a different fibre direction can be prescribed at each point of the muscle. This is particularly important for biological tissues such as muscles as they often exhibit complex muscle fibre architectures.

To assure continuity of the potential across the muscle–body interface and to make sure that no potential flows out of the body region (figure 1), appropriate interface and boundary conditions are introduced. As only Neumann-type boundary conditions for a Laplacian-type equation yield a singular matrix [20], an additional boundary condition is required. Dirichlet boundary conditions, as proposed in [20], are too restrictive for the EMG. Hence, the following zero-mean condition [22] is applied:
2.4
where d*v* denotes an infinitesimal volume element. Based on the above-mentioned modelling assumptions, the governing equations for the EMG model are summarized in figure 2.

### 2.3. Chemo–electro–mechanical model

Previous computational models predicting the EMG are based on a phenomenological approach to describe the AP, e.g. the Rosenfalck approximation or the impulse response [8,11]. In this contribution, the combination of a biophysical Hodgkin–Huxley-type model of the membrane electrophysiology and a transient diffusion equation describes the generation and propagation of APs along the muscle fibres. A major advantage of the multiscale formulation over existing phenomenological descriptions is that biophysical features emerge from the model rather than being *a priori* prescribed as a part of the model construction.

Although the underlying chemo–electro–mechanical muscle model is described in detail in [18], key aspects of the model are summarized in the following for the sake of completeness. Section 2.3.1 describes the AP propagation along skeletal muscle fibres. The generation of APs and the excitation–contraction coupling in the sarcomeres is presented in §2.3.2. Section 2.3.3 links the model of the excitation–contraction coupling to a continuum-mechanical framework of whole muscle deformation and force generation.

#### 2.3.1. Propagation of action potentials

Assuming that the intracellular and extracellular conductivity tensors have equal anisotropy ratios, i.e. *σ*_{i} = *kσ*

_{e}for some scalar

*k*> 0, the bidomain equations, (2.1) and (2.2), simplify to the monodomain equation [23,24] 2.5 where

*σ*_{eff}=

*σ*_{i}

*σ*_{e}(

*σ*_{i}+

*σ*_{e})

^{−}^{1}is the effective conductivity. Since in skeletal muscle electrical activation from one muscle fibre to adjacent ones does not occur, the propagation of APs along a muscle fibre can be considered a one-dimensional problem. Thus, in the present work, the monodomain equation (2.5) is solved on individual one-dimensional muscle fibres to predict the membrane potential along the fibres [18,25]. It is noteworthy that in the context of a one-dimensional problem, the assumption of equal anisotropy ratios is always satisfied and that the one-dimensional monodomain equation is identical to the cable equation (e.g. [26]).

The ionic currents crossing the muscle fibre membrane, *I*_{ion}, in (2.5) (and in (2.1)) will be defined in §2.3.2.

#### 2.3.2. Biophysical half-sarcomere model

A Hodgkin–Huxley-type description of the ionic currents is used to describe the membrane electrophysiology [19,27,28]. This class of models employs voltage-dependent ionic currents that activate and inactivate in time to predict the AP shape. A two-compartment model is used that distinguishes between the sarcolemma and the T-tubule membrane. This description of the ionic currents is part of a model of the excitation–contraction coupling in a half-sarcomere [19], which combines the membrane model with a model of calcium release from the sarcoplasmic reticulum [29], a calcium-dynamics model [30], a cross-bridge dynamics model [31] and a description of metabolic fatigue [19]. The mathematical representation of the excitation–contraction coupling leads to 56 ODEs per simulated half-sarcomere. Two sets of parameters account for slow-twitch and fast-twitch muscle fibres [19].

A major advantage of the biophysical approach is that one can investigate the effect of physiological properties, such as fatigue, on the EMG signal. For example, the magnitude and propagation velocity of APs decrease in fatigued conditions [16,17]. Such effects cannot be taken into account, when using the Rosenfalck approximation or the impulse response in conjunction with a constant AP propagation velocity. Neglecting these effects, however, might lead to inaccurate model predictions and misinterpretations of EMG recordings.

A further advantage of the proposed method is the fact that the half-sarcomere model does not only describe the membrane electrophysiology, but also describes the active stress generation within the half-sarcomere. The sarcomere-based active stress is integrated in the macroscopic continuum-mechanical description (see next section).

#### 2.3.3. Continuum-mechanical muscle model

To describe the force-generating capabilities on the whole muscle scale and the resulting muscle deformation, a continuum-mechanical skeletal muscle model based on the theory of finite elasticity is employed. The model has been proposed in [18] and is briefly summarized for completeness. Neglecting inertia and body forces, the momentum balance equation reduces to div* T* =

**0**, where

*denotes the Cauchy stress tensor. The Cauchy stress tensor, which is defined in the actual configuration, is related to the second Piola–Kirchhoff stress tensor,*

**T***, of the reference configuration via a scaled covariant push forward operation:*

**S***= (det*

**T***)*

**F**

^{−}^{1}

**FSF**^{T}. Therein,

*denotes the deformation gradient tensor, which maps referential line elements d*

**F***to line elements d*

**X***in the actual configuration, i.e. d*

**x****=**

*x**d*

**F***. In this work, the second Piola–Kirchhoff stress tensor for the muscle tissue,*

**X**

**S**_{M}, consists of an isotropic part based on the Mooney–Rivlin material,

**S**^{iso}, a term appealing to stretch in the fibre direction,

**S**^{ani}, [32] and a term representing the muscle's ability to actively generate force,

**S**^{act}: 2.6 Therein,

*=*

**C**

**F**^{T}

*denotes the right Cauchy–Green deformation tensor,*

**F***p*is the hydrostatic pressure,

*is the second-order identity tensor and*

**I**

**a**_{0}is a unit vector in fibre direction defined in the reference configuration. Further, is the fibre stretch with

*I*

_{4}=

**a**_{0}·

**Ca**_{0}being the fourth (mixed) invariant of

*, and denotes the force–length relation. Moreover,*

**C***P*

^{max}is the maximum active stress at optimal fibre length and under isometric conditions, and denotes the normalized active stress, which is obtained through homogenization of the sarcomere-based active stress determined as part of the biophysical half-sarcomere model (§2.3.2). For further details, see [18].

The fat and skin tissues are assumed to behave mechanically isotropic and are modelled using the isotropic and incompressible Mooney–Rivlin material. The second Piola–Kirchhoff stress tensor for the fat and skin is given by 2.7

The material parameters for the constitutive models are listed in table 1.

### 2.4. Computational framework and implementation

Following the monodomain approximation, the computation of the membrane potential (2.5) and the equation for solving for the extracellular potential (2.2) are decoupled. Hence, the solution process is divided into two parts. Both parts are implemented into the open-source software library OpenCMISS [37].

The monodomain equation in conjunction with the biophysical cell model and the governing equations of finite elasticity are solved in an integrated fashion. Details on the implementation and solution process of the resulting multiscale model can be found in [18]. In brief, different discretizations are used for the continuum-mechanical subproblem and the monodomain equation. The geometry of the muscle, fat and skin are meshed using three-dimensional finite elements. Within the three-dimensional mesh, fine-spaced one-dimensional meshes representing individual muscle fibres are embedded (figure 3). The three-dimensional elements are used for the solution of the continuum-mechanical subproblem, while the monodomain equation is discretized and solved on each of the one-dimensional muscle fibre meshes.

To solve for the extracellular potential, see (2.2), a fine three-dimensional finite element mesh is necessary. This three-dimensional mesh is created by connecting the nodes of the one-dimensional muscle fibre meshes such that hexahedral linear Lagrange finite elements are obtained (figure 3).

The previously computed membrane potential enters the equation for the extracellular potential (2.2) as a source term. For each time step, the spatial positions of the finite element nodes and the membrane potential are updated using the values obtained from the chemo–electro–mechanical model.

### 2.5. Material parameters

The parameters of the mechanical muscle model and the biophysical half-sarcomere model are taken from [18] and [19], respectively. The same Mooney–Rivlin parameters as for the passive muscle tissue are used to describe the mechanical behaviour of the fat/skin tissue. This choice yields a very compliant behaviour for the fat/skin tissue. As far as the electrical conductivities are concerned, a large variation of parameter values has been reported for biological tissue. For example, the ratio of the longitudinal conductivity to the transversal conductivity in skeletal muscle has been reported to range from 2.78 to 14.5 [8]. In this work, an anisotropy ratio of 10 is assumed for the muscle's intracellular conductivities, while the muscle's extracellular conductivities were assumed to be isotropic. This yields for the conductivity tensor on the left-hand side of (2.2), *σ* =

*σ*_{i}+

*σ*_{e}, an anisotropy ratio of 2.06. The fat and skin tissues are assumed to be electrically isotropic.

The material parameters used in the presented simulations are summarized in table 1.

## 3. Results

While it is beyond the scope of this paper to investigate biophysical changes in the volume conductor that might result, for example from tissue deformation or fatigue, representative simulations are presented demonstrating the capabilities of the multiscale framework. For all presented test cases, monopolar detections of the EMG are reported.

### 3.1. Surface electromyography due to single-fibre action potentials

The first test case aims to provide some validation of the proposed model. To this end, simulation results are compared to existing EMG models [8]. As these models are not capable of considering any deformation, the analysis is carried out on a rigid rectangular cuboid with a length of *l* = 6 cm (*x*-direction), a width of *w* = 2.9 cm (*y*-direction) and an overall height of *h* = 1.4 cm (*z*-direction). The lower 1.2 cm of the cuboid are defined as transversely isotropic muscle region that is covered by a 0.2 cm thick isotropic fat layer. Within the muscle region, a total of 30 × 13 = 390 equally distributed muscle fibres are embedded. To simulate a parallel-fibred muscle, all fibres are chosen to run in parallel with the cuboid's long edge. Further, all muscle fibres run the entire length of the cuboid. The innervation zone is assumed to be located at the midpoint of the muscle fibres, leading to symmetry with respect to the innervation zone.

For the analysis, an AP is evoked in the middle of a single fibre, and the EMG resulting from the propagating signal is observed at the surface. In repetitive simulations, the distance between the activated fibre and the surface is successively increased from 0.2 to 1.4 cm (from the muscle–body interface to the bottom of the muscle).

Figure 4 shows that the absolute surface potential decreases exponentially with increasing distance between the stimulated fibre and the surface. The result coincides with the findings of Lowery *et al*. [8], and also confirms the findings of Holobar *et al*. [38].

### 3.2. Surface and intramuscular electromyography during fixed-length contractions

Based on the same geometrical model as above, surface and intramuscular EMG signals during fixed-length contractions are analysed. To this end, one assumes for the mechanical problem that all degrees of freedom that are associated with the surfaces that exhibit normals in *x*-direction, are fixed in *x*-direction. Furthermore, for one surface with normal in *y*-direction (*z*-direction) all degrees of freedom in the *y*-direction (*z*-direction) are fixed to mimic symmetric conditions across these surfaces.

In contrast to the above-described test case, the 390 embedded fibres are now randomly assigned to 10 MUs (figure 5). Note that the MU size and territories are randomly chosen and do not reflect reality. As all muscle fibres use the fast-twitch parametrization of Shorten *et al*. [19], the AP propagation velocity is initially identical in all fibres. The MU fibres are stimulated at their midpoint according to the stimulation protocol depicted in figure 6.

Figure 7 shows a sequence of sEMG signals at the onset of the simulation. Starting from the middle of the muscle tissue, where the fibres are stimulated, tripoles are propagating towards each end of the muscle. Despite the simplified example set-up, the predicted EMG compares qualitatively well with the experimental sEMG data of Farina *et al*. [1] and Barbero *et al*. [39]. A quantitative comparison is difficult as experimental EMG signals depend heavily on properties that are difficult to control in experiments, for example the MU territories and the thickness of the subcutaneous tissue [40].

Figure 8 displays the membrane potential in the muscle fibres and the resulting sEMG signal at time *t* = 36.5 ms. Although an isometric contraction is considered, the activation-induced deformation of the muscle tissue is clearly visible at the skin surface. The regular pattern observed in the sEMG is due to the activation protocol, which, for the sake of simplicity, considers stimulation times only with a resolution of 5 ms (figure 6).

In addition to sEMG predictions, the EMG signal can also be reported at any position within the volume conductor. Figure 9 shows the evolution of the raw and rectified EMG signals at a point of the surface and at three points within the muscle. The decrease of the amplitude of the potential due to the volume conducting fat layer is clearly visible. Although a simplified stimulation protocol is used (figure 6), the simulated signal compares qualitatively well to experimental EMG recordings (e.g. [5]).

### 3.3. Fatigue simulation

The aim of the third test case is to investigate the effect of fatigue on the sEMG signal. The biophysical half-sarcomere model employed in the presented multiscale framework is capable of predicting potassium accumulation and sodium depletion within the T-tubules as a result of high-frequency stimulation. As a consequence, a decrease in the amplitude of the membrane potential and a reduction of the AP propagation velocity are expected [16,17,19].

For the numerical experiments, the above-described cuboid is extended by a 1 mm thick layer of skin tissue on top of the fat layer. As before, the biophysical half-sarcomere model uses the fast-twitch muscle fibre parametrization of Shorten *et al*. [19]. To clearly demonstrate the effect of fatigue on the sEMG, only a single fibre located in the middle of the cuboid and 0.4 cm below the skin surface is stimulated. The stimulation with a frequency of 100 Hz is applied at the fibre's left end. Again, for the sake of simplicity, the tissue deformation is neglected in this simulation.

Figure 10 shows the first simulated APs propagating along the stimulated fibre and the resulting sEMG signal on top of the skin surface. The decrease in the amplitude of the potential from the first to the second AP is clearly visible in the sEMG signal. Note that this decrease is due to the high stimulation frequency and becomes less pronounced with larger interstimulus intervals (not shown).

Figure 11 illustrates the effect of fatigue on the sEMG signal at a position in the middle of the skin surface. The amplitude decreases from 0.37 mV to 0.22 mV after 500 ms, which corresponds to a decrease of 40%. This compares well to an experimentally determined mean amplitude reduction of 32% [16].

Furthermore, to investigate changes in the propagation velocity due to fatigue, figure 12 depicts the potential distribution in fibre direction at the surface along the cuboid's centre line at times 30.4 and 500.4 ms. These time instants are chosen such that, in both cases, the next stimulation occurs at the left end of the fibre.

By comparing the sEMG of the non-fatigued state (30.4 ms) with the fatigued state (500.4 ms), one observes that both the propagation velocity and the amplitude of the signal decrease with time: stimulating the fibre for 500 ms at a frequency of 100 Hz, the conduction velocity decreases from 2.0 m s^{−1} to 1.7 m s^{−1} or by 15%. This decrease within 500 ms can be compared to the experimentally determined decrease of approximately 30–50% obtained after 2 min of stimulation [17].

### 3.4. Tibialis anterior

The fourth test case demonstrates that the proposed model is capable of simulating the EMG signal under realistic conditions. To this end, a realistic geometry, a realistic stimulation protocol, different muscle fibre twitch-types and a distributed innervation zone are included in the model.

The geometry of the tibialis anterior muscle (TA) and the fat/skin tissue is based on the visual human male dataset [41]. The bipennate muscle fibre directions are based on DT-MRI data [42]. The geometry of the muscle without the fat/skin tissue has previously been used in [18,36,43]. To reduce complexity and computational time, only 2700 muscle fibres are considered that are grouped into 10 MUs. The general assumption that a muscle typically consists of many slow-twitch (type-I) MUs and few fast-twitch (type-II) MUs [44] is reflected within the model by choosing the first six MUs to use the slow-twitch parametrization, while the other four MUs (MUs 7–10) use the fast-twitch parametrization of Shorten *et al*. [19]. Note that the different parametrizations of the biophysical half-sarcomere model lead to different AP propagation velocities, which are also observed in real muscles and significantly influence the sEMG. Further, an exponential distribution of the innervation number was assumed, where the largest MU had 100 times as many fibres as the smallest MU [44]. The MU fibre distribution is depicted in figure 13, where a different colour is chosen for each MU.

While the MU territories in real muscles are spatially confined to small regions within the muscle's cross-sectional area, this model contains only a random distribution of the fibres. This assumption is made to simplify the model set-up and is not due to any framework-inherent limitations. An algorithm to distribute the MU fibres within the muscle volume is provided in [36].

The MU discharge times are determined using the biophysical model of the *α* motor neurons of Negro & Farina [45], which is based on the description of the motor neurons of Cisi & Kohn [46]. A feature of the motor neuron model of Negro & Farina [45] is that it intrinsically accounts for Henneman's size-principle of recruitment [47,48]. The input to the motor neuron model is a constant mean synaptic current of 0.005 μA cm^{−2} superimposed by two Gaussian-distributed high-frequency oscillating signals, one common to the entire motor neuron pool and one independent for each motor neuron (see [45] for details). Owing to the total synaptic input signal, eight out of the 10 MUs are recruited. The mean discharge frequencies of the recruited MUs are approximately 10–15 Hz. The interstimulus interval exhibits a coefficient of variation (ratio between standard deviation and mean) of approximately 20%. Figure 14 shows the MU discharge times as predicted from the motor neuron model [45].

To model the innervation zone, a Gaussian distribution around the middle of the fibres with a standard deviation of two nodes was assigned to the stimulation point. The maximum deviation was seven nodes, which corresponds to a 0.56 cm wide spread of the innervation zone.

To simulate a fixed-length contraction, all finite element nodes at the proximal and distal ends of the muscle are fixed in the mechanical model.

While the activation-induced contraction of the TA is based on the entire muscle and the surrounding fat/skin tissue, the subsequent EMG computation is only based on the superficial part of the TA covered by a layer of fat and a layer of skin tissue (black lines in figure 13). The TA's superficial part contains 900 embedded fibres. Owing to its larger distance to the skin surface, the deep TA is considered to have a minor contribution to the sEMG signal. It is, however, noteworthy that the presented modelling approach is not limited to the superficial part of the TA. Along the superficial TA, the thickness of the fat and the skin layer varies slightly, i.e. the fat layer has a thickness of approximately 6 mm and the skin layer thickness is approximately 1.5 mm. The finite element mesh for the computation of the EMG consists of 47 908 elements.

The resulting sEMG signal at the skin covering the superficial TA is visualized in figure 15 for three different time steps. To illustrate the activation-induced deformation resulting from the fixed-length contraction of the muscle, the deflection is enhanced by a factor of 5.

As the TA is a pennate muscle, the generated EMG signal is different from the one shown in §3.2, where the fibres are parallel to the skin. Owing to the fibre angle, no pronounced propagation of the potential at the surface can be detected. Rather, the signals propagate along the fibres towards the surface leading to an increase in the surface potential. As reported in [39], this example clearly demonstrates that the end-of-fibre effect and the extinction of the AP are the main contributions to the sEMG of pennate muscles. While these experimental findings are qualitatively reflected by the model results, they make a quantitative comparison difficult. This likewise holds for the mechanical deformation. Experiments using ultrasound together with EMG measurements are performed in [49,50] to study the relation between architectural parameters, such as the pennation angle, and the EMG. As expected from the experimental findings of Hodges *et al*. [49], the pennation angle increases during the isometric contraction. While the model simulates the separated TA, the experiments are performed *in vivo* including interactions of the TA with adjacent tissue. This makes further comparisons within the scope of this work impossible.

## 4. Discussion

The presented multiscale framework for the simulation of EMG signals combines several processes ranging from the subcellular level up to the whole muscle level. Unlike previous, phenomenological models of the EMG, which prescribe the shape and propagation velocity of the AP as part of the model construction, the presented description of the AP generation and propagation is based on a transient diffusion equation in conjunction with a biophysical Hodgkin–Huxley-type model of the membrane electrophysiology. One of the major advantages of the biophysical description is that it intrinsically accounts for physiological effects such as, for example membrane fatigue that causes a change in the amplitude and the AP conduction velocity during sustained contractions. Thus, using the presented multiscale model, one can analyse their effect on the EMG signal, which might improve their interpretation and lead to a better understanding of recorded signals.

The model of the membrane electrophysiology is part of a biophysical model of the excitation–contraction coupling in the half-sarcomeres [19]. As the half-sarcomere model describes the entire signalling pathway from AP generation via calcium release and calcium dynamics to stress generation, the effect of changes in the signalling pathway on the EMG, as might occur, for example in pathological conditions, can potentially be investigated using the multiscale framework. A further potential application of the model is the *in silico* testing of drugs. For example knowing the effect of a certain medication on the conduction of a species of ion channels, its effect on the AP shape and conduction velocity can be studied, and the simulation results can be validated using EMG measurements.

Furthermore, the description of the excitation–contraction coupling in the half-sarcomeres yields a measure for the locally generated active stresses. These sarcomere-based active stresses are incorporated in a continuum-mechanical framework leading to a description of the muscle force generation and the contraction-induced muscle deformation. Hence, the proposed multiscale model is not restricted to isometric contractions like previous models of the EMG, but the chemo–electro–mechanical model can simulate any kind of contraction. Therefore, the model can be used to study the EMG during non-isometric contractions. Furthermore, the effect of changes in the local fibre orientation or a shift of the innervation zone with respect to the skin surface during isometric and non-isometric contractions on the sEMG signal can be investigated [51–53]. These effects influence the EMG, and thus neglecting their influence might lead to erroneous signal interpretations. Using the presented model, one can study these effects in order to take their influence into account when interpreting recorded EMG signals, which might yield to more accurate interpretations. No previous model of the EMG could combine all these processes within one framework.

Moreover, lengthening contractions of skeletal muscles require a neural activation strategy that is different (and currently unknown) from that during isometric contractions [54]. Using the presented framework, activation strategies based on neurophysiological hypotheses can be tested and the simulation results can be compared to EMG recordings during lengthening contractions. This might help to reveal the neural strategy underlying lengthening contractions.

Before a computational model can be applied to biophysical questions, however, it needs to be validated. The validation of the biophysical half-sarcomere model is presented in [19]. A partial validation of the multiscale muscle model is provided in [55], where the entire active behaviour was defined at the microscopic half-sarcomere level, and yet the experimental force–length and force–velocity behaviours of skeletal muscles were recovered at the macroscale. Furthermore, the convergence behaviour of the multiscale approach was demonstrated in [25] by maintaining a fixed number of embedded fibre meshes while successively refining the number of three-dimensional mechanical elements until homogenization is no longer required. The multiscale model had convergence properties similar to those of the mechanical-only problem [25].

The presented model for the simulation of the intramuscular and surface EMG is based on the bidomain model, which is the continuum approximation of the electrophysiological behaviour of electro-active biological tissue [20]. A monolithic coupling of the bidomain model with a detailed biophysical cell model in the muscle and the generalized Laplace equation in the surrounding tissues leads to computationally expensive simulations. One common approach to decrease the computational complexity of the bidomain model is to use the less complex monodomain model [20,23,24]. This approach yields an exact method if the assumption of equally anisotropic extracellular and intracellular conductivity tensors is satisfied, and an approximation if not [23]. By successively solving the monodomain equation and the extracellular bidomain equation, the effect that the extracellular potential might have on the membrane potential is neglected. While the authors believe that this effect is rather small, this limitation of the current model should be further addressed in future works.

Using the monodomain model, the equations for the membrane potential and for the extracellular potential are decoupled from each other, and hence can be solved successively. The successive solution of the monodomain model and the extracellular bidomain equation essentially corresponds to the approach used in previous models of the EMG [8,9,12,21]. In these models, the membrane potential distribution is prescribed and a Poisson equation is solved for the extracellular potential.

Previous models of the EMG, however, do not allow the consideration of arbitrary fibre architectures. Further, existing models are not capable of predicting changes in the AP shape and propagation velocity due to mechanical deformation of the tissue or due to underlying biophysical processes such as membrane fatigue. Other drawbacks of the presented model are the same as for existing models, e.g. they are limited by a lack of reliable experimental data such as accurate descriptions of the fibre directions, material parameters, heterogeneities and MU territories. Because of this and the fact that recorded EMG signals vary at lot, for example due to differences in the thickness of the subcutaneous tissue [40], a quantitative validation of the EMG computation is difficult. The bidomain model, and its simplified version, the monodomain model, however, are well established within the field of biosignal modelling, in particular for simulating the electrical activity of the heart [20]. In this context, Vigmond *et al*. [56] demonstrated that a coarser mesh can be used for the (elliptic) extracellular bidomain equation than for the (parabolic) monodomain equation while still maintaining reasonable accuracy. A convergence analysis for the AP propagation equation in skeletal muscle is presented in [57]. These results also influenced our choice of the mesh size for the monodomain equation. Choosing the same mesh size for the extracellular bidomain equation would not have been necessary [56] but appeared to be, from an implementational point of view, the simplest choice.

Owing to the complexity of the multiscale model and its associated computational cost, the number of muscle fibres that can currently be included in a simulation on a normal desktop PC using a single processor is limited to a few thousand. The human TA consists of about 250 000 fibres [58], and thus only about 1% of the fibres can currently be simulated. However, there exists a large potential for parallelization [18] and model reduction techniques to decrease the computational load and/or allow to include a more realistic number of fibres in the simulations.

Although strategies have been employed to reduce the computational complexity of the model, the computing time is still high for the presented model. On a normal desktop PC (Intel Core i5–3470 CPU, single processor, 3.2 GHz and 32 GB of memory), the solution of the chemo–electro–mechanical model on the TA (2700 fibre meshes, each of which required the solution of 50 biophysical half-sarcomere models, yielding more than 7.5 × 10^{6} degrees of freedom for the bioelectrical problem) required about 3 h to simulate 0.2 ms. The computation of the corresponding EMG was about 15 min. It is, however, noteworthy that about 90% of the EMG computing time is spent on file I/O and only 10% is required for the actual computation of the EMG signal. This is due to the fact that the membrane potential and the extracellular potential distributions are currently stored in ASCII files, which have to be read in and written out, for each time step. A more sophisticated implementation that eliminates the cumbersome I/O of the membrane potential is feasible, as the entire model is solved within a single framework (OpenCMISS). This, however, has not been realized yet. As the application of the model is currently limited by the computation time, this issue has to be addressed in the future.

## 5. Conclusion

A finite element model for the prediction of the intramuscular and surface EMG has been presented. The multiscale model contains (i) a biophysical description of the excitation–contraction coupling at the cellular level, (ii) a model of the AP propagation along muscle fibres, (iii) a continuum-mechanical formulation of the force generation and deformation of the muscle, and (iv) a model of the intramuscular and surface EMG. The overall framework provides novel contributions and new applications to the field of modelling EMG signals through an integrated approach that couples a microscopic Hodgkin–Huxley-type model of membrane electrophysiology to the macroscopic monodomain/bidomain model. In particular, the inherent link from subcellular properties to macroscopic quantities is a major advantage over any other currently existing modelling framework. The coupled formulation intrinsically accounts for physiological properties of muscle such as membrane fatigue, and hence allows to investigate their influence on the EMG signal. Furthermore, previous methods for simulating the EMG did not take into account the contraction-induced deformation of the muscle, which is included in the proposed framework. Further, employing the finite element method for the discretization of the geometry and the governing equations, the model can easily be applied to realistic geometries and can include different types of heterogeneities.

## Funding statement

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007–2013)/ERC Grant Agreement no. 306757 (LEAD).

## Footnotes

One contribution of 11 to a theme issue ‘Multiscale modelling in biomechanics: theoretical, computational and translational challenges’.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.