## Abstract

Optical diffraction tomography (ODT) reconstructs a sample’s volumetric refractive index (RI) to create high-contrast, quantitative 3D visualizations of biological samples. However, standard implementations of ODT use interferometric systems, and so are sensitive to phase instabilities, complex mechanical design, and coherent noise. Furthermore, their reconstruction framework is typically limited to weakly scattering samples, and thus excludes a whole class of multiple-scattering samples. Here, we implement a new 3D RI microscopy technique that utilizes a computational multi-slice beam propagation method to invert the optical scattering process and reconstruct high-resolution ($\mathrm{NA}>1.0$) 3D RI distributions of multiple-scattering samples. The method acquires intensity-only measurements from different illumination angles and then solves a nonlinear optimization problem to recover the sample’s 3D RI distribution. We experimentally demonstrate the reconstruction of samples with varying amounts of multiple-scattering: a 3T3 fibroblast cell, a cluster of *C. elegans* embryos, and a whole *C. elegans* worm, with lateral and axial resolutions of $\le 240\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and $\le 900\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$, respectively. The results of this work lays groundwork for future studies into using optical wavelengths to probe 3D RI distributions of highly scattering biological organisms.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

Fluorescent imaging has enabled stunning visualizations of biological processes at a variety of size scales and resolutions, for studies of gene expression, protein interactions, intracellular dynamics, etc. [1–4]. However, the fluorescent techniques require exogenous biological labels, and thus do not directly give endogenous information about a sample’s biological structure.

Optical diffraction tomography (ODT) also targets 3D biological imaging. In contrast to fluorescent methods, ODT avoids the use of exogenous biological labels and instead utilizes the intrinsic optical variation within a sample to reconstruct its 3D refractive-index (RI) distribution [5–11]. Hence, ODT avoids some of fluorescent imaging’s main drawbacks, such as photobleaching, slow acquisition speed, low signal-to-noise (SNR) ratio, and complex sample-preparation protocol. Furthermore, RI imaging enables the examination of the structural, mechanical, and biochemical properties of a sample, which are important for studies in morphology, mass, shear stiffness, and spectroscopy [9,12–15].

Standard implementations of ODT use either a rotating sample or a scanning laser beam to capture the angle-specific scattering arising from the sample [5,7,16–18]. Under the assumption of weak scattering (i.e., 1st Born or Rytov approximations), 2D electric-field measurements directly yield information about the sample’s 3D scattering potential [19–21]. Standard ODT reconstruction algorithms utilize the Fourier diffraction theorem to project the information contained in each electric-field measurement onto spherical shells (i.e., Ewald surfaces) in the 3D Fourier space of the sample’s scattering potential [22,23]. After sufficient 3D coverage of the sample’s scattering potential has been achieved, 3D RI distributions can be reconstructed [7,24,25]. This reconstruction framework is efficient and fast and has found great success in visualizing individual unlabeled cells with quantitative refractive indexes.

Unfortunately, standard ODT has two main drawbacks that limit its utility in biological research: (1) it requires sensitive electric-field measurements via laser illumination and interferometry, which associate with coherent noise, phase instabilities, and complex system-alignment protocols, and (2) the Fourier diffraction theorem relies on the sample being weakly scattering, which limits applications to predominantly individual cells. This excludes a whole class of biologically-relevant samples, such as densely packed clusters of cells or multicellular organisms, which are optically transparent but multiple-scattering.

Advances in computational phase retrieval have introduced alternative techniques that utilize optimization [26–34] for 3D sample reconstruction. Unlike the standard ODT framework, which relies on an analytical inversion of the scattering process via the Fourier diffraction theorem [7,22], or weak-phase approximation [35,36], such computational techniques iteratively refine an initial estimate of the 3D sample towards a solution that minimizes the difference between the observed measurements and the expected measurements, as predicted by a physical forward model. While analytical inversion for multiple-scattering is often not feasible, forward models based on multiple-scattering are feasible. Hence, for scenarios where the sample is not weakly scattering and there exist no robust analytical inverter, optimization-based iterative solvers are a promising alternative.

Recent work by Tian and Waller [34] introduced a computational technique to invert multiple-scattering. Angled illuminations of the sample, incorporated into the multi-slice beam-propagation (MSBP) forward model, were used to reconstruct the sample’s 3D phase distribution from intensity measurements, without physical sample scanning. To do so, however, the reconstruction framework utilized a 2D gradient-descent protocol to iteratively invert 3D scattering. This simplification accumulates error as the axial sampling density is increased, making the framework untenable for high numerical aperture (NA) 3D RI imaging. Subsequent work by Kamilov *et al.* [33] introduced a framework, based on MSBP and sparse regularization, to rigorously invert high-NA 3D multiple-scattering from electric-field measurements captured by a laser-based interferometer. However, such interferometer systems typically associate with mechanical instabilities and coherent artifacts, which cannot be easily accounted for in the reconstruction model. Thus, this can give rise to instabilities or over-regularization in the nonlinear iterative solver, degrading effective resolution.

