Multiple-scattering suppressive refractive index tomography for label-free quantitative assessment of multicellular spheroids

MCF-7 further

Multiple-scattering suppressive refractive index tomography for label-free quantitative assessment of multicellular spheroids: supplemental document

Conventional formulation of Rytov approximation
The conventional formulation of ODT using Rytov approximation is described here. We first derive the Born approximation and then the conventional Rytov approximation using perturbation theory. The formula provided here is based on a previous study [1]. We introduce a dimensionless parameter and expand ; in as the power series in : ( ; in ) = 0 ( ; in ) + 1 ( ; in ) + 2 2 ( ; in ) + ⋯ = ; in . (S1) We treat the scattering potential ( ) as a perturbation: ( )→ ( ). In this case, on substituting Eq. (S1) in Eq.
(S2) by comparing the coefficients of each power of .
Order 0 : At order 0 , we have 0 ; in = in ; in .
We then equate Eq.
(S2) and (S6), solve it by comparing coefficients of each power of .
Order 1 : At order 1 , we have The Rytov approximation of Eq. (1) is given by the first-order approximation of Eq. (S6). We set →1 and obtain In this conventional formulation, the Rytov approximation is applied to coherent fields ( ; in ). The validity of Rytov approximation depends on the smoothness of the fields. Further, as field ( ; in ) is vulnerable to MS owing to its high coherence, this conventional formulation becomes invalid when the sample shows MS.

Formulation of cMSS-Rytov
Here, we formulate cMSS-Rytov, such that the theory of ODT is compatible with the coherent MSS technique used as in CASS. We expand ; in as a power series in : ( ; in ) = 0 ( ; in ) + 1 ( ; in ) + 2 2 ( ; in ) + ⋯ = ; in . (S10) From Eq. (S3) and (S4), using the relationships ; in ≔ ; in exp -i in ⋅ and in ; in ≔ 0 exp (i in ⋅ ), we obtain the following equations: We expand cMSS ( ) as a power series in : Substituting Eq. (S13) and (S10) in Eq. (2), we obtain We can now solve Eq.
(S14) by comparing coefficients of each power of .
The Born approximation of Eq. (2) is given by the first-order approximation of Eq.(S13). We set →1 and obtain cMSS Born ( ): Next, we derive the Rytov approximation of Eq. (2). We expand cMSS ( ) as a power series in , as follows: We then equate Eq. (S13) and (S18) and solve it by comparing coefficients of each power of .
Order 0 : At order 0 , we have Order 1 : At order 1 , we have The Rytov approximation of Eq. (2) is given by the first-order approximation of Eq. (S18). We set →1 and obtain cMSS Rytov ( ): cMSS Rytov ( )≔ exp( 0 ( ) + 1 ( )) In this formulation, we applied the Rytov approximation to coherently MS-suppressed fields cMSS ( ). As cMMS techniques suppress MS through coherent summation of fields ; in , cMSS ( ) is spatially smoother than ; in . Based on this characteristic, Rytov approximation can be applied to cMSS ( ) even if some MS exists. Therefore, the validity of the Rytov approximation expands from weak-scattering samples to MS samples owing to this formulation.

Formulation of iMSS-Rytov
We formulate iMSS-Rytov such that the theory of ODT is compatible with the iMSS techniques by using QPGI. QPGI is given by the phase of iMSS ( ) in Eq. (4). We expand iMSS ( ) as a power series in : Substituting Eq. (S10) and (S22) in Eq. (4), we obtain We can now solve Eq. (S23) by comparing coefficients of each power of .
This relationship between cMSS ( ) and ( ) is the same as that between iMSS ( ) and ( ). pMSS ( ) lies in between cMSS ( ) and iMSS ( ), such that pMSS ( ) also has the same relationship with ( ):

Source-free Green's function
The Green's function must be modified because a digitally propagated field is different from the solution of the LS equation. The solution of the LS equation, ; in , has some sources in the field: ∇ 2 + 2 b ; in = -( ) ( ; in ). However, we use backpropagated fields from a detector to reconstruct the RI map, and the backpropagated field is source-free. Therefore, we need to determine the modified Green's function, which gives the relationship between the backpropagated field and the scattering potential.
We define the field as B ; in , which is the backpropagated field of all the outgoing components of ; in . Note that B ; in must be a source-free field: To derive the modified LS equation for B ; in , we consider the time-reversed situation in Eq.
(1) (Figs. S1(a), (b)). In this time-reversed situation, the resulting field has the same spatial distribution as ; in but the opposite time-development direction because of the time symmetry of the wave equation. In other words, when we set the incident wave as B * ; in , the total field resulting from scattering potential * ( ) is equal to * ; in : The complex conjugate operation corresponds to a time-reversal operation. Taking the complex conjugate of this equation and equating it with Eq.
(1), we obtain where we assume that ( ) is a real value. This assumption generally holds when the sample is a biological specimen. Here,

