1. INTRODUCTION

Imaging and three-dimensional (3D) visualization of structural information in cellular and subcellular environments will profoundly influence the way we perceive underlying mechanisms of the complex biological world. Many 3D-based imaging techniques have been established for biological studies, from classical confocal microscopy to subdiffraction-limited superresolution techniques ^{[1–4" target="_self" style="display: inline;">–4]}. Most of these 3D imaging techniques require exogenous fluorescent dyes for staining. While fluorescent dyes label specific cellular components with a high degree of reliability, it is also particularly desirable to rapidly acquire high-resolution volumetric data with intrinsic contrast simultaneously concerning the overall morphology ^[5]. In addition, it is sometimes undesirable to stain cells for certain important biomedical research topics, such as stem cell and fertility studies, and the process can become challenging when multiple fluorophores are involved. Therefore, label-free imaging by intrinsic contrast has received significant interest in recent years, spanning various wavelength regimes and modes ^{[6–11" target="_self" style="display: inline;">–11]}. Besides, some alternative bottom-up imaging methods based on various kinds of nanostructures have been reported ^{[12–15" target="_self" style="display: inline;">–15]}. Optical diffraction tomography (ODT) synthesizes information from different angular viewpoints and has been proposed for detecting localized refractive index (RI) variations within a volume ^{[16–20" target="_self" style="display: inline;">–20]}. Many features, including nanofabrication ^[21], partially coherent light microscopy ^[22], dielectrophoretic force, and microfluidics ^[23] have been introduced to ODT to improve the imaging quality and get more control of the samples. Different computational reconstruction techniques like filtered backpropagation algorithms ^[24], numerically extended depth of focus ^[25], and frequency domain interpolation ^[26,27] have been developed and implemented. But the computational complexity from some of these procedures puts a constraint on certain applications in fast and dynamic bioimaging events.

Here, the classical optical tomographic imaging principle is revisited and the 3D relationship between the angular spatial frequencies and the object scattering potential is introduced in a more direct way, leading to the quantitative spatial mapping of subcellular structures with intrinsic values. In addition, we have implemented these newly derived principles in a tomographic microscope with simultaneous fluorescence contrast in a multimodal mechanism. The proposed method unveils interesting subcellular morphologies and brings additional perspectives between function and structure beyond what is available in the single contrast imaging conditions such as either fluorescence or RI-only situations. Finally, to illustrate the accuracy and to demonstrate the quantitative nature of this proposed technique, the Monte Carlo method is employed to simulate the scattered field under the Born diffraction approximation condition. Following the proposed method, the 3D RI distributions of the phantom are well reconstructed, indicating the quantitativeness and accuracy of this method.

2. 3D TOMOGRAPHIC IMAGING PRINCIPLE AND EXPERIMENTAL SETUP

In classical ODT theory, with the Weyl expansion of spherical wave to give two half-spaces, Z+ and Z−, the Fourier diffraction theorem is introduced as ^[28] F˜(kx−ks0x,ky−ks0y,±kz−ks0z)=(kz/2πi)U˜1(s)(kx,ky;z±;s0)exp(∓ikzz±), where F˜ and U˜1(s) denote the Fourier inverse of the scattering potential and scattered field, respectively, s0=(s0x,s0y,s0z) and s=(sx,sy,sz) are the unit vectors of the incident and scattered waves, respectively, and kx=ksx, ky=ksykz=ksz=k(1−sx2−sy2)1/2. In the existing ODT theoretical framework, the analytical object space is broken into two half-spaces, Z+ and Z−. The 3D scattering potential of a specific scatterer, therefore, may be found using the 2D Fourier components of the scattered field, provided that the scattering distance z+ (or z− in Z− half-space) is known ^[28]. However, in an amorphous sample such as a cell, the true scattering distance of an object is often difficult to measure or define precisely.

In our method, instead of expanding the outgoing free-space Green’s function G(r−r′)=exp(ik|r−r′|)/(|r−r′|) in spatial domain, we differ by considering the 3D Fourier transform of the Green’s function.

