2.1.

Theoretical Model

Based on Schmitt et al.⁴ and Izatt et al.,¹¹ the derived single scattering model of the OCT signal can be written as

This occurs because the oblique rays with longer path lengths are attenuated more than the on-axis rays, which effectively reduces the NA of the focusing optics. In our set up with $\mathrm{NA}=0.28$, this effect is negligible. In our fitting model, $z_r$ is assumed to be constant over the entire image.

The focus depth, on the other hand, is not constant over the image because of two confounding effects. First, we must correct for the field curvature effect that arises from the different path lengths that light travels to the image focus for different positions in the image. We calibrated this field curvature correction for our system by imaging a mirror placed at the image focus.²⁶ The second confounding effect is the likelihood that the surface of the brain sample being imaged is not parallel with the image plane. To correct for this effect, we parameterize $z_f$ within the imaging field of view (FOV) as

Finally, the back-scattering coefficient ${\mu }_b$ can be quantified with known $\eta$ and $P_0$. However, this would require knowledge of optical losses of all optical components in the OCT system. A relative back-scattering coefficient was previously introduced as

This relative back-scattering coefficient is system dependent and could be used to estimate the relative spatial variations in ${\mu }_b$.²

2.2.

Spectral-Domain OCT System

We used a commercial spectral-domain OCT system (Telesto TEL320C1, Thorlabs, New Jersey) to measure the optical properties of human brain tissue. The light source is a broadband superluminescent diode with center wavelength of 1300 nm and a full width half maximum bandwidth of 150 nm, yielding an axial resolution of $4.2\mu \mathrm{m}$ in tissue. The spectrometer has a 2048-pixel InGaAs line scan camera operating at an A-line rate of 76 kHz. The total imaging depth is 2.6 mm in tissue. A $10\times$ air objective (Mitutoyo) was used in the sample arm, which yields a lateral resolution of $3.5\mu \mathrm{m}$ with a theoretical Rayleigh range of $\text{40}\mu \mathrm{m}$ in a nonscattering medium. The maximum sensitivity of the system is 109 dB.

In spectral-domain OCT, the sensitivity drops with imaging depth due to the spectrometer resolution and discrete spectral sampling by the line-scan camera. The sensitivity roll-off function is given as²⁷

2.3.

Model Coefficient Fitting

A nonlinear least-squares solver was used for estimating the ${\mu }_s$ and ${\mu }_b^{\prime }$ parameters in the nonlinear theoretical model [Eqs. (17)]. The “trusted region reflective” algorithm was used for the optimization process. The trusted region reflective algorithm is an algorithm used in solving large-scale nonlinear equations. Specifically, this algorithm is a subspace adaptation of the Coleman-Li trust region and interior method.²⁸ Trust region methods form a respected class of algorithms for solving unconstrained minimization problems mostly due to their strong convergence properties. The algorithm first defines a trust region around the current best solution, then forces the size of the region to be updated iteratively only based on a sufficient decrease in the cost function. All the curve fitting was performed in MATLAB^®. Before estimating ${\mu }_s$ and ${\mu }_b^{\prime }$ for each position in an image, we first estimate the values of the constant $z_r$ and the $z_f$ tilt parameterized by Eqs. (2) and (3) by performing a simultaneous fit of OCT depth profiles from multiple positions within a region of the sample where the optical properties are known to be constant. This condition is easily met in phantoms with spatially uniform optical properties. For brain samples where the optical properties are spatially varying, we utilize the agarose material in which our brain tissue is embedded for creating the condition of spatially uniform optical properties for estimating $z_r$ and the parameterization of $z_f$. The value of $z_r$ is then taken to be constant through the entire brain sample, while the value of $z_f$ will vary linearly with $x$ and $y$ as prescribed by Eq. (6). Specifically, we fit the depth profiles from multiple locations within the agarose material around the brain tissue simultaneously to extract the value of $A$, $B$, and $C$ in Eq. (6). The linear coefficients ($B$ and $C$) are always constant for different spatial locations due to flattened surface, whereas the constant term ($A$) is spatially varying due to the angle of the vibratome blade varying from exactly horizontal. Since it is a linear variation, we interpolate the constant term along the direction of surface flattening to obtain the value of $A$ for each image tile.