In this work, we present a new 3D RI microscopy technique, also based on the multi-slice beam-propagation (MSBP) forward model, that uses intensity-only measurements to iteratively reconstruct 3D RI of multiple-scattering biological samples. Similar to the works described above, our technique illuminates the sample from different angles in order to encode 3D information into 2D measurements and then recovers the sample’s 3D RI distribution. Beyond the work by Tian and Waller [34], we present a rigorous nonlinear optimization protocol that performs robustly at high resolutions ($\mathrm{NA}>1.0$) and enables full axial sampling of the system’s diffraction limited point-spread-function. Even better resolution may be achieved due to the multiple-scattering within the sample [34,37,38], but is difficult to specifically characterize. To the best of our knowledge, this work presents the first demonstration of high-resolution (250 nm lateral and 900 nm axial resolution) inversion of 3D biological multiple-scattering from intensity-only measurements.

## 2. THEORY

#### A. Multi-Slice Beam-Propagation Forward Model

Beam-propagations methods generally use a set of initial conditions to calculate the electric-field at a separate location in space (or time) [39,40]. Multi-slice beam-propagation (MSBP) methods specifically approximate 3D objects as a series of thin layers, where light propagation through the sample is modeled via sequential layer-to-layer propagation of the electric-field [33,34,41]. Mathematically, this can be recursively written as

The exit electric-field, ${y}_{N}(\mathit{r})$, accounts for the accumulation of the diffraction and multiple-scattering processes that occurred during optical propagation through the sample and contains information about the sample’s 3D structure. The final electric-field and intensity distributions at the image plane are

#### B. Inverse Problem Formulation

Multiple measurements of the object are captured at varying illumination angles. We denote ${y}_{0}^{\ell}(\mathit{r})=\mathrm{exp}(j{\mathit{k}}_{0}^{\ell}\xb7\mathit{r})$ for $\ell =1,2,\dots ,L$ as a set of $L$ planar electric-fields incident onto the first layer of the sample, at various angles set by their wave-vectors ${\mathit{k}}_{0}^{\ell}$. The corresponding intensity measurements are denoted by ${I}^{\ell}(\mathit{r})$.

The reconstruction framework can be formulated as a least-squares minimization that seeks to estimate the sample’s 3D RI by minimizing the difference between the measured amplitude (square-root of intensity [42]) and those expected via the forward model,

#### C. Reconstruction Framework

Our proposed framework to minimize Eq. (4) and reconstruct $\hat{n}({\mathit{r}}_{3\mathrm{D}})$ draws from previous MSBP works by Kamilov *et al.* [33] and Tian and Waller [34], and utilizes an iterative scheme to interleave back-propagation and gradient-descent steps into each iteration. We describe below the steps taken for the iterative reconstruction.

- (1) Initialize $N$-layer reconstruction volume with constant refractive index, ${n}_{m}$. This will serve as the initial estimate of $n({\mathit{r}}_{3\mathrm{D}})$. Then, initialize the iteration index to $d=0$.
- (2) To start a new iteration, increment the iteration index, $d\leftarrow d+1$, and initialize the per-iteration cost function, $c(d)=0$.
- (3) Randomly choose (without replacement) an illumination angle, with electric-field ${y}_{0}^{\ell}(\mathit{r})=\mathrm{exp}(j{\mathit{k}}_{0}^{\ell}\xb7\mathit{r})$ and raw intensity measurement ${I}^{\ell}(\mathit{r})$, from the complete set of $\ell =1,2,\dots ,L$.
- (4) For the chosen illumination field ${y}_{0}^{\ell}(\mathit{r})$, use Eqs. (1) and (4) to recursively compute and store the electric-field at each layer of the reconstruction volume $\{{y}_{k}^{\ell}(\mathit{r})|k=1,2,\dots ,N\}$ and the estimate of the final intensity measurement as predicted by the forward model, ${\mathcal{G}}^{\ell}\{n({\mathit{r}}_{3\mathrm{D}})\}$.
- (5) Increment the cost function for the current iteration, $c(d)\leftarrow c(d)+\sum _{\mathit{r}}|\sqrt{{I}^{\ell}(\mathit{r})}-|{\mathcal{G}}^{\ell}\{n({\mathit{r}}_{3\mathrm{D}})\}{||}^{2}$
- (6) Initialize a residual term denoted by ${q}_{N+1}^{\ell}(\mathit{r})$. The variable ${q}_{0}(\mathit{r})$ below is used only for notational simplicity$${q}_{0}(\mathit{r})=\mathrm{exp}(j\angle {\mathcal{G}}^{\ell}\{n({\mathit{r}}_{3\mathrm{D}})\})\xb7(|{\mathcal{G}}^{\ell}\{n({\mathit{r}}_{3\mathrm{D}})\}|-\sqrt{{I}^{\ell}(\mathit{r})})\phantom{\rule{0ex}{0ex}}{q}_{N+1}^{\ell}={\mathcal{P}}_{\hat{z}}\{{\mathcal{F}}^{-1}\{\overline{p(\mathit{k})}\xb7\mathcal{F}\{{q}_{0}(\mathit{r})\}\}\}.$$
- (7) For each layer of the reconstruction volume occupied by the sample, compute the back-propagation term ${s}_{k}^{\ell}(\mathit{r})$ by recursively propagating
*backwards*(i.e., $k=N,(N-1),\dots ,2,1$),$${s}_{k}^{\ell}(\mathit{r})=(-j\frac{2\pi \mathrm{\Delta}z}{\lambda})\xb7\overline{{t}_{k}(\mathit{r})}\xb7\overline{{\mathcal{P}}_{\mathrm{\Delta}z}\{{y}_{k-1}^{\ell}(\mathit{r})\}}\xb7{q}_{k+1}^{\ell}(\mathit{r}),$$$${q}_{k}^{\ell}(\mathit{r})={\mathcal{P}}_{-\mathrm{\Delta}z}\{\overline{{t}_{k}(\mathit{r})}\xb7{q}_{k+1}^{\ell}(\mathit{r})\},$$where the notation $\overline{g(\mathit{r})}$ designates the complex conjugate of a 3D complex-valued variable $g(\mathit{r})$. Note that Eq. (8) is the gradient-descent step for the back-propagation process, where $\alpha $ is a manually-tuned parameter to adjust the step-size [33]. - (8) Repeats steps (3)–(7) for each illumination angle $\ell =1,2,\dots ,L$, to incrementally refine $\{{n}_{k}(\mathit{r})|k=1,2,\dots N\}$. After one round through all the illumination angles is complete, consolidate all of the RI layers into a single 3D RI volume, $\{{n}_{k}(\mathit{r})|k=1,2,\dots N\}\Rightarrow n({\mathit{r}}_{3\mathrm{D}})$.
- (9) Implement 3D total-variation (TV) regularization on $n({\mathit{r}}_{3\mathrm{D}})$ to stabilize the iterative convergence in the presence of noise and poor conditioning (see Supplement 1), set by the parameter $\beta $, where the operator $\mathrm{prox}\{f({\mathit{r}}_{3\mathrm{D}}),\gamma \}$ is generally defined for some 3D function $f({\mathit{r}}_{3\mathrm{D}})$ and parameter $\gamma $ as$$\mathrm{prox}\{f({\mathit{r}}_{3\mathrm{D}}),\gamma \}=\underset{g({\mathit{r}}_{3\mathrm{D}})}{\mathrm{arg}\mathrm{min}}\{\frac{1}{2}\Vert f({\mathit{r}}_{3\mathrm{D}})-g{({\mathit{r}}_{3\mathrm{D}})\Vert}_{{\ell}^{2}}^{2}+\gamma \mathrm{TV}[g({\mathit{r}}_{3\mathrm{D}})]\},$$where $\mathrm{TV}[\xb7]$ denotes the 3D TV norm. The parameter $\beta $ is manually tuned to optimize the strength of the regularization. The updated $n({\mathit{r}}_{3\mathrm{D}})$ is the current iterative estimate of the sample’s 3D RI.
- (10) Repeat steps (2)–(9) to continue the iterative process until convergence is reached (when the iterative cost function $c(d)$ levels out with respect to iteration index $d$).

