Fluorescence microscopy point spread function model accounting for aberrations due to refractive index variability within a specimen

Abstract. A three-dimensional (3-D) point spread function (PSF) model for wide-field fluorescence microscopy, suitable for imaging samples with variable refractive index (RI) in multilayered media, is presented. This PSF model is a key component for accurate 3-D image restoration of thick biological samples, such as lung tissue. Microscope- and specimen-derived parameters are combined with a rigorous vectorial formulation to obtain a new PSF model that accounts for additional aberrations due to specimen RI variability. Experimental evaluation and verification of the PSF model was accomplished using images from 175-nm fluorescent beads in a controlled test sample. Fundamental experimental validation of the advantage of using improved PSFs in depth-variant restoration was accomplished by restoring experimental data from beads (6  μm in diameter) mounted in a sample with RI variation. In the investigated study, improvement in restoration accuracy in the range of 18 to 35% was observed when PSFs from the proposed model were used over restoration using PSFs from an existing model. The new PSF model was further validated by showing that its prediction compares to an experimental PSF (determined from 175-nm beads located below a thick rat lung slice) with a 42% improved accuracy over the current PSF model prediction.


Introduction
The three-dimensional (3-D) image of a point is affected by the limiting aperture of the wide-field microscopy system, which prevents the formation of a perfect image. The aperture fails to collect all the light emitted by the point source and also introduces a diffraction pattern. Imaging through stratified media adds an additional challenge due to refraction at media interfaces. 1 The 3-D image formation process can be modeled as an integral operation of the true fluorescence distribution with the point spread functions (PSFs) of the system, which is the impulse response of the system. 2 Computational optical sectioning microscopy (COSM) 3 algorithms restore the true object intensities from the acquired image using prior knowledge about the imaging conditions, including the PSF of the system. 4 The accuracy of the restoration is contingent on the accuracy of the PSF model, which this study aims to improve by incorporating the impact of a specimen on the overall system response.
Biological specimens are thicker (>5 μm) than the depth of field of a high numerical aperture (NA) lens (<0.2 μm) and have a variable refractive index (RI), which is different from the RI of the immersion medium and of the coverslip. Spatial variability in the RI of samples is attributed to the variability of their cellular components. 5 Lung tissue, for example, demonstrates RI variability because there is air (RI ¼ 1.00) in the alveoli, which are surrounded by cellular material with RI between 1.30 and 1.41. 6,7 3-D imaging of such specimen with wide-field microscopy suffers from both defocus and spherical aberration (SA). 3-D imaging with confocal or structured illumination microscopy, although not affected by defocus, is still prone to SA. 8,9 Light emitted by points located deep in the specimen layer must pass through layers whose thicknesses and RI can vary significantly from the layers in the ideally designed system. Refraction between the nondesign layers affects the spherical wave-front of the propagated light resulting in a larger spread of the focus, thus reducing both resolution and peak intensity. 10 This phenomenon poses a challenge to restoration algorithms in COSM, creating the need for a new paradigm.
The variability of the PSF due to change in SA throughout the specimen introduces the challenge of space-variant (SV) imaging. Depth-variant (DV) imaging, i.e., variance along Z has been investigated by many researchers. [11][12][13] Our work has been focused on developing solutions for the SV imaging problem based on integrating four on-going lines of investigation: (1) development of a methodology to derive the sample RI; 14 (2) development of an accurate but practical PSF model (reported in this paper); (3) a block-based imaging model 15,16 that integrates multiple PSFs (computed with the model proposed here) to approximate the SV forward image; and (4) a restoration algorithm based on the block-based SV imaging model. The new PSF model presented here extends the mathematical formulation in existing models 17,18 and is a key component within our broader goal that aims to restore images of thick biological samples using practical model-based computational algorithms. Experimental PSFs, obtained by imaging small point-like sources, such as subresolved fluorescent spheres within a specimen, can provide an accurate representation of the system's response for specific imaging conditions. 19 However, experimental PSFs are not easy to measure, 20 and they tend to have reduced signal-to-noise ratio and poor sampling due to the efficiency of the detector. 21 Analytical PSFs, on the other hand, are limited by model assumptions based on our knowledge of the system and specimen. Improved theoretical PSF models contribute to the development of COSM by enabling computation of more accurate imaging models, and by avoiding the acquisition of multiple experimental PSFs for each restoration.
Initial models 2,22 developed for 3-D PSFs assumed that the specimen is thin with a uniform RI and that it is placed directly under the coverslip; none of these assumptions hold for most biological samples. Prevalent models 17,18 account for SA due to thickness, but the RI of the sample is assumed to be uniform throughout. In the study presented in this paper, a specimen is modeled to have multiple layers of variant RI, which is accounted within the mathematical formulation of the extended PSF model.
The widely used Gibson and Lanni (G&L) PSF model 17 is optimized for a three-interface stratification assumed in microscopy systems, comprising the oil-immersion lens (RI ¼ 1.512), the coverslip (RI ¼ 1.522), and the specimen with an average RI close to water (RI ¼ 1.33). The use of design and actual imaging conditions and efficient calculations render the G&L model prevalent, but it is challenged for high-NA lenses because the extremal rays impinge with large angles of incidence at interfaces. This challenge can be overcome by a vectorial PSF model, such as the one proposed by Török and Varga. 23 The Török and Varga (T&V) model accounts for RI variance within the sample by extending the Born and Wolf model 1 of diffraction for a stratified medium by application of the Fresnel refraction law. The model is rigorous but uses parameters such as the initial aberration function, which cannot be directly determined from the microscope imaging conditions.
Haeberlé 18 combined the optical path difference (OPD) described in the G&L model with the vectorial representation of the PSF provided by the T&V model, to give a more accurate PSF using the design and actual acquisition parameters. The Haeberlé representation of the PSF model has been proven effective in restoring microscope acquired images. 24,25 The PSF model discussed in this paper was developed using an approach similar to the one presented by Haeberlé, and it represents an extension of his model, which assumes that the sample has a uniform RI. The new PSF model presented here accounts for a spatially variant sample RI by assuming that different layers within the sample can have a different RI, thereby approximating the imaging conditions more accurately than previous models while keeping the computations tractable.
In our overall approach, the RI variability within the specimen is represented by a finite arrangement of nonoverlapping blocks [ Fig. 1(a)], where the RI is assumed to be a constant within each block. 15 The block size is determined by the smallest volume within the sample, where the RI can be assumed to be constant, since very small variations in RI do not appreciably affect the PSF. This approximation models specimen RI more Fig. 1 An N-interface point spread function (PSF) model: (a) specimen with space-variant refractive index (RI) divided into blocks with each block having a uniform RI. Light rays from a unique point source (indicated by the dot at the bottom of the axially stacked blocks) travel through blocks, with different RIs, directly above it. The axially stacked blocks (which look like a column) are isolated by the model to calculate the PSF. (b) Schematic diagram of the light path from a point source on the optical axis to the objective lens within a unique column. The light travels through the specimen with multiple layers of varying RIs in addition to the layer of the immersion medium and the coverglass. n i , n g , and n sk are the refractive indices of the immersion medium, coverglass, and sample, respectively, while, t i , t g , and t sk are their corresponding thicknesses. h 1 , h 2 , and h N are the locations of the interfaces. In the design system, the point source is located just under the coverslip at P, but for thicker samples, its location is at A, which is deeper within the sample. Variables marked with an asterisk "*" represent design parameters as opposed to the nondesign (actual) parameters. The design optical path is PQRS, while ABCD is the nondesign (actual) optical path through the stratified medium; their difference is the optical path difference. accurately by allowing it to vary not only axially but also laterally. Light from a point source travels through all the blocks directly above it, and each interface between two blocks introduces a change in the optical path, influencing the formation of its PSF. The PSF model presented here is implemented on a stratified column formed by isolating a series of axially stacked blocks directly above the point source [ Fig. 1(a)]. Each point source is associated with its own stratified column, which allows PSF computation using the new model. This arrangement is further elucidated in Fig. 1. In the rest of the paper, we refer to our proposed model as the N-interface PSF model.
The computation of the OPD in the imaging system using only three layers (i.e., the lens' immersion medium, the coverslip, and the specimen) has been understood to be an approximation, which is not easy to relax because variability within the specimen is hard to model uniquely. Advances in understanding the properties of samples are currently being studied using quantitative phase microscopy approaches. Specimens have the property of changing the phase of light that passes through them due to spatial variations in their RI and thickness; this change can be exploited by specialized microscope setups to generate image intensity variations (contrast). 14 3-D RI maps of objects have been investigated using phase contrast microscopy 26 and tomographic phase microscopy. 27 The RI of lung cells has been studied using a combination of confocal and quantitative phase amplitude microscopy. 5 Kam et al. 28 modeled the change in RI using specimen phase information and a ray tracing method to obtain the aberrant PSF. However, ray tracing based on ray optics is not computationally practical in iterative restoration algorithms because ray tracing attempts to restore one pixel at a time. Hiware et al. 29 extended the OPD of the G&L model to include the variation of turbulent but small RI variances (0.001 to 0.005) within the sample as suggested by Schmitt. 26 The Schmitt sample approximation is not suitable for specimens with large RI variances (1.00 to 1.49), such as those encountered in ex vivo lung tissue 30 (discussed previously). Our investigation with a 63x/1.4 NA oilimmersion lens has shown that for a watery sample (a common assumption for many biological samples), a 0.5% change in the RI (i.e., a change from 1.333 to 1.338) results in >1% change in the PSF (quantified by a mean square error) at a 10 to 50 μm depth range within the sample, which is negligible with respect to imaging and restoration. The proposed PSF model is not limited to such small RI variations.
The paper is organized as follows. Section 2 discusses the mathematical formulation of the N-interface PSF model. The validity of the new PSF model for specimen with RI variance is investigated using experimental results and it is discussed in Sec. 3. The effect of RI variance within the specimen on the response of a typical system is studied using simulated N-interface PSFs in Sec. 4. The use and advantage of N-interface PSFs in the restoration of 3-D images of 6-μm-diameter beads is presented in Sec. 5. The model is further evaluated using images of point sources placed below rat lung slices in Sec. 6. A summary of the findings from this study and further developments are discussed in Sec. 7.

Mathematical Models
This section describes the mathematical formulation developed for the N-interface PSF model. The new PSF model introduced in this paper accounts for RI variation within the specimen modeled over N (≥3) interfaces between layers, such as the lens' immersion medium, the coverslip, and one or more interfaces within the specimen.
The PSF model presented here requires a 3-D sample RI map, nðx; y; zÞ, for the determination of the optical path length (OPL). The OPL is the product of RI and the physical distance that light travels when it propagates through the sample. The determination of the sample RI map is a challenging problem that has been tackled by several groups 5,7,27 including ours. 14 Our overarching method described in Sec. 1 utilizes an automated multimodal microscope capable of 3-D wide-field fluorescence microscopy and differential interference contrast (DIC) microscopy, which facilitates the determination of sample RI and thickness from quantitative phase derived from DIC images. 14 In our formulation, the N-interface PSF model assumes that a specimen 3-D RI map is available.
Our model assumes that the imaging system is comprised of N − 1 stratified layers with distinct RIs that would not cause total internal reflection or generation of evanescent waves, as presented in the T&V model. 23 OPD is the change in the OPL between corresponding rays from the design and nondesign system that pass through the same point in the back focal plane as defined in Ref. 17. The OPD used in the G&L model 17 is modified here to account for a stratified specimen with variable RI, to consequently model the influence of the sample on the amount of SA more accurately.
The OPD (Fig. 1) is based on Snell's law of refraction and the Abbe sine condition 31 as where OPD i is the OPD introduced by the immersion medium, OPD g is the OPD introduced by the coverglass, and P K k¼1 OPD s k is the sum of all the OPDs introduced by the K layers within the specimen. Equation (1) based on derivations elaborated in Ref. 17 can be rewritten as where n i , n g , and n s k are the RI of the immersion medium, coverglass, and sample, accordingly, while t i , t g , and t s k are their corresponding thicknesses, and K is the number of specimen layers. The overall system has N (N ¼ K þ 2) layers with the additional layers belonging to the immersion medium and coverglass. Terms with an asterisk "*" indicate objective lens design conditions. θ is the angle both rays (ABCD and PQRS in Fig. 1) entering the frontal lens of the objective make with the horizontal axis. As assumed in the G&L model, the normalized radius in the back focal plane (where ρ ¼ 0 at the optical axis and ρ ¼ 1 at the edge of the projection of the microscope's limiting aperture onto the back focal plane) can be expressed as ρ ≈ n a sin θ a ∕NA, where NA is the numerical aperture of the objective lens, while the origin of the (x; y; z) coordinate system is at the Gaussian focus point, and thus, Δz is the defocus of the system. θ a is the angle the ray makes with the horizontal axis when entering in the ath layer (i.e., the angle of incidence), and n a is the corresponding RI, where (a ¼ 1; : : : ; N).
In the following equations, we show that the OPD in Eq. (1) and the aberration function (ψ) in Ref. 23 are approximately equivalent, justifying the use of the OPD in the vectorial model proposed by Török and Varga. 23 The entire system is assumed to be composed of individual interfaces of the N stratified layers, which are located axially at z ¼ −h 1 ; −h 2 ; : : : ; −h N−1 as defined in Ref. 23. The layers are assumed to be separate regions of homogenous and isotropic material, with RI of n 1 ; n 2 ; : : : ; n N , correspondingly. The angle of incidence at the first interface is denoted by θ 1 and the angle in the Nth medium by θ N . The aberration function [Eq. (28) from Ref. 23] can be written as In Eq. (4), we map the notation of the modified G&L OPD parameters [Eq. (2)] to the notation of the terms in Eq. (3). The specimen is comprised of K stratified layers as indicated above, with each layer having RI n s k ¼ n l , where (k ¼ 1;2; : : : ; K and l ¼ 3; : : : ; N). The RIs of the immersion medium and coverglass are defined as n i ¼ n 1 and n g ¼ n 2 , respectively. The thickness of the immersion medium and the coverglass are defined as t i and t g ¼ h 3 − h 2 , respectively. The thickness of the Kth (where the point source is located) layer is defined as t s K ¼ h N−1 , while the thickness of all the other specimen layers are defined as t s k ¼ h l−1 − h l , where (k ¼ 1;2; : : : ; K − 1; l ¼ 3;4; : : : ; N − 1). The normalized radius (ρ) in the back focal plane can be used to deter- Modern objective lenses are compensated for design parameters like the RI and thickness of the coverslip and the immersion medium, which are not present in the T&V PSF model. 23 Like Haeberlé, 18 we borrow the compensation terms t Ã i and t Ã g for the design thickness of the immersion medium and coverglass, respectively, and n Ã i and n Ã g for their corresponding RIs, from the G&L model. These compensation factors are subtracted from the aberration function (ψ) defined in Eq. (3) and are the fourth and fifth terms in Eq. (4). The second and last terms in Eq. (3) are combined to form the third term in Eq. (4), which together describe the properties of the specimen. Replacing the variables in the aberration function (ψ) with notations from the G&L model and adding the compensation terms, we can update Eq. (3) to generate Eq. (4).
In Eq. (4), all the parameters can be determined from the microscope apart from h 1 , which varies as the position of the lens varies when it is focused deeper within the sample. We use Fig. 1 to rewrite h 1 in terms of the other described variables such that the mathematical formulation is independent of it. It can be written as Hence, the first term in Eq. (4) can be written as Replacing h 1 [Eq. (5)] in Eq. (4), we observe that all the terms in Eq. (2) are present in Eq. (4). This substitution allows us to show that OPD ≈ ψ for a system with RI variability within the specimen.
For completion, we rewrite the PSF using the OPD in Eq. (2). The PSF at any arbitrary point ðx; y; zÞ is calculated from the electrical field E as in Ref. 23.
where E Nx , E Ny , and E Nz are the directional electric fields (with ϕ in spherical polar coordinates) of a typical ray in the Nth medium. The illumination integrals I of the electrical fields are where J n is the Bessel function of the first kind and of order n, and k N is the wavenumber in the Nth medium. T ðN−1Þ m is the transmission coefficient of the stratified layers, where m ¼ s, p indicates the directionality of the polarized light traversing through the (N − 1)th medium. The derivation of Eq. (7) and further explanation of all the parameters are described in detail in Ref. 23.
This N-interface PSF model is a reformulation of the T&V model 23 using easily accessible parameters as in the G&L model 17 for specimens with nonuniform RI. The new N-interface PSF model was implemented in MATLAB by extending the code originally developed for the Haeberlé PSF model. 32 The experimental evaluation of the N-interface model and its use in restoration is demonstrated in the next sections.

