1. Introduction

Compared with two-dimensional (2D) culture models, three-dimensional (3D) culture models, such as spheroids and organoids, provide a more effective ex vivo approach to replicate the growth of and interaction between cells in in vivo environments [1]. Therefore, 3D culture models have become increasingly important in various biomedical applications [2]. Fluorescence-based tomographic imaging, including multiphoton [3], confocal [4], and light-sheet microscopy [5], has played a central role in visualizing intracellular organelles in 3D culture models or native tissues with high molecular specificity. However, fluorescence imaging could lead to undesired effects, such as phototoxicity [6] and photobleaching [7]. In contrast, label-free microscopy is almost free from such limitations and provides morphological information regarding samples in a label-free manner at the expense of molecular specificity. However, the complex structure of 3D samples induces multiple scattering (MS) of light, which disturbs image formation and limits the imaging ability of label-free microscopy to shallow depths.

Typically, MS degrades the optical imaging performance by primarily causing two problems: speckle noise [8] and degradation of the single-to-multiple-scattering ratio (SMR) [9,10]. Speckle noise is caused by the coherent interference by MS light; the SMR defines the depth limit of the imaging system that develops images using single-scattering light. Numerous approaches have been proposed to address these problems [11,12]. We collectively call such approaches “multiple-scattering suppressive (MSS) techniques.” Among the existing MSS techniques [11,12], the coherent or incoherent accumulations of fields are commonly considered to address these problems. The former approach exhibits strong MS tolerance by selectively extracting coherently interfered photons in an optical [4,13] or a computational [14] manner. The extraction of single-scattering light improves the SMR and suppresses speckle noise. A typical example of the spatially coherent accumulation of fields is confocal microscopy [4], where spatially accumulated fields are selectively extracted by a pinhole. We refer to this approach as “coherent MSS (cMSS) techniques.” The latter approach, i.e., the incoherent accumulations of fields, involves averaging the intensity variations of fields by optically realizing or computationally emulating a spatially incoherent light source [15], resulting in the suppression of speckle noise. Phase contrast and differential interference contrast (DIC) microscopy are representative examples of the latter approach. In addition, quantitative phase gradient imaging (QPGI), which is a quantitative version of DIC microscopy, has been receiving attention [16,17]. We refer to this approach as “incoherent MSS (iMSS) technique.” cMSS techniques provide better MS tolerance than iMSS techniques because cMSS techniques improve SMR whereas iMSS techniques do not. However, the superior imaging performance of cMSS techniques collapses when sample-induced aberration (SIA) is strong.

SIA is predominantly caused by the heterogeneity of the sample and distorts the wavefront of light. Such aberration significantly degrades the imaging performance of spatial cMSS techniques by deteriorating their spatial coherency [18], because the effective size of the spatial gating is broadened by SIA. Conversely, spatial iMSS techniques are more tolerant to SIA than spatial cMSS techniques because the former involves the incoherent accumulation of fields. Therefore, a trade-off exists between the MS and SIA tolerances of MSS techniques. Biological specimens often have a complex structure that induces MS and a spherical surface that induces SIA; thus, SIA and MS are equally important when imaging the specimens. In confocal and multiphoton microscopy, adaptive optics [19] or refractive index (RI) matching media [20] have been utilized to overcome the trade-off problem by removing SIA.

As an established label-free tomography based on MSS techniques, optical coherence tomography (OCT) [13] has been extensively applied in vivo, such as in ophthalmology [21] and cardiology [22]. OCT realizes deep imaging using cMSS techniques: scanning confocal gating and coherent heterodyne detection of spatially and temporally accumulated fields. Owing to the advantages of these techniques, OCT achieves a higher SMR than several other optical microscopy techniques [9]. To improve image quality, the incoherent accumulation of fields has also been utilized to suppress speckle noise in OCT images [15,23]. Limitations owing to SIA have also been addressed by estimating the RI of samples [23]. Recently, several studies have applied OCT to the imaging of the 3D culture models of tissue [24,25]. Unfortunately, OCT signals cannot always visualize single cells and organelles because structures associated with strong backscattering responses conceal weakly backscattering structures [24]. To address this problem, a temporal fluctuation analysis of OCT signals enabled the visualization of subcellular structures and the metabolic contrast of organoids [24,25]. However, this technique requires the acquisition of hundreds of frames at each depth.

