Optical dispersions through intracellular inhomogeneities

Transport of intensity equation (TIE) exhibits a non-interferometric correlation between intensity and phase variations of intermediate fields (e.g., light and electron) in biological imaging. Previous TIE formulations have generally assumed a free space propagation of monochromatic coherent field functions crossing phase distributions along a longitudinal direction. Here, we modify the TIE with fractal (or self-similar) organization models based on intracellular refractive index turbulence. We then implement the TIE simulation over a broad range of fractal dimensions and wavelengths. Simulation results show how the intensity propagation through the spatial fluctuation of intracellular refractive index interconnects fractal-dimensionality with intensity dispersion (or transmissivity) within the picometer to micrometer wavelength range. In addition, we provide a spatial-autocorrelation of phase derivatives which allows the direct measurement and reconstruction of intracellular fractal profiles from optical and electron microscopy imaging.

Introduction. Biological cells are the basic structural and functional unit of life. All cells consist of a cytoplasm enclosed within a membrane, which includes biomolecules (e.g., proteins and nucleic acids) and intracellular organelles (e.g., the endoplasmic reticulum, Golgi apparatus, and mitochondria). These complex intracellular structures can be visualized using optical imaging systems and electron microscopes. As visible light has a wavelength ranging from 360 nm to 760 nm, which are up to five orders of magnitude longer than wavelength of electrons, optical imaging techniques tend to filter out the inner structures of smaller intracellular objects otherwise observable through electron microscopy.
The complex amplitude of a wave function describing optical fields is generally composed of two major physical quantities: amplitude and phase functions [1]. Variations in amplitudes (or intensities) can be captured directly by photosensitive detectors such as the human eye and scientific CMOS cameras. Variations in phases, however, are barely visible using optical microscopy systems and are sensitive to light scattering and fluctuations due to invisible spatial distributions of intracellular optical properties, e.g., refractive index and optical thickness. These phase variations have been formalized into a partial differential equation, termed transport of intensity equation (TIE) [2][3][4], which is dependent on wavelength and the intensity variations of intermediate fields (i.e., photon and electron). This equation represents, in particular, the free space propagation of a monochromatic coherent wave function to a given phase distribution along a longitudinal direction.
A key challenge to bridging the large resolution gap * Corresponding author: masaki@riken.jp that lies between optical and electron microscopy techniques consists of finding meaningful variations in both intensity and phase functions. Over a broad range of wavelengths, intensity and phase variations of intermediate fields are inextricably linked via the TIE. Previous TIE formalizations have generally assumed that a complex plane wave function describing intermediate fields can be propagated through an arbitrary phase distribution in free space. Intensity transport of the wave functions interacting with the complex refractive aspects of cellular interiors, however, remains elusive. In this article, we modify the TIE with fractal (or self-similar) organization models for intracellular refractive index turbulence. We then implement the TIE simulation in the picometer to micrometer wavelength range, computing intensity dispersions and reductions caused by inhomogeneous spatial distributions of intracellular refractive index. Crucially, we show how the a priori lateral intensity distribution of a cell can be described as a function of its fractal dimensions. Such fractal media is optically thin and transparent through the visible wavelengths but exhibits the compression of fractal-dimensionality and the increase of intensity dispersions in the wavelength range of electron microscopy imaging. We also derive a spatialautocorrelation function that interconnects phase derivatives with intracellular fractal profiles and demonstrate the direct measurements and reconstruction of optical properties from fluorescent cell imaging. Our work facilitates further extensions to bioimage simulation modules [5][6][7][8][9][10][11][12][13][14], incorporating, in particular, light scattering and fluctuations caused by nonuniform spatial distributions of intracellular optical properties. Such modification of the simulation is of particular relevance to the numerical evaluation and verification of observational invariance (or symmetry) as postulated in data science, likely leading to more realistic simulations of biological imaging. Transport of intensity equation. For a specific intermediate field wavelength λ (e.g., photon or electron) we reformulate the TIE for intensity propagation by using the laws of conservation of energy (i.e., Helmholtz equation), parameterizing, in particular, a complex refractive index for the intracellular media, where n 0 and ∆n(r ⊥ , z) are the mean and spatial fluctuation of the refractive index, respectively. κ(r ⊥ , z) in the imaginary part denotes the spatial distribution of the attenuation index. In this parameterization, we assume that a monochromatic coherent wave function propagating through the refractive index distribution along the axial axis (see Figure 1) can be written in the form of where the square root of the intensity distribution √ I denotes the amplitude of the wave function; k is the wave number in free space k = 2π/λ. Substituting Eq. (1) and (2) into the Helmholtz equation, we can derive TIE as follows.
where ∇ ⊥ is the two-dimensional Nabla operator in the transverse direction, i.e., ∇ ⊥ = (∂/∂x, ∂/∂y). α(r ⊥ , z) and β(r ⊥ , z) represent amplitude and phase variations arising from the spatial fluctuations of the refractive and attenuation index, and where β-factor vanishes in the TIE modification. This factor may be more relevant to phase changes in transport of phase equation (TPE) [3]. Furthermore, the derivation details are discussed in the section A.1 of Supplementary Material (SM). The left-hand side of Eq. (3) represents the axial intensity differentiation that guides the intensity propagation along the z-axis. The right-hand side of the equation exhibits the total energy variation in a lateral intensity distribution. Such energy variations are composed of two parts. (1) The first part −∇ ⊥ · [I∇ ⊥ φ]/(kn 0 ) represents the intensity and phase variations arising from the spatial fluctuation of the refractive and attenuation index. There are two terms in these variations: I∇ 2 ⊥ φ and ∇ ⊥ I · ∇ ⊥ φ. A scalar product of the intensity distribution and the Laplacian operator acting on the phase distribution can FIG. 1. For a given wavelength λ, the initial intensity distribution of a complex plane wave I(r ⊥ , z initial ) propagates through intracellular refractive index n(r ⊥ , z), exhibiting the final intensity distribution I(r ⊥ , z f inial ).
be interpreted as the intensity variation caused by convergence or divergence of the intensity the radius of local phase curvature of which is inversely proportional to ∇ 2 ⊥ φ. The dot product of the gradient of intensity and phase distributions represents a measure of the translational effects of the phase gradient in the direction of the intensity gradient. (2) The second part −2(∂α/∂z)I relates to Beer-Lambert's law of intensity reduction in an inhomogeneous medium, exhibiting, in particular, the energy absorption and scattering of the intermediate fields traveling through an intracellular medium.
Fractal cell modeling. A wave function of intermediate fields propagating through a biological medium can encounter structures with dimensions ranging from the size of a protein macromolecule (50-100 nm) and intracellular compartments: for example, membranes (10 nm), nucleus (5-10 µm) and mitochondria (0.2-2.0 µm). Cells also contain a cytoskeleton made of filaments (7-25 nm), as well as nucleoli (0.5-1.0 µm) and DNA in the form of chromatin in nucleus. These complex structures are often modeled under the assumption that the constituents of cellular interiors are tightly filled with various discrete particles, with their surfaces pushed together to form contiguous biological components. Such cell modeling, when combined with Mie scattering theory, provides a homogeneous medium in a limited range of wavelengths and voxel size of the cell simulations [15][16][17]. Because of these constraints, the spatial fluctuations arising from the interactions between incoherent wave functions and discrete particles cannot be made to line up at the microscopic levels of intensity transport theory without breaking modeling assumptions.
Interactions of a wave function with intracellular inhomogeneities are associated with a turbulence in the refractive index distribution for which the spatial fluctuations are fractal in nature [18][19][20]. For sake of model simplicity, we consider intracellular constituents as a continuous fractal media rather than aggregation of various discrete particles. For a given fractal dimension D f , which usually represents a mathematical index for char- The WM covariance as a function of distance between two different pixel positions ρ. Each colored lines denote the four representative models of fractal media: medium I (blue), medium II (green), medium III (magenta) and medium IV (orange). (b) An example of the mean differential scattering cross-section per unit volume plotted in spherical coordinates (lc = 30 nm). An incident wave function polarized in the vertical plane can propagate from left to right along the longitudinal direction. The dimple is located at the origin.
acterizing fractal patterns or regularities by quantifying their complexity, continuous random fields of refractive index fluctuation ∆n can be modeled by the Whittle-Matérn (WM) covariance (or correlation) function, where l c are σ 2 n are correlation length and the variance of excessive refractive index, respectively. K ν (·) denotes the ν-th index of the modified Bessel function of second kind ν = (D f − 3)/2. The shape of the WM covariance function can be exhibited in a wide range of plausible fractal dimensions including: power law for D f < 3, Henyey-Greenstein for D f = 3, stretched exponential for 3 < D f < 4, Kolmogorov/von Karman for D f = 3.67, exponential for D f = 4 and Gaussian as D f → ∞. Figure  2(a) shows the WM covariance function for four different fractals models (see the SM section A.2.1 for model parameterizations).
While the attenuation function in the fractal cell model can be represented by the intensity absorption coefficient µ a and scattering coefficient µ s , most cells are transparent, weakly absorbing intermediate fields passing through their interiors. In a limiting case where µ a µ s , we consider four following assumptions: (1) The attenuation function is independent of the spatial variation of scattering coefficient. (2) No backward scattering of the wave function.
where Ω and kn 0 are solid angle and the effective wave number of intermediate fields, respectively. µ s (kn 0 ) represents the mean scattering coefficient derived from an integration of the mean differential scattering crosssection σ(k o ,k i ) over all angles. Figure 2(b) shows an example of the scattering cross-section in spherical coordinates, exhibiting, in particular, forward-directed scattering along the longitudinal direction (see the SM section A.2.2 for the σ definition and more figures for the four representative fractal media). TIE simulations via fractal cell modeling. Fundamental characteristics of the intensity transport through continuous fractal media can be extracted from numerical simulations and analyses of the TIE. Here we assume that phase functions in the TIE are given by Eq. (5), i.e., φ(r ⊥ , z) = −β(r ⊥ , z), and then implement TIE simulations over a wide range of fractal dimension and wavelength, computing the propagation and variation of lateral intensity distribution through the fractal cell models along the axial direction (see Figure 1). The finite difference approximation of the intensity propagation for the n-th image frame can be represented in the form where i and j represent the index of pixel position for a given image frame; with convergence condition for TIE simulations being δz < kn 0 δxδy. To analyze biophysical effects arising from intensity propagation via fractal cell modeling, we run the TIE simulations for the following two initial intensity conditions I(r ⊥ , z initial ): a uniform intensity distribution 1, 000 counts/pixel; and standard image in an intensity range from 500 to 1, 000 counts/pixel. Our simulation results for a thicker cell sample (∼ 10 µm) are summarized in Figure 3 (see the SM section B for more simulation results). For λ = 507 nm, the intensity propagation of the standard image through the fractal medium III (D f = 4.00) exhibits intensity attenuation and intensity dispersion in the final intensity image [see Figure 3(b)]. Figure 3(d) shows a comparison between initial and final intensity histograms, thereby dispersing and reducing the intensity distribution of the initial standard image. However, for the fractal medium I (D f = 3.25), there is significant transmissivity in the final intensity distribution [see Figure 3(a)]; with almost complete overlap between initial and final intensity histograms is shown in Figure  3(a).
Convergence and stability of the TIE simulation can be seen in the size and fractional error of the intensity dispersion which varies with axial distance. Figure 3(e) exhibits relatively fast convergence of the variation of root-mean-square (RMS) values along the axial direction (z > 2-3 µm). The TIE simulation for fractal medium I (D f = 3.25) and medium II (D f = 3.67) converges to steady state at relatively lower RMS values, roughly 50 (green line) and 11 (red line). However, at higher fractal dimensions, i.e., D f = 4.00 and 5.00, the RMS value that gives rise to convergence of intensity dispersion is about 250 (magenta line) and 500 (orange line). In addition, Figures 3(f) shows propagation of averaged fractional errors along z-axis, (z p ) = |I(z p ) − I(z p−1 )|/I(z p ), thus implying no divergence of numerical errors. Electron microscopy is capable of imaging thinner cell sample sizes less than 150 nm [21]. In the wavelength range of electron beams, we consider optical depth (or thickness) to quantify a specific level of transparency of the four fractal media. Optical depth is dimensionless factor that generally represents a measure of scattering and absorption up to a specific "depth" of intracellular media when intermediate fields (e.g., light and electron) travel through the structure inside the sample. In fractal cell modeling, optical depth can be defined as follows.
where L is physical thickness of biological samples. For L = 110 nm, Figure 4 shows optical depth over a broad range of wavelength from electron beams (blue region, λ = 1.00 pm ∼ 1.00 nm) to visible light (green region, λ = 360 ∼ 760 nm). The fractal media are optically thin (τ 1) and transparent through the visible light wavelength, weakly dispersing the light intensities passing through their interiors. Such intensity dispersion thus provides approximately-invariant transport between initial and final intensity distributions, I(z f inal ) ∼ I(z initial ) [see Figures S10(a) in the SM section B2]. However, in a wavelength range of electron microscopy imaging, there is a significant tradeoff between fractal-dimensionality and optical depth. Fractal medium I (blue), II (green) and III (magenta) are optically thick (τ 1), exhibiting, in particular, intensity attenuation as well as intensity dispersion in the final standard image [see Figure S10(d)-(f)]. Mathematically speaking, as an increase of fractal dimensions, optical depth becomes thinner and the transmissivity thereof can be increased in the fractal medium IV (orange).
Applications. A further motivation to our work is the direct measurement and reconstruction of the refractive index profiles of cellular interiors. Recent experimental studies using the TIE via bioimaging [4,22], for example, have reconstructed intracellular phase distributions without modeling the refractive index turbulence. Likewise, J. D. Rogers and his colleagues have verified the apparent similarity in power spectral densities between the WM covariance function and scanning electron microscopy images [19]. These qualitative verifications, however, were made without applying intensity transport theory and so cannot be valid as direct measurements of the refractive index profiles. For a proper comparison and measurement in bioimaging, a mathematical relation between the phase factor and the WM covariance model can be derived from the spatial-autocorrelation of Eq. (5) and then written in the form when satisfying the four assumptions declared in the cell modeling section (see the SM section C.1 for derivation details). In this formula, the WM covariance function B n (ρ) and the mean scattering coefficient µ(kn 0 ) represent the shape and offset of the phase-differential autocorrelation along axial directions. Figure S11 shows no significant shape differences over a broad wavelength range from electron beams (dashed lines) to visible lights (solid lines); a large variation of the offset is found at ρ > 10 µm due to the dependence on the effective wave number. Moreover, we use this mathematical relation to reconstruct intracellular optical properties (e.g., fractal dimension and scattering coefficient) from fluorescent cell imaging. In our demonstration, the covariance function is directly compared and fitted to the observed spatialautocorrelation curves of phase derivatives. Table S4 and Figures S14 not only represent the best estimates in fractal model parameters but also their statistical uncertainties which gives a numerical indication of the level of validity and confidence in our fitting results (see the SM section C.2 for analysis details).
Conclusion. The TIE plays a key role in bridging the resolution gap between optical and electron microscopy techniques, providing, in particular, a noninterferometric correlation between intensity and phase functions of intermediate fields (e.g., electrons and light) over a broad range of wavelengths. These microscopy techniques are capable of directly imaging the intensity variations of intermediate fields passing through the interior of cells. Variations in phase functions, however, are barely visible with bioimaging systems, and sensitive to light scattering and fluctuation through invisible (or unobservable) spatial distributions of intracellular refractive and attenuation indices. In this paper, we reformulated the TIE through fractal modeling of these two indices. Our results from these TIE simulations revealed that intensity propagations through the refractive index fluctuation can lead to nonintuitive interconnections of fractaldimensionality and intensity dispersion (or transmissivity) in the picometer to micrometer wavelength range. We also derived the spatial-correlation of phase derivatives that enables the direct measurements of the WM covariance model parameters from optical and electron microscopy imaging.
Of further significance to our work is a numerical evaluation and verification of observational invariance (or symmetry) as postulated in data science. Most subcellular observations made via optical imaging are dedicated to an empirical (or data-driven) approach, its main function being to allow biophysicists to extract the regularities and patterns of the apparent intracellular properties captured with photosensitive devices [23-27]. These regularities are often formalized into network models, mainly as a function of the observable intracellular components. However, in actual intracellular imaging, propagations and variations of light through inhomogeneous intracellular structures are dependent on the wavelength of light, as well as intracellular optical properties such as the refractive index, optical thickness, and attenuation coefficients. These optical dependencies raise questions regarding whether the subcellular model representations can be conserved or violated through observational processes. Our work sheds light on these questions from the perspective of fundamental optics, and also suggests new developments and extensions to bioimage simulation modules [5][6][7][8][9][10][11][12][13][14], including, in particular, light scattering and fluctuations caused by the nonuniform spatial distributions of intracellular optical properties. Such an implementation is of broad relevance beyond just the TIE simulations presented here, and will likely lead to more realistic simulations of biological imaging. elastic light scattering from continuous media in the A. Intensity transport via fractal cell modeling