Evaluation of the N-Interface PSF Model
To establish the accuracy of the theoretical model developed in Sec. 2, PSFs experimentally acquired from spatially variant samples are compared to simulated N-interface PSFs. This section discusses the details for acquiring experimental PSFs (Sec. 3.1) and simulating corresponding theoretical PSFs (Sec. 3.2), and provides results from a comparison between theoretical and experimental PSFs (Sec. 3.3).

Experimental Acquisition of PSFs from Test Samples with Depth-Variant Refractive Index
Experimental PSFs were imaged from samples with subresolved beads using a Zeiss Axio Imager (Carl Zeiss Microscopy GmbH, Germany) with optical sectioning capability and an AxioCam MRm camera. Molecular Probes FluoSpheres, 0.175 μm in diameter stained with Alexa Fluor, were dried on a microscope slide and sealed with UV-cured optical cement with RI ¼ 1.56 [Norland Optical Adhesive 63 (NOA63), Norland Products, New Jersey] . Another layer of FluoSpheres, stained with 4', 6-diamidino-2phenylindole (DAPI), was dried on top the cured optical cement. This preparation was submerged in ProLong Gold antifade reagent (RI ¼ 1.46) (a mountant used extensively in fluorescence microscopy to prevent photo-bleaching and preserve the fluorescence labeling for long-term storage and analysis). A layer of Rhodamine dyed FluoSpheres was dried on the coverslip, which was carefully placed onto the sample and sealed.
The distance between the best focus of the different color beads present at each layer of the test sample was used to obtain the thickness of each layer. This method has been proven effective to obtain the approximate depth of the sample 24 and has also been verified using guides etched on the coverslip. 33 A schematic diagram of this system is shown in Fig. 2(a). The distance between the coverslip-mountant interface and the mountant-optical cement interface was determined to be ∼45 μm, while the distance from the mountant-optical cement interface to the bead was determined to be ∼43 μm [ Fig. 2(a)]. A 63×∕1.4 NA oil lens was used for imaging. The captured 3-D images have a cubic voxel with size of 0.10 μm in X, Y, and Z directions.