## 3. EXPERIMENTAL METHODS

#### A. Optical System Design

Our experimental setup is shown in Fig. 1. Green ($\lambda =532\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ center wavelength) LED light (Thorlabs M530F2) is fiber-coupled into multimode fiber with a 50 μm core diameter (Thorlabs M14L). Although the LED source is not natively coherent, the fiber acts as a spatial filter to sufficiently increase the LED source’s spatial coherence, while still avoiding speckle and other coherent artifacts. The output light from the fiber is then collimated and directed to a mirror mounted on a motorized kinematic mount (M, Thorlabs KS1-Z8), for programmable tip/tilt capabilities. The plane of the mirror is conjugated to the biological sample through a 4f-system consisting of an achromatic doublet and a condenser objective ($\mathrm{L}1\to \mathrm{OBJ}$), such that tip/tilt of M directly results in scanning the illumination angle at the sample. The light exiting the sample is then collected through a second 4f-system, consisting of an imaging objective and an achromatic doublet ($\mathrm{OBJ}\to \mathrm{L}2$). In the system, the condenser and imaging objectives (OBJ) were identical high-NA microscope objectives (Nikon, CFI Plan Apo Lambda $100\times $, NA 1.45). Finally, a third 4f-system was used to de-magnify the transmission output, before imaging onto a 20-megapixel sensor (CMOS, FLIR BFS-U3-200S6M). An adjustable iris (Ir, Thorlabs ID15) in the Fourier plane enables variable tuning of the system’s NA. To avoid aberrations, we limit the NA of the detection system to $\mathrm{NA}=1.2$. This corresponds to a theoretical lateral and axial resolution of around 240 nm and 900 nm, respectively, in the case of a weakly scattering sample.

#### B. Data Acquisition

The motorized kinematic mount for angle scanning was controlled via MATLAB and was programmed to scan the illumination through 120 angles along a spiral trajectory that fills the pupil [Fig. 1(a)]. As each scan point is reached, an Arduino board sends a trigger pulse to the sensor. Each image was captured with $\sim 40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ms}$ of integration time, resulting in a total acquisition time of 4.8 s. Reconstruction was performed with MATLAB running on a desktop computer equipped with Intel(R) Core(TM) i7-5960X CPU at 3.00GHz 64GB RAM CPU and NVidia’s GeForce GTX TitanX GPU. The total reconstruction time to complete 50 iterations of a $1200\times 1200\times 120$ voxel volume was 20.3 h.

