Emphysema quantified: mapping regional airway dimensions using 2D phase contrast X-ray imaging

: We have developed an analyser-based phase contrast X-ray imaging technique to measure the mean length scale of pores or particles that cannot be resolved directly by the system. By combining attenuation, phase and ultra-small angle X-ray scattering information, the technique was capable of measuring diﬀerences in airway dimension between lungs of healthy mice and those with mild and severe emphysema. Our measurements of airway dimensions from 2D images showed a 1:1 relationship to the actual airway dimensions measured using micro-CT. Using 80 images, the sensitivity and speciﬁcity were measured to be 0.80 and 0.89, respectively, with the area under the ROC curve close to ideal at 0.96. Reducing the number of images to 11 slightly decreased the sensitivity to 0.75 and the ROC curve area to 0.90, whilst the speciﬁcity remained high at 0.89.

. Schematic of the analyser based x-ray phase contrast imaging system (not to scale).
Every material has a wavelength-dependent complex refractive index, n = 1 − δ(λ) + iβ(λ). The wavelength-dependent refractive index decrement, δ(λ), determines how X-rays refract through the medium. If an δ(λ) object is comprised of multiple interfaces that refract the beam many times within an area smaller than a detector pixel, the resultant angular dispersion is often called USAXS [26]. Attenuation of the beam is governed by the linear attenuation coefficient µ(λ) = 2kβ(λ), where k = 2π/λ. Taking multiple images as the crystal is rotated (rocked) through the Bragg condition produces so-called rocking curves (see Fig. 2). From these it is possible to determine how the X-rays were attenuated, refracted and scattered; a procedure sometimes called Multiple-Image Radiography (MIR) [26][27][28][29]. Here we follow the MIR reconstruction process developed by Kitchen et al. [30].
Measuring rocking curves with and without the sample alters the rocking curves, as shown in Fig. 2. Kitchen et al. [30] showed that a Pearson type VII function provides accurate fitting to such curves: Here I(r ⊥ , z, λ, θ) is the intensity recorded by the detector. We define the intensity at the exit surface of the object as I(r ⊥ , z 0 , λ, θ), where r ⊥ denotes the position vector in the (x, y) plane (see Fig. 1). All variables on the right side also depend on r ⊥ , z and λ, but this dependence is dropped for simplicity here and below. θ is the angle of the x-ray beam relative to the Bragg angle θ B . Variables c and m respectively parameterise the amplitude and slope of the rocking curves. Variables a and m define the full-width half-maximum (FWHM) of the curve as ω(r ⊥ , z, λ) = 2a m(2 1/m − 1).
Scattering from the sample causes the width (ω) of the rocking curve to broaden [31], as seen in Fig. 2(a). Deconvolving the fitted sample rocking curve (ω samp ) from the reference rocking curve (ω 0 ) isolates the amount of scatter produced by the sample, giving ω scat (see Kitchen et al. [30]). Deconvolution with Pearson VII distributions is not straightforward [32], hence we instead subtract the widths in quadrature, which is strictly only exact for true Gaussian distributions, but is a good approximation for the experimental curves seen in Fig. 2(b), giving: The attenuation map of the sample I(r ⊥ , z 0 ), normalised to the incident intensity I 0 (r ⊥ , z 0 ), is calculated from the integral of the rocking curve [R sample (r ⊥ , z, λ, θ)] normalized by that of the reference rocking curve with no sample [R 0 (r ⊥ , z, λ, θ)], measured at every pixel: Normalization corrects for variation in the incident intensity and beam divergence across the detector, as well as imperfections in the Laue crystal reflectivity. T is the total projected thickness of the object in the z-direction, and θ min,max are the practical extents of the rocking curve at negligible reflectivity.
Refraction is calculated from the angular difference between the rocking curve peak measured with the sample (θ y s ) and that of the reference rocking curve (θ y 0 ) (see Fig. 2(a)): Finally, the phase information can be calculated from ∆θ y by integration along the direction parallel to the diffraction plane of the analyser crystal (y-direction in Fig. 1): Here y min and y max are the physical extents of the image in the y-direction.