Simulation of N-Interface PSFs for Comparison with Experimental Data
To establish the validity of the new N-interface PSF model introduced in Sec. 2, theoretical PSFs were computed with parameters from the experimental setup described in Sec. 3.1. PSFs from the existing Haeberlé model 18 were calculated for comparison assuming that the specimen has a uniform RI, because the model cannot account for RI variance within the sample. A common practice is to compute an average RI from the different sample constituents. A weighted average RI for a stratified sample can be determined by where t tot is the total thickness of the specimen, while t s k and n s k are the thickness and RI of the kth sample layer, respectively. For this study, four PSFs were generated as follows: (1) an Ninterface PSF using the exact sample parameters of the experimental system, i.e., a specimen with two different layers of media [RI ¼ 1.46 with thickness of 45 μm and RI ¼ 1.56 with thickness 43 μm as in Fig. 2(a)]; and (2) three PSFs using the existing Haeberlé model assuming specimens with total thickness of 88ð¼ 45 þ 43Þ μm and a uniform RI equal to: RI ¼ 1.56 (for optical cement NOA63), RI ¼ 1.46 (for mountant ProLong ® Gold antifade reagent) and RI avg ¼ 1.51 [calculated using Eq. (9) and the experimental information].
The numerical PSFs were convolved with a simulated 175nm-diameter sphere to accurately represent the experimentally determined PSF. The values of the simulated PSFs were scaled by 10 4 AU to match the intensity range of the experimental data for faithful normalization. The intensity values of the PSFs were normalized between 0 and 1, and the peaks were aligned for comparison with the experimental data. Simulated and experimental PSFs were compared to evaluate the accuracy of the new PSF model and its accuracy over the existing model. The results are presented in Sec. 3.3.