For comparison with the full model, we fit the four parameters (${\mu }_s$, ${\mu }_b^{\prime }$, $z_f$, and $z_r$) independently using Eqs. (2) and (3) for each position within a region. Due to the nonlinearity of the model and the high inter-parameter dependency, a two-step procedure was used as previously described in Ref. 2. We fit the four-parameter model to the averaged experimental data and obtained an average $z_r$ for each brain sample. We repeat the fitting with $z_r$ fixed to this estimated value and then estimate the spatially varying optical properties and $z_f$.

2.4.

Phantom Experiment Set Up

A set of eight liquid optical phantoms was prepared by diluting different volume concentrations (4%, 8%, 12%, 16%, 20%, 40%, 60%, and 80%) of Intralipid 20% (batch 10LI4282, Fresenius Kabi, Germany) in deionized water to create phantoms with eight different scattering coefficients. Taking into account the properties of the different Intralipid components, pure Intralipid contains 22.7% (v/v) scattering particles (soybean oil and egg lecithin),^29–32 the concentration-dependent scattering coefficient can be calculated based on the Twersky quotation and corrected by packing factor³³

The eight liquid phantoms were imaged with a $10\times$ air objective. Each phantom was measured at 13 different focus depths, starting with $z_f=0$ and incrementing by $50\mu \mathrm{m}$ up to a depth of $600\mu \mathrm{m}$ in the phantom. The three shallow depths ($z_f=0$, 50, $100\mu \mathrm{m}$) were discarded from further analysis due to the specular reflection from the phantom surface. The deep depths ($z_f>300\mu \mathrm{m}$) were also discarded due to the poor signal-to-noise ratio causing a constant noise offset in the signal, systematically deviating the experimental data from the theoretical model. Each measurement produces an image that spans over a $1\times 1\times 2.5{\mathrm{mm}}^3$ volume with a $2.5\text{-}\mu \mathrm{m}$ isotropic pixel size, but with a lateral resolution of $4\mu \mathrm{m}$ and an axial resolution of $3.5\mu \mathrm{m}$. During acquisition, A-lines were repeated 10 times to improve the signal to noise ratio. The total acquisition time for each phantom was 13.4 s.

The raw data were resampled in $k$-space, compensated for dispersion, and Fourier transformed into the spatial domain following standard procedures.³⁵ Then the volume data were corrected for sensitivity roll off as described in Ref. 2. Finally, the noise floor, defined as the average intensity at an imaging depth between 2.4 and 2.6 mm where no signal is expected, was subtracted.

2.5.

Tissue Preparation and Imaging

A human brain tissue was obtained from the Massachusetts General Hospital Autopsy Suite. The brain was from a neurologically normal subject without previous diagnosis of neurological deficits. The sample was fixed with 10% formalin for at least two months. The postmortem interval did not exceed 24 h. The brain was cut into smaller blocks. Then the blocks of interest were embedded in oxidized agarose and covalently cross-linked with the agarose using borohydride borate solution. Embedded human brain blocks were washed for one month in phosphate buffer saline solution (PBS) 0.01 M at room temperature (RT) while gently shaking. Index matching was performed with serial incubations in 100 ml of 20%, 40%, and 60% (v/v) TDE in 0.01 M PBS (TDE/PBS) each for 24 h at RT while gently shaking.

One human brain block from the primary motor cortex was imaged with a $10\times$ air objective (Mitutoyo) and a $z$-spacer. The surface of the embedded sample was flattened using a vibratome slicer before OCT imaging.³⁶ The volumetric imaging spanned a FOV of $1\times 1\times 2.5{\mathrm{mm}}^3$, with a voxel size of $2.5\text{-}\mu \mathrm{m}$ isotropic. The focus of the light beam was located at roughly $150\mu \mathrm{m}$ below the surface. Motorized $xyz$ stages were integrated under the sample arm to automatically translate the sample between image tiles. The $12\times 12{\mathrm{mm}}^2$ section of the tissue block was imaged with a total of 576 tiles ($24\times 24$) that were acquired with 50% overlap between adjacent tiles to ensure reliable registration and to uniformly reduce speckle noise. Fiji software was used to stitch the image tiles and create the mosaic images.³⁷

