Sonomyographic responses during voluntary isometric ramp contraction of the human rectus femoris muscle

Abstract

This paper aims to investigate the relationship between torque and muscle morphological change, which is derived from ultrasound image sequence and termed as sonomyography (SMG), during isometric ramp contraction of the rectus femoris (RF) muscle, and to further compare SMG with the electromyography (EMG) and mechanomyography (MMG), which represent the electrical and mechanical activities of the muscle. Nine subjects performed isometric ramp contraction of knee up to 90% of the maximal voluntary contraction (MVC) at speeds of 45, 22.5 and 15% MVC/s, and EMG, MMG and ultrasonography were simultaneously recorded from the RF muscle. Cross-sectional area, which was referred to as SMG, was automatically extracted from continuously captured ultrasound images using a newly developed image tracking algorithm. Polynomial regression analyses were applied to fit the EMG/MMG/SMG-to-torque relationships, and the regression coefficients of EMG, MMG, and SMG were compared. Moreover, the effect of contraction speed on SMG/EMG/MMG-to-torque relationships was tested by pair-wise comparisons of the mean relationship curves at different speeds for EMG, MMG and SMG. The results show that continuous SMG could provide important morphological parameters of continuous muscle contraction. Compared with EMG and MMG, SMG exhibits different changing patterns with the increase of torque during voluntary isometric ramp contraction, and it is less influenced by the contraction speed.

Introduction

The relationship between muscle activities and generated force during muscle contraction is an important and challenging topic in neuromuscular research fields because it can provide insight into the muscle contraction mechanism (Beck et al. 2005; Disselhorst-Klug et al. 2009). Various recording techniques have been developed to measure muscle activities from different aspects and some commonly used techniques include electromyography (EMG) (Basmajian and DeLuca 1985), mechanomyography (MMG) (Orizio 1993), and ultrasonography (Narici et al. 1996). EMG is composed of electrical contributions made by the active motor units (MUs) during muscle contraction. Meanwhile, MMG is a recording of the low-frequency lateral oscillations of active MUs and is considered to be the mechanical counterpart of the MU electrical activity as measured by EMG (Orizio 1993). Both EMG-to-force and MMG-to-force relationships have been studied in numerous reports (Akataki et al. 2001; Orizio et al. 2010; Watanabe and Akima 2009; Zhang et al. 2010). Moreover, it has been demonstrated that a combination of EMG and MMG could be used to provide complimentary electrical and mechanical information on muscle contraction (Beck et al. 2004; Youn and Kim 2010).

Muscle force is mainly determined by not only the number of recruited MUs and their firing rate, but also their size (e.g. cross-sectional area). Therefore, the morphological change of muscle is an alternative measurement of muscle activity. Ultrasonography has been effectively employed to estimate the morphological changes of muscles in both static and dynamic conditions (Fukunaga et al. 2001; Narici et al. 1996). Different muscle morphological parameters, including thickness (Shi et al. 2008), fascicle length (Oda et al. 2007), pennation angle (Fukunaga et al. 2001), and cross-sectional area (CSA) (Delaney et al. 2010; Narici et al. 1996), could be extracted from ultrasonography to evaluate the muscle function and activity. We have recently proposed the sonomyography (SMG), which denotes the real-time muscle morphological change detected by ultrasonography, for assessment of isometric muscle contraction (Shi et al. 2008) and for prosthetic control (Chen et al. 2010). The findings of the previous studies demonstrated that ultrasonography could provide a noninvasive method of recording activities from deep muscles without the complication of the cross-talk from adjacent muscles suffered in the EMG analysis (Hodges et al. 2003). In particular, ultrasonography could offer unique insight into the mechanical behavior of fibers and interaction mechanisms between fibers and tendons (Fukunaga et al. 2001).

