Phase-contrast zoom tomography reveals precise locations of macrophages in mouse lungs

We have performed x-ray phase-contrast tomography on mouse lung tissue. Using a divergent x-ray beam generated by nanoscale focusing, we used zoom tomography to produce three-dimensional reconstructions with selectable magnification, resolution, and field of view. Thus, macroscopic tissue samples extending over several mm can be studied in sub-cellular-level structural detail. The zoom capability and, in particular, the high dose efficiency are enabled by the near-perfect exit wavefront of an optimized x-ray waveguide channel. In combination with suitable phase-retrieval algorithms, challenging radiation-sensitive and low-contrast samples can be reconstructed with minimal artefacts. The dose efficiency of the method is demonstrated by the reconstruction of living macrophages both with and without phagocytized contrast agents. We also used zoom tomography to visualize barium-labelled macrophages in the context of morphological structures in asthmatic and healthy mouse lung tissue one day after intratracheal application. The three-dimensional reconstructions showed that the macrophages predominantly localized to the alveoli, but they were also found in bronchial walls, indicating that these cells might be able to migrate from the lumen of the bronchi through the epithelium.


Phase-retrieval based on the TIE
The most commonly used phase-retrieval techniques are based on the Transport of Intensity Equation (TIE) which is well suited for phase-contrast images recorded at short propagation distance or correspondingly large Fresnel numbers. If only a single (defocus) image is recorded, additional assumptions on the objects are necessary for reconstruction. The assumption of a fixed proportionality between phase shift and attenuation coefficient, also known as phase attenuation duality [3] leads to the following reconstruction formula ID 19 setting P10, large FOV P10, zoom single cells TIE α 1.3 · 10 −3 6 · 10 −8 7 · 10 −9 -CTF α0 -6 · 10 −5 5 · 10 −4 -CTF α1 0.01 0.1 0.05 10 −5 CTF δ/β 200 --10 Supplementary Table I. Regularization parameters used for the phase-retrieval. For the CTF-based reconstruction, either the assumption of a pure phase object is used together with a regularization α0 for the low frequencies, or the assumption of a phase attenuation duality is made for which a suitable δ/β ratio is chosen based on visual inspection.
Here |ξ| = ξ 2 x + ξ 2 y is the modulus of the coordinate vector in reciprocal space, κ = δ/β is the ratio of decrement and imaginary part of the refractive index, ∆z the propagation distance, k = 2π/λ the wavenumber and I the measured intensity divided by the (uniformly) illuminating intensity I 0 . Similar to the above equation, but easier to derive is the reconstruction formula for a pure phase object in combination with a regularization for the low spatial frequencies [4] The regularization parameter α can physically be justified by the (neglected) residual absorption which manifests itself at low frequencies. This absorption contradicts the assumption of a pure phase object creating a necessity and providing a physical basis for the regularization parameter α. Since α can thus be identified with the additive term in the denominator of Eq. (2), the structure in Fourier space of the two reconstruction formula is in fact the same. For the comparison between TIE and CTF based reconstructions following below, we used a Tikhonov type regularization for the pure phase-TIE For Eq. (3) and Eq. (4) the regularization parameter was chosen as follows: If the parameter α is too small, reconstruction results in blurred images as the low frequencies are amplified too much. Choosing α too large, the reconstruction qualitatively reproduces the measured phase-contrast image as it becomes a simple global scaling for all spatial frequencies. Based on these two limiting cases the regularization parameter was chosen based on visual inspection. The resulting parameters are shown in Tab. I.

CTF based phase-retrieval
Based on the propagation formula in Eq. (1) and the assumption of a slowly varying phase and weak absorption, the so called contrast transfer function (CTF) can be derived [5] This equation is an expression for the Fourier transform of the intensity I, where δ D (k x , k y ) is the Dirac delta distribution, χ = z/(2k) ξ 2 x + ξ 2 y , φ and µ the projected phase and amplitude distribution and thẽ denotes the Fourier transform. By introducing an error metric for n images S = nĨ exp,n −Ĩ(k x , k y , z), using Eq. (5) forĨ and minimizing S with respect toφ one can derive a reconstruction formula, e.g. for a pure phase object (µ = 0)φ = nĨ exp,n sin χ n n sin 2 χ n + α .
Here, a frequency dependent regularization parameter was added to the denominator to compensate zeros of the sin 2 χ n expression. The regularization parameter was chosen such, that for high and low frequencies Supplementary Figure 1. Phase-retrieval in the direct contrast regime (ID22 experiment): (a) empty beam-corrected intensity distribution (b) "conventional" TIE phase-retrieval (c) CTF phase-retrieval showing slightly improved resolution. Scalebars denote 100 µm a value of α 1 and α 0 is applied, respectively, with an error function transition from α 0 to α 1 at the first maximum of the CTF. α 0 corrects the singularities at low frequencies that occur mainly due to residual absorption and beam inhomogeneities, while α 1 corrects the zero crossings of the CTF. The parameters were chosen based on optical inspection of the resulting image and the resulting values can be seen in Tab. I. The use of multiple distances can improve the reconstruction, as for different Fresnel numbers zero crossings occur at different spatial frequencies. An easy way to achieve different Fresnel numbers is to move the sample in the defocused beam, which leads to slightly different magnifications in a cone beam geometry. Here, we rescaled the images based on the experimentally determined distances z 01 and z 02 by using an bicubic interpolation with antialiasing implemented in Matlab. To register images in different planes, each image was reconstructed assuming that only one distance was measured and these rough estimates were used to determine the shifts of the individual images to each other by using an algorithm for sub-pixel registration based on correlations [6]. The shifts were used to correct the rescaled intensity distributions with the different defocus distance by applying a phase factor in fourier space. Subsequently, the final phase was retrieved by utilizing Eq. (6).
The use of the CTF based phase-retrieval improves the resolution compared to the TIE phase-retrieval, especially for small Fresnel numbers as can be seen in Fig. 1 and Fig. 2. As CTF-reconstructions yield better results and contrary to current practice, are also well suited for data recorded at a single defocus distance (single image acquisitions), just as the TIE approach, we found no reason to use TIE based reconstruction, unless the assumption of a weak object is violated.