To further reduce the effect of noise on fitting depth-profiles and estimating the optical properties, the volumetric images were averaged over $25\mu \mathrm{m}$ in $x$ and $y$ (i.e., $10\times 10\text{pixels}$) and down-sampled in $x$ and $y$ to a $2.5\text{-}\mu \mathrm{m}$ voxel size. The location of the tissue surface was determined on cross-section images using an image processing edge detector.³⁸ A total depth range of $250\mu \mathrm{m}$ was selected for depth-profile fitting, starting from $35\mu \mathrm{m}$ below the tissue surface to reduce the impact of specular reflection from the surface. To compare with conventional OCT contrast, the AIP and maximum intensity projection maps within the same depth range were computed as well.

To further validate our parametrized model, the quantified scattering coefficient of the human brain sample was compared with a conventional low-NA measurement with a negligible axial PSF effect where the OCT signal is modeled as

We imaged the same section of the brain sample with a $2\times$ air objective (Thorlabs) and a $z$ spacer, which yields a lateral resolution of $15\mu \mathrm{m}$ and a theoretical Rayleigh range of $\text{280}\mu \mathrm{m}$ in tissue. The volumetric image tiles have a FOV of $2\times 2{\mathrm{mm}}^2$ with a $5\times 5\times 2.5\mu {\mathrm{m}}^3$ voxel size and 50% overlap between adjacent tiles. A linear function was fit to the logarithmic OCT signal depth profile and the slope was used to extract the ${\mu }_s$ value.

3.1.

Optical Properties Characterization of Liquid Phantoms

We first fit our parametrized model to the OCT signals acquired from the eight Intralipid phantoms measured at different focus depths. Example fits of the OCT signal depth dependence are shown in Figs. 1(a) and 1(b). In Fig. 1(a), we show example fits of the OCT signal depth profile of liquid phantoms with different volume concentrations (0.8%, 2.4%, 4%, and 8% v/v Intralipid). In Fig. 1(b), we show example fits of the OCT signal depth profile of 0.8% v/v Intralipid phantoms with different focus depths ($z_f=150$, 250, 350, and $450\mu \mathrm{m}$). In Fig. 1(c), the estimated scattering coefficient (${\mu }_s$) values are presented as a function of the volume concentration of Intralipid measured at two different focus depths ($z_f$). Error bars indicate the standard deviation of the estimated ${\mu }_s$ across the entire FOV. The error bars are rather small, indicating the high repeatability of the measurements and the coefficient estimation procedure. At high concentration, ${\mu }_s$ tends to be underestimated compared with a prior model of the concentration dependence of ${\mu }_s$ for Intralipid,³² probably due to the impact of multiply scattered photons becoming important for higher scattering coefficients.

Fig. 1

Optical properties of Intralipid liquid phantoms. (a) Representative OCT signal depth profiles and model fitted curves of Intralipid phantoms with different volume concentrations (0.8%, 2.4%, 4%, and 8% v/v) at focus depth of $300\mu \mathrm{m}$. (b) Representative OCT signal depth profiles and model fitted curves of 0.8% v/v Intralipid phantom at different focus depths (150, 250, 350, and $450\mu \mathrm{m}$). (c) Estimation of scattering coefficient (${\mu }_s$) of Intralipid phantoms with respect to Intralipid volume concentration at two different focus depths. Error bars: standard deviation. The solid line indicates a previous theoretical prediction for the scattering coefficient. (d) Effective Rayleigh range ($z_r$) as a function of depth of focus ($z_f$) with respect to Intralipid volume concentration. Deep depths of focus were discarded due to the non-negligible effect of multiple scattering. (e) Shape of $z_f$ of 0.8% v/v Intralipid phantom within a FOV with fixed $z_r$ in model fitting. (f) Shape of $z_f$ of 0.8 % v/v Intralipid phantom within a FOV with freely estimated $z_r$ in model fitting.