Measuring projected lung air volumes
Measuring the airway diameter requires determination of the projected thickness of air in the lungs. This can be achieved by approximating the animal to be composed of two materials, namely water and bone, and first solving for the contact intensity I(r ⊥ , z 0 , λ) and phase ϕ(r ⊥ , z 0 , λ). This is justifiable since the X-ray properties of soft tissues are very similar to those of water. As demonstrated by Kitchen et al. [33], the projected thicknesses of two materials, here water and bone, can be calculated using [33,34]: In these equations, the dependence on λ and z has been dropped for simplicity. By immersing the animal in a water bath of constant thickness T bath , the projected thickness of air in the lungs is given by 2.3. Determining the number of scattering objects in projection von Nardroff [35] theorised that broadening of the rocking curve due to microscopic spheres should be proportional to √ N (where N is the number of spheres), dependent on δ(λ), and independent of their diameter, D. Here we present experimental verification of these facts and determine the relationship between scattering and thickness using microspheres of known diameters. von Nardroff [35] also suggested that if the projected thickness T of the sample is known, then the diameter of the spheres can be readily computed. Our aim was to measure the number of airways from the USAXS signal and determine their diameter. This can be achieved using the relation:

Mouse model
Animal experiments were approved by the Animal Ethics Committees at SPring-8, Japan, and Monash University, Australia. A total of 28 adult male BALB/c mice were studied. Emphysema was induced using elastase seven days before imaging, as described in [36]. Mice were lightly sedated with isoflurane for the orotracheal administration of the elastase solution. Six, three or zero units of porcine pancreatic elastase (six units, n=12; three units, n=8; and control group, n=8) were diluted in phosphate buffered saline (PBS) and administered as a 50 µL bolus per 30 g mouse weight. Immediately before imaging, mice were killed with an intraperitoneal injection of pentobarbitone sodium (>100 mg/kg), then shaved and any remaining fur was covered in glycerol to reduce unwanted phase contrast effects arising from the fur. An endotracheal tube was inserted via tracheotomy into the mid-cervical trachea and connected to a custom designed small animal ventilator [37]. The lungs were inflated with nitrogen to a pressure of 20 cmH 2 O.

Imaging experiments
X-ray images were acquired at the SPring-8 synchrotron in Japan, in Hutch 3 of beamline 20B2 at the Biomedical Imaging Centre, 210 m from the bending magnet x-ray source. A double bounce Si(111) monochromator was used to select a beam energy of 26 keV. Two imaging modalities were employed: ABI and propagation-based phase contrast computed tomography (PC-CT).

Tomographic imaging
To confirm the degree of emphysema induced via the elastase treatment, high-resolution PC-CT data was collected for all mice. To prevent movement artifacts, mouse lungs were ventilated with nitrogen to a set pressure of 20 cmH 2 O and the trachea ligated before the mice were set in a 2% agarose solution. For each mouse, a total of 1800 projections were recorded over 180°r otation. The exposure time was 100 ms per projection and the object-to-detector propagation distance set at 1 m. At this distance, the high spatial coherence of the X-ray source produced strong Fresnel fringes at the air-tissue interfaces, enhancing their visibility [38]. Reconstructing the data using filtered back-projection combined with the phase retrieval algorithm of Paganin et al. [39] provided low noise, high resolution reconstructions of the lungs in 3D [11]. The detector was a Hamamatsu ORCA flash C12849-101U with a straight fiber optic element coupling the sCMOS chip to the 10 µm Gadox (Gd 2 O 2 S; P43) phosphor with 6.5×6.5 µm 2 pixel size. However, this limited the field of view to an area 13.3×13.3 mm 2 , hence not all of the lung tissue was reconstructed. In hindsight, we could have used 360°CT with twice as many projections to capture the whole chest cavity. Regardless, the FOV was sufficiently large to estimate the average airway dimensions in various regions throughout the lungs. Example images are shown in Fig.  5(a)-(c) below for mice with varying levels of lung disease. Alveolar dimensions were measured from the 3D CT data using custom-written greyscale granulometry software (Interactive Data Language; IDL v8.1) [40,13]. Spherical structuring elements of various sizes were created to survey central regions of the lung for airways of similar size using the morphological opening operator. The CT volumes were first magnified by a factor of four and bilinearly interpolated to increase the spatial sampling rate. The 3D granulometry showed that the degree of emphysema present was not readily predictable by the initial elastase dose. That is, some mice were treated with a high dose of elastase, but showed little to no evidence of induced emphysema, while some mice in the low dose group showed evidence of severe emphysema. No mice in the control group showed evidence of significant alveolar enlargement.

