Noninvasive muscle activity imaging using magnetography.

A spectroscopic paradigm has been developed that allows the magnetic field emissions generated by the electrical activity in the human body to be imaged in real time. The growing significance of imaging modalities in biology is evident by the almost exponential increase of their use in research, from the molecular to the ecological level. The method of analysis described here allows totally noninvasive imaging of muscular activity (heart, somatic musculature). Such imaging can be obtained without additional methodological steps such as the use of contrast media.

Exquisite research methodologies have been developed for the study of the structural imaging of brain function, from single-cell level with electron and multiphoton microscopy, to whole organ imaging with magnetic resonance imaging (MRI) (1). Functionally, brain electrical activity has been addressed, over a similar range, from single ionic channel and single-cell activity (patch, microelectrode recording, and dye imaging) to whole brain activity (electroencephalography, magnetoencephalography). Structure and function are combined in fMRI. While fMRI offers a good spatial resolution (∼1 mm), the temporal resolution is poor (seconds to minutes). By contrast, electroencephalography (EEG) and magnetoencephalography (MEG) offer better temporal resolution (ms), and can be combined with structure via coregistration. This particular combination requires the use of mathematical modeling for coregistration and thus improvement has been of recent interest and relevance (2, 3).

Magnetic signals have been recorded from areas other than the brain (4), to include peripheral nerve (5), spinal cord (6), the somatic (7–11), and cardiac musculature (12–14), and see ref. 15. Indeed, recent advances have allowed magnetocardiography (MCG) recordings in the unshielded hospital setting (16).

We have developed a methodology for the analysis and localization of brain activity based on MEG recordings addressed as functional tomography (17–20). This methodology generates a highly spatiotemporal (1 mm/ms) three-dimensional (3D) map of the energy produced by all magnetic sources located at a given point in space. We call the resulting image of such activity a “functional tomogram.” With this augmented analysis technology, we noted that in addition to brain activity, other structures near the head could also be imaged (19). This finding led us to design a study of the magnetic fields associated with the electrical activity in cardiac and somatic musculature in healthy adults. This recording paradigm allowed us to ascertain whether our analysis methods would provide an image of this activity not seen previously. In both nerve and muscle tissue the physiological basis for the generation of magnetic fields is the flow of transmembrane current. The functional tomograms are a representation of such magnetic fields (Figs. 1, 2, 3, Bottom, 4, Center, and 5, Upper Left).

Fig. 1.

Average functional tomogram of resting activity reconstructed from the frequency spectra from 10 subjects (). The tomogram is embedded in 3D space and can be viewed from any axis. The outline of a head has been superimposed on the tomogram as a guide to the orientation selected for this figure. Note that the spatial distribution of activity includes signals from the brain (yellow-range signal) and muscles of the head and neck region (red-range signal). (Inset) provides the position of the head and shows the nasion and left preauricular fiducial markers. The positions of these markers are superimposed on the tomogram as an aid to orientation. (Subject’s head inside the MEG sensor helmet, resolution = 3 mm.)

Fig. 2.

MCG and heart structure reconstructed from the frequency spectra of one subject. (A) Three cycles of the MCG. (B) Multichannel MCG spectrum. (C, Lower) Functional tomogram of the heart; (C, Upper), image of heart for reference. (Left side of chest under the sensor helmet, resolution = 2 mm.)

Fig. 3.

MCG of a single cardiac cycle. (Upper) Single cardiac cycle with stages of the cardiac cycle marked. (Lower) Current dipoles, calculated at every stage, are plotted over the functional tomogram, reconstructed from the spectra such as that shown in Fig. 2. (Left side of chest under the sensor helmet, resolution = 2 mm.)

Fig. 4.

Magnetic structure of the activated hand musculature, reconstructed from the frequency spectra. (Left) The hand showing the fingers bound by tape and its location within the sensor helmet. Note that the back of the hand and the index finder touched the sensor helmet. Recordings were made with the hand relaxed (no signal) and during isometric extension (when the subject attempted to separate the fingers). (Center) Dorsal view of the magnetic image reconstructed from the frequency spectra of the contracting muscles during isometric extension. The dot above represents the tip of the index finger touching the upper part of the sensor helmet (resolution = 1 mm). (Right) Anatomical drawing of a hand illustrating the dorsal interossei musculature. (Left hand and wrist within the sensor helmet, resolution = 1 mm, color scale as in Fig. 3.)