Another label-free technique, RI tomography, which provides an intrinsic contrast of samples with a subcellular resolution, has been attracting increased interest [26–28] and has been realized in various configurations, including an interferometer with angle-scanning illumination [29–31], a structured illumination microscope [32], a flowcytometry with sample rotation [33], an LED array microscope [34,35], and an add-on module to a commercial microscope [36]. Optical diffraction tomography (ODT) is a representative RI tomography technique, which involves capturing interferometric fields with angle-scanning illuminations [30,31]. The reconstruction algorithms of ODT are generally based on linearized models, the first Born [37] or Rytov [38] approximation, and improved the quality of reconstructed tomograms using regularization techniques [39,40]. The linearized models are computationally efficient and fast. However, they are invalid for MS samples, because they are formulated on the assumption of weak-scattering approximation and do not involve any MSS techniques, which degrades the resolution of the reconstructed images of MS samples and limits the size of reconstruction to a single-cell level [41].

Recent studies on RI tomography have tackled MS using nonlinear scattering models, which calculate the MS process more accurately than linear models [42]. One approach iteratively solves the Lippmann–Schwinger (LS) equation, which describes scalar-wave propagation most accurately, to obtain scattering fields as close to the analytic solution as possible [43–46]. Another approach is the multislice method, which splits the sample into infinitesimally thin layers that sequentially interact with incident light [35,47–49]. These nonlinear approaches can provide accurate RI distributions. However, they involve considerably high computational costs owing to the complex nonlinear models used.

Herein, we present an approach, called MSS-Rytov, to realize RI tomography of MS specimens under weak-scattering approximation by leveraging the concept of MSS techniques using the conventional ODT system (Fig. 1). In MSS-Rytov, MSS techniques allow for the conversion of the distorted fields into the MS-suppressed smoothened fields that validate the weak-scattering assumption. The conventional formulation of Rytov approximation is not formulated against these MS-suppressed fields. Therefore, to reconstruct an RI distribution from these MS-suppressed fields, we formulate MSS-Rytov such that the MSS operations are compatible with the weak-scattering assumption. Although MSS-Rytov is inferior to the nonlinear approaches in terms of reconstruction accuracy, our approach achieves a linear reconstruction of the RI distribution with greater computational efficiency than the nonlinear approaches even in the presence of MS. In this sense, MSS-Rytov prioritizes computational cost over reconstruction accuracy. In addition to MS, to address the diversity of the MS and SIA strengths of biological specimens, we introduce three MSS-Rytov schemes that involve different MSS operations. These MSS operations are as follows: (1) coherent accumulation of fields, (2) incoherent accumulation of fields, and (3) an intermediate method between (1) and (2). The third method is introduced to balance SIA and MS tolerances.

 

Fig. 1. Principle of MSS-Rytov. MSS-Rytov utilizes a conventional ODT setup. A plane wave is illuminated on an MS sample, resulting in a scattered field. The resulting field was measured using an interferometer. MSS-Rytov first suppresses MS using a coherent or incoherent accumulation of fields and yields MS-suppressed fields. Then, the RI distribution is reconstructed from these MS-suppressed fields under weak-scattering approximation.

Download Full Size | PDF

2. Multiple scattering suppressive Rytov approximation for RI tomography of multiple-scattering specimens

2.1 Conventional linear model

