A computational model to investigate the effect of pennation angle on surface electromyogram of Tibialis Anterior

This study has described and experimentally validated the differential electrodes surface electromyography (sEMG) model for tibialis anterior muscles during isometric contraction. This model has investigated the effect of pennation angle on the simulated sEMG signal. The results show that there is no significant effect of pennation angle in the range 0° to 20° to the single fibre action potential shape recorded on the skin surface. However, the changes with respect to pennation angle are observed in sEMG amplitude, frequency and fractal dimension. It is also observed that at different levels of muscle contractions there is similarity in the relationships with Root Mean Square, Median Frequency, and Fractal Dimension of the recorded and simulated sEMG signals.


Introduction
Tibialis Anterior (TA) is essential for posture and gait stability. Age-associated weakness of TA is known to be a major cause of falls [1][2][3]. As there is no direct way for measuring TA contraction force, sEMG provides a useful tool to measure the electrical activity of the muscle associated with the force of contraction.
Numerical modelling has been considered to improve the interpretation of the relationship between neuromuscular parameters, generated force and sEMG signal [4] and these models have demonstrated the relationship of many parameters with the signal [5][6][7][8][9]. However, their applications are very limited because of number of assumptions involved, such as the parallel orientation of muscle fibres to the surface [9]. These assumptions make the models unsuitable for Tibialis muscles which are pennate. In pennate muscles, fibres run at an angle to the axis of traction. The anatomical cross sectional area (ACSA) does not represent the cross section perpendicular to all fibres in the muscle, i.e., the physiological CSA (PCSA) [10]. The maximum force of a muscle depends on its PCSA rather than its ACSA [10,11]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Modelling sEMG of TA requires consideration of pennation angle of the muscle [11]. SEMG signals are recorded using differential electrodes as it enhances the signal to noise ratio [12]. However, some researchers have investigated pennate muscles [13][14][15][16][17], they have not considered the differential electrodes and made assumptions on the shape of the signals [18]. While the pre-defined shape of MUAP is reasonable for muscles such as biceps, it may lead to erroneous results because of the possible age-associated change in muscle shape due to pennation angle in TA. It is reported that power spectrum depend on shape of MUAP [19,20]. Howard et al [21] reported that properties of action potential changes due to aging. Hence, for the model to be effective for investigating changes to sEMG with ageing or disease, it is essential for the model to provide the examiner the accurate estimate of the shape of SFAP.Zuniga et al [22] reported that there is no significant difference in the EMG amplitude and mean power frequency MPF responses for EMG recorded from the electrode oriented parallel and perpendicular to the muscle fibres of the vastus lateralis muscle (VL). Shi et al [23] reported exponential relationship between the RMS and pennation angle of brachialis. These studies have considered only the amplitude and spectral based features and have not investigated the features describing its chaotic nature. However, muscle contractions are reported to exhibit fractal characteristics [24,25]. We report a computational model that describes the generation of sEMG signals for TA muscle where we have investigated the effect of pennation angle on sEMG recorded using differential electrode. The model has been validated experimentally based on the temporal, spectral and chaotic signal properties by comparing simulated and experimental recordings using the three features; Root Mean Square (RMS), Median Frequency (MDF) and Fractal Dimension (FD).