We employ the Fourier sine transform and write the 3D Fourier transform of outgoing free-space Green’s function as ^[29] where λ is the wavelength, f=|f|=(fx2+fy2+fz2)1/2, G˜H(f) and G˜IN(f) are the 3D Fourier transforms of the imaginary and real parts of the Green’s function, isin(kr)/r and cos(kr)/r, respectively. Note that G˜H(f) and G˜IN(f) are obtained by Fourier transformations of the solutions of homogeneous and inhomogeneous Helmholtz equations; we call them homogeneous and inhomogeneous parts, respectively. Now introducing a specific orientation, fz, the transform of the Green’s function can be written as If we consider the propagating components and ignore the evanescent components, i.e., only counting for the frequencies of fx2+fy2≤1/λ2, which is the case for conventional far-field detection, the inverse Fourier transforms of the two terms on the right-hand side of Eq. (2) are identical for fz>0, while they cancel each other for fz<0 in the transmission configuration (z>0). The transform of the Green’s function becomes a purely out-going, forward-propagating field, represented by a hemispherical shell in reciprocal space. The opposite is true in the reflection configuration (z<0), where the components for fz<0 reinforce. Thus, we can write the Fourier transform of Green’s function in the simple expressions, where M=(1/λ2−fx2−fy2)1/2. It is worth noting that G˜p is a hemisphere in either transmission or reflection configuration, where it includes contributions from both the homogeneous and inhomogeneous parts. But due to the cancellation ^[30], the value of G˜p in either configuration (z>0 or z<0) is twice the value of the homogeneous components.

Return to the classical integral equation of potential scattering. Using the first-order Born approximation and substituting Eq. (3), finally we arrive at the analytical expression of the scattering potential and the scattered field, Essentially, Eq. (4) contains only the propagating waves and relates the 3D light diffraction with a hemispherical cap in the Fourier domain (Fig. 1). We have set up the weakly scattered field U1(s)(r) with the Green’s function in reciprocal space, and this unifies with the theoretical framework of coherent transfer function (CTF) that we presented previously ^[29,30]. Henceforth, the Fourier synthesis of 3D diffraction can be approached using differential geometry solutions ^[31]. Because the entire theoretical framework is derived in 3D, no iterative plane-by-plane “backpropagation” is needed. The effects of diffraction are inherently accounted for when using the spherical support in K space.

Fig. 1. Comparison of the 2D and 3D diffracted field using Fourier relationship. (a) 3D light field of a cell, shown as different planes of wavefront information; (b) 2D light field of a cell extracted from (a); (c) 3D and (d) 2D Fourier spectrum of the light fields. Under Fourier operation, a “slice” of the diffracted field results in a 2D spectrum as shown in (d), whereas the entire volumetric diffracted field results in a 3D “spherical cap,” as shown in (c).

下载图片 查看所有图片

The complex field of the object was first captured and Fourier transform was applied to obtain the object’s 2D spectrum. We then constructed a spherical cap in digital space according to Eq. (4) by “scooping out” from a stack of 2D spectra of the object. The N.A. of the cap is determined by the subtended aperture of the illumination, and the cap thickness is kept to a minimum to reflect the quasi-monochromatic nature of the coherent illumination. For one illuminating condition, we position the “spherical cap” as the zero-angle reference cap. The key is to synthesize the complex spectrum of the scattering potential F˜(f) correctly into K space under various angular conditions. We used our previously derived CTF ^[31] to guide the correct positioning of each angular spherical cap into the reciprocal space, noting that all the zeroth-order diffraction points coincide (Fig. 2). All the angles relative to the initial zero angle were calculated, and the 3D movements of the caps were calculated accordingly. Examples of mapping the three angles are shown in Fig. 2. After the spectra of the scattering potential F˜(f) under different illuminations were constructed in the 3D conjugate K space, the scattering potential F(r) of the object was inversely retrieved by a single Fourier transform. The 3D RI distribution then can be calculated: where km is the wave vector in the background medium of RI nm, U˜j(s)(fj) is the 3D Fourier transform of the scattered field under incident direction sj, δj=±fzjδ[fzj∓(1/λ2−fxj2−fyj2)1/2] is the hemispherical cap for transmission or reflection configurations, and F3−1 represents the 3D inverse Fourier transform. Subscript j denotes the different illumination angles. The proposed tomography is valid under the condition of weak absorption, as is the case for visible light in an unstained biological sample. A difference with computed tomography (CT) ^[32] is that, here, due to diffraction, the information from one monochromatic image must be transcribed to the conjugate K space under the spatial support of a thin spherical cap instead of a plane. In addition, the method here has enhanced accuracy because phase unwrapping is not needed.