Fig. 5.

Magnetic image of pain source in back muscles. (Upper Left) Functional tomogram of resting muscle activity in one subject. (Upper Right) Superposition of MEG activity on an image of the right hip showing psoas and iliacus muscles. (Lower Left) Photo of subject with his lower back under the sensor helmet. (Lower Right) Drawing for orientation showing location of psoas and iliacus muscles. (Lying prone with lower back under the sensor helmet, voxel size, 3 × 3 × 6 mm, color scale as in Fig. 3.)

Results

The results obtained in the present set of functional tomograms correspond to actual anatomical data, very much to those obtained by head MRI, or by direct observation in other cases. Recordings were made from the head (Fig. 1), heart (Figs. 2 and 3), and several muscle groups (Figs. 4 and 5) at rest and during voluntary contraction, or under conditions where muscles were intrinsically activated during muscle pain (Fig. 5).

Recordings with the Head in the Sensor Helmet.

To verify that the magnetic recording profiles obtained from neck and body musculature observed from a single subject (19) were generalizable, MEG recordings were obtained from 10 subjects. The subjects sat quietly with their head inside the MEG sensor helmet during the recording (). The resulting functional tomogram is shown as colored voxels in Fig. 1. In addition to the recording of brain activity (Fig. 1, light voxels), an image of the magnetic profile of the head and neck was obtained that reconstructed their anatomy in three dimensions (Fig. 1, dark voxels). An outline of a head and upper neck is superimposed on the tomogram for orientation. The approximate locations of the nasion and left preauricular fiducial markers have been superimposed on the tomogram for reference. The inset of a head image (Fig. 1, Upper Left) with highlighted fiducial markers is included to show the orientation of the tomogram. Thus, imaging of neck muscles and other musculature, while weaker than that of brain, could be obtained by combining MEG recording with the present analysis technology.

Recordings with the Left Chest under the Sensor Helmet.

Given the large electrical signal of the heart, we proceeded to attempt imaging of the cardiac muscle activity by having the subjects lay supine on the recording bed with their left chest under the sensor helmet. Before beginning the recording, the position of each subject was adjusted until the largest cardiac signal was obtained. As shown in Fig. 2, cardiac muscle activity could be imaged, even when the sensors are arrayed in a helmet shape, rather than in a planar configuration (21). Three cycles of cardiac activity are shown in panel A. The multichannel spectrum of this activity, calculated from a 2-min cardiogram, is shown in panel B. A 3D image of the heart, reconstructed from the magnetic data, is shown in the lower image in panel C (for reference, the upper image in this panel shows the orientation of the heart in the tomogram below).

A more rapid display of the heart activity could also be obtained (Fig. 3). Short time image analysis provided a demonstration of the stages of heart activation, allowing single heartbeats to be decomposed into time fragments. This analysis generated an image of the sequential activation of the cardiac muscle and illustrates the current dipoles during a heartbeat. The P wave (atrial depolarization), the QRS complex (ventricular depolarization), and T wave (ventricular repolarization) were distinguished as shown in Fig. 3.

Recording with Several Muscle Groups under the Sensor Helmet.

Hand.

When a hand was imaged during isometric flexion (squeezing an object) or extension (opening against a load, Fig. 4), the respective activation of the lumbrical or interossei musculature could be imaged. Fig. 4, Left shows the hand in its position in the sensor (an outline of the sensor helmet is provided for reference). Fig. 4, Center is the image of the active dorsal interossei muscles with the index finger touching the top of the sensor helmet. Because this analysis method allows the magnetic signal to be imaged in three dimensions, the tomogram was rotated to show the magnetic signal as viewed from the dorsal aspect of the hand. (For reference, Fig. 4, Right illustrates the dorsal interossei muscles.)

Recording from Spontaneous Back Muscle Pain.