A biological sample with a refractive index $n({\boldsymbol r} )$ can be modeled by its scattering potential $V({\boldsymbol r} ): = k_\textrm{b}^2\; [{{{({n({\boldsymbol r} )/{n_\textrm{b}}} )}^2} - 1} ]$, where ${k_\textrm{b}}: = 2\pi {n_\textrm{b}}/{\lambda _0}$ denotes the wavenumber in the medium, ${\lambda _0}$ is the wavelength in the vacuum, and ${n_\textrm{b}}$ is the refractive index of the medium. In ODT, incident waves are angularly scanned, and the resulting diffracted fields are detected by an interferometer. For an incident plane field ${u_{\textrm{in}}}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} ): = {a_0}\textrm{exp}({\textrm{i}{\boldsymbol k}_{\textrm{in}}^j \cdot \,{\boldsymbol r}} )$ with the jth incident wavenumber ${\boldsymbol k}_{\textrm{in}}^j$ ($j = 1 \ldots {N_{\textrm{in}}}$) $({N_{\textrm{in}}}$ is the total number of incident angles) and the amplitude ${a_0}$, the resulting total field $u({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )$ is described by the following LS equation:

2.2 MSS-Rytov linear model

MSS-Rytov solves this problem by giving MS tolerance to Rytov approximation with the use of spatial MSS techniques. Furthermore, in addition to MS, we need to consider SIA due to the curved surface of spheroids. Considering this issue, we describe three methods that utilize different MSS techniques. we first describe coherent MSS-Rytov (cMSS-Rytov) that utilize the coherent accumulation of fields, incoherent MSS-Rytov (iMSS-Rytov) using the incoherent accumulation of fields, and finally partially coherent MSS-Rytov (pMSS-Rytov) that is formulated in a compromised manner between iMSS-Rytov and cMSS-Rytov to balance the MS and SIA tolerances.

First, we derived cMSS-Rytov, which is a linear reconstruction model involving the coherent accumulation of fields. Recently, collective accumulation of single scattering (CASS) microscopy [14,50], which utilizes angle-scanned data obtained via interferometric measurements, has realized deep tissue imaging by the coherent accumulation of fields. This method selectively enhances single-scattering components by collectively accumulating components with the same momentum changes in all captured data. This operation realizes a computational spatial gating effect equivalent to confocal microscopy. Although CASS was demonstrated in the reflection mode, this method is potentially valid for transmission modes such as in ODT. Applying this to the resulting fields $u({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )$ of ODT, we obtain an MS-suppressed field:

Next, we derived iMSS-Rytov, which involves the incoherent accumulation of fields. To reconstruct quantitative RI distributions, we utilized computational QPGI calculation that involves the incoherent accumulation of fields. We calculate QPGIs from angle-scanned data and relate them to $V({\boldsymbol r} )$. When QPGI microscopy utilizes a spatially incoherent light source, the resulting image, ${\bar{W}^{\textrm{iMSS}}}({\boldsymbol r} )$ is emulated as follows:

Finally, we derive pMSS-Rytov, which addresses the trade-off by balancing the characteristics of cMSS-Rytov and iMSS-Rytov. To this end, we propose an intermediate MSS method between the coherent and incoherent accumulation of fields. This is realized by the subtraction of MS components from ${\bar{W}^{\textrm{iMSS}}}({\boldsymbol r} )$. We split the total field $\bar{u}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )$ into single-scattered ${\bar{u}^\textrm{S}}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )$ and MS component ${\bar{u}^\textrm{M}}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )$: $\bar{u}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )= {\bar{u}^\textrm{S}}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )+ {\bar{u}^\textrm{M}}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )$. From Eq. (4), we obtain the following relationship [when ${N_{\textrm{in}}}$ is sufficiently large]:

Here, we neglected the cross terms of ${\bar{u}^\textrm{S}}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )$ and ${\bar{u}^\textrm{M}}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )$, as they interfere destructively. This equation indicates that ${\bar{W}^{\textrm{iMSS}}}({\boldsymbol r} )$ is the sum of the incoherent accumulation of single-scattered components and MS components. As MS components can be approximately calculated by ${\bar{u}^\textrm{M}}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )\cong \bar{u}({{\boldsymbol r};{\boldsymbol k}_{\textrm{in}}^j} )- {\bar{U}^{\textrm{cMSS}}}({\boldsymbol r} )$ using Eq. (2), we define

2.3 Reconstruction based on MSS-Rytov

Figure 2 illustrates the numerical process of RI tomography based on the three proposed schemes. Interferograms obtained were captured by angle-scanning illumination using the conventional ODT system. After the application of the standard phase reconstruction using a digital holographic technique [51] to the interferograms, complex fields of different angles of incidence were digitally propagated to the entire sample volume via angular spectrum propagation [35,52]. Therefore, we obtained volume data for each incident angle. The z stacked images were processed via one of the three MSS techniques at each z position. cMSS-Rytov calculated confocal QPI with the use of Eq. (2) and then reconstructed RI distributions by deconvolution using the kernel given by Eq. (3). iMSS-Rytov and pMSS-Rytov calculated QPGI or MS-subtracted QPGI. In iMSS-Rytov and pMSS-Rytov, two types of QPGI with orthogonal shear directions were calculated to compensate for missing information at the differentiation operation of QPGI. Then, RI distributions were then calculated by deconvolution using the kernel given by Eq. (5) (Supplementary Section 6). Here, two modifications to these deconvolution kernels are required for accurate reconstruction. First, the limited aperture of illumination and detection must be considered (Supplementary Section 5.2). Second, the Green’s function must be modified because the digitally propagated field is different from the solution of the LS equation. This is because digitally propagated fields are source-free, whereas the solution of the LS equation is generated from the scattering potential (Supplementary Section 5.1). These two modifications to the kernel enable the realization of quantitative RI tomography. We also present the pseudo code for calculating phase maps using the three proposed methods in Supplementary Section 7.

 