Figures 1(b), 1(d), and 1(f) and Fig. 1(c), 1(e), and 1(g) show the raw intensity measurements and associated Fourier magnitude, respectively, of a 3T3 fibroblast cell under different illumination angles. The Fourier transforms demonstrate that the spatial-frequency information is mostly contained within two circular regions in Fourier space, symmetrically positioned around the origin. This is expected for intensity-based imaging of weakly scattering objects with off-axis illumination [35,43] (see Supplement 1 for details).

#### C. Calibration of Illumination Angles

Our reconstruction requires precise knowledge of the illumination angles at the sample. Though illumination angles are outputted by the motorized kinematic mount, we experimentally found that such values did not have sufficient accuracy for satisfactory 3D RI reconstruction (likely due to system imperfections, sample-induced changes and misalignments). Thus, we instead use post-acquisition algorithmic angle-calibration as our means for precise angle measurement. We leverage previous developments in self-calibration to precisely estimate the incident illumination angles from the raw measurements [43]. Figures 1(c), 1(e), and 1(g) show that the illumination wavevectors ${\mathit{k}}_{0}^{\ell}$ can be estimated from the acquisitions’ Fourier transforms, since the center of the circles (illustrated with red arrows) correspond to the illumination angle (and its conjugate). For comparison, Fig. 1(a) shows the illumination angle values outputted by the hardware alongside those estimated by the self-calibration algorithm.

## 4. EXPERIMENTAL RESULTS

We experimentally demonstrate our method for 3D RI reconstruction of calibration objects and biological samples with varying amounts of multiple-scattering. We first image polystyrene beads for validation and then image a fibroblast cell, as an example of a weakly scattering sample. Finally, we look at multiple-scattering samples—*C. elegans* embryos and whole worm. Because all objects considered in this work are highly transparent with essentially no absorption, we show only the real component of the complex-valued RI reconstruction. Supplement 1 provides a comparison between the reconstruction quality obtained by our MSBP framework and the 1st Born reconstruction framework [32].

#### A. Polystyrene Microspheres

Figure 2 shows the 3D RI reconstruction of a conjoined pair of polystyrene spheres, mounted in immersion oil ($n(\lambda )=1.552$ at $\lambda =532\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$). The reconstruction used iterative step-size and regularization parameter set to $\alpha =6\times {10}^{-4}$ and $\beta =4\times {10}^{-4}$, respectively. Figure 2(a) shows a lateral slice through the center of the reconstruction volume. Figure 2(b) shows an axial slice taken across the horizontal dashed-white line in Fig. 2(a), and Fig. 2(c) shows a 3D rendered visualization of the reconstruction volume, clearly capturing the spherical geometry of the microspheres. A line-profile across a single microsphere [Fig. 2(d)] shows that the experimentally measured RI is in good agreement with the expected shape and RI of polystyrene (${n}_{s}(\lambda )=1.598$ at $\lambda =532\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$) [44].

#### B. Fibroblast Cells

We next demonstrate RI reconstruction (Visualization 1, Visualization 2, and Visualization 3, $700\times 700\times 40$ voxels, $\alpha =4\times {10}^{-4}$, $\beta =1\times {10}^{-5}$) of a weakly scattering fibroblast (3T3) cell, using our MSBP model. Figure 3 shows both grayscale and colored (for easier RI visualization after halo-removal [45]) lateral cross-sections at different depths through the reconstruction volume ($z=0.0,1.05$, and 2.10 μm, where $z=0.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$ designates the attachment interface of the fibroblast cell to the coverslip). As adherent, migratory cells, fibroblast cells are known to form dynamic membrane protrusions [46]. These protrusions are clearly visualized in the reconstructed volume, with long filopodial extensions [indicated by yellow arrows in Fig. 3(c)] and broader lamellipodia and membrane ruffling at cell edges [red arrows in Figs. 3(a) and 3(c)]. Intracellular compartments also have strong contrast, notably the cell nucleus with the surrounding nuclear envelope and internal nucleoli. We show a zoom-in of this region in Fig. 3(e) to highlight the nuclear envelope (red arrows), internal nucleoli (yellow arrows), and outer region of optically-dense cellular structure.

The dense fibrous network within the cell is clearly visualized and exhibits a RI of 1.335-1.34 for the individual fiber tendrils. Intracellular lipids and nucleoli consistently demonstrate a $\mathrm{RI}>1.35$ and $\sim 1.335$, respectively, and the internuclear space shows a RI of $\sim 1.33$, matching closely the surrounding media. It has also been observed in previous works that the nuclear RI is lower than that of the general cell cytoplasm for various cell lines [7,47,48]. Figure 3(g) shows a tomographic rendering for easy identification of the 3D cell morphology.

#### C. Caenorhabditis elegans Embryos

Next, we image *Caenorhabditis elegans* (*C. elegans*) embryos captured at varying developmental stages [49]. Unlike the 3T3 cell, embryos are multicellular clusters and are, therefore, more condensed than any one individual cell. Thus, we expect such samples to exhibit more multiple-scattering. In Figs. 4(a)–4(c), we show the raw intensity measurements and associated Fourier amplitudes of the embryo clusters under three different illumination angles. Immediately, we see from Fig. 4(a) that the condensed internal structure within the *C. elegans* embryos results in significant intensity contrast even under orthogonal plane-wave illumination. This contrasts with the observations of the 3T3 fibroblast cell sample, which had virtually no intensity contrast under orthogonal illumination [see Fig. 1(b)]. This suggests that though the fibroblast cell can be modeled as weakly scattering, the *C. elegans* embryos cannot.