Among the muscle morphological parameters, CSA represents the maximal number of cross-bridges that can be activated in parallel during contraction; therefore, it is directly related to the force-generating capacity of a given muscle (Lieber and Friden 2000). However, most studies focused the static relationship between CSA and maximal force (Aagaard et al. 2001; Suetta et al. 2004); and the knowledge about continuous CSA changes during dynamic contraction is still seriously limited. It was hypothesized that since EMG and MMG reflect the muscle activities form the electrical and mechanical aspects, while CSA reflects the anatomical information of muscles, CSA would exhibit different patterns from EMG and MMG regarding continuous muscle behaviors. More recently, we have developed a system to successfully collect EMG, MMG, SMG and torque signals in a simultaneous manner from the rectus femoris (RF) muscle during isometric contraction (Guo et al. 2010). Our purpose was to identify a dynamic SMG pattern reflecting continuous morphological changes of muscle and compare it with EMG/MMG patterns. Unlike EMG and MMG which can easily produce instantaneous parameters (i.e. root mean square, RMS) to represent continuous electrical and mechanical activities of muscles, few techniques have been developed to automatically extract SMG parameters (i.e. CSA) from a large number of continuously recorded ultrasound images (such as 25 frame/s). For instance, in our previous study (Guo et al. 2010), only a few separate contraction conditions were analyzed in the ramp contraction by manually extracting CSA from a very limited number of ultrasound images. Apparently, such “sparse” SMG parameters can hardly represent the continuous muscle behaviors, and it is highly desirable to overcome the technical difficulty of automatically extracting continuous SMG patterns prior to the comparisons between SMG and EMG/MMG. Addressing these questions would considerably advance the developments and applications of the SMG and the simultaneous EMG/MMG/SMG techniques, and, consequently, contribute to reveal the underlying physiological mechanisms of muscle contraction.

In this paper, we employed a system for simultaneously collecting EMG, MMG and SMG from the RF muscle during isometric ramp contraction at different contraction speeds. A novel image tracking algorithm was firstly developed to extract CSA automatically from ultrasound image sequences, making SMG a practical technique to characterize continuous morphological changes of muscles. The aim of the study was to examine the continuous SMG-to-torque relationship during isometric ramp contraction, and compared it with EMG-to-torque and MMG-to-torque relationships.

Methodology

Subjects and experimental procedure

Nine healthy adults, six males and three females, volunteered to participate in this study (mean ± SEM, age = 31.2 ± 1.8 years; body weight = 63.7 ± 5.3 kg; height = 168.1 ± 3.1 cm). No participant had a history of neuromuscular disorders, and all were aware of experimental purposes and procedures. The human subject ethical approval was obtained from the relevant committee in the Hong Kong Polytechnic University before carrying out the experiment. The subjects were briefed about the procedure of experiment and written consents were collected from all subjects prior to the experiment.

During measurement, a subject was seated on a test bench of an isokinetic dynamometer (Humac Norm Testing and Rehabilitation System, Computer Sports Medicine, Inc., MA, USA) with straps secured around the trunk, and the lateral condyle of the femur was aligned with the rotation axis of the dynamometer. All isometric torque assessments were performed with both hip joint angle and knee joint flexed at 90°.

Before the test began, an ultrasound image was collected from the subject on the RF muscle under the relaxed condition as reference. Prior to the isometric testing sessions, the subject was required to perform two maximal, 6-s isometric knee extensions to determine the maximal voluntary contraction (MVC), which was defined as the highest value of torque recorded during the entire isometric contraction. Then the subject performed several warm-up trials to get familiar with the experiment protocol. After a rest of 5 min, the subject was instructed to perform ramp contractions which produced torques increasing from 0 to 90% of his/her MVC. During each contraction, a torque production template, serving as a target, and the subject’s real-time torque signal were displayed together on a computer screen placed in front of him/her. The subject was required to adjust his/her torque production to track the target torque which consisted of a linear progression from 0 to 90% MVC. Three kinds of ramp time durations were tested: 2, 4 and 6 s; the corresponding speeds were 45, 22.5, and 15% MVC/s. For each speed, three trials were repeated with a rest of 5 min between adjacent trials.

Data acquisition and synchronization

A real-time B-mode ultrasonic scanner (EUB-8500, Hitachi Medical Corporation, Tokyo, Japan) with an electronic linear array probe (L53L, Hitachi Medical Corporation, Tokyo, Japan) was used to obtain transactional ultrasound images of the RF muscle. The ultrasound probe was fixed by a custom-designed multi-degree adjustable bracket. The long axis of the probe was arranged perpendicularly to the long axis of the thigh on its superior aspect, 40% distally from the knee (measured from the anterior superior iliac spine to the superior patellar border). The video output of the scanner was digitized by a video capture card (NI PCI-1411, National Instruments Corporation, Austin, TX, USA) with a frame rate of 25 Hz and resolution of 0.15 mm; the instant when every frame was captured was recorded as timestamp for off-line signal processing.