Fig. 2. MSS-Rytov flowcharts including three different schemes. Interferograms with different angles of incidence were converted to complex fields by holographic processing and then digitally propagated to obtain focus-scanned volumetric data for each angle of incidence. cMSS-Rytov sums up these complex images coherently according to Eq. (2) and calculates confocal QPI images. iMSS-Rytov incoherently accumulates the images through Eq. (4) and calculates QPGI images. pMSS-Rytov calculates MS-subtracted QPGI images according to Eq. (7). Finally, deconvolution yields RI distributions.

Download Full Size | PDF

3. Simulations

We evaluated the MS and SIA tolerances of our methods via numerical simulation and compared them with the results of conventional Rytov approximation. We used different simulation models to assess MS and SIA tolerances. In each simulation, the sample was impinged by 441 incident waves with different angles in a regular grid pattern in k-space, such that the maximum illumination $\textrm{NA}$ was 1.0, and the diffracted light was detected by an objective lens of $\textrm{NA} = 1.0$. In this study, to improve the reconstruction accuracy and to mitigate the missing cone problem [40], we imposed a non-negativity constraint on the reconstructed scattering potentials of all methods, and set the tuning parameter $\gamma $ in pMSS-Rytov as 0.5. All processing was performed using an NVIDIA QUADRO RTX 8000 GPU and an Intel i9-9980XE CPU installed on a desktop computer. We also evaluated the resolving power of our approach to visualize subcellular structures of cells using a cell phantom in Supplementary Section 8.

3.1 Evaluation of multiple-scattering tolerance

For evaluating the MS tolerance, we applied our methods on a multiple-scattering simulated sample consisting of 591 spheres (Figs. 3(a) and (b)) within a volume of $144 \times 144 \times 144$ voxels with ${\lambda _0}/5$ resolution. The RI and diameter of each sphere were given by $1.390\,\textrm{and}$ 2${\lambda _0}$. The spheres are positioned slightly off the body-centered cubic lattice with a lattice constant of $4{\lambda _0}$. The displacements were given by Gaussian distribution with a standard deviation of $0.4{\lambda _0}$ in each direction to suppress the periodicity of the sample. We calculated the resulting field from the cell using the forward model of sparsity-driven image reconstruction under multiple scattering (SEAGLE) [44] in the same manner as an existing study [49]. This method iteratively solves the LS equation and provides accurate diffracted fields. The simulated fields with an illumination angle $({\theta _x},\; {\theta _y}) = ({0^\circ ,\,0^\circ } ),\,({18.2^\circ ,\,18.2^\circ } )$ are shown in Fig. 3(c), where the wavefront of simulated fields are distorted because of the MS. Figure 3(d) shows the reconstructed RI and reconstruction error of each method, where the reconstruction error is given by ${n_{\textrm{true}}} - {n_{\textrm{recon}}}$, where ${n_{\textrm{true}}}\; \textrm{and}\,{n_{\textrm{recon}}}$ denote the true and estimated RI distribution. We also quantitatively evaluated the performance of each method with the normalized reconstruction error (NRE) defined as

 

Fig. 3. Simulations for MS tolerance in conventional Rytov, cMSS-Rytov, iMSS-Rytov, and pMSS-Rytov. (a) Simulation condition. The incident fields first interact with the 3D cell phantom and then the diffusers. The diffuser works as a phase mask, and the MS strength depends on the number of diffusers. (b) x-y and x-z cross-sections of the ground truth of the RI map of the cell along the dotted line in (b). (c) Simulated fields with an illumination angle $({\theta _x},\; {\theta _y}) = ({0^\circ ,\; 0^\circ } )\; \textrm{and}\; ({18.2^\circ ,\; 18.2^\circ } )$. (d), x-y and x-z cross-sections of RI distributions and reconstruction error obtained by each method. Scale bars: 5${\lambda _0}$.