More rigorous evidence of multiple-scattering comes from the Fourier transforms of the intensity acquisitions. As mentioned earlier and described in more detail in the Supplement 1, off-axis illumination of a weakly scattering sample results in intensity acquisitions where the spatial-frequencies lie within two circular regions in Fourier space. While the intensity spatial-frequency content of the 3T3 fibroblast cell [Figs. 1(c), 1(e), and 1(g)] was mainly contained in these circular regions, the *C. elegans* embryos show significant spatial-frequency content outside the two circles, as indicated by yellow arrows in Figs. 4(a)–4(c). This provides further evidence of multiple-scattering.

We show lateral slices from the 3D RI reconstruction (Visualization 4, Visualization 5, and Visualization 6, $1100\times 1100\times 120$ voxels, $\alpha =6\times {10}^{-4}$ and $\beta =3.5\times {10}^{-4}$) of the embryo sample in both gray-scale [Figs. 4(d)–4(f)] and color-scale Figs. 4(h)–4(j), for three different depths ($z=-2.6,0$, and $+2.6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$). Figures 4(g) and 4(k) show axial cross-sections, and Fig. 4(l) shows a 3D rendering. The features visualized in this 3D RI reconstruction fit well with their developmental stages, and individual cellular compartments can be identified. Yellow arrows in Fig. 4(d) and Fig. 4(e) indicate embryos visualized during their comma stage and general 8-cell to 26- cell stages, respectively. The axial cross-section in Fig. 4(g) demonstrates in 3D the embryos’ globular nature and heterogeneous composition, with RI values ranging from 1.33 to above 1.37, and shows a bulk RI increase towards the periphery of the embryos. Distinct eggshells are visualized encapsulating all embryos and demonstrate a moderately high RI $\sim 1.35$, which matches their known lipid-rich composition. Moreover, the lipoprotein complexes containing lipid and yolk proteins, localized between the embryo and eggshell, and show the highest RI as $>1.375$ [indicated with yellow arrows in Figs. 4(i) and 4(j)].

#### D. Whole Caenorhabditis elegans Worms

Our final experiment generates a high-resolution 3D RI reconstruction of a whole adult hermaphrodite *C. elegans* worm (Visualization 7, $1914\times 10408\times 118$ voxels, $\alpha =4\times {10}^{-4}$, $\beta =4\times {10}^{-4}$). Because the imaging objective’s field-of-view did not fit the whole worm, the final reconstruction volume was stitched together from 14 individually-reconstructed patches (details in Supplement 1). We note that unlike the embryos sample, the spatial-frequencies contained within the intensity measurements of the worm (shown in Supplement 1) do not contain any discernable symmetric circular regions. This indicates that the multiple-scattering in the worm sample dominates over the single-scattering and that the whole *C. elegans* worm is even more multiply-scattering than the embryos, as expected.

In Figs. 5(a) and 5(b), we show two lateral slices through the reconstruction volume at axial positions of $z=0.0$ and $-4.9\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, respectively. The major components of the reproductive and digestive systems are clearly identified and include the pharynx, intestine, gonads, spermathecal, and fertilized eggs [50]. Figure 5(c) presents an axial cross-section of the pharynx and clearly demonstrates the interior three-fold rotationally symmetric organization of the muscles and marginal cells around the pharyngeal lumen. Visualizing this morphology typically requires a heavily invasive preparation protocol (cross-sectional slicing through the *C. elegans* pharynx followed by electron microscopy) [51]. Figures 5(d) and 5(e) present a zoom-in of the worm’s mouth and pharyngeal tip at axial positions of $z=+2.42$ and $+5.19\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$. Yellow arrows indicate *E. coli* bacteria, a food-source for the worm, present on the coverslip at the time of sample preparation. Red arrows indicate the cuticle surrounding the pharyngeal lumen and the buccal cavity, respectively.

We outline three regions-of-interest (ROI) within Fig. 5(a), labeled as ①, ②, and ③, which highlight prominent structures within the *C. elegans* worm. Figure 5(f) shows a zoom of the ROI ① and qualitatively highlights the features pertaining to the pharynx-bulb, intestinal cavity, and the start of the intestinal lumen, as indicated by the yellow, red, and green arrows, respectively. Figure 5(j) shows a zoom-in of the ROI ② and depicts fertilized eggs and the distal gonad, as indicated by the yellow and green arrows, respectively. Figure 5(n) shows a zoom of the ROI ③ at the tail-end of the worm and highlights the distal and proximal gonads. To quantitatively visualize the worm’s 3D biological RI, we present RBG-colored cross-sectional images of all ROIs at various axial positions. Figures 5(g)–5(i) show quantitative RI cross-sections of the ROI ① at $z=0.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, $+2.4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, and $+4.9\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, respectively. Figures 5(k)–5(m) and Fig. 5(o)–5(q) show quantitative RI cross-sections of the ROI ② and ③, respectively, at $z=0.0,-2.4$, and $-4.9\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$. The bulk tissue of the pharynx and gonads are of relatively consistent RI $\sim 1.35$. The fertilized eggs, however, demonstrate high heterogeneity (corroborated by DIC imaging) with RI values ranging from 1.335 to 1.35, and regions with high density of lipid droplets, such as the intestine or the proximal gonads, consist of the most variable scattering. Individual lipid droplets can demonstrate RI values as high as 1.40 and are often directly surrounded by features of low RI, $\sim 1.335$. Figures 5(r)–5(t) show 3D tomographic visualizations of the ROIs ①, ②, and ③, respectively.