Description of the model
The model is based on the works reported by Wheeler et al [9] and Siddqi et al [18] with two significant enhancements; i) differential electrodes and ii) pennation angle. The use of differential electrodes overcomes the need for pre-defined shape of Motor Unit Action Potentials (MUAP), while the second allows the investigation of muscles such as TA.
Volume conduction. The motor-unit territory cross-section was considered to be circular and described using Cartesian coordinate system located at the centre of muscle. The crosssection area was obtained based on the type of muscle fibres and number of fibres in the muscle. The electrodes were separated from the muscle by two layers; skin and subcutaneous [7]. The innervation zone was located at the centre of the fibre length. The model is described in Fig 1. Muscle fibre action potential was calculated by Rosenfalk Eq 1 [8,26]. The current distribution between the inside and outside of the muscle cell was modelled as the second spatial derivative of action potential voltage V m [27] given in Eq 2.
where, I m (t) is trans membrane current, C is a constant and is given as , σ i is internal conductivity, σ e is external conductivity, d is fiber diameter, v is the conduction velocity, A is the amplitude of action potential, B is resting membrane potential, λ is a scaling factor, z e is the distance between the electrode and neuromuscular junction and σ is assumed to be direction independent. Earlier models [9,18] have considered a single electrode with a pre-defined tri-phasic motor unit action potential. However, such models are limited because these cannot be used to investigate the effect of electrode placement, inter-electrode distance and pennation angle. To overcome these limitations, the model described in this paper has considered differential electrodes and the TA muscle has been modelled as uni-pennate structure [18] with angle θ to the surface. Thus, Eq 3 which describes the volume conduction effect in the muscle has been modified to Eq 4.
Motor unit recruitment and rate coding. In the model described in this paper, the motor-units (MU) are recruited according to the size principle [28]. All the motor units are recruited at a given force recruitment range (RR) according to Fuglevand model [29]. RR is expressed in percentage of MVC. The i th MU is activated once the neural drive crosses its recruitment threshold (RT i ) which is based on the force recruitment range (RR) based on Fuglevand model [29].
where TMu is the total number of motor units and 'i' is the motor unit index based on its size rank. RR is the recruitment range and is assumed to be the voluntary input based on fraction of MVC. The firing rate is controlled by input force RR shown in Fig 2 and given in Eq 6 (3). Firing rate of active motor-units has been considered to increase linearly from Minimum Firing Rate (MFR) to Peak Firing Rate (PFR) as a function of the fraction of maximum voluntary contraction of the muscle. It is assumed that all motor units reach their PFR at 100% MVC. MFR ranges from 7 to 23 Hz while PFR ranges from 14 to 50 Hz for human skeletal muscles (5). This has been summed in Eqs 6 and 7; where FR i,RR represents the firing of i th motor unit at input force level RR and α i represents the slope associated with i th motor unit;

Model simulation
The values of the parameter were taken from experimental studies reported in the literature and are listed in Table 1 [30]. The model was simulated 20 times with values randomly obtained from the range shown in Table 1 such that it represented more realistic values and was comparable with the experimental data.

Experiments
Five young volunteers with age-range: 23-30 years Mean (±SD): 26.1 ±2.9 years and body mass index-range: Mean (±SD): 22.3 Kg/ m 2 +2.9 with no clinical or self-reported history of neuromuscular disease or ankle injury participated in this study. Seven participants were approached, and two subjects did not participate in the study. All research involving human participants has been approved by Royal Melbourne Institute of Technology (RMIT) University College Human Ethics Advisory Network (CHEAN) committee and all clinical investigations have been conducted according to the principles expressed in the Declaration of Helsinki. Prior to the participation in the experiment, written informed consent was obtained from the participants. The experimental set-up and protocol for recording surface Electromyography used in this study have been reported in Siddiqi et al [18]. SEMG activity was recorded from TA (Fig 1) muscle using SENIAM configuration. The ground electrode was placed at the patella. The locations for electrodes were shaved, abraded and cleansed with an alcohol swipe to ensure good connectivity. The signals were recorded using the Delsys Myomonitor 4 (Delsys, Boston) with fixed gain of 1000, CMRR of 92 dB and bandwidth of 2-450 Hz, with 12 dB/ octave rolloff. The sampling frequency was 1000 Hz with a resolution of 16 bits/ sample. The Delsys single-channel active differential surface electrodes with an embedded preamplifier and interelectrode distance of 10 mm were used.
Participants sat with their right leg strapped, hip at 90˚, knee at 140˚and ankle at 90˚. SM100-type strain gauge force sensor (Interface S type) was used to immobilize the foot-plate and measure the isometric force produced during dorsiflexion or plantar flexion. Its output was displayed to give visual feedback to the subject for maintaining steady contraction. Participants performed an isometric contraction at 25%, 50%, 75% of their maximum voluntary contraction (MVC) for 15 seconds each. Throughout the experiment, the left leg was planted firmly on the ground. Absence of heel lift, and foot or toe movement during plantar flexion and dorsiflexion was ensured with the foot and ankle secured to the footplate.