Download Full Size | PDF

3.2 Evaluation of sample-induced aberration tolerance

Subsequently, for the evaluation of SIA tolerance, we simulated scattering using a homogeneous sphere with various RIs. The diameter of the sphere is 10${\lambda _0}$. The spheres were within a volume of $96 \times 96 \times 96$ voxels with ${\lambda _0}/5$ resolution. We changed the RI of the sphere (${n_\textrm{s}}$) from 1.370 to 1.450; this range approximately represents cellular contents. The diffracted fields were calculated using SEAGLE. The ground truth and reconstructed results are shown in Figs. 4(a) and 4(b), respectively. We further show the cross-section of reconstructed results for ${n_\textrm{s}} = 1.370,\; 1.390,\; \textrm{and}\,1.410$ in Fig. 4(c). To evaluate the SIA tolerance of each method, we calculated the NRE for each method at various values of ${n_\textrm{s}}$ (Fig. 4(d)). Among all methods, cMSS-Rytov failed to reconstruct RI distribution for high ${n_\textrm{s}}$ and thus shows the largest error at high SIA. iMSS-Rytov and pMSS-Rytov successfully reconstructed RI distributions with all tested values of ${n_\textrm{s}}$, demonstrating higher reconstruction accuracy than conventional Rytov and cMSS-Rytov.

 

Fig. 4. Simulations for SIA tolerance in conventional Rytov, cMSS-Rytov, iMSS-Rytov, and pMSS-Rytov. (a) x-y and x-z cross-sections of the ground truth of a homogeneous sphere. (b) x-y, and x-z cross sections of the reconstructed RI distributions obtained by each method. From the left columns, the RIs of the spheres, ${n_\textrm{s}},$ are 1.370, 1.390, and 1.410. (c) 1D cross sections along the white dotted lines in (b). (d) NRE plot for each method at various values of ${n_\textrm{s}}$. Scale bars: 5${\lambda _0}$.

Download Full Size | PDF

4. Experiments

Since simulation results suggest that pMSS-Rytov has well-balanced MS- and SIA tolerances, the method is suitable for the sample with MS and SIA. In the following experiment, we utilized pMSS-Rytov to visualize multicellular spheroids, because they induce both MS and SIA. We constructed a standard ODT setup with a helium-neon laser (${\lambda _0} = $ 632.8 nm) as the light source. A plane wave was incident on the sample via a condenser lens (×60, $\textrm{NA} = 1.0$, water immersion) and angularly scanned with 441 angles in a regular grid pattern in k-space, such that the maximum illumination NA was 0.9. The diffracted fields were collected using an objective lens (×60, $\textrm{NA} = 1.0$, water immersion). The complete experimental setup and the procedure for biological sample preparation is provided in Supplementary Sections 10 and 11. To reduce the aberration in the optical system, numerical aberration estimation and correction were also performed, and the procedure is provided in Supplementary Section 12 in detail. All the experimental data presented in this manuscript were generated using aberration correction with an estimated phase map, as presented in Supplementary Section 12.

4.1 3D imaging of weak-scattering sample

First, we imaged a weak-scattering sample, the single MCF-7 cell. We prepared floated MCF-7 cells by trypsinization, not adherent ones, to better compare the results of its 3D cultured version (Section 4.2). Figures 5(a) and (b) show the reconstructed RI distribution by conventional Rytov and pMSS-Rytov, respectively. Visualization 1 is a video showing the whole cell. The reconstruction volume contains $200 \times 200 \times 76$ voxels with a voxel size $0.14 \times 0.14 \times 0.4\,\mathrm{\mu}{\textrm{m}^3}$. We also rendered MCF-7 cells for visualization in a 3D volume (Fig. 5(c)). As this sample is a weak-scattering object, both conventional Rytov and pMSS-Rytov approximations are valid. As can be observed in Figs. 5(a) and (b), both methods successfully visualized the cell. The total reconstruction times were 6.0 and 1.8 s, respectively. This difference is mainly because pMSS-Rytov skips phase unwrapping. It required 4.7 s for conventional Rytov to perform phase unwrapping because phase unwrapping was performed on the CPU.

 