Evaluation of the Theoretical N-Interface PSF Model Using Experimental Data
In this section, simulated PSFs (discussed in Sec. 3.2) generated from the N-interface PSF model are compared to experimental PSFs discussed in Sec. 3.1 for model evaluation and also to simulated PSFs generated using the existing Haeberlé PSF model to highlight differences between the new and existing PSF model. The experimental PSF image acquired using the conditions described in Sec. 3.1 is shown in Fig. 2(b). Figure 2(c) shows the XZ cross-section of the simulated PSF computed using the N-interface model for a four-layer system with parameters that match the experimental system. Qualitative comparison shows that the PSFs in Figs. 2(b) and 2(c) are similar, especially with respect to the direction of the PSF side lobes. PSFs computed using the Haeberlé PSF model for specimen RI equal to that of ProLong Gold [ Fig. 2(d)], optical cement [ Fig. 2(e)], and averaged RI [ Fig. 2(f)] deviate from the experimental PSF [ Fig. 2(b)]. In Fig. 2(c), light experiences negative SA (Ref. 34) because it travels from a lower to a higher RI medium (1.46 < 1.512), while in Fig. 2(e), light experiences positive SA because it travels from a higher RI (1.56) to a lower one.
Quantitative comparison using axial profiles through the PSFs [ Fig. 2(g)] shows that the N-interface-system PSF is in better agreement with the experimental data than the PSFs from the three-interface systems. The central lobes of the experimental and simulated PSF from the N-interface system overlap. The location of the secondary maximum and minimum also colocalize indicating agreement between the experiment and simulation. The axial profile through Fig. 2(f) has the same direction of aberration, but the secondary lobes do not localize with the four-interface PSF [ Fig. 2(c)] and are also higher in intensity, indicating that it experiences less SA. To quantify the deviation of the different theoretical PSFs from the experimental PSF, a normalized mean square error (NMSE) metric 35 was computed using the experimental PSF as a reference. The calculated NMSE values are 0.13 for the N-interface PSF and 0.41, 0.57, and 0.23 for the PSFs predicted by the existing model assuming specimen RI ¼ 1.46, RI ¼ 1.56, and RI ¼ 1.51, respectively. The NMSE reinforces the fact that the PSF generated using the N-interface PSF is closest to the experimental PSF.
The accuracy of the N-interface model PSFs for a specimen with variable RI was shown in this section. The mismatch in PSFs caused by assuming the specimen has uniform RI is also highlighted here.