Fig. 2. Complex spectra synthesis of the scattering potential using three example angles. (a) A cell is subject to a multiangle plane-wave illumination with a maximum subtended half-angle θ; hence, N.A.=sinθ. The 3D angular spectra under different incident angles can be acquired (only three angles shown). (b) 2D side view of the angular compensation method introduced during the synthesis of the 3D angular spectrum (only three angles shown). The green arc represents the center reference cap, the red and blue arcs represent the “compensated” caps from incident angles θ and −θ, respectively. ΔFz is the amount of shift in axial direction for the incident angle θ. (c) 3D view of a “synthesized” angular spectrum under three example angles. In (a), the center illumination and two maximum subtended angles illuminations formed three individual spherical caps in the Fourier domain as the 3D angular spectra. The respective 3D angular spectra are then translated and shifted accordingly, as shown in (b), so that angular information introduced in the incident beams is compensated. A synthesized angular spectrum for the three example angles is shown in (c), where all the 0th-order components coincide at a point.

下载图片 查看所有图片

It is noted that although the Fourier mapping procedure from 2D to 3D has been mentioned before ^[18], we uniquely established the analytical expression in Eq. (4) and elucidated how we arrived at this important relationship by canceling the inhomogeneous parts. In fact, the proposed method applies universally to any type of coherent tomography, such as X-rays or electron waves in the case of weak scattering.

The feasibility and robustness of the method for cellular imaging are tested firsthand with a benchtop cellular tomography system with varied illumination scan patterns. It is based on a scanning Mach–Zehnder interferometer (Fig. 3). The sample is illuminated by a collimated beam with different incident angles introduced by the galvo mirror. The transmitted light that has been collected by the microscope objective interferes with the reference beam in the off-axis configuration to produce the image holograms. The intensities of the holograms under all different illuminations are recorded by a high-speed scientific complementary metal-oxide-semiconductor (sCMOS) camera. In addition, we combine our technique with epifluorescence microscopy (Fig. 3) as a dual-contrast modality. A data acquisition (DAQ) card (USB-6341, National Instrument) is employed to control the whole system by switching between tomography mode and fluorescence mode rapidly for simultaneous in-suit dual-contrast imaging.

Fig. 3. Schematic diagram of the multimodal setup combining tomography with fluorescence. The tomography module is based on an off-axis digital holographic microscope setup in transmission. Different illumination angles are introduced by galvo mirror. The inset illustrates the tilted angle illumination situation. The fluorescence module is based on the epifluorescence microscopy. BS, beam splitter; GM, galvo mirror; SL, scanning lens; C, condenser; O, objective lens; TL, tube lens; Ex, excitation filter; Em, emission filter; DM, dichroic mirror.

下载图片 查看所有图片

3. RESULTS AND DISCUSSION

3.1 A. 3D Mappings of the RI Distributions

Various human cells are cultured and 3D mappings of the RI distributions are recovered using the above procedure. We first imaged human primary buccal epithelial cells, where, in Fig. 4(a), the reconstructed 3D structural information is encoded with pseudo-colors representing different regions of similar RI values. The reconstruction is completely 3D, and we are able to generate a rotating view of the sample object (Visualization 1). The shape of the nucleus can be clearly delineated. Cytoplasmic components are divided into two layers through segmentation of the RIs. This reveals cytoskeletal remodeling at the edge of a cell, which is difficult to observe for unlabeled cells in conventional light microscopes. We also examined neurons. In Fig. 4(c), the reconstructed neuron cell reveals clustering of organelles having similar biomasses (Visualization 2). The higher RI value [represented with red color in Fig. 4(c)] exhibits a hollow shape that cannot be observed in traditional 2D phase imaging [Fig. 4(d)]. This is a good example of how this method isolates variations in thickness and RI values.

Fig. 4. Results for 3D reconstructions. (a) Buccal cell and (c) neuron, with internal structures delineated by groupings of refractive indices. The rotating views of the buccal cells and neuron can be found in Visualization 1 and Visualization 2, respectively. Optical phase distribution of the same (b) buccal cell (in dashed box) and (d) neuron calculated from the interferometric measurements. The units for the color bars in (b) and (d) are radians. Part of the buccal cell in (a) is removed for clear visualization of the subcellular structures. The magnification is 20× with an N.A. of 0.45.

下载图片 查看所有图片

3.2 B. Comparisons of Colocalized Fluorescence Distributions