Limited aperture of objective lens
We must consider the limited aperture of an objective lens. The 3D pupil function of the objective lens, whose numerical aperture is , is given as follows: where ( ) = 1 if ≤ 1 and 0 otherwise, and ( ) = 1 if ≥ 0 and 0 otherwise. From Eqs. (S34) and (S35), we define the modified Green's function m : where ℱ -1 is the inverse 3D Fourier transform operator.

cMSS-Rytov
cMSS-Rytov reconstructs the RI map from confocal QPI images by deconvolution, whose kernel is given by the summation of modified Green's functions (Supplementary Section 5).
Here, let G, v, and ϕ be the kernel given by Eq. (S41), vector representation of the scattering potential, and vector representation of the 3D confocal QPI image, respectively. We also introduce a cropping operator [2] D to avoid aliasing, as the kernel is heavy-tailed. The estimation of v from ϕ is given by the following optimization problem: This problem was solved using the alternating direction method of multipliers (ADMM) [3].
The program was stopped after 50 iterations.

iMSS-Rytov and pMSS-Rytov
iMSS-Rytov and pMSS-Rytov reconstruct the RI distribution from QPGI images by deconvolution. Here, let G d , G d , ψ , and ψ be the discretized derivative of the 3D Green's function with respect to x, the discretized derivative of the 3D Green's function with respect to y (Eq. (S41)), vector representation of QPGI differentiated with respect to x, and vector representation of QPGI differentiated with respect to y, respectively. Unlike QPI, QPGI does not have a constant component owing to its differentiation operation; hence, the estimated v is also uncertain about its constant value. To address this problem, the estimation of v was performed using the following optimization problem: The L1-norm term was introduced to limit the range of v to near-zero values. This term and the nonnegativity constraint help the solution converge to the accurate solution. If is too large, the estimated RI distribution can be different from the exact values. We chose = 0.1 for simulation and 0.2 for experiments. The shear amount of QPGI is set as 0 /7.5 in simulations and 420 nm in experiments, respectively. The optimization problem was solved using ADMM, and the program was stopped after 50 iterations.

Pseudo code for the reconstruction of phase maps through the cMSS, iMSS and pMSS techniques
Tables S1, S2 and S3 present the pseudo-code for the reconstructions of phase maps ∠ cMSS ( ) , ∠ iMSS ( ) and ∠ pMSS ( ) from measured complex fields, respectively.    for ←1 layer do 2.