Two surface bipolar Ag–AgCl EMG electrodes (Axon System, Inc., NY, USA; 2 cm diameter; center-to-center inter-electrode distance of 3.2 cm) were placed on the RF muscle belly parallel to the long axis of the muscle on both sides of ultrasound probe, and a reference EMG electrode was placed near the kneecap. The MMG signal was detected by an accelerometer (EGAS-FS-10-/V05, Measurement Specialties, Inc., France) placed near the ultrasound probe and fixed to the skin using double adhesive tape (Fig. 1). The surface EMG and MMG signals were amplified by a custom-designed amplifier with a gain of 2,000, filtered separately by 10–400 and 5–100 Hz band-pass analog filters, and digitized by a 12-bit data acquisition card (NI-DAQ 6024E, National Instruments Corporation, Austin, TX, USA) with a sampling rate of 1 kHz. The isometric torque output from the dynamometer was also sampled by the NI-DAQ card in synchronization with the ultrasound image capture. Ultrasound images, surface EMG, MMG and torque signals were simultaneously collected and stored by custom-developed software for Ultrasonic Measurement of Motion and Elasticity (UMME, http://www.tups.org) developed in VC++ platform for further analysis.

Fig. 1

Experimental setup for collecting EMG, MMG and ultrasound images from the subject’s right rectus femoris (RF) muscle during isometric ramp contraction. The ultrasound probe was aligned perpendicularly to the RF muscle belly using a multi-degree adjustable bracket. The two surface EMG electrodes were located at the two sides of the ultrasound transducer and in parallel with the muscle fibers. The MMG sensor was placed as near as possible to the distal EMG electrode, and they were in parallel with the ultrasound probe

SMG processing

In this study, a novel image processing method, named constrained mutual information (MI)-based free-form deformation (C-MI-FFD) tracking, was developed to automatically extract the continuous CSA changes in the ultrasound image sequence. In the proposed C-MI-FFD method, we aim to determine the transformation function which describes the deformation between two successive images by minimizing a MI-based objective function. For accurate matching, the transformation in the C-MI-FFD method is carried out in two steps: (1) the global transformation and (2) the local transformation. In global transformation, the transformation parameters are determined by matching the two images globally so as to model their relative scale, translation and rotation. In the local transformation, the local deformation is defined by a 2D spline function and the transformation parameters are refined as displacement values at a regular grid to interpolate the spline function. To further improve the tracking performance, we incorporate three image structure-derived constraints, including smoothing constraint, feature point constraint, and edge constraint, into the MI objective function. The derivation and details of the C-MI-FFD method can be found in the “Appendix”.

For each trial, the first image in the sequence was selected as reference and the boundary of the RF muscle was outlined with smooth lines by the investigator using ImageJ software (ImageJ, National Institutes of Health, USA). Then the C-MI-FFD method was applied to track the CSA boundaries in the subsequent images. The CSA value was then normalized as the percentage decrease ∆CSA, with respect to the value at relaxed condition.

EMG and MMG processing

The EMG and MMG signals were segmented as 256-ms epochs. The center of each epoch was aligned in time with the corresponding ultrasound image according to the timestamp, so that the epochs were synchronized with the image sequence in time domain. The root mean square (RMS) values of EMG and MMG (RMS_EMG and RMS_MMG) were calculated for each epoch and expressed as a percentage of their maximal values at 90% MVC.

Statistical analyses

To evaluate the reliability of the automatic image tracking algorithm, a total of 100 randomly selected ultrasound images were measured twice by the image processing algorithm and the investigator. The intra-class correlation coefficient (ICC) based on one-way random model and standard error of the measurement (SEM) were calculated (Weir 2005).

There are totally 9 (subjects) × 3 (speeds) × 3 (trials) sets of data. For each subject at each speed, the mean values of RMS_EMG, RMS_MMG, ∆CSA, and torque were calculated by averaging across the three trials. Then the inter-individual means among the nine subjects were computed at each speed. The relationships for RMS_EMG, RMS_MMG, and ∆CSA versus torque were investigated for each subject and the mean relationships over all the subjects were computed for each contraction speed. The patterns of these relationships were examined using polynomial (linear, quadratic) regressions:

$$ Y = a_{0} + a_{1} X\quad ({\text{linear model)}} $$

and

$$ Y = a_{0} + a_{1} X + a_{2} X^{2} \quad ({\text{quadratic}}\;{\text{model}}) $$

where X is torque, Y is RMS_EMG, RMS_MMG, or ∆CSA, a ₀, a ₁ and a ₂ are regression coefficients.

The best-fit model was selected for each subject according to the coefficients of determination R ² based on the method described in (Beck et al. 2004). One-sample t test was applied to test whether the quadratic coefficient a ₂ was positive for EMG and MMG and negative for SMG (for a linear model, a ₂ was set as 0). Further, to check whether the RMS_EMG/RMS_MMG/∆CSA to torque relationships are influenced by contraction speed, a hypothesis testing based on adaptive Neyman test (Fan and Lin 1998) was applied to compare the relationship curves of RMS_EMG, RMS_MMG, and ∆CSA at different contraction speeds. All the data were analyzed using SPSS (SPSS Inc. Chicago, IL, USA). Statistical significance was set at the 0.05 probability level.

Results

EMG and MMG amplitudes

The torque, ∆CSA, EMG, and MMG signals during one representative trial (the contraction speed is 22.5% MVC/s) are shown in Fig. 2. The relationships for EMG/MMG versus torque for this trial were best-fitted with quadratic models for both RMS_EMG (R ² = 0.954) and RMS_MMG (R ² = 0.921) (Fig. 3a, b). For all the subjects, the inter-individual means of the EMG and MMG amplitudes at three contraction speeds are plotted as functions of torque in Fig. 4. The composite relationships were investigated using polynomial regression and the best-fit models are presented in Fig. 4. All the six composite patterns were best-fitted with quadratic regressions. The quadratic coefficients a ₂ of composite EMG-to-torque and MMG-to-torque curves had positive values at all three contraction speeds (EMG: 0.216, 0.469, 0.594; MMG: 0.435, 0.421, 0.668; for 45, 22.5, 15% MVC/s, respectively). One-sample right-tail t test on the coefficients a ₂ showed that a ₂ of EMG and MMG was significantly greater than 0 at all the contraction speeds, except at 15% MVC/s of MMG (Table 1).

Fig. 2

The time courses of a torque, b ∆CSA, c EMG, and d MMG from the RF muscle during a representative trial. The torque signal was overlaid onto the ramp template as it appeared for the subject during the trial. The EMG and MMG values were the captured voltages from the acquisition card

Fig. 3

The relationships between torque and a RMS_EMG, b RMS_MMG, and c ∆CSA during one trial. The best-fit polynomial regression curves are shown as solid lines

Fig. 4

The inter-individual mean RMS_EMG values at different speeds: a 45% MVC/s, b 22.5% MVC/s, c 15% MVC/s; and RMS_MMG values at different speeds: d 45% MVC/s, e 22.5% MVC/s, f 15% MVC/s. RMS values are shown as mean (circles) ± SEM (black horizontal bars). The best-fit polynomial regression curves are shown as solid lines

Table 1 P values of the one-sample one-tailed t test on the quadratic coefficients a ₂ of quadratic regressions of EMG/MMG/SMG-to-torque relationships. The test for EMG and MMG is that the quadratic coefficient a ₂ is positive, and the test for SMG is that the quadratic coefficient a ₂ is negative

∆CSA of SMG

Four images of a typical trial, including the first image, the first image with manually drawn boundary, the image at 50% MVC, and the image at 50% MVC with automatically tracked boundary, are shown in Fig. 5. The ICC for CSA measurements of the RF muscle in this study was 0.987 (P < 0.0001), and the SEM was 0.15 cm². The average CSA for the nine subjects was 6.15 ± 0.48 cm² (mean ± SEM) under the relaxed condition. The example of time course of ∆CSA is shown in Fig. 2b, the relationship between ∆CSA and torque for this trial was best-fitted with a quadratic model (R ² = 0.985) (Fig. 3c). The composite curves are shown in Fig. 6 as means of all subjects at each contraction speed. The quadratic coefficients a ₂ of composite SMG-to-torque curves had negative values at all three contraction speeds (−4.227, −1.477, −9.171 for 45, 22.5, 15% MVC/s, respectively). However, one-sample left-tail t test on the coefficients a ₂ showed that a ₂ of SMG was not significantly less than 0 at any contraction speed (Table 1).

Fig. 5

Ultrasound images of the RF muscle in one trial. a The first image in the image sequence, where the neighboring quadriceps muscles [vastus lateralis (VL), vastus intermedius (VI), and vastus medialis (VM)] are indicated. b The first image in the image sequence, with manually outlined boundary as reference for further image processing. c The image at 50% MVC. d The image at 50% MVC, with automatically outlined boundary by the proposed C-MI-FFD algorithm

Fig. 6

The inter-individual mean ∆CSA values at different speeds: a 45% MVC/s, b 22.5% MVC/s, and c 15% MVC/s. ∆CSA values are mean (circles) ± SEM (black horizontal bars). The best-fit polynomial regression curves are shown as solid lines

Contraction speed

We next compared EMG, MMG, or SMG-to-torque relationships between different contraction speeds. The mean curves of RMS_MMG, RMS_MMG or ∆CSA at three speeds were calculated across all subjects (Fig. 7). Table 2 shows the results of adaptive Neyman test on pairwise comparisons of the RMS_EMG, RMS_MMG or ∆CSA curves among three speeds (45 vs. 22.5%, 45 vs. 15%, 22.5 vs. 15%). The significance level of the adaptive Neyman test is 0.05. For RMS_EMG, the 22.5% MVC/s curve had no significant difference with the 45 and 15% MVC/s curves (P = 0.60 and 0.95, respectively), while a significant difference was found between 45 and 15% MVC/s curves (P < 5e−5). For RMS_MMG, the 22.5% MVC/s curve had no significant difference with the 45 and 15% MVC/s curves (P = 0.96 and 1.00, respectively), while a significant difference was identified between 45 and 15% MVC/s curves (P = 0.046). For ∆CSA, no significant difference was observed between any pair of curves.

Fig. 7

The inter-individual mean curves at different contraction speeds for a RMS_EMG, b RMS_MMG, and c ∆CSA

Table 2 P values of the adaptive Neyman test (Fan and Lin 1998) for pairwise comparisons of the EMG, MMG, or SMG curves among different speeds

Discussion

In this study, we collected EMG, MMG and SMG simultaneously from the RF muscle during isometric ramp contraction at different speeds. The purpose is to examine the continuous SMG-to-torque relationship during isometric ramp contraction, and to compare it with those of EMG and MMG. To our knowledge, this is the first attempt to investigate the continuous muscle behavior with simultaneously measures of EMG, MMG and SMG.

Previous studies showed that muscle functions, such as the amplitude of force production, are closely related to architectural characteristics (e.g. CSA) (Lieber 2002). Many researches have studied the static relationship between CSA and muscle strength under different conditions (Aagaard et al. 2001; Suetta et al. 2004). Narici et al. (1996) applied ultrasonography to measure human gastrocnemius architecture and investigated for the first time the in vivo change of CSA with changing joint angle at rest and during graded isometric contraction. The study concluded that human muscle architecture was significantly affected by both changes of joint angle at rest and isometric contraction intensity. Recent studies expanded this method by simultaneously measuring ultrasonography, EMG, and MMG (Delaney et al. 2010; Guo et al. 2010). However, all the previous ultrasonography studies measured the CSA values at separate contraction conditions by manually outlining the muscle boundary, and the continuous relationship between CSA and force was not clear yet. In this study, we developed a C-MI-FFD image tracking algorithm to automatically extract the continuous CSA changes of RF muscle during isometric ramp contraction, which makes it possible to identify muscle morphology in a high temporal resolution. The reliability of this method was investigated by ICC test in comparison with the manual operation. The excellent results showed that the C-MI-FFD algorithm was highly reliable and its results are consistent with the manually obtained results.

EMG and MMG provide assessments of the electrical and mechanical activities of MUs. In isometric ramp and step contractions, EMG amplitude is typically characterized by linear or curvilinear increase with the increase of force (Disselhorst-Klug et al. 2009). The results of the present study showed quadratic curvilinear increase in EMG amplitude across the force spectrum at each contraction speed (Fig. 4a, b, c). These findings were consistent with previous studies (Ryan et al. 2008b; Watanabe and Akima 2009) which also reported curvilinear isometric torque-related increase in EMG amplitude for the RF muscle. For MMG amplitude, it has been proposed to rely mainly on MU recruitment for increasing isometric torque at low force region, and above which, modulation of MU firing rates becomes more important for further increasing torque (Akataki et al. 2001; Ryan et al. 2008b). Therefore, the MMG amplitude–force relationship may be broadly divided into two contiguous force regions: curvilinear increase in the low force region and decrease or plateau in the high force region (Akataki et al. 2001). During the fastest contraction in the present study (45% MVC/s, Fig. 4d), the RMS_MMG showed a first increase up to about 80% MVC, and thereafter it decreased till the end of the contraction. The deflection point was similar to that reported in the previous study (Ryan et al. 2008b). Such observation was also presented at the other speeds of contraction (Fig. 4e, f), although it was less evident.

In the present study, the SMG of the RF muscle was found to decrease with torque. It was hypothesized that the curvilinear change pattern of CSA was different from those of EMG and MMG. As shown in Fig. 3, the quadratic coefficients a ₂ of the EMG and MMG regression curves were positive and thus EMG/MMG-to-torque relationships were convex downward, implying the change rate of EMG and MMG had an increasing manner. More precisely, EMG and MMG amplitudes increased relatively slowly in the low torque region and fast in the high torque region. In contrast, the quadratic coefficients a ₂ of the SMG regression curve was negative, representing a concave SMG-to-torque relationship. Therefore, SMG amplitude changed fast in the low torque region and slowly in the high torque region. These findings were also supported by the regression analysis and the one-sample t test. As shown in Table 1, the coefficient a ₂ of EMG or MMG was significantly positive at all the contraction speeds, except at 15% MVC/s of MMG. Although significance was not found for SMG coefficient under all conditions, the quadratic coefficients a ₂ of composite SMG-to-torque curves had negative values at all three contraction speeds (Fig. 6: −4.227, −1.477, −9.171 for 45, 22.5, 15% MVC/s, respectively), strongly suggested that the coefficients a ₂ of SMG trended to be negative.

The relationship between CSA and torque observed on RF muscle in this study is similar to that reported in the recent study (Delaney et al. 2010). In both studies, the CSA values showed an initial small increase followed by a curvilinear decrease with the increase of torque (Fig. 6). Many factors during muscle contraction may contribute to this phenomenon. (1) During isometric contraction the series elastic component is under tension and is shortened a finite amount. Therefore, the transverse diameter of individual muscle fiber may bulge and contribute to the overall CSA increase. This could explain the small CSA increment in initial phase of the contraction. (2) As a consequence of the intra-muscular pressure increase during contraction, the intra-muscular fluid may be squeezed outside the space between fibers and therefore decrease dramatically the muscle dimension. (3) As indicated in Fig. 5, the RF muscle is surrounded by three other muscles, vastus intermedius, vastus lateralis, and vastus medialis, which also contribute to the knee extension. The space constraint between these muscles may distort the shapes of RF. It is of interest in the future to study the specific effect of this factor on a muscle without surrounding muscles. On the basis of the factors mentioned above, SMG has the potential to contribute to the study of the muscle contraction by providing unique morphological parameters. In particular, it could help to investigate the relationship between the changes of the output force and of the muscle morphology resulting from the fiber shortening and interaction between fibers.

According to (Akataki et al. 2001), the continuous MMG/force relationship may provide more detailed information on the MU activation strategy. However, it is technically difficult to obtain the continuous morphological change during muscle contraction. To our knowledge, this study is the first attempt to investigate the continuous CSA behavior during isometric ramp contraction. The findings on CSA in the present study are also in accordance with our previous study (Guo et al. 2010), which measured isometric ramp contraction but calculated the CSA values at separate contraction levels. The present study showed great improvement by automatically calculating CSA values and therefore providing continuous CSA change during isometric contraction. Comparing the relationship curves of both studies, the findings of present study demonstrated more information: (1) the initial increase of CSA cannot be observed in previous study because it occurred below the first contraction level; and (2) the reliability of the polynomial regression was improved since there were much more sample points in the present study.

The joint angle has a significant effect on the muscle dimension. The RF muscle is a bi-articular muscle which flexes the thigh at the hip joint and extends the leg at the knee joint. Therefore, the RF dimension will be affected by both hip and knee joint angles which were fixed during isometric knee extension. In our experiment setup, both angles were set at 90°. Although this setup may not be optimal for knee extension, it relates to a natural sitting posture and is adopted by other research on RF muscle during knee extension (Delaney et al. 2010). Some studies have used other combinations of the hip/knee joint angles, such as 125°/90° (Watanabe and Akima 2009), 90°/120° (Pincivero et al. 2003; Ryan et al. 2008a). Since the purpose of our study is not to investigate the influence of joint angles, we just chose these specific angles (90°) and investigated the relationships under this condition. In future, it will be of great value to study the influence of hip angle or initial shortness of RF muscle on the behavior of SMG during knee extension.

Most of the studies on isometric ramp contraction were conducted at a fixed contraction speed. Although previous studies showed that the contraction speed had significant effect on the EMG amplitude and spectral characteristics (Orizio et al. 2010; Sbriccoli et al. 2003), the effect of speed on morphological characteristics of muscles has not been investigated. In the current study, three different contraction speeds, namely 45, 22.5 and 15% MVC/s, were employed in the ramp contraction. As shown in Fig. 7, the relative positions among the three curves of EMG are similar to those of MMG, but different to those of SMG. For EMG and MMG curves, the curve of intermediate speed lies between those with lower and higher speed, and it is more close to the curve of lower speed. These observations were validated with the statistical test results in Table 2 which indicated that EMG/MMG curves of 45 and 15% MVC/s were significantly different; while SMG curves were less influenced by the speed. The different effects of contraction speed on EMG/MMG and SMG may be derived from their different characteristics. According to the previous study (Sbriccoli et al. 2003), MU recruitment strategies are significantly related to the contraction speed. Since EMG/MMG amplitude relies mainly on MU recruitment, it is considerably related to the contraction speed. However, SMG represents the anatomical information of muscle which is directly related to the force (or torque) production. The SMG-to-torque relationship is more direct than that of EMG/MMG, and therefore, is less affected by other factors like contraction speed.

In conclusion, this study demonstrates that the SMG-to-torque relationship is different from those of EMG and MMG during isometric ramp contraction in two aspects. (1) The mean regression curves of SMG-to-torque relationships have negative coefficients a ₂ at all the contraction speeds, while the mean regression curves of EMG/MMG-to-torque relationships have positive coefficients a ₂. (2) The SMG-to-torque relationship is less affected by contraction speeds than the EMG/MMG-to-torque relationships. These results show that continuous SMG could provide important morphological parameters of continuous muscle contraction. Compared with electrical and mechanical characteristics measured by EMG and MMG, SMG-measured morphological characteristics exhibit different changing patterns with the increase of torque during voluntary isometric ramp contraction, and they are less influenced by the contraction speed.

Appendix

The MI-based free-form deformation (MI-FFD) method was originally proposed in (Huang et al. 2006) to do shape registration in 2D space based on the mutual information between two images and it was extended in this paper for object tracking. More precisely, the intensities of two images, say A and B, to be matched are treated as random variables with probability density functions (pdfs), p _A (i _A ) and p _B (i _B ), and joint pdf p _AB (i _A , i _B ). B is then deformed by mean of a transformation function T(B) with parameters to be determined. Ideally, when the two images are registered, the mutual information between A and T(B),

$$ \begin{gathered} {\text{MI}} = I(A;T(B)) = H(p_{A} (i_{A} )) + H(p_{T(B)} (i_{B} )) - H(p_{A,T(B)} (i_{A} ,i_{B} )) \hfill \\ \, = \iint\limits_{I(A),I(B)} {p_{A,T(B)} (i_{A} ,i_{B} )\log \frac{{p_{A,T(B)} (i_{A} ,i_{B} )}}{{p_{A} (i_{A} )p_{T(B)} (i_{B} )}}{\text{d}}i_{A} {\text{d}}i_{B} }, \hfill \\ \end{gathered} $$

(1)

will be maximized, where i _A and i _B are the intensity valuables of A and T(B), H(p _A (i _A )), H(p _T(B)(i _B )), H(p _A,T(B)(i _A , i _B )) are, respectively, the entropies of A, T(B) and their joint entropy. Therefore, by maximizing I(A, T(B)) using the parameters of T(·), the two original images can be registered and the segmented object can be tracked continuously. In practice, the pdfs are approximated by the kernel method, which can be computed directly from the image data.

However, for the highly noisy ultrasound images under study, the conventional MI-FFD method cannot achieve a satisfactory tracking performance. For a more accurate matching, the transformation in the proposed C-MI-FFD method is carried out in two steps: (1) the global transformation and (2) the local transformation. In the first step of global transformation, the parameters are determined by matching the two images globally so as to model their relative scale, translation and rotation. The global deformation function T _G is chosen as an affine transformation T _G  = Ψ[x,y,1]^T, where x and y are the coordinates of the images and Ψ is a 3 × 3 affine transformation matrix in the homogeneous coordinate that can be obtained by maximizing Eq. 1. Let B′ be the transformed image obtained by the affine transformation after the global transformation. Further, B′ will be refined by the subsequent local transformation. In the second step of local transformation, the local deformation T _L (B′) is defined by a 2D spline function. The transformation parameters, which are the displacement values at a regular grid to interpolate the spline function, are, again, determined by maximizing the MI objective function in Eq. 1. The local transformation T _L (B′) is parameterized by the displacement vectors at a uniform grid of control points C, $\user2{P}_{c} (m,n) = [P_{{c,\overline{x} }} (m,n),P_{{c,\overline{y} }} (m,n)]^{T} ,$ which are indexed by m = 1, …, M, n = 1, …, N. If (X, Y) is the resolution of the input image, the spacing of the control points in the $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {x}$ and $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {y}$ directions are $\Updelta_{{\overline{x} }} = X/M$ and $\Updelta_{{\overline{y} }} = Y/N$, respectively. The deformation of any pixel in the image is obtained by spline interpolation of those at the grid points C. Therefore, the deformation of pixel (i, j), $\user2{P}(i,j) = [P_{{\bar{x}}} (i,j),P_{{\bar{y}}} (i,j)]^{T} ,$ where 1 ≤ i ≤ X, 1 ≤ j ≤ Y can be written as:

$$ \user2{P}(i,j) = \sum\limits_{\mu = 0}^{3} {\sum\limits_{\gamma = 0}^{3} {\beta_{\mu } (u)\beta_{\gamma } (v)\user2{P}_{c} (m + \mu ,n + \gamma )} } , $$

(2)

where β _μ (u) and β _γ (v) are, respectively, the cubic B-spline function with $u = i/\Updelta_{{\overline{x}}}-\left\lfloor {i/\Updelta_{{\overline{x}}}} \right\rfloor$ and $v = j/\Updelta_{{\overline{y}}}-\left\lfloor {j/\Updelta_{{\overline{y}}}} \right\rfloor,$ and {(m + μ, n + γ)|(μ, γ) ∈  [0, 3]} are the neighboring control points of (i, j). Then, the local transformation of B′ can be obtained as:

$$ T_{L} (B'(i,j)) = B'(i + P_{{\overline{x} }} (i,j),j + P_{{\overline{y} }} (i,j)). $$

(3)

By substituting (2) into (1), one gets the local matching objective function E _local = −MI to be minimized.

The global transformation only needs to estimate nine model parameters in the 3 × 3 matrix Ψ, but the local transformation has to estimate a large number of parameters (i.e., each pixel has two parameters to be estimated in two directions). In order to reduce the variance of the local transformation parameters, some prior constraints, including the smoothing constraint, the feature point constraint, and the edge constraint, are added into the local matching objective function. In the proposed C-MI-FFD method, we incorporate several constraints, which are derived from structure information of the images, into the MI objective function so as to improve the tracking performance.

1. Smoothing constraint. A popular smoothing constraint is the L ₂ norm of the displacement $E_{\text{smooth}} = \sum\nolimits_{(m,n)} {( {\| {\frac{{\partial P_{c} (m,n)}}{\partial m}}\|^{2} + \| {\frac{{\partial P_{c} (m,n)}}{\partial n}}\|^{2} })} .$ The smoothing term can effectively limit the displacement and thus decrease the estimation variance.

2. Feature point constraint. If the pair of images being registered does have distinct geometric features as correspondences, incorporating this geometric feature information can greatly improve accuracy and efficiency. In our algorithm, the scale invariant features (Lowe 2004) can be used as a feature constraint. When feature points and correspondences are available, these constraints can be conveniently integrated into our registration framework. Assuming that the total number of features is K, and for each feature, there is a pair of corresponding points, I ^(T) _k , k = 1, …, K, on the target image and I ^(S′) _k  = T _L (I ^(S) _k ), k = 1, …, K, on the locally transformed source image, then the following term can be incorporated as feature constraint:

$$ E_{\text{feature}} = \sum\limits_{k = 1}^{K} {D(\user2{I}_{k}^{(S')} - \user2{I}_{k}^{(T)} )} , $$

(4)

where D(I ^(S′) _k  − I ^(T) _k ) is the Euclidean distance between I ^(T) _k and I ^(S′) _k .

3. Edge constraint. The edge term is expressed as:

$$ E_{\text{edge}} = \sum\limits_{i = 1}^{M} {\sum\limits_{j = 1}^{N} {[A_{\text{edge}} (i,j) - B_{\text{edge}}^{(T)} (i,j))]} } , $$

(5)

where A _edge and B ^(T)_edge are, respectively, the edge maps of A and T(B), and M and N are the height and width of the images. A color tensor-based edge detection (van de Weijer et al. 2006) was used in our study to extract the image edges A _edge and B ^(T)_edge .

With all these constraint terms, the final local objective function to be minimized reads:

$$ E_{\text{local}} = - {\text{MI}} + E_{\text{smooth}} + E_{\text{feature}} + E_{\text{edge}} . $$

(6)

To solve the nonlinear optimization problem, an L-BFGS algorithm (Byrd et al. 1995) is used in the proposed C-MI-FFD method. An advantage of the L-BFGS algorithm is that the explicit evaluation of the Hessian matrix is not required, so it can be recursively estimated. Moreover, the L-BFGS algorithm is much faster than the conventional level set method in (Huang et al. 2006) for solving the nonlinear optimization problem.