## 5. DISCUSSION

This work demonstrates a cost-effective and simple optical hardware system that uses computational imaging for biological imaging problems where no analytical solution exists, as is the case with multiple-scattering. This comes at a cost of higher computational requirements than standard ODT. In this work, every angled illumination corresponds to two passes through all layers of the reconstruction volume, which increments the sample’s predicted 3D RI one step along its convergence curve. Factors such as increasing the number of acquisitions (to increase the likelihood of correct convergence for the reconstruction) or the resolution of the reconstruction volume further slow the computation. To address this, future work will leverage the recent advent of big-data processing via GPU acceleration or cloud-computing to significantly reduce reconstruction times.

Another interesting extension would be to explore the maximum resolution achievable using the introduced framework. Previous works have demonstrated that multiple-scattering, if considered appropriately, can significantly improve a microscope’s resolving power to beyond the diffraction limit [34,37,38]. However, this capability for resolution enhancement is difficult to quantify because it depends on the multiple-scattering characteristics of the object, which are generally unknown. Furthermore, strongly scattering objects typically also exhibit strong back-scattering, which is not mathematically treated in our MSBP model. Thus, future work will also explore strategies to incorporate back-scattering into our imaging model as well as to estimate the maximum attainable resolution, based on inferences of the sample’s multiple-scattering characteristics.

## 6. CONCLUSION

In this work, we have introduced a new computational technique that uses the multi-slice beam-propagation model to reconstruct 3D RI with high-resolution, even in multiply-scattering samples. Compared to the standard ODT techniques, our method offers two key advantages: (1) intensity-only acquisitions enable phase-stable and speckle-free measurements with a dramatically simpler and cheaper optical system and (2) we impose no constraints upon the sample to be weakly scattering, such that suitable samples need not be limited to a sparse distribution of individual cells. This opens up 3D RI microscopy to the more general class of biological samples that include dense cell-clusters or multicellular organisms.

We experimentally demonstrated our method by first reconstructing 3D RI of 3T3 fibroblast cells and *C. elegans* embryos, as examples of single-cell (weakly scattering) and multi-cell cluster (multiple-scattering) samples, respectively. In both cases, the MSBP model enabled high-fidelity and robust RI reconstruction. Our final demonstration was the 3D RI reconstruction of a whole *C. elegans* worm as an example of the applicability of this technique to whole multicellular organisms, even in the presence of multiple-scattering. To the best of our knowledge, this is the first demonstration of high-resolution ($\mathrm{NA}>1.0$) intensity-based 3D RI reconstruction of multiple-scattering biological samples. Due to the simplicity of the introduced optical system, easy hardware additions can be implemented to adapt the system to existing fluorescent microscopes and enable 3D multimodal imaging [47,52,53].

## Funding

Gordon and Betty Moore Foundation’s Data-Driven Discovery Initiative (GBMF4562); NIH Ruth L. Kirschstein National Research Service Award (F32GM129966).

## Acknowledgment

The authors would like to thank Dr. Emrah Bostan and other members of the Computational Imaging Lab for helpful discussions.

See Supplement 1 for supporting content.

## REFERENCES

**1. **D. T. Ross, U. Scherf, M. B. Eisen, C. M. Perou, C. Rees, P. Spellman, V. Iyer, S. S. Jeffrey, M. Van de Rijn, and M. Waltham, “Systematic variation in gene expression patterns in human cancer cell lines,” Nat. Genet. **24**, 227–235 (2000). [CrossRef]

**2. **A. Rustom, R. Saffrich, I. Markovic, P. Walther, and H.-H. Gerdes, “Nanotubular highways for intercellular organelle transport,” Science **303**, 1007–1010 (2004). [CrossRef]

**3. **M. Okuda, K. Li, M. R. Beard, L. A. Showalter, F. Scholle, S. M. Lemon, and S. A. Weinman, “Mitochondrial injury, oxidative stress, and antioxidant gene expression are induced by hepatitis C virus core protein,” Gastroenterology **122**, 366–375 (2002). [CrossRef]

**4. **R. Rizzuto, M. Brini, P. Pizzo, M. Murgia, and T. Pozzan, “Chimeric green fluorescent protein as a tool for visualizing subcellular organelles in living cells,” Curr. Biol. **5**, 635–642 (1995). [CrossRef]

**5. **W. Choi, C. Fang-Yen, K. Badizadegan, S. Oh, N. Lue, R. R. Dasari, and M. S. Feld, “Tomographic phase microscopy,” Nat. Methods **4**, 717–719 (2007). [CrossRef]

**6. **V. Lauer, “New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope,” J. Microsc. **205**, 165–176 (2002). [CrossRef]

**7. **Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express **17**, 266–277 (2009). [CrossRef]

**8. **R. Fiolka, K. Wicker, R. Heintzmann, and A. Stemmer, “Simplified approach to diffraction tomography in optical microscopy,” Opt. Express **17**, 12407–12417 (2009). [CrossRef]