We were particularly interested in determining whether images of back muscle activity could be obtained. There is no direct methodology that assesses spontaneous back muscle activity. This is particularly significant as such activity may be pain-related. We proceeded to obtain magnetic recordings while the subjects lay with their lower back beneath the sensor helmet of the MEG instrument as shown in Fig. 5, Lower Left. The subjects were either relaxed or were contracting their back muscles. Several methods for lower back activation were tested. The most effective method consisted of lifting both legs while lying prone. However, while muscle activity was present in the raw magnetic recordings, the imaged signal was weak in all but one subject as described below.

Unexpectedly, strong muscle activity was present in one subject when the back muscles were relaxed. Off-line analysis of this subject revealed a well-defined image of activity localized to the area of the iliacus and psoas muscles (Fig. 5, Upper). Upon further inquiry, the subject mentioned the presence of lower back pain during the magnetic recording. The upper left panel is the tomogram of back muscle activity. The upper right panel superimposes this tomogram on an image of the iliacus and psoas musculature. (The lower right panel shows these muscles in reference to the skeleton.) Thus, an individual with unilateral lower back muscle pain was recorded with the MEG sensor helmet immediately over the lower back. Indeed, the location of muscle magnetic activity was well correlated with the anatomy of the psoas and iliacus muscles, by contrast this activity was not present on the opposite pain-free side.

Representation of the Data and Their Comparison with the Baseline.

The general spectral features of the experimental data were estimated by using the sum of the powers in all channels: , where k is the channel number, the frequency is n, and is the Fourier amplitude (22).

This power spectrum was used to reconstruct the site and timing of electrical activity. Further, note that the inverse problem for each elementary source is solved “in empty space,” without any knowledge about the object under study. This holds correctly for all magnetic studies due to well-known magnetic transparency of the biological tissues.

Fig. 6 is a comparison of the signal power from the contracted hand (blue) and baseline noise power (red). Note that the muscle signal power is 20–30 times higher than the power of the baseline. When the 3D view of the functional tomogram is rendered (as in Figs. 1, 2, and 3, Bottom, 4, Center, and 5, Upper), the transparency of the voxel is taken to be inversely proportional to the power located in this voxel. So, the voxels with low power (baseline) are transparent and cannot be seen when represented together with voxels containing high power (signal).

Fig. 6.

Comparison of the summary power of the magnetic recording for the isometrically extended hand, as described in Fig. 3 (blue) and the summary power for the “empty room” recording or baseline noise (red). Note that the contracted muscle signal is 20–30 times larger than the power of the baseline.

Discussion

A type of human imaging spectroscopy is presented based on magnetic field recordings obtained from different human body sites. This methodology allows a totally noninvasive approach to functional imaging that offers a significant step toward passive functional spectroscopy. A prerequisite for its use requires the generation of somatic electrical fields as occurs during muscle contraction or during the activation of other tissues or organs such as those found in the abdominal cavity. As a possible instance, the activation of uterus and abdominal musculature, in early labor, may also be imaged (11). In addition, the imaging of real-time pathway conduction and the localization of activity in peripheral nerve, spinal cord, brainstem, deep mesencephalic structures, and thalamocortical activity, may also be obtained.

While the utilization of body magnetography is possible using an MEG instrument, our device does not have an optimal configuration, given that it has been optimized for the recording of brain magnetic activity. Thus, a planar recording has been constructed for MCG (see ref. 14) that will provide an optimal configuration. Indeed, Baule and McFee used two coils to record cardiac activity in an unshielded, room-temperature environment (21). The advance of MCG (as all magnetic recording) followed the introduction of the superconducting quantum interference device (SQUID) (4). Indeed, human peripheral muscle activity was recorded in 1972 using a single SQUID magnetometer (7). A summary of the history of MCG, its clinical applications, and comparison with other cardiac recording approaches, has been ably summarized (12–14), and see ref. 15.

Concerning utilization of the obvious advantages, it should be noted that beyond the noninvasive character of magnetic recording is the fact that, in addition to static imaging, this spectroscopy can describe the dynamics of organ function over protracted periods. This could even include functional response to pharmacological and other interventional therapies. At this point the issue, from the engineering perspective, is the optimal configuration of the sensor distribution matrix that would allow for maximum data gathering and optimal dynamic properties.