Data analysis
The experimentally recorded and simulated sEMG signals were analysed to obtain three features; RMS, MDF and FD. While most of the earlier studies have investigated the RMS and MDF of simulated sEMG, these suffer from large inter-subject variation which makes it difficult for validating the model based on experimental recordings. FD is an important feature and is an intrinsic property of a muscle [37,38] and has lower variability compared with other parameters. To determine the relationship of the angle of pennation with the signal property, the model was simulated for different values of θ and single fibre action potential (SFAP) was obtained. Correlation analysis was conducted to identify the similarity between SFAP for fibres with θ pennation angle with the parallel fibres (θ = 0). Consider sEMG with parallel muscle fibres is x and sEMG corresponding to fibres with pennation angle of θ is y; then r xy = Sample correlation coefficient, s xy = Sample covariance, s x = Sample standard deviation of x, s y = Sample standard deviation of y Correlation coefficient, given by The following features of sEMG were computed to compare the simulation and experimental data: x i represents the i th sample of the signal while N is the window size which was taken to be 2500 ms based on [39][33].

Median Distribution Frequency (MDF):
P j is the power spectrum at frequency bin j and M is length of the frequency bin(33).

Higuchi's Fractal Dimension (FD)
: FD was computed using the algorithm described in [40,41] and briefly described below.
Considering y (n) is the signal for epoch n, n = 1 to N, and this is sub-divided in k segments,y k m y k m ¼ fyðmÞ; yðm þ kÞ; yðm þ 2kÞ; . . . :yðm þ where m represents the initial time and k is interval time. b NÀ m k c denotes integer part of NÀ m k The average length L m (k) is calculated as The total length of the entire epoch in k time interval is calculated as L(k) is proportional to k -D for total average length for scale k. D is the FD in Higuchi's method. The slope of the least squares linear best fit of the curve ln(L(k)) versus ln(1/k) estimates the fractal dimension.
Tuning of k parameter is necessary in Higuchi's algorithm because the correct selection of k parameter has an important role to obtain reliable FD values [41] (36). In this work, k = 6 is used to get the best estimate of fractal dimension [40,41]. Correlation analysis was performed to identify the similarity in the shape of the SFAP for different pennation angles.
Statistical testing.
a) Paired sample t-test: RMS, MDF and FD are expressed as mean (standard deviation), and compared using two-sample t-test for sEMG signals simulated and experimental data. The differences were considered significant at p<0.05.

b) Two One Sided Test (TOST):
The similarity of the experimental and simulation data was tested using Two One Sided Test (TOST) [42]. The test was performed for the three levels of MVC and for the three features independently; RMS, MDF and FD. Bland Altman plots are used to select the lower bound and upper bound for TOST.