Fig. 5. 3D RI reconstruction of a single MCF-7 cell using conventional Rytov and pMSS-Rytov. (a) x-z cross-section along the white dotted line in (b). (b) x-y cross-section along the white dotted line in (a). In this case, as the sample is a weak-scattering object, both methods succeeded in imaging RI distributions of the cell (see Visualization 1). (c) 3D rendering of the cell reconstructed by pMSS-Rytov. Scale bars: $5\,\mathrm{\mu}$m.

Download Full Size | PDF

4.2 3D imaging of multiple-scattering sample with sample induced aberration

Next, as an MS sample, we imaged a multicellular MCF-7 spheroid. Figures 6(a) and (b) show the reconstructed RI distribution by conventional Rytov and pMSS-Rytov, respectively. Visualization 2 is a video showing the whole spheroid. The reconstruction volume contains $470 \times 470 \times 150$ voxels with a voxel size $0.14 \times 0.14 \times 0.4\,\mathrm{\mu}{\textrm{m}^3}$. Notably, the total time required for reconstruction was almost the same for the two methods, 27.3 s in conventional Rytov and 18.5 s in pMSS-Rytov. Here, conventional Rytov took 12.5 s to perform phase unwrapping. These results demonstrate that pMSS-Rytov is also a computationally efficient algorithm and almost comparable to conventional Rytov, despite the superior imaging capability of MSS-Rytov. We also rendered the spheroid for visualization in a 3D volume (Fig. 6(c)). Because a multicellular spheroid is not a weak-scattering object, conventional Rytov failed to reconstruct high-frequency content such as particle-like structures as expected. We also provide captured complex amplitudes and reconstructed phase maps in Supplementary Section 13 and 14, respectively. The complex amplitudes show a great deal of discontinuity in its phase. These discontinuities lead to loss of image quality in conventional Rytov. In contrast, pMSS-Rytov succeeded in visualizing organelles [53,54], which appeared similar to nucleoli and mitochondria, even inside the spheroid. Thus, we confirmed the superior imaging performance of pMSS-Rytov for imaging MS objects compared to conventional Rytov-based ODT.

 

Fig. 6. 3D RI reconstruction of a multicellular MCF-7 spheroid using conventional Rytov and pMSS-Rytov. (a) x-z cross-section along the white dotted line in (b). (b) x-y cross section along the white dotted line in (a). pMSS-Rytov succeeded in visualizing some vesicles (blue arrowheads) and organelles, which appeared similar to nucleoli (yellow arrowheads) and mitochondria (red arrowheads). In contrast, conventional Rytov approximation failed to visualize these organelles. (see Visualization 2) (c) 3D rendering of spheroid reconstructed by pMSS-Rytov. Scale bars: $10\,\mathrm{\mu}$m.

Download Full Size | PDF

4.3 Quantification of lipid accumulation in 3D-cultured hepatocytes

Finally, to demonstrate the biomedical application capability of MSS-Rytov, spheroids consisting of human hepatocytes, a major liver component, were imaged. A liver is a central organ for metabolism in vertebrates, with more than 500 vital functions and unrevealed complexities [55,56]. This implies that spheroids and organoids of liver cells intrinsically contain a wide variety of RI distributions and their dynamics, which could be one of the ideal inherent contrast agents for visualization using RI tomography. As a representative of 3D-cultured hepatocytes, we selected human hematoma HepG2 spheroids. The HepG2 cell line was established from the liver biopsy of a 15-year-old Caucasian male patient [57] and these cells have been extensively used in metabolism and hepatotoxicity research [58]. 2D-cultured HepG2 cells can be induced to accumulate lipids by adding free fatty acids to the culture medium and have been widely used in lipid metabolism assays [59,60]. We applied this method to 3D-cultured HepG2 spheroids and imaged them under two different culture conditions: normally cultured and oleic acid treated. Figure 7 shows reconstructed RI distributions of HepG2 spheroids by pMSS-Rytov. More lipid droplets can be observed in the oleic-acid-treated spheroid, and the droplet sizes are larger in the oleic-acid-treated spheroid than in the normally cultured spheroid (Figs. 7(a), (b), (d), and (e)). We also rendered the spheroids for visualization in 3D volume (Figs. 7(c) and 7(f)). We estimated the lipid concentration from the reconstructed RI distributions. We set the region of the spheroid as $n({\boldsymbol r} )> 1.34$ and the region of the lipids as $n({\boldsymbol r} )> 1.385$, and estimated the lipid concentration, ${C_{\textrm{lipid}}}$, from the ratio of these volumes: ${C_{\textrm{lipid}}} = \mathop \sum \nolimits_{n({\boldsymbol r} )> 1.385} \; /\mathop \sum \nolimits_{n({\boldsymbol r} )> 1.34} $. The results were ${C_{\textrm{lipid}}}$ = 2.2% for normally cultured spheroids and 4.9% for oleic-acid-treated spheroids. Therefore, we confirmed an increase of 2.23 times in the lipid concentration in the oleic-acid-treated spheroids compared to those normally cultured. We also provide the RI distributions of these HepG2 spheroids reconstructed by conventional Rytov in Supplementary Section 15. The result demonstrates that conventional Rytov is unable to visualize lipid droplets in the hepatocytes that can be captured by pMSS-Rytov.

 