Analyser based phase contrast imaging experiments
We used a Laue analyser crystal to produce a large enough beam to view the entire thorax of each mouse in a single exposure, as shown in Fig. 1. The beam was collimated to 31 mm × 31 mm to be just larger than the detector field-of-view (FOV). The Laue analyser crystal was made from a 100 µm thick silicon wafer attached at the base to the monolithic silicon slab to minimise distortion of the atomic planes (Sharan Instruments Corporation, Aomori, Japan). The Si(111) planes used to reflect the beam in this experiment were perpendicular to the surface of the crystal. The detector for the ABI experiments was a Hamamatsu ORCA-Flash4.0 (C11440-22C) coupled to a tandem lens system with a 25 µm thick Gadox (Gd 2 O 2 S; P43) phosphor. The effective pixel size was 15.2 µm (full FOV = 30.72×30.72 mm 2 ). Images were acquired using a 300 ms exposure time per frame.
Rocking curve measurements were made with and without the specimens in the beam. Each scan rotated the crystal from -50.9 µrad to 50.9 µrad about the Bragg peak. A total of 141 images were recorded for each scan in increments of 0.73 µrad. Rocking curves were numerically approximated at each pixel by least squares fitting using the Pearson type VII function (Eq. (1)). Since the rocking curves are not perfectly symmetric, we only used some of the images spanning just over half of the data to cover one side of the rocking curve and about 10% down the other side. For the images acquired without the sample, we used 80 of the 141 images. This improves the curve fitting using a symmetric function such as the Pearson VII [30]. With the sample present, we used either 80 or 11 images for image reconstruction, as shown in Fig.  2(b). Using 80 images would give the technique the best chance of success, and using only 11 images shows the potential for significant reduction of dose and acquisition time. In the sparsely sampled case, image data was acquired at 11 angular positions (−38 to + 5 µrad) about the Bragg peak. Since the refraction and absorption information is especially sensitive to the rocking curve angular position and maximum intensity, respectively, points were more finely sampled about the expected peak position (see Fig. 2(b)). The choice for 11 angles stems from the work of Majidi et al. [41] who showed that no more than 11 angular positions are required for quantitative image reconstruction. They also showed that non-uniform angular separations provided increased reconstruction accuracy. All images were first corrected for detector dark current offsets.
In order to quantify the relationship between the number of scattering objects and the amount of USAXS, we imaged polymethyl methacrylate (PMMA) microspheres. Microspheres have previously been used as a model for alveoli with good results [13,40,42,43]. Microspheres employed here ranged in size from 90 to 150 µm, which is similar in size to the typical mean linear intercepts for the alveoli of mice (80 µm), rats (100 µm) and humans (210 µm) [44]. Microspheres were imaged in containers of depth 1, 2, 5, 10, 15 and 20 mm. To determine the number of particles, we divided the container depth by the average particle size, multiplied by the approximate packing density of 70%. The known packing density for polydisperse spheres is known to be between 64% (monodisperse) and 75% [45].
As Fig. 3 shows, we found that scattering is independent of particle size. For each material the scattering width was approximately proportional to √ N (actually 2.2627 √ N). That is, the variance of the scatter distribution was approximately proportional to the number of scatterers. These results are consistent with previous studies of von Nardroff [35], Khromova et al. [46] and Khelashvili et al. [47]. At 26 keV, the refractive index decrement, δ, for soft tissue (3.59×10 −7 ) and PMMA (3.94×10 −7 ) are sufficiently close that we can use the scattering relationship for PMMA microspheres (Fig. 3) to estimate the number of air/tissue boundaries in projection through the lung. The relationship to convert scattering from PMMA microspheres to the number of particles (or airways), was found to be: where ω scat is measured in radians. Substitution of the measured scattering from the lungs into Eq. (11), followed by substitution into the airway diameter equation (Eq. (10)) reveals the average airway dimensions.

Measuring emphysema severity from ABI data
In order to measure the projected airway dimensions of mice lungs, we immersed the mice in a water bath to aid in the calculation of the projected lung air volume (see Eq. (9)). Upon performing the ABI experiment, we extracted the scattering, absorption and refraction information using Eqs. (3), (4) and (5), respectively, as shown in Fig. 4. To calculate the lung air volume, we employed the simplifying assumption that the thorax is comprised only three materials, water (for all soft tissues), air (approximated as vacuum), and bone. At 26 keV, µ and δ of water were calculated to be 47.2 m −1 and 3.40×10 −7 , respectively. Equivalent µ and δ values for bone were 372.5 m −1 and 6.09×10 −7 , respectively. From the refraction angle maps (e.g. Figure 4(b)), phase maps were calculated via Eq. (6), numerically integrating the image parallel to the diffraction plane (horizontal in these images). The projected thickness Eqs. (7), (8) and (9) were then used to determine the projected thickness of air in the lungs (see Fig. 5(j)-(l)). The number of airspaces in the lung tissue were found using the scattering data in combination with Eq. (11) (see Fig. 5(g)-(i)). The average projected airway diameter was then computed via Eq. (10) (see Fig. 5(d)-(f)). Figure 5(a)-(c) shows axial slices of the lungs reconstructed from the high resolution CT scan for comparison of the actual airway structure.