A.1. Modification of the TIE
A plane wave function propagating through the refractive index distribution along the z-axis (see Figure 1 in the main text) obeys Helmholz equation, where ∇ 2 denotes Laplacian operator in three dimensions; A(r ⊥ , z) is a complex amplitude of the wave function. Substituting Eq. (1) in the main text into Eq. (S1), we can deduce a paraxial wave equation as follows, where ∇ ⊥ is the two-dimensional Nabla operator in the transverse direction. The amplitude propagator of the paraxial wave equation can be represented in the finite difference form of where δz denotes small changes in axial axis. The first term in the exponential factor represents the kinetic energy in amplitude propagation. The second term is coupled with an imaginary part of the wave function, representing, in particular, phase variation arising from the spatial distribution of refractive and attenuation index (δβ). The third part exhibits the amplitude variation in a real part of wave function (δα). Using Eq. (4) and (5) in the main text, Eq. (S2) can be rewritten in the form of Here we assume that the scalar complex amplitude is given by where I(r ⊥ , z) and φ(r ⊥ , z) are intensity and phase factor of the wave function, respectively. We then modify Eq. (S4) as follows.
where the β-factor vanishes through this modification.

A.2.1. The index of refraction
We construct four different fractal cell models from focused ion beam scanning electron microscopy (FIBSEM) images of Cyanidioschyzon merolae cells [5,6]. Model parameterizations are shown in Table S1. Our procedure for model construction is given as follows.
(1) 3D cell imaging: C. merolae is a unicellular haploid red alga that inhabits sulphate-rich environments such as acidic hot spring water (pH 1.5, 45 • C) [7]. The cell has a size of about 2-5 µm with relatively simple cellular architecture, containing only 4775 protein-coding genes in its nucleus and the minimum number of membraneenclosed organelles (i.e., a single nucleus, a single mitochondrion and a single chloroplast), lacking a vacuole and cell wall. Furthermore, cellular and organelle divisions can be highly synchronized through cell cycle progression. Therefore, the C. merolae cells is an excellent model system in which to study cellular and organelle division processes as well as structural biology.
FIBSEM imaging techniques involve combining serial gallium-ion beam etching of a resin block-face of cell samples with electron scanning of the exposed block surface. Such combinatorial operation allows accurate visualization of the 3D intracellular structure of C. merolae at the whole cell level. Parameters for the FIBSEM cell imaging are configured as follows: gallium-ion beam current, 1.9 pA; electron wavelength, λ = 4.00 pm; defocus distance, δz = 11.7 nm; image pixel length, δx = δy = 5.6 nm; total number of images, 796 frames; and each image size, 1679 × 1402 pixels. Figure S1(a) shows an example snapshot of the intensity distribution for FIBSEM images of a C. merolae cell block.
(2) In the FIBSEM images, multi-Otsu thresholding method [8] can be applied to split between intracellular region and non-cellular background region. The cell segmentation results for the 100-th image frame are shown in Figures S1(b) and (c).
(3) In binary images, we use a spatial statistical tool such as GSTools [9] to generate spatial random fields of intracellular refractive index, following, in particular, the WM covariance function B n (ρ) [i.e., Eq. (6) in the main text].  In fractal cell modeling, an incident monochromatic wave function polarized in the directionê i propagates along the longitudinal directionk i . The wave function can be scattered through continuous fractal media and then observed in the far field region at the directionk o [see Figure S2]. The mean differential scattering crosssection per unit volume is defined as the averaged power per unit solid angle at directionk o scattered by unit volume of the fractal medium and then can be written in the form of where sin 2 χ = 1 − sin 2 θ cos 2 φ and k s = 2kn 0 sin(θ/2). The power spectral density of refractive index fluctuation Φ n (k s ) can be derived from Fourier transformation of the WM covariance function B n (ρ). Figures S3 show the examples of the mean differential cross-section plotted in a spherical coordinate. In Figure S3(a), the dimple is located right on the middle of an ellipsoidal-like shape of the scattering cross-section, representing, in particular, uniform scattering in all direction. Figures S3(b) and (c), however, exhibit forward scattering along the longitudinal direction (or z-axis). The dimple cannot be clearly seen in Figure S3(c) but definitely located at the origin. Furthermore, the mean differential cross-section for the four representative models of fractal media (see Table  S1 for model parameterization) are compared in the right column of Figures S4.  spatial distributions of refractive index and its histograms for the 100-th image frame, respectively. In the right column figures, the mean differential scattering cross-section per unit volume are plotted in a spherical coordinate and then compared between the four different fractal media: medium I (blue), medium II (green), medium III (magenta) and medium IV (orange). The dimple of the cross section cannot be clearly seen in these column figures but definitely located at the origin. Root-mean-square (RMS) and its relative values of intensity distributions between the p-th and q-th index of image frames are given by where N is the total number of image pixels.     C. Application C.1. Derivation of the autocorrelation function in differential-phase factors Using Eq. (5), (6) and (7) in the main text, the Eq. (10) is derived as follows.
FIG. S11. The spatial-autocorrelation of phase derivatives normalized by (kn0) 2 is shown as a function relative pixel distance ρ. Solid and dashed colored lines represent wavelength ranges for visible light (λ = 507 nm) and electron beams (λ = 4.00 pm), respectively. Each colors denote the four different fractal models.