Numerical simulation of cell phantom
To validate that our approach visualizes the subcellular structure of cell samples, we applied our methods to a cell phantom immersed in water (RI = 1.333). The cell phantom ( Fig. S2(a)

Comparison of forward models
We evaluated the forward model accuracy of each model using the model in Fig. 3. The left column of Fig.S3(a) shows true fields calculated using the SEAGLE forward model.
; in was calculated using the SEAGLE forward model, and ∠ cMSS ( ), ∠ iMSS ( ), and ∠ pMSS ( ) were calculated from ; in using Eqs. (2), (4) and (7), respectively. The right column of Fig.S3(a) Fig.S3(b) shows the error of these fields. We also quantitatively evaluated the performance of each method with the normalized data fit (NDF) defined as where true is the true field calculated using results of SEAGLE and recon is calculated using each forward model. Hence, ( true , recon ) = ; in , Rytov ; in , The NDF of conventional Rytov, cMSS-Rytov, iMSS-Rytov and pMSS-Rytov was 0.152, 0.004, 0.144 and 0.086, respectively. According to the NDF, the three proposed methods have higher forward-model accuracy than conventional Rytov.

Experimental setup
We constructed a transmission-mode interferometric microscope with a Mach-Zehnder configuration (Fig. S4). We used a He-Ne laser (632.8 nm, 5 mW) as the light source. The light was coupled into a single-mode fiber and split into a sample beam and a reference beam using a beam splitter (BS1). The sample beam was illuminated onto the sample at various angles of incidence using a 2D scanning mirror (Optotune). We set the scanning pattern as a grid pattern, such that the numerical aperture of incidence equaled 0.9. The sample was placed between a condenser lens (OB1, Olympus, 60×, 1.0NA, water immersion) and an objective lens (OB2, Olympus, 60×, 1.0NA, water immersion). The scattered light from the sample was projected onto a CMOS camera (GO-2400M-USB, JAI) using a 4f system consisting of OB1 and a tube lens (f = 125 mm) after combining with the reference beam at another beam splitter (BS2). On the camera plane, the sample beam interfered with the reference beam, which was slightly tilted against the propagation direction of the sample beam to generate off-axis interferograms. To reconstruct a 3D RI tomogram, we recorded 441 holograms of the sample at a frame rate of 62 fps, which required 7 s.

Preparation of cellular samples
The human breast cancer cell line, MCF-7, was obtained from the European Collection of Cell Cultures and cultured in an advanced minimum essential medium (Gibco) supplemented with 2.5% (v/v) fetal bovine serum (FBS, Gibco), 100 U/mL of penicillin, 100 µg/mL of streptomycin, and 292 µg/mL of L-glutamine (1× Penicillin-Streptomycin-Glutamine, Gibco). The human hepatoma HepG2 cell line (JCRB1054) [4] was obtained from the Japanese Collection of Research Bioresources and cultured in Dulbecco's modified Eagle's medium (DMEM, Gibco) supplemented with 10% (v/v) FBS (Gibco) and 1× penicillin-streptomycinglutamine (Gibco). Multicellular spheroids were formed using a 3D cell culture container, the EZSPHERE 6-well plate (AGC Techno Glass). MCF-7 cells were seeded at 35,000 cells per well in 2 mL of the culture medium and incubated for 1 d to allow spheroid formation. HepG2 cells were seeded at 30,000 cells per well in 2 mL of the culture medium and incubated for 2 d to form spheroids. To induce lipid-droplet formation in HepG2 spheroids, the medium was replaced with phenol-red-free DMEM (Gibco) containing 1% (w/v) fatty-acid-free bovine serum albumin (BSA, FUJIFILM Wako Pure Chemical) and 0.5 mM sodium oleate (Nacalai Tesque) 1 d after seeding, and the spheroids were further cultured for 1 d. All cell samples were cultured at 37 ℃ in 5% CO2 humidified atmosphere. The spheroids were washed in phosphate-buffered saline (Gibco) or Hanks' balanced salt solution (Gibco) and collected via centrifugation at 150 × g for 3 min.

Numerical aberration estimation
In general, the aberration of an optical system is expressed as a phase map in the pupil plane. In most cases, the phase map is expressed as the sum of orthogonal modes, such as Zernike polynomials. Our strategy to identify the system aberration is to automatically estimate the coefficients of Zernike polynomials that maximize a sharpness metric calculated from aberration corrected images. As the metric to determine the coefficients, we used QPGI image given by Eq. (4), which is given by the incoherent summation of the angle-scanned complex fields. We define the sharpness metric using QPGI as follows: Using this metric, we estimated 4th to 28th-order coefficients of Zernike polynomials. The order of Zernike polynomials follows Noll's sequential indices [5]. To automatically optimize the coefficients, we used the particle swarm optimization (PSO) algorithm [6]. The number of particles of PSO was set to 100, and the particle parameters were initialized with a Gaussian distribution with standard deviation =1. Two-hundred iterations of PSO were performed. Polystyrene beads were used as the target image for the optimization. Figures S5(a) and (b) show the estimated aberration map and coefficients of Zernike polynomials, respectively. Figure S5(c) shows the QPGI image of a polystyrene bead used as the target image before and after aberration correction, respectively. Compared to the former, the latter is sharper and has higher peak phase values. This result indicates that aberration correction was successful. Using the estimated aberration map, we also performed aberration correction to the image of an MCF-7 cell used in Fig

RI distributions of HepG2 reconstructed by conventional Rytov
We present the RI maps of normally cultured and oleic-acid-treated HepG2 spheroids obtained by conventional Rytov approximation (Fig. S8). According to Fig. S8, conventional Rytov approximation failed to visualize lipid-droplet structures that can be observed by pMSS-Rytov, as shown in Fig. 7.

Computational cost of our method and BPM-based approach
First, we consider the complexity of the BPM-based approach. In the BPM case, 2 2D-FFT is required for end-to-end propagation for each angle of incidence. Therefore, for the endto-end propagation for all angles of incidence, it requires a complexity of ( angle 3 log ). The same complexity is necessary for backpropagation. The BPM-based approach repeats this computation iter times for optimization. As a result, it requires a complexity of ( iter angle 3 log ). Next, we consider the complexity of our approach, which first calculates a multiplescattering suppressed image from angle-scanned data. Similar to the forward process of the BPM approach, it requires a complexity of ( angle 3 log ). Subsequently, we calculate the RI distribution for the multiple-scattering suppressed image using ADMM. In our implementation of pMSS-Rytov, 7 3D-FFT are required for each iteration, requiring a complexity of ( iter 3 log ) for all iterations. As a result, the complexity of our approach is angle 3 log + ( iter 3 log ).