Airway measurement accuracy
To quantify the accuracy of the lung structure measurements, we compared small regions (100×100 pixel areas in the centre of the lungs) from ABI reconstructions against CT-based measures of alveolar diameters. For the 80-image ABI data, we plotted the projected airway thickness, number of airways and airway diameter against airway dimensions measured using CT data. We found very poor correlation for the projected airway thickness and number of airways, with R 2 values of just 0.34 and 0.05, respectively (data not shown). Biological variability can greatly affect the number of airways between different animal's lungs. Differences in tissue compliance and airway pressure can also significantly affect the total volume of air in the lungs. Combining these two components to compute the average airway dimensions was expected to greatly reduce this variability. We found that plotting the average airway dimensions measured with ABI against the CT measurements, the correlation significantly increased to R 2 = 0.64, as shown in Fig. 6(a). Although this correlation is not very strong, much of the variability comes from the large uncertainty in the CT granulometry data, as the x-axis error bars reveal. We further discuss the sources of uncertainty in Section 4. Importantly, the gradient of linear best fit in Fig. 6(a) was 1.04, showing a near ideal 1:1 relationship between the ABI size measurements with the CT granulometry results. We see that all of the control animals cluster closely and had small alveoli. However, the treated mice showed a large variety in alveolar dimensions regardless of the treatment amount. This is most likely because the delivery of elastase was inconsistent, with some mice potentially avoiding inhalation through coughing or sneezing.
To test the accuracy of the ABI technique with a reduced amount of input data, we recomputed all parameters using only 11 images (see Fig. 2(b)), as shown in Fig. 6(b). Despite the increased image noise and larger uncertainty bars, for this data we have R 2 = 0.61 with a gradient of 1.18. This shows that the technique is surprisingly robust to noise for quantitative measurement of the airway dimensions, and that relatively few images are required for analysis. The most important consequence of this is the significant reduction in radiation dose and acquisition time, which are both important clinical considerations.

Sensitivity and specificity
To assess the diagnostic accuracy of the ABI data, we measured the sensitivity and the area under the Receiving-Operator Curve (ROC). To do this we used the CT data from control animals to provide an accurate measure of alveolar dimensions in situ, in normal healthy mice. The mean airway dimension of 16 control lung regions was found to be 143.5 µm with a standard deviation (SD) of 27.3 µm, with an upper boundry (mean + 2SD) of 198 µm for healthy airways. Using this value to categorise lungs as healthy or diseased, we determined the number of true (T) and false  Table 1 shows the data for the reconstructions. To compute the ROC curves, we constructed histograms of the average airway diameters measured using ABI (Fig. 7). Setting the airway diameter threshold at 198 µm for emphysema classification, slightly more false positives values were obtained using 11 images than with 80 images, as expected from Table 1. By varying the airway diameter threshold from 75 µm to 325 µm across the histograms and plotting the True Positive Fraction (Sensitivity) against the False Positive Fraction (1-Specificity) gives the ROC curves in Fig. 8. The area under the respective ROC curves is 0.96 and 0.90 for the 80 and 11 ABI image reconstructions, respectively. This shows excellent diagnostic capability in both cases, with diagnostic accuracy getting close to 100% with only a small fraction of the number of images compared to the micro-CT data. Fig. 7. Histograms of the lung airway dimensions measured using ABI for mice categorised as healthy or emphysematous based on CT data airway diameters being smaller or larger than 198 µm, respectively. Measures of the airway dimensions were calculated from images reconstructed using (a) 80 and (b) 11 ABI images.

