Four-dimensional in vivo X-ray microscopy with projection-guided gating.

Visualizing fast micrometer scale internal movements of small animals is a key challenge for functional anatomy, physiology and biomechanics. We combine phase contrast tomographic microscopy (down to 3.3 μm voxel size) with retrospective, projection-based gating (in the order of hundreds of microseconds) to improve the spatiotemporal resolution by an order of magnitude over previous studies. We demonstrate our method by visualizing 20 three-dimensional snapshots through the 150 Hz oscillations of the blowfly flight motor.

Studies of animal movement are founded upon our ability to measure and model the underlying musculoskeletal mechanisms, and have long driven the advances of imaging techniques. In the late 19th century, Muybridge's photographic time series of galloping horses provided the first visualizations of locomotor behaviors too fast to be resolved by the human eye1, presaging the age of cinematography. A century and a half on, modern motion-capture technology has made the imaging of fast external movements almost routine, and has also facilitated radiographic imaging of some internal movements2. X-rays penetrate tissues that are opaque to visible light3, and have been used extensively in visualizing vertebrate motion4. Radiographic imaging of live invertebrates has only become practical more recently, thanks to the very short exposure times and high spatial resolution permitted by the brilliance of X-rays of third-generation synchrotron light sources57. Radiographic time series of living insects have already led to deepened understanding of the mechanics of tracheal respiration7, and of pistoning motions in the gut of crawling caterpillars8. But, because they overlay all of the tissues in the X-ray's path, radiographic projections alone cannot describe complex three-dimensional movements. The four-dimensional (4D) reconstructions that this requires may instead be obtained through computed tomography (CT) using radiographic projections acquired from multiple viewing angles at multiple phases of the movement cycle910.

Non-repetitive motions can only be reconstructed in 3D if they occur on a timescale much longer than the time needed to obtain a complete set of angular projections10. However, for periodic or quasi-periodic movements, a three-dimensional image can be obtained by gating the projections at defined phases of multiple, consecutive movement cycles. The resulting tomograms do not represent a time-continuous description of a single period of movement, but instead represent a series of averages over multiple movement cycles. In cardiac clinical CT, for instance, radiographic acquisition is gated by the electrical activity of the heart as an external indicator of movement phase11. The same technique has been extended to image cardiac and respiratory movements in mice, achieving 12 ms temporal resolution and 150 μm voxel spacing in a clinical CT scanner12, and 17 ms temporal resolution and 20 μm voxel spacing at a synchrotron facility13. Impressive as these results are, each prospective or retrospective gating technique published so far relies on the acquisition of an external signal that is then used for gating. The movements of many animals, including those of most insects, fall beyond the reach of this approach.