Fig. 7. 3D RI reconstruction of a multicellular HepG2 spheroid under normal culture conditions and oleic acid induction reconstructed by pMSS-Rytov. (a-c) MIPs in two different planes, x-y cross section along the white dotted line in (a), and 3D rendered image of the normally cultured HepG2 spheroid. (d-f) MIPs in two different planes, x-y cross section along the white dotted line in (d) and 3D rendered image of the oleic-acid-treated HepG2 spheroid. More lipid droplets can be observed in the oleic-acid-treated spheroid, and the droplet sizes are larger in the oleic-acid-treated spheroid than in the normally treated spheroid. Scale bars: $10\,\mathrm{\mu}$m.

Download Full Size | PDF

5. Discussion

We proposed a RI tomography approach, MSS-Rytov, for the reconstruction of MS objects under weak-scattering approximation. This approach was realized by reformulating the ODT theory under weak-scattering approximation to ensure compatibility with MSS techniques. The conventional Rytov approximation collapses when the field distortions induced by MS are significant, whereas MSS-Rytov overcomes this limitation by suppressing MS using MSS techniques. Therefore, MSS-Rytov expands the applicability of this approximation from weak-scattering objects to MS objects. We proposed three MSS-Rytov schemes to address the diversity of the MS and SIA strengths of biological samples. cMSS-Rytov utilizes the coherent accumulation of fields to produce a confocal QPI before 3D reconstruction and exhibits the highest MS tolerance among the three methods at the expense of SIA tolerance. iMSS-Rytov uses iMSS techniques to generate a QPGI and exhibits the highest SIA albeit lowest MS tolerance. pMSS-Rytov balances this trade-off, and numerical simulations demonstrate that pMSS-Rytov is the most suitable method for observing biological samples with MS and SIA. To demonstrate the proposed method, we visualized a multicellular spheroid using pMSS-Rytov and observed subcellular structures that cannot be observed using conventional Rytov approximation, confirming that pMSS-Rytov has superior imaging performance over conventional Rytov approximation with a comparable computation time (Figs. 6, 7, and S3).

To demonstrate our proposed approach has greater computational efficiency than nonlinear approaches, we consider the computational complexity of each approach. As a nonlinear approach, we consider a method based on the beam propagation method (BPM) [47,49] because it is the most computationally efficient nonlinear approach to the best of our knowledge. We consider a space with $N \times N \times N$ voxels. We also define ${N_{\textrm{angle}}}$ and ${N_{\textrm{iter}}}$ as the number of incident angles and the number of optimization iterations, respectively. The computational complexity of our and BPM-based methods are given $O({{N_{\textrm{angle}}}{N^3}\textrm{log}N} )+ O({{N_{\textrm{iter}}}{N^3}\textrm{log}N} )$ and $O({{N_{\textrm{iter}}}{N_{\textrm{angle}}}{N^3}\textrm{log}N} )$, respectively (detailed derivation of the complexity is included in Supplementary Section 16). As ${N_{\textrm{iter}}}$ and ${N_{\textrm{angle}}}$ typically take values between 50 and 500, our approach is expected to be faster than the nonlinear approach. As explained in Section 4.2, pMSS-ODT required 18.5 s to reconstruct $470 \times 470 \times 150$ voxels, whereas the BPM-based approach required 1.5 h to reconstruct $450 \times 450 \times 150$ [49]. This result demonstrates that the proposed approach can enable computationally efficient RI reconstruction of MS specimens.