N-Interface PSFs Computed for Different Sample Conditions
For comparison, two sets of theoretically predicted PSFs were computed and investigated assuming that the point light source is in (1) a specimen with a uniform RI and (2) a specimen that has two or three different layered RIs. The RIs used in the simulation were chosen to be similar to RIs of materials commonly found in microscopy systems, such as air (RI ¼ 1.00), water (RI ¼ 1.33), and glycerol (RI ¼ 1.47).
Simulations were carried out on a 256 × 256 × 256 grid with a cubic voxel of size 0.1 × 0.1 × 0.1 μm. The emission wavelength was set to 488 nm and a 63×∕1.4 NA oil-immersion lens (RI ¼ 1.512) was used. Five PSFs were generated assuming that the point source of light is at a depth of 10 μm below the coverslip within different specimens (Fig. 3). In the five cases presented here, the specimen is comprised of 1. a single medium with RI ¼ 1.47 with thickness 10 μm [ Fig. 3(b) Fig. 3(c)], To quantify the deviation of the different theoretical PSFs from a reference PSF, an NMSE metric 35 was computed. PSF comparison and discussion of the variability of the PSF with the RI change within the specimen is summarized in the next section.

Quantification of PSF Change with Depth-Variant RI Within the Specimen
The N-interface PSF model was used to generate PSFs for point sources of light located under layers of different RI within different numerical samples, while keeping other imaging conditions constant (including the overall sample thickness) as described in Sec. 4.1.  Fig. 3(b)] indicated by the highest peak intensity and lowest secondary lobe intensities. The highest amount of SA is experienced by the PSF for specimen case 3 [ Fig. 3(d)] due to the highest RI mismatch with the immersion medium. The PSF from specimen case 4 [ Fig. 3(e)] experiences different amount of SA than the PSF computed using an averaged RI [ Fig. 3(f)] of the layers of specimen case 4. This demonstrates that using a single RI based on an average of specimen RI can lead to PSF errors (consistent with the result in Fig. 2). The change in the PSFs is also quantified using the NMSE calculated with respect to Fig. 3(b), which was chosen as the reference because it experiences the lowest SA. The NMSE is equal to 0.43 for the PSF from specimen case 2 [ Fig. 3(c)], 1.02 for specimen case 3 [ Fig. 3(d)], 0.89 for specimen case 4 [ Fig. 3(e)], and 0.92 for specimen case 5 [ Fig. 3(f)]. The NMSE values corroborate with the inference from the intensity plots in Fig. 3(g). These results highlight the PSF variation due to SA at a fixed depth of 10 μm below the coverglass. The SA in this study is only due to the depth variability of the RI within the specimen. Thus, the study demonstrates the impact of RI variability on the PSF and that averaging the RI is not an acceptable solution to address specimen variability, underlining the need for appropriate PSF modeling. This section emphasizes the need to account for specimen properties in order to obtain a more accurate PSF of the imaging conditions.