To address this challenge, we developed a 4D synchrotron imaging technique combining single-exposure phase retrieval with retrospective gating. The resulting tomographic time series enable us to visualize the fast internal movements of the blowfly flight motor on sub-millisecond and micrometre scales. Said to be the most complex hinge in nature1415, and containing more than 40 individual muscles oscillating at 150 Hz, the dipteran flight motor serves as a challenging target on which to demonstrate the potential applications of this technique to a range of biomechanical systems. In a first iteration of this technique16, we gated the projections retrospectively in a traditional way using an external indicator of the phase of the movement cycle (the fly's wingtip trajectory, filmed using a pair of high-speed cameras). Using an external gating signal is cumbersome, however, and may not always be possible. Hence, in order to generalize the method to other systems, we present here a refined technique that uses retrospective gating guided by the radiographic projections themselves, and provide a detailed validation of the performance of the method under various modes of data acquisition and analysis. Two sets of experiments were performed at the TOMCAT beamline17 of the Swiss Light Source. In the first set of experiments, we evaluated two different acquisition protocols using broadband radiation, imaging the anterior part of the flight motor, and the head-neck motor system. In the second set of experiments16, we imaged the complete flight motor, comparing monochromatic and broadband radiation. This allowed us to validate the robustness of retrospectively gated tomography to differences in experimental configuration. Synchrotron-based tomographic imaging requires the sample to be rotated quickly during radiographic acquisition, which could affect the behavior of a living sample in different ways depending upon rotation speed. To allow for an informed choice of this behaviorally significant parameter, we therefore explored the effects on image quality of using different rotation speeds.

Results

We began by acquiring radiographs of blowflies flying tethered to a rotating stage in single slow scan (SSS) mode (rotation speed: 54°s−1), single and multiple fast scan modes (SFS and MFS) (rotation speed: 332°s−1 or 347°s−1), and multiple medium scan mode (MMS) (rotation speed: 180°s−1), where a scan refers to a sample rotation of 180°. In the acquisition of multiple fast scans, we rotated the sample continuously. We compared the results obtained using monochromatic versus broadband (polychromatic) radiation. Single exposures of 0.2–0.4 ms and a total scanning time as short as 0.5 to 3.3 s allowed us to capture the fly spontaneously beating its wings with minimal motion artefacts. We retrieved the spatial distribution of the refractive indices of the tissues from the recorded phase deviation of the partially coherent synchrotron X-rays in the sample18. Compared to classical attenuation contrast, phase contrast enhances the visibility of tissue boundaries at a reduced radiation dose deposited in the specimen (Refs. 19,20, Methods). The projections from a single sequence were sorted according to the phase of the wingbeat cycle, using the spatial cross-correlation of each projection with the next to estimate the anteroposterior movement of the thorax between projections, which served as a retrospective gating signal (Fig.1b–c). Finally, for each of the 10 to 20 evenly spaced phase steps that we identified, we reconstructed tomographic sections through the blowfly thorax from approximately 500 projections. The wingbeat frequency of the blowflies was approximately 150 Hz, so the resulting tomographic time series have an effective temporal resolution of 1.5 to 3 kHz, depending upon scan mode.

Figure 1

Experimental arrangement and retrospective gating.

(a) The X-ray micro-tomography acquisition system (b) Normalized projections at four stages of one wingbeat cycle, with local image motion vectors (yellow arrows). The translation vectors were found by spatial cross-correlation of consecutive projections tracking the movements of twelve evenly spaced regions each containing 24 × 24 pixels. Orientation: head-neck joint at the top, posterior end of the thorax at the lower end of the images. Sequence of wingbeat phases from left to right: dorsal stroke reversal, downstroke, ventral stroke reversal, upstroke. (c) The gating signal (blue) derived from spatial cross-correlation of successive projections (after smoothing). The peaks of the signal (red points) identified the beginning of each new wing beat cycle, which was divided into uniformly spaced wingbeat phases. The green, pink, and orange points correspond to the phases of the wing beat shown in (b). (d–g) Rendering of the fly (d,e) thorax and (f,g) wing hinge () during (d,f) dorsal and (e,g) ventral stroke reversal (orientation: anterior towards the top, posterior towards the bottom). A frontal clipping plane has been introduced to reveal the large dorsal longitudinal and dorsoventral flight muscles, DLMs () and DVMs (*), respectively. The ladder-like structure spanning the flight muscles is the tracheal network (→). At the level of the wing hinge, clear deformations of the exoskeleton are visible.

We validated the performance of our method under each acquisition condition using standard image quality metrics (see Methods). The contrast between air and tissue was more than one order of magnitude higher in the reconstructed phase maps than it was in the attenuation-based reconstructions (Supplementary Fig. 1), although it remained difficult to distinguish between exoskeleton and soft tissue. Nevertheless, important functional elements of the fly's anatomy, such as the individual muscle units of the main flight muscles, were readily resolved at a spatial resolution of 9 to 31 μm, depending upon acquisition protocol (see Methods). The best overall scan quality was achieved using the ‘slow’ scan protocol with broadband radiation (spatial resolution: 10 μm; signal-to-noise ratio: 24 as detailed in the Methods). We validated our method by comparing its results to those of standard (i.e. non-gated) tomography of a fly at rest, acquired using monochromatic X-rays20. The spatial resolution of 15 μm and signal-to-noise ratio of 26 that we achieved using standard tomography were close to that of the retrospectively gated tomography (see Supplementary material). In conclusion, the two acquisition protocols, one involving radiographic acquisition during a single slow scan, and the other involving multiple ‘fast’ or ‘medium’ scans both have their own advantages and disadvantages. The lower sample rotation speed (54°s−1) that we used in the ‘slow’ scan diminished the perturbation experienced by the animal, but perhaps more importantly, it decreased the motion blur due to aperiodic movements and resulted in a more uniform sampling of the projection space. The fast rotation protocol, on the other hand, is likely to exhibit increased motion blur in this case, as well as strong streaking artefacts due to the uneven filling of the projection space12. Nevertheless, a pragmatic advantage of using multiple fast scans is that reconstructions will still be available, albeit of lower quality, if the animals stops moving before all scans are completed.

As a visual demonstration of the quality of the data obtained, we used a reconstruction obtained using monochromatic X-rays at 18 and 15 keV in the fast and medium scan modes, respectively, to create volume renderings of the fly's dorsoventral and dorsal longitudinal flight muscles (DVM and DLM, respectively; Fig. 1 d–g). The six DVMs attach to the top and the bottom of the thorax, and contract alternately with the DLMs, which run along the anteroposterior axis of the thorax. Together the DVMs and DLMs power the wingbeat by causing the thorax to deform in a way that transmits force to the wing hinge. These muscles are prominent in the three-dimensional renderings of the thorax shown at ventral and dorsal stroke reversal (Fig. 1 d,e, Supplementary Video 2). Volume renderings of the wing hinge at the points of stroke reversal reveal the strong deformations of the cuticular exoskeleton that these muscles produce (Fig. 1 f,g). As a further demonstration, Fig. 2 shows the first ever 4D visualization of the tracheal network of a living insect. The different orders of tracheae in a blowfly cover diameters from hundreds of micrometers down to less than a micrometer; Fig. 2 shows branches with diameters between 20 and 50 μm. The major branches of the tracheal network separate the individual subunits of the DLMs, and can clearly be seen spanning the flight muscles. Animated renderings of the tracheal network (Supplementary Video 2,3) allow us to visualize the heaving motion of this complex structure at wingbeat frequency.

Figure 2

Visualization of the air-filled tracheal network spanning the dorsal longitudinal muscles (DLMs), and rendered in isolation from the rest of the thorax at the time of dorsal stroke reversal.

The inset panel provides a similar rendering of the tracheal network within the rest of the thorax for context (see also Supplementary Movie 2). The large opaque green structures in the inset panel are air spaces inside the DLM.

In Fig. 3(a,b) we show one frame from the volume renderings of the left, anterior DVM for ten equally-spaced phase steps through the wingbeat cycle. The section of cuticle covering the flank of the thorax is clipped away in the volume renderings to allow for an unobstructed view of the flight muscles. Thoracic deformations are clearly visible in the rendered time series, as are the contractions of the DVMs that drive these movements. The length of the left, anterior DVM was calculated at each time step for 3 pairs of virtual anatomical landmarks at the dorsal and ventral ends of the left DVM, corresponding to the intersections of the four subunits of the muscles (Fig. 3a,b). As a demonstration of the practicality of making automated measurements on these data, template matching was applied to measure the translational movement of each marker by maximizing the three-dimensional cross-correlation between a 25 × 25 × 25 voxel region around the marker between consecutive time steps (Fig. 3c). The length oscillations normalized by the average distance between a pair of markers were consistent between pairs, with the ventral markers moving up to 46.8 μm per time step. This movement mainly occurred in a plane almost normal to the muscle's dorsoventral axis. The maximum muscle strain was calculated by normalizing the peak-to-peak amplitude by the mean distance between the dorsal and ventral markers. The average muscle strain for the three marker pairs was 1.62% (standard deviation: 0.30%), which is consistent with videographic measurements of the displacements of external markers painted on the exoskeleton of fruit flies and observed during tethered flight21.

Figure 3

3D visualization of the flight muscles.

(a) Rendering of dorsoventral muscle (DVM) with markers (wing beat stage: ventral stroke reversal). The dorsal longitudinal muscles (DLMs) are hidden behind other muscles in these visualizations. Three marker pairs were manually chosen in the region of the attachment site of the first unit of the DVM. The frontal (l-f), medial (l-m) and hind (l-h) legs are visible as is the attachment point of the fly holder (+). (b) Magnification of the ventral area of (a) with paths of the three markers at the ventral attachment site of DVM1. (c) Length oscillation and strain estimates of the first unit of the dorsoventral muscle (DVM1).

Effects of radiation dose

The technique that we have developed here can be implemented using either broadband or monochromatic radiation. Flies exposed to either spectrum of X-rays (broadband or monochromatic) would typically stop flying within ten seconds of the start of the exposure, and later died, suggesting that lethal doses were incurred within seconds of the start of imaging. The estimated total dose16 was between 350 Gy and 1300 Gy: levels similar to those reported by Socha et al.5 to have no observable or permanent physiological effects, respectively, at lower incident beam fluxes. However, in our experiments, the dose was applied on a timescale of seconds rather than minutes, and the applied dose rates (between 90 and 325 Gy s−1) were one to two orders of magnitude higher than those applied in Ref. 5. Furthermore, whereas the insects in Ref. 5 were irradiated in a relatively inactive state, the insects in our experiments were performing the metabolically expensive action of flight. The flies' sensitivity to these high radiation rates is not unexpected5, but certainly merits further investigation.

Discussion

In summary, using a projection-guided retrospective gating method, we were able to resolve and quantify the contractions of individual muscle sub-units powering the 150 Hz wing beat cycle of the blowfly, and were also able to visualize the volumetric changes of the tracheal network that ventilates and cools the flight motor. These results complement our earlier measurements of steering muscle strains, made using a similar technique where the wingtip kinematics were used to provide an external gating signal16. We achieved 4D reconstructions of the flight motor in vivo at 1.5–3 kHz effective temporal resolution with a voxel size of 3.3–5.0 μm, and corresponding true spatial resolution of 10–15 μm. This represents an order of magnitude improvement in spatiotemporal resolution as compared to previous studies. The technique that we have presented here is applicable in principle to any periodic biomechanical system that can be imaged successfully in a synchrotron. Considering its unmatched spatiotemporal resolution, the method that we have presented has the potential to revolutionize in vivo studies of functional anatomy, biomechanics, and physiology in small animals.

Methods

Preparation of specimens

Blowfly pupae (Calliphora vicina) were taken from the lab culture at the Department of Bioengineering, Imperial College London and kept on a 24 h (12:12) light-dark cycle. Specimens were used within two weeks of eclosion. They were briefly anesthetized for 8–10 minutes by cooling to 4°C, and fixed to a wooden tether at their scutum with a mixture of beeswax and collophonium. Care was taken to avoid the wax coming into contact with the scutellum. We fixed the wooden tether to a custom-made holder attached to the rotation stage of the beamline end station, and aligned the anteroposterior axis of the animal with the tomographic rotation axis. The results presented here were obtained from three specimens per acquisition mode (SFS reconstructions were based on a subset of projections from the MFS acquisition). Each specimen was scanned once.

X-ray instrumentation and acquisition

The X-ray source was a superbending magnet located 25 m from the sample. For experiments with monochromatic light, a double crystal multilayer monochromator was used 7 m downstream of the source to extract monochromatic X-rays with a bandwidth of 2% at 18 and 15 keV photon energy (wavelength = 0.07 and 0.08 nm). The beam was 4.1 mm high and 10 mm wide at the sample position, and the flux at the sample was of order 1012 ph/s1/mm2 in this configuration. We selected a region of interest of 4.1 × 3.32 mm (1242 × 1008 pixels) to accommodate the whole thorax. The monochromator was omitted for experiments with broadband light, giving about 50 times greater total photon flux at the sample position and increasing the beam height to 5.7 mm. The broadband X-ray beam was filtered to 5% of its original power and accordingly the bandwidth was narrowed down around a peak value of about 35 keV (0.03 nm).

Radiographic projections were recorded at either 1840 Hz (0.4 ms exposure time) with broadband radiation, or at 2500 Hz (0.3 ms exposure time) or 4432 Hz (0.2 ms exposure time) with monochromatic radiation of 18 and 15 keV respectively, using a pco.Dimax 12 bit CMOS detector system. The detector was mounted on a custom made high numerical aperture microscope (Elya solutions, s.r.o.) featuring a continuous zoom-in option and a 100 μm thick LuAG:Ce scintillator crystal in its object plane for converting the X-ray photons to visible light photons that were collected by the CMOS sensor. The distance between sample and scintillator was 150 mm in broadband configuration, and 350 mm in monochromatic configuration. Isotropic voxel edge lengths were 3.3 μm (SSS, SFS, MFS) and 5.0 μm (MMS). These differences are the result of optimizing the sample to scintillator distance separately in each configuration: because we used a shorter single frame exposure time in the monochromatic configuration, a greater propagation distance provided better SNR as defined by Equation (1) after phase retrieval, with some accompanying loss in resolution.

Phase retrieval

Previous studies67 conducted at partially coherent X-ray sources relied on the ‘edge-enhancement’ effect caused by interference of X-rays at the places where the refractive index has edges, and did not consider retrieving the actual phase shift of the probing X-rays in the sample. In phase retrieval experiments, ‘edge’-type contrast is transformed into ‘area’-type contrast, which provides segmentable images in 3D. We performed phase retrieval in a qualitative manner using our implementation of the single-image phase retrieval algorithm described previously1822 under the assumption that the object consisted of a homogeneous soft tissue material. For the monochromatic beam (18 or 15 keV), we assumed that the refractive index decrements and absorption indices of the material were δ18 = 7 × 10−7 or δ15 = 9.9 × 10−7 and β18 = 5 × 10−10 or β15 = 9.8 × 10−10, respectively. For the broadband beam, we assumed that the mean x-ray energy was 35 keV and used values of δ35 = 2 × 10−7 and β35 = 10−10.

Projection-guided retrospective gating and tomographic reconstruction

Each complete set of projections comprised 8350 radiographs taken for the angular range of 180° for the ‘slow’ rotation case (SSS), 4432 for the medium rotation speed (MMS) and 1000 or 1300 projections for the ‘fast’ rotation case (SFS, MFS). The first scan (i.e. 180° rotation) of an MFS acquisition was used for each SFS reconstruction. The fly beat its wings approximately 480 times during each 180° scan in the ‘slow’ rotation case, and approximately 80 times during each 180° scan in the ‘fast’ rotation case. Each 180° scan therefore sampled multiple wingbeats, and multiple stages of the wingbeat cycle. We developed projection-guided retrospective gating to group projections taken from different angles but at the same phase during the wing beat cycle. This is important in keeping experimental complexity to a minimum, by avoiding the need to identify and record a separate external gating signal.

The main movement of the animal's thorax during tethered flight was along the anteroposterior body axis. We therefore used the spatial cross-correlation of each projection with the next to estimate the vertical movement of the thorax between projections (Fig.1 b,c). The time course of these vertical translation estimates was smoothed by convolution with a discrete Gaussian filter (1.4 frames). We used the peaks in the resulting periodic signal to identify the beginning of each wingbeat cycle, and sorted projections taken at the same phases on this basis (Fig.1c). Visual inspection of the projection groups was used to verify the success of the algorithm in sorting the projections. Tomographic reconstruction was performed using a Fourier transform-based algorithm23.

Estimation of scan quality using signal-to-noise ratio and resolution

The local spatial resolution was quantified using the edge response function at naturally occurring sharp edges in reconstructed slices of the object. Eight such edges from muscle tissue to air were chosen manually in the region of the dorsal longitudinal muscles. We fitted an error function to a short line profile of two pixels width, and identified the point spread function with a Gaussian function as described in Ref. 24. The full width at half maximum of the calculated point spread function served as a measure of the spatial resolution. This metric indicates the minimal distance between two point sources at which their images can still be distinguished.

A second criterion was employed to measure the global spatial resolution in twenty equally spaced transverse slices at the level of the anterior thorax. We applied the technique described by Modregger et al.25 to divide the power spectrum of the images into object image signal and process noise. The resolution was then defined as the critical wavelength at which the contribution of the noise-free image was equal to the contribution of the noise.

The signal-to-noise ratio (SNR) describes the distance of the object intensity distribution from the background intensity distribution and is closely related to the ability to automatically segment the image. SNR is defined as:μo and μb are the mean intensity of patches of 12 × 12 object and background pixels, respectively, and where σo and σb are the standard deviation of pixel intensity in the same object and background patch. Twelve pairs of patches were chosen pseudo-randomly across the field of view in each of four slices, and for each of the data sets. We assumed that the true refractive index of the object represented by these patches was uniform.

Tomograms from 480 projections in the slow scan setup (Fig. 1) had a SNR of just 0.43 (±0.18 SD) if only attenuation and edge-enhancing phase contrast was reconstructed (Supplementary Fig. 1a). This was just sufficient to separate object and background voxels. Phasing fringes around the tissue-air-interfaces are visible as a halo-like edge enhancement effect in the attenuation-based reconstructions, indicating that there is additional information contained in the data that can be recovered through phase retrieval. Reconstructions obtained using phase retrieval (Supplementary Fig. 1b, Supplementary Video 1) had a much higher SNR, ranging from 16.2 (±5.7 SD) at worst (near the middle of the upstroke, when movement is fastest), to 17.8 (±4.4 SD) at best (at the end of the upstroke, when movement is slowest). The spatial resolution, defined as the full width at half maximum of the point spread function, was estimated at 16.2 (±2.4 SD) μm at worst (near the middle of the upstroke), and 14.3 (±2.0 SD) μm at best (at the end of the upstroke). The lowest spatial frequency that could still be resolved corresponded to 10.2 (±0.4 SD) μm at the end of the upstroke.

Tomographic sections reconstructed from 83 phase contrast projections acquired across a single fast scan with filtered broadband radiation (Supplementary Fig. 2a,b) exhibited a best-case SNR of 6.8 (±1.6). Combining the data acquired across six consecutive ‘fast’ scans of this fly gave a total of 507 projections (Supplementary Fig.2c,d) and improved the SNR to 15.7 (±2.9). The local spatial resolution in both cases (11.4 ± 1.3 μm and 15.4 ± 3.4 μm, respectively) was similar to that obtained with a single slow scan, as was the global spatial resolution (9.3 ± 0.2 μm and 10.3 ± 0.3 μm, respectively) (Supplementary Fig. 3). Visual inspection of the tomographic data revealed streaking artefacts in all of the data obtained from single and multiple fast scans (Supplementary Fig. 2a–d). These artefacts were not present in the data obtained from a single slow scan (Supplementary Fig. 1). Tomograms of the fly thorax acquired using monochromatic radiation at 18 keV (Supplementary Fig. 2e,f) had a higher SNR (23.5 ± 6.8) compared to the broadband case, but had somewhat decreased spatial resolution (31.4 ± 9.7 μm) as measured by the point spread function. These results are summarized in Supplementary Fig.3. The spatial resolution of gated tomographic imaging depends not only on the properties of the instrumentation used, but also on the consistency of the structure captured in different projections. This, in turn, is affected by the uniformity of the behavior or process under investigation, and by the movement velocity of parts of the specimen and the frame rate of the detector26. It is therefore not surprising that we obtained both good spatial resolution (13.5 ± 2) and SNR (16.3 ± 2.5) for the MMS (Supplementary Fig. 2g,h) with only 0.2 ms single exposure time and 15 keV monochromatic X-rays at a larger pixel size (5.0 μm instead of 3.3 μm).

Author Contributions

D.A.S. and R.M. conceived this study. R.M. and M.S. designed the fast, synchrotron-based tomographic microscope. D.A.S., M.W., H.G.K., G.K.T., S.M.W., M.S. and R.M. designed the experiment. D.A.S., G.K.T., M.D., M.W., R.M., S.M.W. and T.M. performed research. D.A.S., M.D., R.M. and S.M.W. processed and analyzed data. D.A.S., R.M., H.G.K., M.S., S.M.W. and G.K.T. wrote the paper.