Result
The effect of pennation angle on the shape of simulated single fibre action potential (SFAP) for θ ranging from 0 to 20 degrees is shown in Fig 3. The pennation angle is increased up to 20 degree because maximum pennation angle for TA is reported 20 degree in literature [22] (22). It is observed that the shape of SFAP is similar and triphasic for all the four values of θ. Table 2 represents the correlation coefficients and P values between shape of SFAP of parallel muscle and pennate muscle with different pennation angle. It is observed that the correlation is high for all the angles in the range (r>0.9). The relationship of pennation angle with amplitude, frequency and fractal dimension (FD) of sEMG is shown in Figs 4-6 respectively. From these figures, it is observed that amplitude, frequency and FD of sEMG increase with increased pennation angle, while the shape remains unchanged. Table 3 shows the mean values of RMS, MDF and FD for the three different levels of %MVC for simulated and experimental data. Two sample t-test is performed between  experimental and simulated data for each level of MVCS and corresponding p-values are given. From the Box-whisker plots in Figs 7-9, it is observed that while the three features of the experimental and simulated data follow similar trends, one common difference is that the experimental data has larger range of the quartiles compared with the simulated data.
The results of the similarity test performed by TOST are shown in Table 4. To select the lower and upper boundaries for equivalent test, Bland Altman plots are given for each feature in different MVCs as shown in Figs 10-12. This test shows that there is similarity between the two data sets; simulated and experimental, within the intervals (±2 SD) at 90% confidence interval as observed in Bland Altman plots. The corresponding 'p' values shown in the Table 4 indicate that the two data are significantly similar within the range.

Discussion
We have investigated the effect of pennation angle on the shape of SFAP and observed that there is no significant difference between the signals for θ ranging from 0 to 20˚. This supports the approximation which was proposed by Siddiqi et al [18]where TA had been approximated to be parallel fibre muscles. Amplitude of sEMG increases as pennation angle increases from 0 degree to 30 degree and this is in line with literature [23]. Studying the effect of pennation angle is important when comparing sEMG recorded from young and older cohort since it is known that pennation angle changes with age for pennate muscle [11]. However, this analysis suggests that there is no change in the shape of the SFAP with pennation angle, and this would indicate that the age-associated change in sEMG is not due to this factor. From Figs 7-9 and Table 3, it is observed that the trend of RMS, MDF, FD of simulated and experimentally recorded sEMG are similar. While the RMS results are in agreement with previous studies [43,44] for biceps brachii muscle and as has been discussed by Moritani et al [45], no earlier study has examined the FD of simulated sEMG. The results show that while the trends of MDF of the simulated and experimental recordings with respect to the fractional MVC are similar, the absolute values are different. There may be number of reasons that could contribute to this difference. For example, the spectrum of sEMG is highly dependent on the electrode placement and inter-electrode distance [46] and hence comparing the absolute values of the MDF may be difficult.
One common difference between the simulated and experimental data is that the quartile range for experimental is much higher than the simulated. Moreover, this study has considered the parameters to belong to a range of values rather than single value and repeated simulations have been initialized randomly within the range, the results indicate that the actual variation appears to be larger. This may be because we have assumed normal distribution which would be require very large number of experiments while this study reports the data from five participants. This study has investigated the FD of simulated and experimentally recorded sEMG. While the trends of the two are similar, there is a difference in the absolute values. The FD corresponding to the experimental data is higher than that of the simulated data. It is known that the Higuchi's FD is dependent on the complexity of the signal sources and this difference may be attributable to the necessity of models being simplified versions of the real system.
The study has some limitations as it has experimentally tested only with small number of healthy participants. The model reported in this study is two-dimensional, and future studies will investigate the influence of physiological conditions such as aging and disorder on the skeletal muscle system with large number of participants.

Conclusions
This study has developed a computational model that describes sEMG suitable for pennate muscles. The model was validated for the TA muscle by comparing the experimental with the simulated sEMG for different fractions of MVC. Three features were considered; RMS, MDF and FD. The results show that while the relationship between these features and fraction of MVC are similar between simulated and experimental data, there is difference in the absolute values. This indicates that the simulated values are subjected to normalization and while it can be used to predict the trends, but is unsuitable for predicting the absolute values of the recording.  Surface EMG model to investigate the effect of pennation angle The model has investigated sEMG for different muscle fibre types, pennation angle and the results show that there is no significant difference in the shape of SFAP because of change in the pennation angle in the range of 0 to 20˚. However, the changes were observed in amplitude, frequency, and chaotic nature of sEMG signals. This reveals that parallel model of the TA can be used to approximate the TA, and this eliminates the need of computationally intensive volume conductor model with angle of muscle fibres to represent the muscle.