Restoration of Experimental and Simulated
Images of Samples with Nonuniform RI Using N-Interface PSFs An important role of the N-interface PSF model is its use in restoring images of samples with nonuniform RI. This section provides results from a restoration study using the proposed PSF. Data acquisition from a controlled sample with RI variance is first described in Sec. 5.1. Simulation of the experimental system to serve as a control set is discussed next in Sec. 5.2. The restoration results obtained from the experimental and simulated data are discussed in detail in Sec. 5.3 to highlight the advantage of using N-interface PSFs for spatially variant samples.

Experimental Data Acquisition from a Sample with Variable RI
FocalCheck microspheres, 6 μm in diameter labeled with DAPI throughout and a 1 μm fluorescent outer shell labeled with Alexa Fluor, were imaged with a 63×∕1.4 NA oil-immersion lens (RI ¼ 1.515) using a Zeiss Axio Imager. The microspheres were air-dried on the glass slide, and then they were sealed with UV-cured optical cement (NOA63, RI ¼ 1.56). A drop of mountant [ProLong Gold antifade reagent cured for 24 h (RI ¼ 1.42)] was placed on top of the optical cement layer to introduce RI variability within the specimen. A layer of FluoSpheres (0.175 μm in diameter), dyed with Rhodamine, was air-dried on top of the cured layer to demarcate the interface between the optical cement and ProLong Gold. A layer of FluoSpheres dyed with DAPI was air-dried on the coverslip to locate the exact location of the coverslip. A 256 × 256 × 300 size image with a cubic voxel size equal to 0.10 μm in the X, Y, and Z directions was acquired. The location of the microscope stage recorded at the best focus of the various beads in the sample was used to obtain the thickness of each layer. In future studies, the accurate thickness and RI of the specimen will be derived using phase-microscopybased methods, described in Sec. 1. The distance from the coverslip-mountant interface to the glycerol-mountant interface was determined to be ∼44.5 μm using the DAPI-and Rhodaminestained marker beads. The distance from the Rhodamine-stained marker beads to the 6-μm microspheres was found to be equal to ∼73 μm [ Fig. 4(a)].
The FocalCheck microsphere test sample was also imaged using the Zeiss ApoTome structured illumination attachment to assess more information about the underlying sample needed to generate a numerical object that resembles the sample accurately. The generation of simulated images comparable to the experimental image is discussed in the next section.

Simulation of Images from a Sample with Variable RI
Simulated images were computed based on a strata approximation model that accounts for SA introduced due to changes in imaging depth within the sample. 11,36 PSFs were computed using the N-interface model (N ¼ 4) for a point light source in a specimen with two different RIs equal to 1.42 and 1.56. The thickness of the mountant layer (RI ¼ 1.42) was set to 44.5 μm and the thickness of the optical cement (RI ¼ 1.56) was assumed to be between depths of 70 and 76 μm, respectively, covering the diameter of the 6-μm bead. This simulated system was completely matched to the experimental system. The numerical object used in the simulations was generated on a 256 × 256 × 300 grid and a high intensity (10 4 AU) shell with an outer diameter of 6 μm and inner diameter of 4 μm. The same simulated image was restored using four different sets of PSFs (A to D, described in Sec. 5.3) to evaluate the effect of PSF model mismatch on the restoration result. Noise was not introduced in the simulated data to isolate the effect of the PSF variability on the restoration. The restored images computed from the experimental and simulated data were studied and compared to demonstrate the effect of the N-interface PSF on the restoration.