Iterative phase-retrieval
The derivation of the CTF assumes a weakly absorbing object and a slowly varying phase. Especially in the case of a single cell fed with barium particles, the assumptions imposed in the CTF based phaseretrieval are invalid and result in halo-like artifacts around the dense particles which can be seen in the reconstruction in Fig. 3 (c) and (d) which is based on measuring the intensity distribution in eight different planes. An alternative reconstruction method is based on iterative propagation between object and detection planes, with sequential projection onto the different constraint sets (modulus constraint corresponding to the square root of the measured intensities, object constraints). Since the empty beam-corrected intensities still showed some residual low frequency variations which do not stem from the object, these artifacts had to be suppressed as they cause inconsistencies between the individual measurement planes. A suitable correction was implemented by calculating image profiles in areas devoid of Barium labeled cells and dividing the image by a combination of these profiles along with a Fourier filter which suppresses the very low frequencies. An exemplified empty beam-corrected image is shown in Fig. 3 (a) before correction and in Fig. 3  Starting with a random initial guess, the wave function ψ was propagated to a measurement plane and the wavefunction was updated such that it satisfied the measured intensities and kept the current guess for the phases. The updated wave function was subsequently propagated back to the object plane, where positivity of the refractive index decrement δ and β are enforced by setting |ψ| ≤ 1 and arg (ψ) < 0 as in [7]. This step was repeated for each measurement plane, all of which have slightly different Fresnel numbers, and the total outer loop was repeated 5000 times. The scheme can be considered a multi-plane Gerchberg-Saxton algorithm, but with the modification that the modulus constraint in the object plane is relaxed to a positivity constraint, and is thus applicable to a much wider range of samples. The resulting phase is shown in Fig. 3 (e) and (f) for different contrast settings.
As the results show, the halo-like artifacts introduced by the CTF reconstruction are not present in the iterative reconstruction. However, the need of a very homogeneous background makes it difficult to apply the iterative reconstruction for a complete tomography data set (in which empty beam images are not available for all projection angles). This difficulty was particularly prominent for region-of-interest (ROI) tomography.

The empty beam-correction
Each recorded image was divided by an empty beam image to reduce the influence of the illumination profile. It is known, that this empty beam-correction is mathematically incorrect and thus small artifacts may be introduced in the intensity distributions [8] which are amplified by phase-retrieval techniques. In the experiments this can be especially seen in the case of a pure KB beam illumination, which has many high frequency artifacts (compare Fig. (1) (c) to Fig. (2) (e) and (f)). Approaches to circumvent the flawed empty beam-division by simultaneously reconstructing the illuminating beam from the measurement along with the object have been proposed [9,10], but require additional data recorded by scanning. Here we used the fact, that x-ray waveguides provide a very clean illuminating beam, consisting only of a few optical modes. Thus, the illuminating beam is much more uniform and higher resolutions can be achieved although empty beam-correction is performed.

Tomographic reconstruction
A standard filtered backprojection algorithm implemented in Matlab was used for tomographic reconstruction of the P10 data. Cone beam effects typically encountered for divergent beam were negligible due to the small divergence angle. During a tomographic scan, several empty images were recorded to compensate beam instabilities. Projections for the same angle were aligned before and after recording a new empty beam image using the algorithm presented in [6]. As all samples were larger than the FOV the 3D reconstructions are results of local tomography. Thus, no additional alignment of different projections was performed for the P10 data as especially for the high zoom setting only a small region of interest (ROI) was regarded and features outside this ROI would have disturbed the alignment. The projections obtained at the ID22 beamline were processed with the online tools available during the beamtime and reconstructed using PyHST.

SUPPLEMENTARY METHODS 3: FOURIER SHELL CORRELATION
The Fourier shell correlation (FSC) is a measure for the resolution of a 3D dataset. To calculate the FSC the reconstruction is splitted into two datasets by taking every second angle. The correlation is then calculated as follows where F 1 and F 2 denote the splitted datasets in Fourier space and < · > denotes the 3D angular average. The resolution can then be estimated by an intersection of the FSC with a threshold curve. Different threshold curves are discussed in detail in [11]. The 3σ threshold indicates whether significant information above the noise level has been recorded and the 1/2 bit threshold curve indicates the level which is sufficient for interpretation of the data. The resulting FSC curves together with the used threshold curves are shown in Fig. 4. With the effective voxel size of the measurements the resolution of the large FOV dataset could be estimated at 648 nm half period and at 173 nm for the zoomed dataset respectively. Figure 5 shows an overview of the different 3D structures obtained from samples measured at the ID 22 beamline. The tomograms were obtained from five different mice where blank denotes a mouse without any instilled macrophages and without asthma. Asthma and control denotes that barium labeled macrophages were given to a mouse with and without asthma, respectively. The renderings show orthoslices through the structure, segmented barium particles which indicate the position of macrophages (green) and selected labeled structures like blood vessels (purple), bronchial tubes (red) and alveolarepithelium (yellow). To label different structures in the reconstructed 3D volumes, automatic threshold based segmentation was used for barium sulfate clusters. Alveolar walls were segmented with a "magic wand" tool followed by application of a watershed algorithm, both available in Avizo (VSG studios). Bronchial tubes and blood vessels were labeled with the contrast based "blow tool" in many slices with subsequent interpolation between the slices.