As important as the configuration of the sensors is a method to analyze the data they obtain. Functional magnetic tomography, as illustrated here and previously (18–20), does provide precise and reversible decomposition of the initial time series into the sum of elementary coherent oscillations that are distributed in space. Given this fact, it becomes possible to select partial spectra of the biological source and to restore the time course for the electrical activity of such a source (19).

The main advantage of the method described here relates to the algorithms utilized. Indeed, applying such algorithms, we transform the set of multichannel time series into the spatial power distribution generated by elementary sources. Furthermore, the resultant time series are totally utilized; thus no experimental information is lost or arbitrarily selected. The number of elementary sources is determined by the time of measurement and frequency band, Nsources ≥ T(νmax−νmin). The precision of each source localization is defined by the grid resolution, which is specified by the researcher (). The number of channels can vary from several probes to hundreds; more important is their sufficiently wide distribution over or around an object.

In short, this methodology for frequency-pattern analysis can be used to investigate the dynamic properties of systems. The images obtained incorporate the size and shape of the active source and at the first level are independent of coregistration with anatomical imaging such as MRI. It can also be applied to other than biological systems. Indeed, the method has been recently applied to the problem of the localization of magnetic nanoparticles (23) that would allow description of nonelectromagnetic explicit biological morphology.

Methods

The magnetic recordings of brain and muscle activity were carried out in a mu-metal shielded room. A CTF Systems, 275-channel instrument was used for all recordings. The sampling rate was 600 or 1,200 Hz and distant noise and artifacts were reduced with a third-order gradientometer (24). Before each recording session, the instrument and distant noises were recorded under conditions identical to those of the recording session. Recordings were made at the New York University (NYU) Center for Neuromagnetism in Bellevue Hospital Center.

The NYU Institutional Review Board approved the study and an informed written consent was obtained from all subjects before the MEG recordings. Participants were recruited from the NYU Medical Center and the local community.

MEG Recordings.

Participants.

The database for this set of recordings comprised 10 healthy adults (5 men and 5 women) 28–76 y of age (mean 41.8 ± 5.4 y; median age 33.5 y).

MEG data acquisition.

The magnetic activity of the brain was recorded following the protocol used in our previous studies (25, 26). The subjects sat upright in the chair of the MEG instrument with their head in the sensor helmet. A set of three MEG recordings, each lasting 7 min, was made from each subject. Two sets of recordings were made with the eyes closed (EC) and one with the eyes open (EO). The subjects were asked to relax, but to stay awake. The location of the head was monitored during the recordings using electrodes placed at three standard fiducial marker points (left and right preauricular points, and the nasion). The head shape, including the location of the three fiducial markers, was obtained for each subject using a 3D tracking system by moving a stylus to each fiducial point and over the surface of the head (Fastrak, Polhemus).

MCG and Magnetic Myographic Recordings.

The database for this study comprised 10 healthy adults (2 women and 8 men; 19–75 y of age; mean, 37.9 ± 6.6 y; median 27.6 y).

For the cardiographic recordings each subject lay on the bed of the MEG instrument with the left side of the chest directly under the sensor helmet. Before the recording, the position of the subject was adjusted until the largest signal was obtained. One 2-min recording was then obtained. For the myographic recordings, each subject lay on the MEG bed with the area of interest (hand, or hand and arm) within and touching the side of the sensor helmet or directly under the sensor helmet (lower back). Two 2-min recordings were obtained for each area. During the first recording the muscles in the area of interest were relaxed. During the second recording the subject performed isometric contraction of these muscles. Fiducial markers were placed on the sensor helmet to provide an orientation signal.

MEG Data Analysis.

The experimental data were analyzed using precise frequency-pattern analysis as detailed in previous publications (18–20). The 275 channels (sensors), K, of the CTF instrument provided a set of experimental vectors, k = 1,…,K. These vectors represent the discrete sampling of continuous functions (20).

Consider the Fourier transform