Discussion
We have developed a technique for measuring the mean length scale of pores, or particles, that cannot be directly resolved by the imaging system. For many applications, the water bath will not be required if the projected thickness of the scattering object is known and there will be no need for the phase retrieval process outlined in Eqs. (4)-(8). Alternatively, the need for a water bath can also be removed by measuring the shape of the object by other methods, such as laser profilometry.
Emphysema results in the destruction of alveoli, leading to increased airway dimensions and an increase in lung gas volumes due to a reduction in elastic recoil [48]. Therefore, we expect that animals with emphysema would have large lung volumes, expanded thoraxes, reduced numbers of airways in projection and enlarged airspaces. Surprisingly, in the lung regions we studied, we only found a weak correlation with lung volume (R 2 = 0.34) and essentially no correlation with airway number (R 2 = 0.05). The lack of correlation with number may arise if the variability in animal size and lung maturity creates a larger variation in airway numbers than the treatment alters. This finding is consistent with other publications using scatter ('dark field') signals for studying emphysema [19][20][21]. Those studies showed weak correlation between disease and scatter, but stronger variation in the transmittance (related to lung air volume) in excised mice lungs. Schleede et al. [19] imaged mice in a water bath and found that taking the ratio of the natural log of both the visibility (V) (dark field or USAXS signal) and the transmittance (I) (i.e., ln[V]/ln[I]), gave a very clear measure of emphysematous lesions in projection. There, the transmittance directly related to lung air volume since their excised lungs were in a water bath. Equations (9) through (11) above show that since their lungs were excised, hence free of bones, this ratio is closely related to the airway dimensions, which explains their strong correlation with airway dimensions. In that case their 'normalised visibility' (ln[V]/ln[I]) is closely related to the inverse of our airway diameter equation (Eq. (10)) above, but is missing the direct conversion between scatter and the number of scattering objects, as per Eq. (11) (airway number). However, for lungs imaged in situ (with the ribcage), and without knowledge of the total projected thickness, the 'normalised visibility' ratio cannot directly measure airway dimensions. Regardless, Yaroshenko et al. [49] showed that the 'normalised visibility' did provide a negative proportional relationship to airway dimensions in a study of lung injury in mice measured in situ. Interestingly, in contrast to their early work, this group later showed direct correlation between scattering (visibility) information and airway enlargement [22,24,49]. Importantly, those studies showed that the 'normalised visibility' always provided a stronger correlation with airway dimensions than did the dark field signal alone.
Our technique for measuring the average projected airway diameter has shown clear correlation with CT data and can accurately diagnose emphysema, as seen in Figs. 6-8. However, the accuracy of the airway dimensions measured from the micro-CT data limits the diagnostic accuracy of the ABI data. Figure 6 shows the large uncertainty of the CT measurements, which are mainly due to the variability of airway dimensions across the lung and to airway pressure variations between animals. It is possible that some lungs lost considerable airway pressure prior to CT acquisition. For example, Fig. 6 shows that two of the control animals had surprisingly small airway diameters according to the CT data. Such variability would have led to incorrect categorisation of disease and reduced the measured sensitivity and specificity of the ABI data.
In future studies, it will be interesting to determine how reducing the spatial resolution will affect the accuracy of the reconstructions. Employing a larger pixel size will enable faster image acquisition with reduced radiation exposure. Although the USAXS signal is robust against resolution, measurements of the phase gradients are less forgiving, so the accuracy of the airway dimension measurements may be reduced.
ABI requires a highly stable mounting system and a monochromatic beam that, when using a beam large enough to image the entire mouse thorax without scanning, must also have very low divergence. These challenges are difficult to overcome without a synchrotron radiation source. However, there are alternative phase contrast techniques that can separate absorption, phase and scatter information that are more suitable for non-synchrotron sources. These include grating interferometry [50], edge-illumination [51], speckle tracking [52] and propagation based phase contrast imaging [53]. In future, we will employ such a technique in an effort to translate this work for potential pre-clinical/clinical application.

Conclusions
This paper demonstrates that ultra-small angle x-ray scattering information can quantitatively measure the mean length scale of porous media. Here we combined phase and absorption information with scatter information to measure the spatial distribution of airway dimensions in projection across the lungs. Employing a mouse model of emphysema, the technique was capable of accurately distinguishing healthy lungs from those with mild and severe emphysema. Using 80 images, the sensitivity and specificity were measured to be 0.80 and 0.89, respectively, with the area under the ROC curve close to ideal at 0.96. Reducing the number of images to 11 slightly decreased the sensitivity to 0.75 and the ROC curve area to 0.90, whilst the specificity remained high at 0.89. Using alternative phase contrast imaging systems that are more stable and noise robust will likely see these results further improve in future studies.

Disclosures
The authors declare no conflicts of interest.