The other schemes, i.e., cMSS-Rytov and iMSS-Rytov, could also be valuable for biomedical research although their experimental data were not included in this study. Simulations indicated that cMSS-Rytov has a strong MS tolerance but weak SIA tolerance. This trade-off can be overcome by adjusting the RI of the sample environment [20], which could assist in leveraging the strongest MS tolerance for cMSS-Rytov. iMSS-Rytov is expected to be beneficial for visualizing densely spread samples throughout the field of view of objective lenses with no background, such as cell sheets. In standard ODT, the absolute or relative phases between captured fields with different incident angles must be corrected using the background. In contrast, iMSS-Rytov can reconstruct RI distributions without the phase corrections, because these phases vanish in the incoherent summation given by Eq. (4). This advantage would help expand the potential of RI tomography to cell sheets and samples confluent in the field of view, although resulting RI values may be relative owning to the loss of the indefinite constant during differential operations.

To observe thicker and more complex specimens than the multicellular spheroids we visualized here, we must overcome the trade-off between MS and SIA because we solely provide a compromised solution as pMSS-ODT in this study. This is a limitation of our approach. Although adjusting the RI of the sample environment mentioned above is a conceivable solution, this type of adjustment may limit the applicability of RI tomography. In the deep optical imaging field, several label-free methods significantly improve their imaging depths through numerical compensations up to higher-order modes of aberration [50,61]. Integrating these techniques and ODT theory may help fundamentally overcome the trade-off.

Another limitation of the proposed approach is its image fidelity. The reconstructed tomograms show elongation along the z-axis (Fig. 3). In particular, the results of pMSS-Rytov and iMSS-Rytov show greater elongation than those of conventional Rytov and cMSS-Rytov. This effect is known as the missing cone problem, and has been addressed by numerous studies [40]. The reason iMSS-Rytov and pMSS-Rytov are more susceptible to the missing cone problem is that QPGI images are used in their reconstruction, as explained in Section 2.3. The optimization in the reconstruction minimizes the squared error between QPGI form experimental data and that calculated from the RI distribution. As the squared error is susceptible to a large difference (outliers) and QPGI is given by the derivative of QPI, the fitting error is integrated and has a greater effect on the resulting RI distribution, leading to elongated distribution. Previous studies using conventional Rytov approximation have demonstrated that this problem can be alleviated by regularization techniques using a priori information [39,40,62]. Therefore, incorporating these regularization techniques in our methods would improve their fidelity, even for MS specimens.

Further, we discuss how our method is superior and inferior to recent nonlinear model-based approaches. Previous approaches [48] based on nonlinear models show better reconstruction accuracy than the conventional Rytov approach because these nonlinear approaches significantly alleviate the missing cone problem. As MSS-Rytov also suffers from this problem, our approach is inferior to nonlinear approaches in this respect. However, MSS-Rytov can provide quantitative RI distributions of MS specimens with subcellular resolution in a computationally efficient manner using a linearized model. This is one of the advantages of our approach.

As a biomedical application of MSS-Rytov, we demonstrated the visualization of HepG2 spheroids as an example of the human hepatocyte model (Fig. 7). We successfully quantified lipid droplets even inside the spheroids, thereby proving that MSS-Rytov can be applied to evaluate 3D fatty liver models [63]. The demand for 3D cell cultures has recently increased in various biomedical fields, such as drug development, and regenerative medicine [64–66]. MSS-Rytov has the potential to conduct label-free quantitative analyses of multicellular spheroids (Figs. 6 and 7), and its visualizing ability and computational efficiency would benefit a wide range of spheroid- and organoid-based assays.

6. Conclusion

We proposed an approach termed MSS-Rytov. MSS-Rytov was formulated such that Rytov approximation is compatible with MSS techniques. Thanks to this formulation, MSS realized computationally efficient reconstructions of MS samples, overcoming the conventional limit of Rytov approximation. In addition, to visualize multicellular spheroids, which caused MS and SIA, we introduced a compromised solution, pMSS-Rytov, to balance MS and SIA tolerance. Our method requires a conventional ODT setup without any modification, hence, prevalently used ODT systems can immediately acquire the imaging capability of MS objects by applying MSS-Rytov. Moreover, the proposed approach can be applied to other wave diffraction tomography techniques such as acoustic and microwave tomography [67]. We demonstrated substantial improvements in MSS-Rytov over conventional Rytov-based ODT through simulations and experiments. We believe that MSS-Rytov can open a new direction for RI tomography of MS specimens and become a standard imaging modality in observing multicellular samples in life science and biomedicine.