where are Fourier coefficients for the frequency in the channel array k, and , where is the highest desirable frequency. The coefficient is not considered further since the constant field component has no meaning in SQUID measurements (20, 22). The high sampling rate of 600 or 1,200 Hz represents a continuous function and has sufficient precision to allow the integrals in Eq. to be effectively calculated with discrete Fourier transform (20).

The frequency resolution is determined by T: . This implies that to reveal the detailed frequency structure of the system, it is necessary that: 1) Data are recorded for a sufficient length of time; and 2) All spectra are calculated for the whole duration of the recording procedure time T. In our experiments T was equal to 100–400 s, thus providing frequency resolution 0.01–0.0025 Hz.

The next step of the analysis requires that the multichannel signal be restored at every frequency and that the time functions be analyzed. In all of the channels, the multichannel signal is restored at frequency .

where , and are Fourier coefficients, found in Eq. .

If , then Eq. describes a coherent multichannel oscillation and can be written as

where is the amplitude, and is the normalized pattern of oscillation (20, 22).

Thus, Eq. provides the separation of time and space. Furthermore, it is possible to determine the spatial structure of the source from the solution to the inverse problem. This is the case since the normalized pattern has been obtained and the structure is assumed to be stable throughout the period of oscillation.

The time course of the field is determined by , which is common for all channels, i.e., this source is oscillating as a whole at the frequency (20, 22).

The development of the theoretical bases for the reconstruction of static entities such as neuronal circuits has been detailed in previous papers (18–20). The reconstruction is based on a detailed frequency analysis and extraction of the frequencies that have high coherence and similar patterns.

The algorithm for the mass precise frequency-pattern analysis was formulated in four steps (20, 22):

Discrete Fourier transform of the multichannel signal for the whole recording time T.

Inverse Fourier transform–restoration of the signal at each frequency.

If the coherence at the particular frequency is close to 1, then the pattern and frequency were used as an elementary coherent oscillation (see Eq. ).

If the restored signal consists of several phase-shifted coherent oscillations, then extract those oscillations using independent component analysis (see ref. 27).

After the fourth step of this analysis, the initial multichannel signal is represented as a sum of elementary coherent oscillations:

where is maximal number of coherent oscillations, extracted at the frequency .

Each elementary oscillation is characterized by frequency , phase , amplitude , and normalized pattern and is produced by the functional entity having a constant spatial structure.

We, thus, define the functional tomogram as the electrical functional structure of the system, reconstructed from the analysis of the set of normalized patterns (20). The functional tomogram displays a 3D map of the energy emitted by all of the sources located at a given point in space. In order to build a functional tomogram, the space under study is divided into elementary cubicles with centers in . The cubicle’s edge is determined by the desired precision; in this study, we selected 1.0 mm for the hand (Fig. 4), 2.0 mm for the heart data (Figs. 2 and 3), and 3.0 mm for the whole-head images (Fig. 1). For the back muscles the voxel size was 3 × 3 × 6 mm (Fig. 5). The energy produced by all of the sources at the center of a given cubicle is calculated from the set of trial dipoles .

The magnetic induction at point , for dipole , located at , is calculated as a current dipole in a spherical conductor (20, 28).

All trial dipoles are orthogonal to , since the vector product is nonzero only for those dipoles. Trial dipoles cover the circle in directions with 360/ degrees step (In this study =8). The set of normalized trial patterns is then calculated as

More than 2.5 million trial patterns were used for each object, calculated at the grid, covering the space of experiment (20).

For each normalized pattern from Eq. , the following function was calculated, giving the difference between this pattern and one of the actual trial patterns:

The position and direction of the source producing the experimental pattern were determined by numbers, providing the minimum to the function over the variables .

The energy of this source is added to the energy produced from the cubicle with the center at .

Following this procedure for all normalized patterns: , it is possible to determine the spatial distribution of energy, produced by all oscillations from Eq. . The result of such distribution is then the functional tomogram of the system under study (20, 22).

Data Availability.

The raw data (converted to the GNU Octave format) and computational scripts are available from the Open Science Framework repository (29).

The lists of data and scripts, and instructions on how to use them are presented in .