In Fig. 1(d), the estimated Rayleigh range $z_r$ is plotted as a function of both the volume concentration of Intralipid and the focus depth. We found that the estimated $z_r$ differed by $<1\mu \mathrm{m}$ across different focus depths for each volume concentration. This indicated that we could assume $z_r$ to be constant across the FOV. We next evaluated the spatial dependence of $z_f$ within the FOV both with $z_r$ set as a constant and as a freely fitted parameter that was allowed to vary over the FOV, as shown in Figs. 1(e) and 1(f), respectively. We found that the spatial dependence of $z_f$ was almost identical with $z_r$ set as a constant or as a freely fitted parameter. Moreover, when $z_r$ was set as a freely fitted parameter, the coefficient of correlation for the estimated $z_f$ and $z_r$ was calculated to be $-0.043$, indicating a strong independence of these two parameters.

To quantify the improvement in fitting quality using the parametrized model compared with the full fitting model, we calculated the coefficient of variance (CV), which is defined as the standard deviation divided by the mean of extracted ${\mu }_s$ and ${\mu }_b^{\prime }$ over the imaging field, as shown in Fig. 2. We first estimated both ${\mu }_s$ and ${\mu }_b^{\prime }$ in an 0.8% volume concentration Intralipid phantom with varying focus depths ($z_f$), as shown in Figs. 2(a) and 2(b). We found that our parametrized model showed a smaller CV in both ${\mu }_s$ and ${\mu }_b^{\prime }$ for all focus depths. We then estimated both ${\mu }_s$ and ${\mu }_b^{\prime }$ in a set of liquid phantoms with different concentrations (0.8%, 1.6%, 2.4%, 3.2%, 4%, and 8%) with a fixed focus depth of $150\mu \mathrm{m}$ below the liquid surface, as shown in Figs. 2(c) and 2(d). We found that our parametrized model showed a smaller CV in both ${\mu }_s$ and ${\mu }_b^{\prime }$ for all scattering levels.

Fig. 2

Coefficient of variance (CV) of estimated parameters of liquid phantoms shows improvement of fitting quality. (a) CV of estimated ${\mu }_s$ of 0.8% volume concentration Intralipid phantom over imaging field with various focus depths of both parametrized model and full model. (b) CV of estimated ${\mu }_b$ of 0.8% volume concentration Intralipid phantom over imaging field with various focus depths of both parametrized model and full model. (c) CV of estimated ${\mu }_s$ of Intralipid phantoms at $150\text{-}\mu \mathrm{m}$ focus depth over imaging field with various scattering levels of both parametrized model and full model. (d) CV of estimated ${\mu }_b$ of Intralipid phantoms at $150\text{-}\mu \mathrm{m}$ focus depth over imaging field with various scattering levels of both parametrized model and full model.

3.2.

Optical Properties Characterization of Human Brain

3.2.1.

Validation with low-NA model