**9. **S. Lee, K. Kim, A. Mubarok, A. Panduwirawan, K. Lee, S. Lee, H. Park, and Y. Park, “High-resolution 3-D refractive index tomography and 2-D synthetic aperture imaging of live phytoplankton,” J. Opt. Soc. Korea **18**, 691–697 (2014). [CrossRef]

**10. **T. Kim, R. Zhou, M. Mir, S. D. Babacan, P. S. Carney, L. L. Goddard, and G. Popescu, “White-light diffraction tomography of unlabelled live cells,” Nat. Photonics **8**, 256–263 (2014). [CrossRef]

**11. **T. H. Nguyen, M. E. Kandel, M. Rubessa, M. B. Wheeler, and G. Popescu, “Gradient light interference microscopy for 3D imaging of unlabeled specimens,” Nat. Commun. **8**, 210 (2017). [CrossRef]

**12. **Y. Park, M. Diez-Silva, G. Popescu, G. Lykotrafitis, W. Choi, M. S. Feld, and S. Suresh, “Refractive index maps and membrane dynamics of human red blood cells parasitized by Plasmodium falciparum,” Proc. Natl. Acad. Sci. USA **105**, 13730–13735 (2008). [CrossRef]

**13. **W. J. Eldridge, A. Sheinfeld, M. T. Rinehart, and A. Wax, “Imaging deformation of adherent cells due to shear stress using quantitative phase imaging,” Opt. Lett. **41**, 352–355 (2016). [CrossRef]

**14. **N. T. Shaked, L. L. Satterwhite, N. Bursac, and A. Wax, “Whole-cell-analysis of live cardiomyocytes using wide-field interferometric phase microscopy,” Biomed. Opt. Express **1**, 706–719 (2010). [CrossRef]

**15. **J. Jung, K. Kim, J. Yoon, and Y. Park, “Hyperspectral optical diffraction tomography,” Opt. Express **24**, 2006–2012 (2016). [CrossRef]

**16. **S. S. Kou and C. J. Sheppard, “Image formation in holographic tomography: high-aperture imaging conditions,” Appl. Opt. **48**, H168–H175 (2009). [CrossRef]

**17. **M. Habaza, B. Gilboa, Y. Roichman, and N. T. Shaked, “Tomographic phase microscopy with 180 rotation of live cells in suspension by holographic optical tweezers,” Opt. Lett. **40**, 1881–1884 (2015). [CrossRef]

**18. **B. Simon, M. Debailleul, M. Houkal, C. Ecoffet, J. Bailleul, J. Lambert, A. Spangenberg, H. Liu, O. Soppera, and O. Haeberlé, “Tomographic diffractive microscopy with isotropic resolution,” Optica **4**, 460–463 (2017). [CrossRef]

**19. **A. Devaney, “Inverse-scattering theory within the Rytov approximation,” Opt. Lett. **6**, 374–376 (1981). [CrossRef]

**20. **F. Lin and M. Fiddy, “The Born-Rytov controversy: I. Comparing analytical and approximate expressions for the one-dimensional deterministic case,” J. Opt. Soc. Am. A **9**, 1102–1110 (1992). [CrossRef]

**21. **B. Chen and J. J. Stamnes, “Validity of diffraction tomography based on the first Born and the first Rytov approximations,” Appl. Opt. **37**, 2996–3006 (1998). [CrossRef]

**22. **E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. **1**, 153–156 (1969). [CrossRef]

**23. **Y. Cotte, F. Toy, P. Jourdain, N. Pavillon, D. Boss, P. Magistretti, P. Marquet, and C. Depeursinge, “Marker-free phase nanoscopy,” Nat. Photonics **7**, 113–117 (2013). [CrossRef]

**24. **K. Kim, H. Yoon, M. Diez-Silva, M. Dao, R. R. Dasari, and Y. Park, “High-resolution three-dimensional imaging of red blood cells parasitized by Plasmodium falciparum and in situ hemozoin crystals using optical diffraction tomography,” J. Biomed. Opt. **19**, 011005 (2014). [CrossRef]

**25. **S. Chowdhury, W. J. Eldridge, A. Wax, and J. Izatt, “Refractive index tomography with structured illumination,” Optica **4**, 537–545 (2017). [CrossRef]

**26. **M. H. Maleki and A. J. Devaney, “Phase-retrieval and intensity-only reconstruction algorithms for optical diffraction tomography,” J. Opt. Soc. Am. A **10**, 1086–1092 (1993). [CrossRef]

**27. **G.-Z. Yang, B.-Z. Dong, B.-Y. Gu, J.-Y. Zhuang, and O. K. Ersoy, “Gerchberg-Saxton and Yang-Gu algorithms for phase retrieval in a nonunitary transform system: a comparison,” Appl. Opt. **33**, 209–218 (1994). [CrossRef]

**28. **U. S. Kamilov, D. Liu, H. Mansour, and P. T. Boufounos, “A recursive Born approach to nonlinear inverse scattering,” IEEE Signal Process. Lett. **23**, 1052–1056 (2016). [CrossRef]

**29. **D. Ancora, D. Di Battista, G. Giasafaki, S. E. Psycharakis, E. Liapis, J. Ripoll, and G. Zacharakis, “Optical projection tomography via phase retrieval algorithms,” Methods **136**, 81–89 (2018). [CrossRef]