Because fluorescent dyes can label specific cellular components with a high degree of reliability, the fluorescence imaging can be employed as a reference for label-free imaging techniques ^{[33–37" target="_self" style="display: inline;">–37]}. We demonstrate the dual-contrast modality (Fig. 3) by measuring and visualizing the volumetric structural information with its in-suit fluorescence signals simultaneously in live human immune cells (T cells). The colocalized RI and fluorescence-labeled distributions of nuclear components during the apoptosis (a form of programmed cell death) procedure are obtained. Specifically, human Jurkat T cells were stained with Hoechst 33342 dye to monitor the localization of nuclear contents before imaging or induction of apoptosis. 3D mappings of the RI distributions were recovered [Fig. 5(a)] and distributions of wide-field fluorescence images were obtained simultaneously [Fig. 5(c)]. The 3D RI image shows more clearly than its fluorescent counterpart that the nuclear contents are condensed and have the hollow shape common to apoptosis ^[38] (Visualization 3). Chromatin condensation can be easily observed through a 2D slice of RI distributions along the xy plane [Fig. 5(b)]. The nuclear contents highlighted with blue contour lines in Fig. 5(b) are based on the intensity distributions of fluorescence [Fig. 5(c)]. By overlapping the fluorescence image with the reconstructed RI values, the histogram distance is calculated, which provides the similarity of the two distribution areas. The histogram distance for the overlapping parts in Figs. 5(b) and 5(c) is 0.8548. Our method is not only compatible with conventional fluorescence microscopy, but also enhances the latter by revealing previously unseen morphological features and mechanisms. Like Z-stack fluorescence imaging, we can show detailed 3D interior cell structures [Fig. 5(a)]. Furthermore, the method can potentially provide the entire spatial information inside the cell and its subcellular structures in real time. Additional comparisons between RI values and fluorescence intensities for Jurkat T cells under different conditions are shown in Fig. 6. The histogram distances calculated for the overlapped parts in the first, the second, and the third column of Fig. 6 are 0.8815, 0.8772, and 0.8003, respectively, indicating pretty good similarity. The materials and methods for T cells used here can be found in Appendix A.

Fig. 5. Structures reconstructed using RI spatial mappings in apoptotic T cells are shown to be compatible with wide-field fluorescence counterparts. (a) Cells with internal structures delineated from groupings of RI values (top view). Inset gives an opened-up view of the cell through the central z plane of the nucleus. The rotating view of the cells can be found in Visualization 3. (b) 2D slice of RI distributions along xy plane; (c) wide-field fluorescence image of the same cells. The nuclear contents highlighted with blue contour lines in (b) are based on the intensity distributions of fluorescence in (c). The cell nuclei were stained with Hoechst dye. The value of RI in the range of 1.388–1.419 (pseudo-colored red orange) matches very well with what the fluorescent signal is showing, indicating that our technique is able to differentiate functional substructures within the cell. The magnification is 40× with an N.A. of 0.6.

下载图片 查看所有图片

Fig. 6. Distributions of RI values (row 1) and fluorescence intensities (row 2) for Jurkat T cells under different conditions (scale bar, 10 μm). (a), (b) Viable cell; (c), (d) UV light-induced apoptotic cells with trovafloxacin treatment; (e), (f) necrotic cell. The positions of stained DNA, where many proteins are localized, are highlighted with blue contour lines in RI mappings based on fluorescence imaging. The magnification is 40× with an N.A. of 0.6.

下载图片 查看所有图片

3.3 C. Validate the Method Using Monte Carlo Simulation of the Volumetric Scattered Field

For the investigations of naturally 3D complex objects like cells, the validation and accuracy of these methods are crucial. Monte Carlo methods are a broad class of computational algorithms that are mainly used for optimization, numerical integration, and generating draws from a probability distribution ^[39] and have been employed in many bioimaging studies such as CT ^[40], bioluminescence tomography ^[41], optical coherence tomography ^[42], and fluorescence molecular tomography ^[43]. As important aids, Monte Carlo methods are used to simulate the optical models that can make the comparison with measurements or reconstructed results or can optimize the instrumental parameters.

To validate the accuracy and to demonstrate the quantitative nature of this proposed technique, we first simulated tomographic images of a 3D scattering object using Monte Carlo methods ^[39] and then reconstructed the phantom using the above principles. The configuration of the simulation is shown in Fig. 7(a), a 3D semispherical shaped object immersed in a 3D environment resembling typical cellular physiological conditions. One-dimensional (1D) scanning illumination is applied on the object to introduce different incident angles. Volumetric diffraction within the computational boundary is derived from a numerical implementation of Born’s diffraction integral. Once diffraction conditions are established, a detected 2D image of the phantom was constructed by taking the intensity of each simulated optical field of the plane just exiting the sample. Then a 2D Fourier transform was performed on the simulated image before the spatial frequencies being mapped to the 3D support, as guided by Eq. (4). After a series of 2D images under different illumination angles were simulated and the corresponding spherical Fourier caps were obtained through the same processes, the complex spectra were synthesized into K space under various angular conditions (as shown in Fig. 2 for just three angles). Finally, an inverse 3D Fourier transform was performed, and the correct RI values were recovered as well as the shape of the phantom [Fig. 7(c)], indicating the accuracy and the quantitativeness of the proposed method.

Fig. 7. Monte Carlo simulation of a phantom object volumetric diffraction and its quantitative reconstruction. (a) The configuration of the Monte Carlo simulation; (b) 3D view of the semispherical-shaped phantom object; (c) reconstructed RI distribution of the object in central yz plane; (d) cross section of the CTF (central fxfz plane), which is synthesized by 21 spherical Fourier caps.

下载图片 查看所有图片

To calculate the 3D scattered field of a weak scattering object, a numerical approach was devised using Monte Carlo principles ^[39]. Then the first-order Born diffraction integral ^[28] is rewritten as Us(r)=1nfv∑i=1nF(ri)eiks0·rieik|r−ri||r−ri|,(6)where fv=1/V if r∈V and 0 if r∉V. V is the scattering volume, and n is the number of random events in the numerical procedure. For higher-order 3D Born scattered fields, Ums(r)=1nfv∑i=1nF(ri)Um−1(ri)eik|r−ri||r−ri|,(7)where m is the order number of 3D Born scattered fields. For calculation here, the dimension of the object is 2.5 μm in diameter in the top flat region and the depth is 2 μm [see 3D view of the object in Fig. 7(b)]. The immersion medium is set with an RI of 1.3336 and an object region of 1.3736 to satisfy Born’s first-order approximation condition (the maximum optical path is approximately 0.8 rad for a 633 nm wavelength). The illumination angle ranges from 0° to 30°, corresponding to an objective with an N.A. up to 0.5, which is similar to the realistic experimental conditions. The steps for Monte Carlo computation of Born approximation in 3D scattered fields can be found in Appendix A. We calculated the first-order Born scattered fields using a graphics processing unit (GPU); the MATLAB code is included in Appendix A as well. For higher-order calculations, one can follow the steps, but longer computation time is required. Following our proposed methodology, the conjugated K space was filled with 21 individual spherical caps under 1D scan configuration. Figure 7(d) shows the central fxfz plane of K space, which indicates the cross section of the CTF. The RI distribution was retrieved, and the RI value distribution in the central yz plane of the object is shown in Fig. 7(c). We count the reconstructed RI value higher than 1.3450 in the central yz plane as the reconstructed object, which is in good accordance with the shape of the phantom. The average RI value of the reconstructed phantom object is 1.3577±0.0082, which is compatible with the setting value, considering the low sampling rate of the object. The correctly reconstructed RI values and the shape indicate the accuracy and the quantitativeness of the proposed method. It is worth noting that here only 21 scan angles are used for reconstruction, which means the proposed technique can be used when only a limited amount of data is acquired for quick evaluation of the subcellular conditions, particularly for rapid live cell studies.

4. CONCLUSION

In summary, we have demonstrated an accurate tomographic cellular imaging method that elucidates 3D cellular and subcellular structures with both intrinsic and exogenous contrast. This is significant because both specificity and background context are available now as a combined modality, more accurately revealing the relationship between the function and structural information, as shown in the Jurcat T cell example. The reconstruction of Monte Carlo simulations of the 3D phantom validated the quantitativeness and accuracy of this proposed method. A cell is a complex 3D object containing many components and organelles with RI values that vary across these minute structures. Each RI quantifies the speed of light in the media and is related to the regional concentration of protein (sometimes called biomass). Encoding similar refractive indices is, therefore, an intuitive way to map out regions sharing similar biophysical content—a useful perspective that could greatly complement fluorescence-based imaging. Label-free and fluorescence information can be cross-linked so that cellular events can be more accurately visualized. Based entirely on Fourier principles, the proposed technique can be rapidly implemented in Fourier-enabled signal processing tools or GPUs for observation of quick and dynamic events in cells. Requiring minimum sample preparations and working both with adherent and nonadherent cells, the proposed dual-contrast imaging modality is a powerful and versatile tool that can be applied in all kinds of bioimaging conditions for new discoveries.

5 Acknowledgment

Acknowledgment. This work was largely supported by the start-up grant of S. S. K. from La Trobe Institute for Molecular Science (LIMS). S. S. K. acknowledges the Partnership Grant from La Trobe Research Focus Area (RFA) of Understanding Diseases. The authors declare that there are no conflicts of interest related to this paper.