To validate the quantitative estimates of the optical properties of the brain tissue derived from the proposed parametrized model with the conventional method using a single exponential model, the same section of a human brain sample was imaged with both a high NA (10×, Mitutoyo) and a low-NA ($2\times$, Thorlabs) objective. The impact of axial PSF was negligible for the low-NA objective since the Rayleigh range was longer than the depth range analyzed. We first plot the back-scattering (${\mu }_b^{\prime }$) and scattering (${\mu }_s$) maps along with AIPs for both the low-NA and high-NA objectives, as shown in Fig. 3. Moreover, for the high-NA objective, we plot the ${\mu }_b^{\prime }$ and ${\mu }_s$ maps using both the parametrized model [Fig. 3(b)] and the full model [Fig. 3(c)]. To demonstrate the ability of our parametrized model to quantitatively extract the optical properties of human brain tissue, two regions of interest (ROIs) were manually selected, including a $2\times 2{\mathrm{mm}}^2$ GM region and a $2\times 2{\mathrm{mm}}^2$ WM region, as shown in Figs. 3(a) and 3(b). The averaged scattering coefficient (${\mu }_s$) within these two ROIs are plotted in Fig. 3(d). The estimated ${\mu }_s$ of WM were $6.43\pm 0.96{\mathrm{mm}}^{-1}$ and $6.12\pm 0.99{\mathrm{mm}}^{-1}$ for the low-NA and high-NA objectives, respectively. The estimated ${\mu }_s$ of GM were $3.32\pm 0.89{\mathrm{mm}}^{-1}$ and $3.58\pm 0.64{\mathrm{mm}}^{-1}$ for the low-NA and high-NA objectives, respectively. We then performed a pairwise t-test on the extracted ${\mu }_s$ from both ROIs using the full and the parametrized models to determine if there is a significant difference in the mean. We found no significant difference in the mean of the estimated ${\mu }_s$($p=0.49$ for the GM ROI and $p=0.25$ for the WM ROI). We further performed a pairwise t-test on the extracted ${\mu }_s$ from both ROIs for the low-NA and high-NA objectives to determine if there is a significant difference in the mean. We found no significant difference in the mean of the estimated ${\mu }_s$($p=0.26$ for the GM ROI and $p=0.09$ for the WM ROI). Although the scattering coefficient is comparable for the two objectives, as expected, we see that structures such as vessels and axonal fiber tracts are more detailed in the high-NA images. As shown in Figs. 3(b) and 3(c), the full four-parameter model overestimated the scattering coefficient of GM and had higher standard deviation within both ROIs compared to the parametrized and the low-NA exponential model. Recall that we take the low-NA measurement as the ground truth result since the model fitting is a simple exponential with no dependence on the confocal parameters.

Fig. 3

Comparing the scattering coefficient map of primary motor cortex imaged by a low-NA $2\times$ and higher-NA $10\times$ objective. (a) The relative back-scattering (left), scattering coefficient map (middle), and AIP (right) of the brain section with a $2\times$ objective. Some but not all vessels and fiber tracts are labeled with blue arrows. Scale bar: 1 mm. (b) The relative back-scattering (left), scattering coefficient map (middle), and AIP (right) of the brain section with a $10\times$ objective using the parametrized model. Some but not all vessels and fiber tracts are labeled with blue arrows. Scale bar: 1 mm. (c) The relative back-scattering (left) and scattering map (middle) of the brain section with a $10\times$ objective using the full model. Some but not all vessels and fiber tracts are labeled with blue arrows. Scale bar: 1 mm. (d) Averaged ${\mu }_s$ within selected ROIs for both high NA and low-NA objectives. Both parametrized model and full model were used to extract ${\mu }_s$ under high-NA objective.

3.2.2.

Comparison with the full fitting model

We compare the results from our parametrized model with results from the previously established full four-parameter estimation model by imaging the same section of a human motor cortex sample. The same depth range data was used to fit both the full four-parameter model and our proposed parametrized model. To show the improvement for our parametrized model, we generated the CV images for both the back-scattering (${\mu }_b^{\prime }$) and scattering (${\mu }_s$) maps for both the parametrized model and the full model using a $5\times 5\text{pixel}$ sliding window, as shown in Figs. 4(a) and 4(b). The value of each pixel in the CV images is defined as the standard deviation of values over the sliding window divided by the mean of values over the sliding window. We observed from the CV images that our proposed parametrized model showed a significantly smaller CV compared to the full model. The histograms of the CV values are shown in Fig. 4(d). The median of the CV for ${\mu }_s$ was 0.28 for the full model and was reduced to 0.13 for the parametrized model. Similarly, the median CV of ${\mu }_b^{\prime }$ was 0.57 for the full model and was reduced to 0.26 for the parametrized model. This confirms that our proposed parametrized model is effective at reducing inter-parameter covariation and can provide more precise estimates of the optical properties while maintaining accuracy.

Fig. 4

Comparison between the full model and the proposed parametrized model. (a) CV image of the extracted ${\mu }_s$ using the parametrized model. (b) CV image of the extracted ${\mu }_s$ using the full model. (c) CV image of the extracted ${\mu }_b$ using the parametrized model. (d) CV image of the extracted ${\mu }_b$ using the full model. (e) Histograms of CV for the estimated optical properties. Left: CV of the extracted ${\mu }_s$ using both the parametrized and full models. Right: CV of the extracted ${\mu }_b$ using both the parametrized and full models.