Comparison of Restoration Results from Experimental and Simulated Images
The use of N-interface PSFs for restoration of experimental and simulated data (presented in Secs. 5.1 and 5.2) is discussed in this section. The restoration results were obtained after 1000 iterations (a number determined experimentally to yield acceptable results in practice) of the DV expectation maximization algorithm 11 implemented in the COSMOS software. 37 For the restoration, four different sets of DV PSFs were used. Each set consisted of seven DV PSFs each computed at different depths (1 μm apart) covering a range of 6 μm within the bead, using the following PSF models and assumptions: A. The N-interface model (N ¼ 4) for a point light source in a specimen with two different RIs equal to 1.42 and 1.56. The thickness of the mountant layer (RI ¼ 1.42) was set to 44.5 μm, while the thickness of the optical cement (RI ¼ 1.56) was set between depths 70 and 76 μm, covering the diameter of the 6-μm bead. The total thickness of the specimen, composed of two distinct layers, was ∼117ð¼ 44.5 þ 73Þ μm. This simulated system was completely matched to the experimental system (Secs. 5.1 and 5.2). B. The existing PSF model 18 for a point light source in a specimen with uniform RI equal to 1.56. For this case, it was assumed that the axial extent of the specimen (composed only of optical cement) was 117 μm, while the PSF depths were taken in the range of 114 to 120 μm. C. The existing PSF model 18 for a point light source in a specimen with uniform RI ¼ 1.42 for a thickness of 117 μm. PSF depths were set as in B. D. The existing PSF model for a point light source in the specimen described in A approximated by a single averaged RI ¼ 1.50, calculated using the experimental conditions and Eq. (9), for the same thickness range and PSF depths as in B.
Experimental images are degraded by aberrations other than depth-induced SA, which could introduce restoration artifacts. Because these aberrations are not accounted by the PSF models, in order to show how the restoration deteriorates due to SA, a simulated image of a numerical spherical shell was also computed and restored with the same PSFs used in processing the experimental data. Comparison of simulated and experimental results shows a similar trend in restoration artifacts and validates that the artifacts observed in the restorations of the experimental data could be attributed to changes in the PSFs used for the restoration.
Restorations with these PSFs are summarized in Fig. 4. The experimentally acquired image [ Fig. 4(c)] of a 6-μm spherical shell exhibits spreading of intensities and loss of structural integrity, which is typical of wide-field images experiencing SA. The axial asymmetry of the intensity and the direction of the spread indicate the presence of positive SA. The experimental image restored with the PSF computed from the N-interface PSF model [ Fig. 4(d)] appears qualitatively to be the closest to a spherical shell, when compared to the restored objects obtained with PSFs based on a model that assumes the specimen has uniform RI [Figs. 4(e), 4(f), and 4(g) for the mountant, the optical cement, and the averaged RI, respectively]. Figures 4(e), 4(f), and 4(g) show artifacts because the PSFs used in each case for the restoration do not account for the RI variability within the sample. Figure 4(g) shows improved restoration when compared to Figs. 4(e) and 4(f), but it shows more artifacts compared to Fig. 4(d), indicating the need to account for specimen variability in the PSF.
Results from simulation (bottom row of Fig. 4) show consistency with restoration results from the experimental data (top row of Fig. 4) Fig. 4(c)] with respect to spreading and elongation of the spherical ring structure. The simulated image was restored using the matched N-interface PSFs, which were also used to generate the image (hence, there is no PSF mismatch in the restoration of the simulated image in [Fig. 4(j)], and consequently, the restored object [ Fig. 4(j)] converges to the true object. The artifacts observed in the restored object when the three-interface PSFs from the existing model are used instead [Figs. 4(k), 4(l), and 4(m), respectively] follow the trend of the artifacts observed in the reconstruction of the experimental data [Figs. 4(e), 4(f), and 4(g)]. The difference in the experimental restoration [ Fig. 4(d)] from the ideal simulated result [ Fig. 4(j)] suggests a residual model mismatch that could be attributed to inaccuracies in the experimental parameters used for the computation of the four-interface PSFs or the presence of aberrations (other than depth-induced SA) not accounted in the PSF model.
To emphasize the advantage of using the N-interface PSF model in restoration over the existing three-interface PSF model, intensity profiles from the restoration results with different PSFs are compared in Fig. 4(h). The profile through the restored object computed with the N-interface PSF [ Fig. 4

(d)]
is not only different from the restoration with the averaged RI [ Fig. 4(g)], but it also comes closer to the numerical object [ Fig. 4(b)], which approximates the actual spherical shell. This result highlights that accounting the effect of sample properties in the PSF model results in improved restoration of the sample's fluorescence intensity. To further quantify the quality of restoration, NMSE results are reported, which were computed using the numerical object as the reference. For the restorations obtained from the experimental data using the four different PSFs considered, the NMSE value is as follows: 1.80 with the PSF set A [ Fig. 4 Fig. 4(g)]. The smallest NMSE value is exhibited by the restoration with the PSF set A, from the new Ninterface model [ Fig. 4(d)], indicating that the restored object is closest to the simulated object. The discussion in this section allows us to conclude that the N-interface PSF model is advantageous for the restoration of samples with variable RI.

N-Interface PSF Model Evaluation Using a Biological Test Sample
In this section, experimental PSFs imaged from a biological test sample were compared to simulated PSFs using the N-interface model to evaluate the effectiveness of the model to predict PSF changes due to a biological sample.

Experimental PSF Acquisition from a Lung-Tissue Slice Test Sample
In this experiment, PSFs were acquired from images of 0.175μm-diameter microspheres, dyed with Alexa Fluor, placed below a rat lung slice dyed with calcein, and prepared to observe the effect of methacholine. For this purpose, a catheter was used to insert agarose in the tracheae and the sample was sectioned using the Vibratome into ∼150-μm sections. Because the air tubules were filled with agarose, the RI of the lung slice was approximated by an average RI ¼ 1.336. This value was derived from available information about the lung tissue (RI ¼ 1.342) 38 and agarose (RI ¼ 1.33). The lung slice submerged in ProLong Gold antifade reagent (RI ¼ 1.46) introduced additional RI variability within the sample. A layer of FluoSpheres (0.175 μm in diameter), dyed with Rhodamine, was air-dried on the coverslip to locate it. Fluo-Spheres, dyed with Alexa Fluor, were air-dried on the glass slide and the rat lung slice was carefully placed on top and left overnight in a cool and dark location. The top of the lung slice was demarcated using FluoSpheres, dyed with DAPI. The mountant was then placed on the lung slice along with the coverslip and was left to cure for two days allowing the mountant's RI to stabilize to 1.46.
The sample was imaged with a 40×∕1.3 NA oil-immersion lens (RI ¼ 1.515) using a Zeiss Axio Imager. A 677 × 765 × 300 size image with a voxel size equal to 0.16 × 0.16 × 0.1 μm 3 in the X, Y, and Z directions was acquired. The location of the microscope stage recorded at the best focus of the various beads was used to obtain the thickness of each layer (as described Sec. 3.1). The thickness of the mountant and the lung slice was derived to be 69 and 160 μm, respectively [ Fig. 5(a)]. From this information, the depth of the FluoSpheres below the coverslip was determined to be 229 μm. Four bright beads were chosen for the study [shown enclosed by the red box in Fig. 5(b)]. The averaged intensity distribution of the four beads was used to represent a single experimental PSF [ Fig. 5(c)].

Simulation of PSFs Matching Imaging
Conditions of the Lung-Tissue Slice Sample Simulated PSFs were calculated for two of the cases for comparison with the experimental PSF in Fig. 5(c). PSFs were generated as follows: (1) an N-interface PSF based on the experimental imaging parameters, i.e., a specimen with two different layers of media (one with RI ¼ 1.46 and thickness of 69 μm and another with RI ¼ 1.336 and thickness of 160 μm) and (2) a PSF using the existing Haeberlé model assuming a specimen with uniform RI ¼ 1.37 [the RI average obtained using Eq. (9) of the two actual RIs in the sample] and thickness equal to 229ð¼ 69 þ 160Þ μm. Since the size of the FluoSpheres is greater than the pixel size, a 2-pixel-diameter sphere was convolved with all the simulated PSFs for an improved comparison with the experimental PSF.

Comparison Between Experimental and Simulated PSFs
The experimental PSF image acquired using the conditions described in Sec. 6.2 is shown in Fig. 5(c), and the corresponding simulated PSF using the N-interface model is shown in Fig. 5(d). For comparison, the simulated PSF using a uniform averaged RI ¼ 1.37 is also shown in Fig. 5(e). Figure 5(f) compares axial intensity profiles through the simulated and experimental PSFs. There is better agreement between the simulated PSF generated using the N-interface PSF model and the experimental PSF than between the simulated PSF, which assumes the entire sample has a uniform RI ¼ 1.37, and the experimental one. There is disagreement between the simulated and experimental PSFs beyond the secondary lobes, which can be possibly attributed to residual intensity from surrounding structures. PSF comparison is quantified by the NMSE computed between each simulated PSF and the experimental PSF (which was taken as the reference). The NMSE was found to have values equal to 0.0182 and 0.0312 for the N-interface PSF and the existing model PSF (assuming specimen RI ¼ 1.37), respectively. This section further validates that the N-interface PSF model is able to predict PSF changes, due to thick biological sample imaging, more accurately than the existing three-interface PSF model.

Discussion and Conclusion
The new N-interface model can be used to obtain PSFs for specimens with variable RI. The model allows the use of microscope-derived parameters in a high precision vectorial model to obtain a PSF that allows the existence of more than three interfaces within the system. Our results emphasize the need of a model that can account for sample variance (such as layers with variable RI), which changes the amount of SA experienced by the light traveling through it. The N-interface PSF model was experimentally verified for a four-interface system by the demonstrated agreement between simulated and experimental data. Its accuracy over PSFs from an existing model, which can only account for a three-layer stratification, was also highlighted. PSFs from the N-interface model were also used to solve the inverse imaging problem using experimental data acquired from a bead test sample. The use of accurate PSFs generated using the N-interface model in the restoration as opposed to the PSFs generated using the existing PSF model reduced the NMSE between a numerical object (that approximates the physical test object) and the restored object at least 18% and at most 35% for the cases compared, indicating a definite improvement in the restoration. In the experimental case, the PSF mismatch is not exactly the same to the one predicted by the simulated case; thereby, the restoration artifacts are not exactly alike either. However, the overall consistency between the simulated and experimental results is very promising, and it suggests that our PSF model is able to capture the main features of the experimental system to the extent that our model-based restoration method can produce an improved restoration. Experimental PSFs from point light sources below a thick lung-tissue slice were also imaged and their average was shown to be predicted well by the N-interface PSF model, which represents the system better than the existing PSF model, quantified by a 42% reduction in the NMSE computed between the simulated and experimental PSF.
The proposed N-interface PSF model requires information about the RI at different layers within the sample. In experimental samples discussed in this work, we used subresolved spheres as markers in between layered media with known RI and we assumed that the specimen can be modeled by stratified layers. In biological samples, the RI information varies spatially and different phase imaging techniques have been under development over the recent years to enable determination of this variability. The new PSF model presented here will be used in a block-based forward imaging model to ensure accurate representation of the sample. This representation of the sample will improve restoration of the image, allowing COSM algorithms to effectively restore thick biological data. Future studies involve integration of a biological sample RI map in PSF determination suitable for SV restoration fluorescence microscopy.