C.2. Direct measurement of intracellular fractal profiles
The Eq. (10) in the main text is applicable to optical as well as electron microscopy imaging. A direct measurement of the WM covariance model parameters in fluorescent cell imaging is demonstrated as follows.
(1) Fluorescent cell imaging: Fluorescence microscopy system has been configured to visualize intracellular distribution of α-tubulin and histone H2B in HeLa cells. The HeLa cells stably expressing histone H2B-GFP were a kind gift from Dr. Hiroshi Kimura (Tokyo Institute of Technology, Japan) [10]. The HeLa cells were fixed with 4% paraformaldehyde in 250 mM Hepes-NaOH (pH 7.2) for 10 min at room temperature (RT) (approximately 26 • C). After permeabilization and blocking by blocking buffer (10% Blocking One (Nacalai Tesque, Kyoto, Japan; Cat No. 03953-95) and 0. (2) Phase retrieval based on original TIE [11,12]: We use the geometric multigrid (GMG) method to reconstruct the phase distribution from fluorescent cell images [13, 15? , 16]. The lateral phase distribution of the p-th axial index φ(r ⊥ , z p ) can be reconstructed from three different images along the axial direction [see Figure S13(a)]. In numerical analysis, the GMG method aims to solve original TIE by using a hierarchy of discretization as well as approximating derivatives with finite differences. This hierarchical discretization enables acceleration of the convergence of a basic iterative relaxation, efficiently reducing errors in high spatial-frequency components by global correction of the fine grid solution approximation in a step-by-step manner. The original equation can be rewritten in the approximate form of finite differences, following, in particular, the assumption that intensity differences between the two images I(z p+5 ) − I(z p−5 ) is approximately linear with an equal axial distance 5δz (1.00 µm). The finite difference approximation of the phase distribution for a specific axial position of the p-th index is given by FIG. S12. Example snapshots of fluorescent cell images for the 5-th, 10-th and 15-th index. where i and j represent the index of lateral intensity and phase distributions along the x-and y-axes; δx and δy are lateral resolutions; δz is axial distance.
The entire procedure for the GMG algorithm [see Figure S13(b)] can be summarized in the following five steps: presmoothing, reducing errors in high frequency components by iterative solvers such as Jacobi, Gauss-Seidel and conjugate gradient methods; residual estimations, computing residual errors after the iterative solvers; restriction, downsampling the residual errors to a coarser grid; prolongation, interpolating a correction estimated on a coarser grid into a finer grid; and, correction, modifying the interpolated coarser grid solution onto the finer grid. These five steps are implemented into the simplest multigrid process, termed the V-cycle, reconstructing the phase distributions from the 31 fluorescent cell images. Figures S13(c)-(e) show an example snapshot of the reconstructed phase distribution of the 10-th axial index, as well as error convergence of intensity distributions in the V-cycle process.
(3) Parameter fitting: For the sake of analysis simplicity, we first omit the offset of Eq. (S12) and then rewrite the spatialautocorrelation equation in the form of where phase derivative ∂φ/∂z is given by taking a linear difference of the two phase distributions in axial defocus distance δz, ∂φ(r ⊥ , z) ∂z ≈ φ(r ⊥ , z p+1 ) − φ(r ⊥ , z p ) δz (S15) Figure S14(a) shows an example output of the phase derivative in fluorescent cell imaging.
Fractal properties of cellular interiors can be found in the shape and offset of the spatial-autocorrelation curve of phase derivatives. In our analysis, Eq. (S12) can be fitted to the observed spatial-autocorrelation curves of the derivatives obtained from the 5-th to 25-th index of phase distributions in fluorescent cell imaging. The best fit parameters are these which minimize χ 2 = −2 ln L where L is the likelihood function, where E i and O i are the expected and the observed spatial-autocorrelation at the i-th input. σ i is statistical error in O i . Furthermore, no penalty term for nuisance parameters is introduced in this minimization function. Table  S4 and Figure S14 show the results from model parameter fitting without analyzing systematic uncertainties (or errors) arising from various biophysical sources in fluorescent cell imaging. In Figure S14(b), the best fit (red solid line) of spatial-autocorrelation function is directly compared with the observed spatial-autocorrelation curve (blue bars). We also estimate the statistical uncertainties in the model parameters to indicate numerically the validity and confidence of our fitting results [see Figures S14(c)-(f)].

Model parameters
Fitting results  FIG. S14. Fitting results. (a) An example output of differential-phases along axial direction. (b) The spatial-autocorrelation curves of differential-phases in axial direction: blue bars and red line represent the observed spatial-autocorrelation and the best fit curves. (c-e) ∆χ 2 -contour plots for fractal dimension D f , fluctuation (n0σn) 2 and the correlation length lc. Blue points indicate the best estimates. Black solid lines around the best fit values are allows at the 1σ (68%), 2σ (95%) and 3σ (99%) confidence levels of the fitting results. (f) The hyperbolic curve of ∆χ 2 is represented as a function the mean scattering coefficients µs converted from the fitting results (assuming n0 = 1.38).