3.2.3.

Optical properties of human brain before and after index matching

A section of a superior frontal cortex brain sample was first surface flattened using a vibratome and imaged with the high-NA objective ($10\times$, Mitutoyo). Then the sample underwent incremental index matching with 20%, 40%, and 60% v/v TDE/PBS. The tissue section was imaged with the high-NA objective at each TDE concentration after it reached equilibrium. Due to nonuniform distortion introduced by index matching, the surface of the tissue block was reflattened at each TDE volume concentration before imaging with the high-NA objective. As shown in Figs. 5(a)–5(d) and 5(j), the AIP indicates a reduced GM/WM contrast after index matching. The GM/WM contrast is defined as the mean value of the WM intensity divided by the mean value of the GM intensity. Prior to index matching, the averaged intensity of WM over the depth range is smaller compared to GM due to the higher scattering coefficient of WM. After index matching, we see that the WM AIP becomes brighter as the ratio between WM and GM intensity reached 0.95 at 60% v/v TDE/PBS equilibrium, suggesting that the index matching process reduced the scattering coefficient of WM more than GM and resulted in a more homogenous scattering level across the sample, as indicated in the extracted ${\mu }_s$ map, as shown in Figs. 5 (e)–5(h).

Fig. 5

Estimating optical property change during index matching. (a) AIP of the brain section before index matching. (b) AIP of the brain section after equilibrium in 20% v/v TDE/PBS. (c) AIP of the brain section after equilibrium in 40% v/v TDE/PBS. (d) AIP of the brain section after equilibrium in 60% v/v TDE/PBS. (e) Scattering map of the brain section before index matching. (f) Scattering map of the brain section after equilibrium in 20% v/v TDE/PBS. (g) Scattering map of the brain section after equilibrium in 40% v/v TDE/PBS. (h) Scattering map of the brain section after equilibrium in 60% v/v TDE/PBS. (i) Contrast between WM and GM during each TDE concentration. (j) Extracted scattering coefficient from WM and GM ROI during each TDE concentration. Error bar: standard deviation.

Two $3.75\times 3.75{\mathrm{mm}}^2$ ($150\times 150\text{pixel}$) ROIs were selected, one from GM and one from WM, as shown in Figs. 5(e)–5(h). The average ${\mu }_s$ of each ROI was calculated to quantify the scattering coefficient change with index matching at each TDE concentration, as shown in Fig. 5(k). We found that after index matching to 60% v/v TDE/PBS, ${\mu }_s$ of GM decreased from $2.29\pm 0.79{\mathrm{mm}}^{-1}$ to $0.45\pm 0.27{\mathrm{mm}}^{-1}$ and that ${\mu }_s$ of WM decreased from $5.61\pm 0.86{\mathrm{mm}}^{-1}$ to $0.77\pm 0.65{\mathrm{mm}}^{-1}$. We found that the scattering coefficient of WM decreased to $3.06\pm 0.73{\mathrm{mm}}^{-1}$, $1.97\pm 0.89{\mathrm{mm}}^{-1}$, and $0.77\pm 0.65{\mathrm{mm}}^{-1}$ at 20%, 40%, and 60% v/v TDE/PBS equilibrium, respectively, whereas the scattering coefficient of GM decreased to $0.25\pm 0.17{\mathrm{mm}}^{-1}$ and $0.57\pm 0.37{\mathrm{mm}}^{-1}$ at 20% and 40% v/v TDE/PBS equilibrium, respectively, and then reached $0.45\pm 0.27{\mathrm{mm}}^{-1}$ at 60% v/v TDE/PBS equilibrium. The different behavior of GM and WM under different TDE concentration index matching is a result of the expected lower average index of refraction for GM relative to WM.³⁹ To better justify the use of the parametrized model, we also fitted the scattering coefficient at each TDE concentration with the full model, as shown in Fig. 5(k). We found that, in general, both models agreed with each other on the extracted scattering coefficient whereas the standard deviation of estimated ${\mu }_s$ was smaller for the parametrized model at each TDE concentration.