High-resolution 3D refractive index microscopy of multiple-scattering samples from intensity images

Optical diffraction tomography (ODT) reconstructs a samples 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 (NA>1.0) 3D RI distributions of multiple-scattering samples. The method acquires intensity-only measurements from different illumination angles, and then solves a non-linear optimization problem to recover the sample 3D RI distribution. We experimentally demonstrate 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 250 nm and 900 nm, respectively.

(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, then recovers the sample's 3D RI distribution.Beyond the work by Tian et al [34], we present a rigorous nonlinear optimization protocol that performs robustly at high resolutions (NA > 1.0) and enables full axial sampling of the system's diffraction limited pointspread-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.

A. Multi-slice beam-propagation forward model
Beam-propagation 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 modelled via sequential layer-to-layer propagation of the electric-field [33,34,41].Mathematically, this can be recursively written as: where  denotes the 2D spatial position vector, {  () |  = 1,2, … } denotes the electric-field at the th layer of an object that is approximated with  equally-spaced layers separated by distance ) ℱ{. }}, with ℱ{. } and ℱ −1 {. } denoting the 2D Fourier and inverse Fourier transforms, respectively, and  denting the 2D spatial frequency vector.The boundary condition to initialize Eq. ( 1) is simply the incident planar electric field illuminating the sample at a particular angle, i.e.,  0 () = exp(   • ), where   is the 2D illumination wave-vector.
The exit electric-field,   (), 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, where () denotes the system's pupil function, and ̂ refers to the distance between the plane of   () and the plane within the object's volume at conjugate focus to the imaging plane of the system.When the imaging system is focused at the center of the object, ̂= Δ /2.

B. Inverse problem formulation
Multiple measurements of the object are captured at varying illumination angles.We denote  0 ℓ () = exp(    • ) for ℓ = 1, 2, … ,  as a set of  planar electric-fields incident onto the first layer of the sample, at various angles set by their wave-vectors    .The corresponding intensity measurements are denoted by  ℓ ().
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, Here,  3 = 〈, 〉 is a 3D spatial position vector, such that ( 3 ) =   (),  = 1,2,3, … , .The non-linear operator  ℓ {. } denotes the ℂ 3 → ℂ 2 forward model operation that predicts the 2D electric-field distribution measured at the camera after illuminating ( 3 ) with the incident electric-field  0 ℓ (), as described by Eqs (1) and (2).