**30. **E. Froustey, E. Bostan, S. Lefkimmiatis, and M. Unser, “Digital phase reconstruction via iterative solutions of transport-of-intensity equation,” in *13th Workshop on Information Optics (WIO)* (IEEE, 2014), pp. 1–3.

**31. **H.-Y. Liu, D. Liu, H. Mansour, P. T. Boufounos, L. Waller, and U. S. Kamilov, “SEAGLE: sparsity-driven image reconstruction under multiple scattering,” IEEE Trans. Comput. Imag. **4**, 73–86 (2018). [CrossRef]

**32. **R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang, “Diffraction tomography with Fourier ptychography,” Optica **3**, 827–835 (2016). [CrossRef]

**33. **U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comput. Imag. **2**, 59–70 (2016). [CrossRef]

**34. **L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica **2**, 104–111 (2015). [CrossRef]

**35. **R. Ling, W. Tahir, H.-Y. Lin, H. Lee, and L. Tian, “High-throughput intensity diffraction tomography with a computational microscope,” Biomed. Opt. Express **9**, 2130–2141 (2018). [CrossRef]

**36. **M. Chen, L. Tian, and L. Waller, “3D differential phase contrast microscopy,” Biomed. Opt. Express **7**, 3940–3950 (2016). [CrossRef]

**37. **J.-H. Park, C. Park, H. Yu, J. Park, S. Han, J. Shin, S. H. Ko, K. T. Nam, Y.-H. Cho, and Y. Park, “Subwavelength light focusing using random nanoparticles,” Nat. Photonics **7**, 454–458 (2013). [CrossRef]

**38. **K. Belkebir, P. C. Chaumet, and A. Sentenac, “Influence of multiple scattering on three-dimensional imaging with optical diffraction tomography,” J. Opt. Soc. Am. A **23**, 586–595 (2006). [CrossRef]

**39. **J. Van Roey, J. Van der Donk, and P. Lagasse, “Beam-propagation method: analysis and assessment,” J. Opt. Soc. Am. **71**, 803–810 (1981). [CrossRef]

**40. **Y. Chung and N. Dagli, “An assessment of finite difference beam propagation method,” IEEE J. Quantum Electron. **26**, 1335–1339 (1990). [CrossRef]

**41. **A. M. Maiden, M. J. Humphry, and J. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A **29**, 1606–1614 (2012). [CrossRef]

**42. **L.-H. Yeh, J. Dong, J. Zhong, L. Tian, M. Chen, G. Tang, M. Soltanolkotabi, and L. Waller, “Experimental robustness of Fourier ptychography phase retrieval algorithms,” Opt. Express **23**, 33214–33240 (2015). [CrossRef]

**43. **R. Eckert, Z. F. Phillips, and L. Waller, “Efficient illumination angle self-calibration in Fourier ptychography,” Appl. Opt. **57**, 5434–5442 (2018). [CrossRef]

**44. **N. Sultanova, S. Kasarova, and I. Nikolov, “Dispersion properties of optical polymers,” Acta Phys. Pol. A **116**, 585–587 (2009). [CrossRef]

**45. **M. E. Kandel, M. Fanous, C. Best-Popescu, and G. Popescu, “Real-time halo correction in phase contrast imaging,” Biomed. Opt. Express **9**, 623–635 (2018). [CrossRef]

**46. **J. V. Small, T. Stradal, E. Vignal, and K. Rottner, “The lamellipodium: where motility begins,” Trends Cell Biol. **12**, 112–120 (2002). [CrossRef]

**47. **S. Chowdhury, W. J. Eldridge, A. Wax, and J. A. Izatt, “Structured illumination microscopy for dual-modality 3D sub-diffraction resolution fluorescence and refractive-index reconstruction,” Biomed. Opt. Express **8**, 5776–5793 (2017). [CrossRef]

**48. **M. Schürmann, J. Scholze, P. Müller, J. Guck, and C. J. Chan, “Cell nuclei have lower refractive index and mass density than cytoplasm,” J. Biophoton. **9**, 1068–1076 (2016). [CrossRef]

**49. **J. E. Sulston, E. Schierenberg, J. G. White, and J. Thomson, “The embryonic cell lineage of the nematode Caenorhabditis elegans,” Dev. Biol. **100**, 64–119 (1983). [CrossRef]

**50. **J. G. White, E. Southgate, J. N. Thomson, and S. Brenner, “The structure of the nervous system of the nematode Caenorhabditis elegans,” Philos. Trans. R. Soc. B **314**, 1–340 (1986). [CrossRef]

**51. **S. E. Mango, *The C. Elegans Pharynx: a Model for Organogenesis* (2007).

**52. **Y. Park, G. Popescu, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Diffraction phase and fluorescence microscopy,” Opt. Express **14**, 8263–8268 (2006). [CrossRef]

**53. **K. Kim, W. S. Park, S. Na, S. Kim, T. Kim, W. D. Heo, and Y. Park, “Correlative three-dimensional fluorescence and refractive index tomography: bridging the gap between molecular specificity and quantitative bioimaging,” Biomed. Opt. Express **8**, 5688–5697 (2017). [CrossRef]