Skeletal muscles not only act as force transducers allowing for macroscopic displacements but also play an important role for the stabilization of vertebra and joints. Neuromuscular disorders such as those due to peripheral nerve injury can thus severely affect the overall systemic function and health condition. Fine-wire electromyography (EMG), the standard diagnostic technique in the clinic, allows to monitor the electrical activity of single motor units. However, the pain associated with electrode introduction and the sensitivity to spasm-induced artifacts prevent the method from being used for chronical monitoring of patients, e.g., during pharmaceutical treatment or rehabilitation. Tissue Doppler imaging, on the other hand, is a noninvasive method sensitive to tissue velocity, but it suffers from susceptibility to artifacts from global muscular displacements and relatively low sensitivity.1 Early ex vivo experiments on frog muscle have suggested that quasi-elastic light scattering (QELS) is a sensitive probe for muscle contraction.2 Attempts at using the laser Doppler method, the frequency-domain analog of QELS, for in vivo measurements have been shown to differentiate denervated muscles in the human hand.3 Unfortunately, the small Doppler widths of typically several that are dictated by the short source–receiver distance of a few mm do not allow to resolve the contraction process.4
Here, we show that the contraction of skeletal muscles can be measured noninvasively and with millisecond temporal resolution using near-infrared diffusing-wave spectroscopy (DWS). DWS (also called diffuse correlation spectroscopy, or DCS) is the extension of QELS to the regime of optically dense, opaque media.5, 6 Being highly sensitive to subwavelength displacements of scatterers within tissue, e.g., from erythrocytes moving within the vasculature, DWS has recently been introduced as a method for the noninvasive monitoring of blood flow in tumors7, 8, 9 and for the noninvasive measurement of functional brain activity.10, 11, 12, 13 DWS signals correlate strongly with tissue blood flow rates, as was demonstrated for skeletal muscle microvasculature.14 Its high sensitivity to tissue perfusion changes and the possibility for bedside measurements makes DWS an attractive method for monitoring blood flow, in particular in situations where conventional methods such as PET or Xe CT cannot be used, such as in neurointensive care, neonatology, or rehabilitation. Indeed, a recent study demonstrated that skeletal muscle perfusion and oxygenation can be monitored simultaneously with a compact hybrid instrument combining DWS and near-infrared spectroscopy.15 While it has recently been recognized that characterization of muscular blood flow by DWS during exercise does require correction for artifacts arising from muscle fiber motion,16 the potential of DWS for the explicit investigation of muscle contraction has so far not been explored.
Our results show that DWS signals from contracting skeletal muscles may be dominated by tissue shearing due to sarcomere shortening and advection of erythrocytes within the microvasculature. Using a multispeckle detection setup,17 we measure DWS autocorrelation functions with an integration time of , which allows us to resolve the initial stages of the contraction and to differentiate between fast- and slow-twitch muscles. The analysis of DWS autocorrelation functions using tissue optical properties measured with time-resolved reflectance spectroscopy yields quantitative values for the evolution of the muscle strain .
All experiments were performed on a single, healthy, male volunteer (age ). Muscle contraction was induced by applying bipolar voltage pulses ( width per phase) generated by a function generator, at a frequency (see Fig. 1 ), via two electrode pads placed symmetrically about the muscle belly at a distance of about . Measurements on the M. biceps brachii were carried out with the volunteer’s arm placed horizontally on a table at the height of the heart. Stimulation blocks lasted typically . Forces were measured with a piezoelectric transducer (Kistler 9321B) placed on the wrist. The force sensor was preloaded in order to guarantee good mechanical coupling during both contraction and relaxation. Recording the stimulus trigger signal allowed to synchronize the stimulation with the optical measurements. Optical measurements were started about after the start of a stimulation block, allowing the muscles to reach a physiological steady state. Heart beat rate was monitored with a pulse oxymeter at the index finger tip.
Time-Resolved Reflectance Spectroscopy
In order to quantify the absorption and reduced scattering coefficient, and , during a contraction–relaxation cycle, we used a dual-wavelength time-resolved reflectance spectroscopy system consisting of two diode lasers operating with pulse width of at a repetition rate of and wavelengths and for illumination, and a time-correlated single-photon counting system for the detection of diffusely transmitted photons (for details, see Ref. 18). The light from each of the two lasers is coupled into a multimode fiber that is connected to an optical switch; output channel 1 injects light into the muscle, while channel 2 is put over a solid phantom (made of epoxy resin, particles as scatterers, and black toner powder as absorber) to continuously monitor system performances and possible time and power drifts. Time-multiplexing is accomplished by the different lengths of the multimode fibers. Diffusely transmitted light is collected by a fiber optic bundle placed at a distance from the source and detected by photomultiplier tubes whose output signals are processed by two time-correlated single-photon counting boards. The optical sensor was placed along the direction of the muscle fibers determined by the vector between Fossa cubitalis and the shoulder. Using an integration time of ( for each of the two wavelengths), we obtained time-resolved reflectances with about events per data point. The values of the optical properties as a function of time after the stimulus, and , were determined by fitting the solution of the diffusion equation19, 20 convoluted with the experimental instrument response function, to the measured time-resolved reflectance. In Eq. 1, is the speed of light in the medium with refractive index , is the photon diffusion coefficient, is the depth of diffuse source, and is the extrapolation length that accounts for reflections at the tissue surface. 2 To reduce dispersion of the fitted absorption coefficient values, we use the method described by Nomura known as the modified Lambert-Beer law.21 First, for each wavelength, a reference time-resolved reflectance curve is derived by averaging the data within an initial resting period of typically . Fitting Eq. 1 to yields the reference absorption value . For each further data point, the variation with respect to the reference value is computed using , yielding the value . The combined thickness of the skin and adipose tissue layers covering the muscle, as measured with a standard caliper, was about .
Diffusing-wave spectroscopy measurements were carried out using a fiber-based multispeckle setup.17 Light from a diode laser operating at in a single longitudinal mode was coupled into a multimode optical fiber whose other end was placed on the skin, resulting in an illumination spot of about in diameter. A bundle of 32 single-mode optical fibers ( SMC-780; nominal cut-off wavelength , mode field diameter , numerical aperture 0.13) placed at a distance from the source, was used to collect light that had been diffusely transmitted through the tissue. In order to perform experiments in the presence of ambient light, we equipped the fiber receiver with a bandpass filter (Semrock FF01-800/12-25) with spectral width and 93% transmission at the center wavelength . The light from each fiber was guided to an avalanche photodiode (APD; Perkin-Elmer SPCM-AQ4C). The normalized photon count (intensity) autocorrelation function as a function of lag time was computed by a 32-channel multitau autocorrelator ( www.correlator.com) from the output signal of the APDs. Averaging the over the fibers within the bundle followed by a stimulus-synchronized average over typically 200 trials allowed us to reduce the integration time per data point to at an average photon count rate of about . Field autocorrelation functions were calculated using the Siegert relation , with a constant coherence factor accounting for the two orthogonal polarizations guided by the receiver. In order to quantify the contraction-induced changes of the DWS signal in a model-independent way, we computed the average decay time (Ref. 22).
To separate tissue dynamics from changes in optical properties that also affect , we computed the mean square phase fluctuation for a single scattering event as a function of lag time and latency time from the field autocorrelation function with a numerical zero-finding routine, using Eqs. 6, 7 and the experimentally determined values of and . Assuming that scatterer displacements arise from a combination of diffusion and shearing, we fitted the modelaccounts for reductions of the coherence factor from its theoretical value due to detector afterpulsing and dead-time as well as slight normalization errors in due to the restricted integration time. [From the inversion of simulated intensity autocorrelation functions, we estimate that the errors in due to a reduced intercept are at most 4.6%.]
Force Measurements and Physiology
The forces measured on the wrist during contraction of the M. biceps brachii show, for a stimulation amplitude between and , a rapid increase to a peak and a slower decay during relaxation (see Fig. 2 ). The minimal stimulation voltage at which an increase in the force after stimulation can be observed varies between about and , depending on the mechanical contact between force transducer and the wrist. The latencies for the peak force vary between about for low and for the highest stimulation voltages, in good agreement with the literature value (Ref. 23. The slight undershoot of the force during relaxation at high stimulation amplitude is due to slight rotations of the wrist. Within experimental error, the heart beat rate measured with the pulse oxymeter did not show any difference between baseline and stimulation periods.
Tissue Optical Properties
Figure 3 shows that representative single-trial data for are well described by the theory for the homogeneous semi-infinite geometry. Upon electrical stimulation, the reduced scattering coefficient, , increases from its baseline value at about to a maximum whose height increases with stimulation voltage (see Fig. 4 ). For the lowest voltage at which a change in can be observed , reaches its maximum at , about 17% above the baseline value, within about . The peak value of increases at higher stimulation amplitudes; the latency of the peak is, to within the temporal resolution of , independent of stimulation amplitude.
In contrast to the behavior of , the absorption coefficient, , is observed to decrease rapidly from its baseline value of about to a minimum whose depth increases with the stimulation voltage (see Fig. 4). For the lowest stimulation voltage at which a change in is observed , the minimum value is by about 5% lower than the baseline value. The latency of the peak in is about and decreases slightly with increasing stimulation amplitude.
The variations of the baseline values and of about 5.6% and 5.4%, respectively, over the range of stimulation amplitudes are likely due to variations in the optical coupling between the time-resolved reflectance spectroscopy (TRS) sensor and the tissue. We estimate that the average number of scattering events for and varies between 45 and 60, which is consistent with the good agreement between measured and the diffusion theory Eq. 1.
Figure 5 shows that the average decay time of the field autocorrelation function follows, in contrast to the evolution of and , a biphasic pattern. For the lowest stimulation voltage showing a functional signal, drops within about from the baseline value of about to about . At later times, increases to a maximum whose height decreases with stimulation voltage. After this maximum, recovers the baseline value via a second, more shallow minimum. The average decay times in the first and in the second minimum reach plateau values of about and for stimulation voltages exceeding . The position of the intermediate maximum in at about shows very weak dependence on the stimulation amplitude.
The reduced intensity autocorrelation functions show a pronounced nonexponential decay that, in accordance with the variations of , strongly varies during the contraction–relaxation cycle (see Fig. 6 ).
Mean square phase fluctuations extracted from measured cover the range of to for lag times . Clearly, is not described by a single power law in , indicating that different mechanisms are responsible for the scatterer displacements (see Fig. 7 ).
By fitting the mixed diffusion-shear model Eq. 2 to the measured mean square phase fluctuations, we obtain, using measured optical properties and in Eqs. 3, 4, time-dependent diffusion coefficients , and the modulus of the strain rate (see Fig. 8 ). Both quantities follow a biphasic pattern similar to the one of , indicating that the temporal evolution of the latter is dominated by temporal variations of the scatterer dynamics, while variations of the optical properties have a minor impact on the DWS signal. The two maxima of reflect the points of maximal elongation rate during contraction and relaxation, respectively, while the minimum separating the peaks reflects maximal contraction. The minimum of is reached at for the lowest stimulation voltage showing a functional DWS signal , in very good agreement with the peak time of . The peak amplitude of is observed to increase with stimulation voltage, likely due to an increased number of activated motor units.
The tissue strain is obtained usingis Heaviside's unit step function. The minus sign in the second term accounts for the fact that the DWS signal does not allow us to determine the sign of the strain rate. This procedure results in effective strains with a single, asymmetric peak at the position of the minimum of (see Fig. 9 ). The height increases with stimulation amplitude, indicating the recruitment of an increasing number of muscle fibers. For the lowest supra-threshold stimulation amplitude of , we find that returns, to within experimental error, back to the baseline level after relaxation. The peak latencies of are practically independent of stimulation amplitude and agree very well with the ones of and .
DWS and TRS were measured to estimate tissue strains associated with contraction and relaxation of the biceps muscle. For the same source–receiver distance, stimulation-induced changes of and are observed for voltages beyond and , respectively. At a stimulation voltage of where both and first show a change with stimulation, the one of (17%) is considerably higher than the one of . The increasing functional reduction of at higher stimulation amplitudes suggests reduction of vascular volume during the repetitive stimulation. Our results contrast the TRS data from plantar flexion of the M. gastrocnemius,24 where and change in the opposite direction. The differences are likely due to the different orientation of the optodes in Ref. 24, and due to the muscle that was able to move under the loosely attached TRS probe.
In our experiments, we observed that the peaks of and reverse their sign when the optical sensor is oriented perpendicular to the muscle fiber. While the sign reversal of might reflect the anisotropy of light diffusion arising the structural anisotropy of the muscle fiber,25 the sign reversal of could have various origins: (1) cross talk with related to our analysis of the TRS data with an optically isotropic model or (2) orientation-dependent weighting of contracting and noncontracting muscle regions characterized by reduced and increased vascular volume, respectively. However, a detailed elucidation of these issues is beyond the scope of this work. 3
Interestingly, the analysis of the mean square phase fluctuation including the measured optical properties is able to resolve the first peak of for a stimulation voltage of (see inset in Fig. 8). At this voltage, and are constant, and no functional DWS signal can be seen in . This shows that the decomposition of into the diffusive and shear modes provides a sensitivity superior to that of . The reason for this is that strongly weights the contributions of long lag times to that arise from dynamics within superficial layers. Indeed, a functional signal can be seen already at by computing , which, for , provides a stronger weighting for the long photon paths. This increase of sensitivity provides direct evidence that DWS and TRS are indeed probing the muscle.
The observation that the mean square phase fluctuations can be explained only by including a term scaling quadratically with lag time during contraction and relaxation indicates that a substantial part of the DWS signal has its origin in deterministic scatterer displacements. In perfused striated muscle, spatial fluctuations of the dielectric constant leading to light scattering arise mainly from the periodic arrangement of sarcomeres and from the erythrocytes within the microvasculature. At , blood at a physiological hematocrit has (Ref. 26). At a typical vascular volume fraction of 3% (Ref. 27) and the value for muscle from our experiments, we estimate that the erythrocytes account for about 12% of the total number of scattering events. We expect that during the rapid contraction shortly after stimulation, the erythrocytes will follow the strain field within the muscle. By virtue of their higher scattering cross section, the erythrocytes within the microvasculature could thus enhance the contrast for the strain rate signal from the muscle.
The origin of the diffusive mode in that closely parallels the strain mode is not entirely clear. It has been noted earlier11 that DWS autocorrelation functions from perfused tissue are described by diffusion rather than by random shear flow for which . Whether the diffusion coefficient measured for muscle tissue is related to residual microvascular blood flow or whether it reflects other diffusive processes, such as crossbridge dynamics, remains to be clarified.
Recently, Yu reported DWS experiments on the M. gastrocnemius during a plantar flexion exercise where reduced DWS decay times were interpreted in terms of enhanced blood flow.28 Our experiments using the same protocol, but with high temporal resolution show that the reduced decay time coincides with the transition from toe-up to toe-down, while the pulsatile background of due to the muscular blood flow is unaffected by the exercise (data not shown).
The question arises whether the increase of the peak effective strain with stimulation voltage shown in Fig. 9 is due to increasing recruitment of motor units or due to increasing peak contraction, i.e., decreasing sarcomere length at maximal contraction. Looking at the force recordings in Fig. 2, we observe that the peak force increases with stimulation voltage, indicating enhanced motor unit recruitment resulting in increased number of contracting fibers within the volume probed by the DWS experiment. Note that in order to extract true volume-averaged muscle strains, the effective strains measured with DWS need to be normalized with the fraction of fibers that are actually contracting. For stimulation voltages above , we observe a saturation of the peak effective strain at , indicating the recruitment of all excitable fibers in the muscle. We thus estimate a value for the true volume-averaged strain of about 0.18, which is in good agreement with fiber length measurements on the cat gastrocnemius.29
In order to explore whether DWS is able to discriminate between fast- and slow-twitch fibers, we performed experiments on the M. gastrocnemius of the same subject, using the same electrical stimulation protocol as in the experiments on the biceps muscle. The subject was in prone position on a table, and the DWS probe was positioned on the muscle belly. In contrast to the biceps muscle, which consists mainly of fast-twitch fibers, the gastrocnemius muscle predominantly contains slow-twitch fibers, which is reflected by a larger latency of the peak force of (Ref. 30). The temporal evolution of the DWS decay time is qualitatively similar to the one in the M. biceps brachii; in contrast, the latency of of the peak in is significantly larger (see Fig. 10 ). Clearly, this difference is not explained by possible differences in excited muscle volume between biceps and gastrocnemius muscle, but rather suggests that it reflects the physiological difference between fast- and slow-twitch fibers. More detailed experiments are required to verify this hypothesis.
In conclusion, we show on a single subject that DWS sensitively and noninvasively measures skeletal muscle contraction with millisecond temporal resolution. Single-photon counting and multispeckle detection allow to place source and receiver at a distance that is sufficient to overcome adipose tissue layers several millimeters thick. Further improvement of signal-to-noise ratio and spatial resolution could thus open the way for quantitative, noninvasive characterization of muscle function in sports medicine and clinical applications.
We thank J. Brader, H. Sprott, and H. Jung for helpful discussions, and G. Maret for continuous support. This work was funded by the Center for Applied Photonics (CAP) Konstanz and partially by Grant Agreement No. 228334 LASERLAB-EUROPE II (EU FP7-INFRASTRUCTURES-2008-1).
Appendix: DWS from a Uniaxially Contracting Cylinder
For a point source at the origin and a point receiver a distance away, the normalized field autocorrelation function for a semi-infinite, optically homogeneous and isotropic medium is given by11. The quantities and are related to the depth of the effective diffuse light source and the extrapolation length . The dynamic absorption coefficient is related to the mean square phase fluctuation of a single scattering event by , , resulting in a linear increase of with lag time .
In order to compute the mean square phase fluctuation resulting from the contraction of the muscle, we use the formalism of Bicout and Maynard for DWS from inhomogeneous flows.31 For the strain rate tensorof scatterers within the tissue, the mean square phase fluctuation is is the probability density that a photon that is injected at the origin and detected at a distance on the surface is scattered at the position ; the index denotes the number of scattering events. The quantity is the wave vector in the medium.
We assume that the volume probed by the diffuse photon cloud undergoes a homogeneous uniaxial deformation characterized by a time-dependent Cauchy strain . The Cauchy strain for a uniaxial deformation from an initial length to a final length is defined as . Within the integration time of the DWS experiment, the rate of change of the Cauchy strain, , is assumed to be a constant. The time-dependent strain rate tensor is given by
The mean square phase fluctuation arising from this deformation is9 involves the quadratic invariant of the strain rate tensor , the mean square phase fluctuation does not depend on overall rotations of the muscle during contraction. Using , we obtain the dynamic absorption coefficient