C. Reconstruction framework
Our proposed framework to minimize Eq (4) and reconstruct  ̂( 3 ) draws from previous MSBP works by Kamilov et al [33] and Tian et al [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 a -layer reconstruction volume with constant value   .This will serve as the initial estimate of ( 3 ).Then, initialize the iteration index at this step to  = 0.
Note that Eq. ( 8) illustrates the gradient-descent step for the back-propagation process, where  is often a manually-tuned parameter to adjust the step-size.
9. Conduct a 3D total-variation (TV) regularization process on ( 3 ) to stabilize the iterative convergence in the presence of physical effects in the imaging system that are unaccounted for in the MSBP forward model, such as camera noise, coherent source noise, optical aberrations, etc.The strength of regularization is set by the parameter .
where the operator prox{ ( 3 ), } is generally defined for some 3D function ( 3 ) and parameter  as, where TV[.] denotes the 3D total-variation norm.The final ( 3 ), after being output by the regularizer, now becomes the current iterative estimate of the object's 3D RI.
10. Repeat steps (2)-( 9) to continue the iterative process, until convergence is reached.This point can be identified when the iterative cost function () levels out with respect to iteration index .

A. Optical system design
Our experimental setup is shown in Fig. 1. Green (λ = 532 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 (L1 → 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 (OBJ → L2).In the system, the condenser and imaging objectives (OBJ) were identical high-NA microscope objectives (Nikon, CFI Plan Apo Lambda 100x, 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 NA = 1.1.This corresponds to a theoretical lateral and axial resolution of 240 nm and 890 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 ~40 ms of integration time, resulting in a total acquisition time of 4.8 seconds.Reconstruction was performed with MATLAB running on a desktop computer equipped with Intel(R) Core(TM) i7-5960X CPU @ 3.00GHz 64GB RAM CPU and NVidia's GeForce GTX TitanX GPU.Total reconstruction time to complete 50 iterations of a 1200⨯1200⨯120 voxel volume was 20.3 hours.Fig. 1(b,d,f) and Fig. 1(c,e,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 Supplementary Information for details).

C. Calibration of angular illuminations
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 postacquisition 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].Fig. 1(c,e,g) shows that the illumination wavevectors    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.

IV. 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.The Supplementary Information provides a comparison between the reconstruction quality obtained by our MSBP framework and the 1st Born reconstruction framework [32].

A. Polystyrene microspheres
Fig. 2 shows the 3D RI reconstruction of a conjoined pair of polystyrene spheres, mounted in immersion oil (() = 1.552 at  = 532 nm).The reconstruction used iterative step-size and regularization parameter set to  = 6 × 10 −4 and  = 4 × 10 −4 , respectively.Fig. 2(a) shows a lateral slice through the center of the reconstruction volume.Fig. 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 (  () = 1.598 at  = 532 nm) [44].

B. Fibroblast cells
We next demonstrate RI reconstruction (Supplementary videos 1-3, 700⨯700⨯40 voxels,  = 4 × 10 −4 ,  = 1 × 10 −5 ) of a weakly-scattering fibroblast (3T3) cell, using our MSBP model.Fig. 3 shows both grayscale and colored (for easier RI visualization after halo-removal [45]) lateral cross-sections at different depths through the reconstruction volume ( = 0.0, 1.05, and 2.10 um, where  = 0.0 um 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 Fig. 3(a,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 RI of 1.335-1.34for the individual fiber tendrils.Intracellular lipids and nucleoli consistently demonstrate RI > 1.35 and ~1.335, respectively, and the internuclear space shows a RI of ~1.33, matching closely the surrounding media.Nuclear RI being lower than that of the general cell cytoplasm has also been observed in previous works for various cell lines [7,47,48].Fig. 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 thus more condensed than any one individual cell.Thus, we expect such samples to exhibit more multiple scattering.In Fig. 4(a-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 modelled as weakly-scattering, 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 Supplementary Information, 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 (Fig. 1(c,e,g)) was mainly contained in these circular regions, the C. elegans embryos show significant spatial-frequency content outside these regions, as indicated by yellow arrows in Fig. 4(a-c).This provides further evidence of multiple scattering.
We show lateral slices from the 3D RI reconstruction (Supplementary videos 4-6, 1100⨯1100⨯120 voxels,  = 6 × 10 −4 and  = 3.5 × 10 −4 ) of the embryo sample in both grayscale (Fig. 4(d-f)) and color-scale Fig. 4(h-j), for three different depths ( = −2.6,0, and +2.6 um).Fig. 4(g,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 heterogenous 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 ~1.35, which matches their known lipid-rich composition.Moreover, the lipoprotein complexes containing lipid and yolk proteins, localized between the embryo and eggshell, show the highest RI > 1.375 (indicated with yellow arrows in Fig. 4(i,j)).

D. Whole Caenorhabditis elegans worms
Our final experiment generates a high-resolution 3D RI reconstruction of a whole adult hermaphrodite C. elegans worm (Supplementary videos 7-10, 1914⨯10408⨯118 voxels,  = 4 × 10 −4 ,  = 4 × 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 Supplementary Information).We note that unlike the embryos sample, the spatial-frequencies contained within the intensity measurements of the worm (shown in Supplementary Information) 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 Fig. 5(a,b), we show two lateral slices through the reconstruction volume at axial positions of  = 0.0 and −4.9 um, respectively.The major components of the reproductive and digestive systems are clearly identified, and include the pharynx, intestine, gonads, spermathecal, and fertilized eggs [50].Fig. 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].Fig. 5(d,e) present a zoom-in of the worm's mouth and pharyngeal tip at axial positions of  = +2.42 and +5.19 um.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), labelled ①, ②, and ③, that highlight prominent structures within the C. elegans worm.Fig. 5(f) shows a zoom of 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.Fig. 5(j) shows a zoom-in of ROI ②, and depicts fertilized eggs and the distal gonad, as indicated by the yellow and green arrows, respectively.Fig. 5(n) shows a zoom of 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.Fig. 5(g,h,i) show quantitative RI cross-sections of ROI ① at  = 0.0, +2.4 and +4.9 um, respectively.Fig. 5(k,l,m) and Fig. 5(o,p,q) show quantitative RI cross-sections of ROI ② and ③, respectively, at  = 0.0, −2.4 and −4.9 um.The bulk tissue of the pharynx and gonads are of relatively consistent RI ~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 ~1.335.Fig. 5(r,s,t) show 3D tomographic visualizations of ROIs ①, ②, and ③, respectively.

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 down 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 sample, which are generally unknown.Thus, future work will also explore strategies to estimate the maximum attainable resolution of a sample, based on inferences of its multiple-scattering characteristics.

Conclusion
In this work, we have introduced a new computational technique that uses the multi-slice beampropagation model to reconstruct 3D refractive index with high-resolution, even in multiplyscattering samples.Compared to standard ODT techniques, our introduced method offers two key advantages, 1) intensity-only acquisitions enable phase-stable and speckle-free measurements with a dramatically simpler and cheaper optical system.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 multicell 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 (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); Ruth L. Kirschstein National Research Service Award (F32GM129966) spectra shows virtually no distinct "brightfield" signal.This signifies that the C. elegans worm is an even stronger multiple-scattering sample than the C. elegans embryo.Fig. S3(g,h,i) and Fig. S3(j,k,l) show the lateral cross-section through the center of the 3D RI volumes of the three samples, as reconstructed using the 1st Born and MSBP inversion models.Because the 3T3 fibroblast cell is a weakly-scattering sample, we expect the 1st Born assumption to be valid and to enable high-fidelity reconstruction.As expected, the reconstructed RI via 1st Born and MSBP match well (Fig. S3(g,j)).In the case of the multiple-scattering C. elegans embryo and worm samples, however, the reconstructions via the 1st Born (Fig. S3(h,i)) and MSBP Fig. S3(k,l)) scattering models show distinct and fundamental differences.1st Born demonstrates an inability to reconstruct the lower spatial-frequencies in the presence of multiple-scattering, and essentially high-pass filters the RI content.This observation has been affirmed by previous works as well [7].The degree of high-pass filtering is dependent on the degree of multiple-scattering within the sample, and the 1st Born RI reconstruction of the C. elegans worm (Fig. S3(i)) shows greater high-pass filtering than in that of the embryo (Fig. S3(h)).In both cases, important biological features that are easily visualized in the MSBP reconstruction, such as the individual cells within the embryo and the pharynx muscles and buccal cavity within the worm's head, cannot be visualized in the 1st Borne reconstructions.

Fig. 1 .
Fig. 1.(a) Optical system design with illumination angle trajectories, comparing the raw position values outputted from the kinematic tilt mirror to that after algorithmic self-calibration.(b,d,f) Raw intensity acquisitions and (c,e,g) amplitudes of associated Fourier transforms, after illuminating the sample with varying angles.

Fig. 2 .
Fig. 2. 3D reconstruction of refractive index (RI) for two 3 um diameter polystyrene microspheres (n=1.598) in index-matched oil (n=1.552).(a) Lateral (x-y) slice, (b) axial (x-z) slice, taken along the horizontal dashed line in (a), and (c) 3D rendering of the reconstructed microspheres.(d) Cross-cut of the refractive index along the vertical dashed line in (a) matches well with expected values from the sample.

Fig. 4 .
Fig. 4. 3D RI reconstruction of a C. elegans embryo cluster is presented.(a-c) Raw intensity-acquisitions at varying illumination angles are shown with corresponding Fourier transforms.(d,e,f) Lateral cross-sectional planes through the 3D reconstruction volume are shown, at axial positions  = −2.6,0, +2.6 um, respectively.(g) Axial cross-sectional plane, taken from the reconstruction volume at the location indicated by dashed yellow in (d-f) is shown.(h,i,j,k) Lateral and axial cross-sectional views shown in (d,e,f,g) are redisplayed in RGB-color to facilitate easy visual inspection of refractive index.Color-bar to decode RI from color is located in (h).(l) 3D tomographic rendering of the embryo cluster is shown.

Fig. 5 .
Fig. 5. 3D RI reconstruction of a whole C. elegans worm is presented.(a,b) Lateral cross-sectional planes through the 3D reconstruction volume are shown, at axial positions  = 0, −4.9 um, respectively.Major components of the reproductive and digestive systems are labelled.(c) Axial cross-sectional plane through the pharynx is shown.(d,e) Lateral cross-sectional zooms of the worm's head region at axial positions  = +2.42,+5.19 um, respectively, are shown.Zooms of ROIs ①, ②, and ③ are shown in (f,j,n), respectively.Lateral RGB-colored (for easy visual inspection of RI) cross-sections of ROIs ①, ②, and ③, are shown with defocus in (g,h,i), (k,l,m), and (o,p,q), respectively.(r,s,t) 3D tomographic renderings of ROIs ①, ②, and ③, respectively, are shown.