Inverse spectroscopic optical coherence tomography (ISOCT) for characterization of particle size and concentration

Inverse spectroscopic optical coherence tomography (IS-OCT) methods apply inverse problem formulations to acquired spectra to estimate depth-resolved sample properties. In the current study, we modelled the time-frequency-distributions using Lambert-Beer law and implemented IS-OCT using backscattering spectra calculated from Mie theory, then demonstrated the algorithm on polystyrene microspheres under idealized conditions. The results are significant because the method is compatible with high concentrations and broad diameter distributions of scattering particles such as organelles in biological cells.


Introduction
In medicine and life sciences, the characterisation of light scattering, absorption and polarisation in biological tissues reveals information on their underlying structure and composition, allowing for diagnosis and monitoring in diverse applications such as arterial plaque in atherosclerosis, the size and nature of tumours, blood perfusion and oxygen saturation levels, and drug and nanoparticle delivery [1]. In particles, at the cellular and subcellular size-scale, light scattering and absorption are wavelength dependent (spectral) properties that can be described using Mie theory [2], a well-established stochastic model that forms the basis for the Monte Carlo model of multiple light-particle interactions [3][4][5][6]; alternatively, more recent scattering models also exist that include electromagnetic wave descriptions of light [7,8].
In optical coherence tomography (OCT), low coherence spectroscopy (LCS) and low coherence interferometry (LCI), light from a low coherence source interact with a sample and a Michelson interferometer is used to extract depth-dependent information about the sample from the backscattered light [9]. In spectroscopic OCT (S-OCT) and LCS, broadband illumination and spectrally sensitive signal detection allows for characterisation of spectrally resolved absorption and scattering properties [6] which have been used to quantify media composition from absorption [10][11][12][13][14] and, when combined with Mie scattering models, to quantify particle size [4,5,[14][15][16][17][18][19][20][21]. The use of inverse problem approaches to analyse particle size distribution in a sample is prevalent in systems using forward-scattering of light (see for example: [22,23]) but is less common in backscattered light techniques such as OCT where it is known as inverse S-OCT (IS-OCT) [19,20,[24][25][26] In this paper, we briefly review the theory of S-OCT to show where each of the existing IS-OCT techniques are derived from [19,20,[24][25][26] before extending the technique presented in [26], that used Lambert-Beer's law to describe the time frequency distributions, to backscattering spectra from Mie theory.
Where these earlier techniques have demonstrated capability in estimating particle size [19,20,[24][25][26], the technique presented here provides a quantitative estimate of particle size and concentration. Results are presented from data acquired on solutions of polystyrene microspheres with mean diameters of 0.17, 1.04 and 8.91 µm at concentrations of 0.156, 0.625 and 2.5 volume %. The results demonstrate the capability of IS-OCT to separate experimentally acquired spectra into several Mie scattering components. The results are significant because biological tissues are largely heterogenous on the micro-scale of organelles within cells and so there is a need for IS-OCT analysis methods which are can operate in the presence of high particle concentrations and broad diameter distributions.

S-OCT
In OCT, the fields originating from the sample arm and the reference arm interfere with one another to produce the interferogram detected by the photodetector, : where = is the wavenumber, the wavelength, 2 the difference in optical path length between the two arms, and angled brackets indicate an average over the detector acquisition time which is assumed to be much longer than the time of an optical cycle. The reference and sample arm fields are functions of the source spectral density ( ) and the optical properties of the reference mirror and sample which modify the spectrum, such that [27]: where ( , ) is the spectral backscattering profile of the light reflection within the sample, ( , ) the spectral modification by media before the scatterer, and ( , ) the spectral modification by optical components in the system.
In a swept source OCT system, a laser with narrow linewidth, , spectrally sweeps across the light source bandwidth while the photodetector acquires the interferogram ( ), which is then Fourier transformed, ℱ, to produce the sample depth dependent OCT signal ( ) [6]: Assuming that only light from single scattering events are present in the interferogram (first order Born approximation), that light attenuation is due to scattering and absorption (Lambert-Beer's law), and plane wave illumination from the light source, produces the general form for ( ) of [28][29][30]: where is a system constant which describes the influence of system components such as the photodetector, ℎ( ) is a function which describes the depth-dependent influence of system components such as beam divergence, the backscattering coefficient, and = + is the attenuation coefficient which is the sum of the scattering and absorption coefficients, and .
With swept source OCT, the spectral variation in and is encoded within the interferogram ( ) as the laser sweeps across the source bandwidth [6]. Extraction and analysis of the depth-resolved spectral variation from the spectral interferogram is known as spectroscopic OCT (S-OCT), and is performed using a timefrequency-transformation (TFT) method such as a short-time Fourier transform (STFT), double window (DW), or the wavelet transform [6,31]. In the case of STFT, a window, ( ), is scanned across the interferogram and where ( ) is a window function with spectral width ∆ . The process of windowing introduces a trade-off between the spectral and axial resolution in the resultant TFD due to the uncertainty principle, which states that sweeping a narrow window allows high spectral resolution but poor axial resolution and vice versa.

Parameters
Researchers have employed a range of experimental configurations using OCT to extract , , and from Eqn. 4 for tissues and tissue phantoms [12,13,15,29,31]. Once extracted, the spectral variation of and allows characterisation of particle diameter using Mie theory [5,16,32], while the spectral variation of aids Figure 1: In conventional OCT, an A-scan is produced by taking the Fourier transform of the entire bandwidth of the interferogram, producing high depth resolution but no spectral resolution. In S-OCT, as the window slides across the interferogram, each STFT produces an A-scan for a narrow spectral band, shown here as different coloured lines. The depth-resolved spectral profiles are then generated by extracting depth indexed values from this set of A-scans. The wideband and narrowband A-scans can all be described using Eqn. 4.
in tissue classification [31]. The parameters , and are local material properties which, for identical scattering particles, can be calculated from the product of the number density of particles, , the cross-sectional area of the particle, , and either the backscattering or scattering efficiency, [2]: where the subscript indicates either backscattering or scattering, is the effective cross section, ( , , , ) are functions of the particle diameter, , light wavelength, , and the refractive indices of the particle and the medium which can be found using Mie theory [2,33].
While it is not the focus of the current study, other methods have resolved the depth dependent spectra by TFT of the A-scan [31], ( ), and then estimated the particle size from autocorrelation of the resultant spectra [ , ; ( )] [4,19,34]. A further, alternative method for characterising and is to describe spatial fluctuations in refractive index in the sample using the refractive index (RI) spatial auto-correlation function. In Refs. [24][25][26], the authors model this RI function using the Whittle-Matern family of functions that includes a shape factor determining the type of function, and then derive expressions for and using the shape factor among other parameters.

IS-OCT to estimate particle size distribution
The linear inverse problem = can be used to estimate a sample's properties, , from experiment or modelling data, , via inversion of the coefficient matrix, , which relates the two. Researchers have applied this linear problem formulation to IS-OCT using several different parameters: in Ref. [20], the diameters of microspheres were estimated from spectral interferogram data, Eqn. 2, and a coefficient matrix that was populated with Mie spectra for different diameter particles each evaluated across different depths; in Ref. [19], the diameters of microspheres were estimated from frequencies in the TFT of the depth interferogram ( ) and a coefficient matrix that was populated using principal component analysis on training data; and, lastly, in Refs. [24][25][26], and the contributions of and to were estimated from the depth dependent spectra, [ , ; ( )], Eqn. 5, using the absorption spectra of chlorophyll and the Whittle-Matérn function derivations of scattering and backscattering.
In this section, we apply the linear inverse problem formulation to estimate particle size and concentration from depth dependent spectra, [ , ; ( )], Eqn. 5, using a coefficient matrix that is populated with Mie backscattering efficiencies. Using the depth dependent spectra, as in Ref. [26], permits incoherent averaging of multiple data sets to reduce speckle [35] without loss of signal amplitude -which can occur in the spectral interferogram method [20] if phase instability is present in the spectra, as is common in some wavelength-swept lasers.
The depth dependent spectra can be modelled in the same format as Eqn. 4 because they are constructed from a set of intensity profiles of this form, but with each parameter now a 1 × vector of values corresponding to each of the spectral ranges analysed in the TFT. For an arbitrary depth, : where the operator ⨀ indicates Hadamard multiplication and is equivalent to converting the vector of values produced in each function into a diagonal matrix and then performing matrix multiplication. When there are multiple scattering particles in the imaging volume, the voxel intensity of the imaging volume is the coherent sum of the backscattered spectra from each of the particles in the imaging volume [32], albeit with speckle noise.
Therefore, for an imaging volume containing a set of particles with different diameters, the spectral modification by sample backscattering can be expressed as the sum of the backscattering coefficients from the particles in each diameter range: where is the number density of diameter range , , is the effective backscattering cross section for a single particle of diameter range . Eqn. 8 can be expressed in matrix notation as = , where is a × matrix of backscattering cross sections for a single particle analysed for different diameters which can be populated using Mie theory, and is a × 1 vector of number densities for each particle diameter. Substituting = into Eqn. 7 and rearranging to solve for yields: where ϯ indicates the pseudoinverse. If of the sample is known, then Eqn. 9 provides a means to estimate the diameter distribution of scattering particles, , in the imaging volume at . In a homogeneous sample, can be estimated by fitting a straight line to the logarithm of the set of intensity profiles (or A-scans) [ , ; ( )]; while in a heterogeneous sample, the local attenuation coefficient can estimated by dividing the intensity in the pixel of interest by the sum of the intensities of the pixels located at deeper depths [36].

Coefficient matrix inversion
For the linear problem = , the Moore-Penrose pseudoinverse method for inversion of the coefficient matrix is equivalent to the least squares method that minimizes the expression ‖ − ‖ : where superscript is the transpose, is the singular value matrix and and are orthogonal matrices produced by singular value decomposition of , is the number of non-zero singular values and the subscript indicates the first columns of the corresponding matrix [37]. In the current study, where represents number densities, negative values are not physically possible, and so a non-negative constraint can be imposed by changing negative values in to 0. Tikhonov regularisation can stabilize the inverse problem by introducing a penalty term, ‖ ‖ , with regularisation parameter, , and then minimising ‖ − ‖ − ‖ ‖ . Where in Tikhonov regularisation the L2 norm in the penalty term smooths the output by penalising large and sparsely distributed values in , using the L1 norm on the penalty term promotes sparse solutions -i.e. those with several 0's in [38].
In IS-OCT studies, researchers have used the Moore-Penrose pseudoinverse [19] and Nelder-Mead algorithm [26] to minimise the expression ‖ − ‖ , and an iterative thresholding algorithm which minimises the expression ‖ − ‖ − ‖ ‖ with the L1 norm implemented on the penalty term. The IS-OCT problem becomes non-linear when cannot be extracted from the data or if non-linear effects such as multiple scattering events are included. In these scenarios, non-linear formulation and inversion techniques are required, for which there are several techniques that are commonly applied to forward scattering problems [23].

OCT system
The OCT system is polarisation sensitive with separate reference arms for horizontal and vertical polarisation and a quarter waveplate in the sample arm to circularly polarize the incident beam, and has been described in detail in earlier studies [39]. In the current study, this polarisation sensitivity was not utilised because the samples were not birefringent; the total back-reflected intensity was calculated as the sum of horizontal and vertical polarisation intensities. A-scans were acquired at the laser sweep rate. This scanning approach produced uncorrelated speckle noise across the A-scans that was used to reduce speckle noise during data processing by averaging the resultant TFDs from a total of 80,000 A-scans [35].  Table 1. The concentrations were selected to be representative of organelles in biological cells, which for example can vary from 0.01 to 20 % in yeast [40].

Data processing
TFDs were generated for each of the 80,000 interferograms by stepping a 20 nm Gaussian window at 0.  Backscattering efficiencies were calculated using Mie theory equations implemented in MATLAB [33] across the wavelength range of the source, of 1.259 to 1.370 µm (from product test sheet supplied by manufacturer), and for diameters from 0.1 to 10 in 0.01 µm increments. The refractive indices used for water and polystyrene microspheres were 1.3225 and 1.59, respectively. The wavelength range of each backscattering efficiency array was then clipped to 1.2771 to 1.3594 µm to account for the window size of the TFT and the data acquisition method that acquired the first 1024 points of the total 1120 points generated by the laser source.
The backscattering efficiency arrays spanned two orders of magnitude, and so each was normalised against its maximum value to reduce sensitivity to noise and improve matrix conditioning, and then populated into columns of a matrix, Fig. 3a. The distance matrix, which calculates the Euclidean distance between each column, revealed the backscattering efficiencies had minimal variation across the range of particle diameters when observed across the narrow bandwidth of the laser source, particularly across 0.1 to 2 µm range that encompassed the 0.17 and 1.04 µm diameter samples, Fig. 3b. Therefore, to construct the coefficient matrix, a reduced diameter range was implemented which encompassed the diameter ranges sample: 0.10, 0.14, 0.18, 0.22; 0.9, 1.0, 1.1, 1.2; and, 8, 8.4, 8.8, 9.2, 9.6, 10.0, Fig. 3c. Columns with a Euclidean distance less than 1 were removed to improve the condition number of the coefficient matrix to a value comparable to the SNR of the TFDs, Fig. 3d.
This targeted and restricted diameter range essentially limits the reconstruction to a classification problem but may be improved with broader bandwidth sources.
The particle diameter distribution was estimated by solving Eqn. 9 using the Lawson-Hanson algorithm [41] that iteratively finds a solution to the least squares expression with a non-negative constraint using an active-set method and is implemented in the MATLAB function lsqnonneg.

Results
For each microsphere sample, the extracted depth-dependent spectra were used to estimate depthdependent diameter estimates by finding the least-squares fit that minimized the residual error ‖ − ‖ , The spectra are very similar across this narrow bandwidth. The coefficient matrix used in the current study (c) employed a limited number of microsphere diameters (c) with similar spectra removed by implementing a threshold on the Euclidean distance (d).
column). The concentrations were consistently underestimated by a factor of 4, which was likely due to the contribution of the large spectral patterns at 1.1 µm. The depth reconstructions of the 1.04 µm sample contained clear spectral patterns corresponding to a diameter of 1.1 µm, but also some spectral patterns at smaller at larger diameters, Fig. 4d (centre column). The concentration was underestimated in the 0.15 % sample, at around 0.04 % total, but was reasonably close in the other two samples, at 0.4 % in the 0.625 % sample and 2 % in the 2.5 % sample. The depth reconstructions of the 8.9 µm sample contained spectral patterns across the Figure 4: Diameter distribution estimates for each microsphere sample shown graphically with feature labels for interpretation (a). Each sample graphic contains a plot of wavelength versus normalized magnitude for an examplespectra (solid line) and its least squares fit (dashed line) taken at the depth indicated by the black diamondtipped line connecting the plot to the adjacent colour-map. In the colour-map, the colour-scale indicates microsphere concentration in units of volume % for each microsphere diameter (x-axis) and sample depth (y-axis). Along the colourmap's x-axis, a three-colour bar indicates the diameter ranges for the three microsphere samples used in the current study. Along the colour-map's y-axis, a heat map indicates the residual error of the least squares fit ‖ − ‖ at the corresponding depth. The spectral patterns contained in the coefficient matrix are shown with three Mie scattering spectra that were dominant in the sample results highlighted in bold to aid interpretation (b). In the current IS-OCT study, the coefficient matrix and sample spectra are used to estimate the diameter distribution by rearranging the equation = . Nine graphics present results for the nine samples listed in Table 1 diameter range, Fig. 4e (right column). The spectral patterns in the 8 to 10 µm range were typically near the correct concentration in the 0.15 % and 2.5 % samples.
In the 0.17 µm and 1.04 µm samples, oscillations in the TFDs are thought to be caused by speckle and phase instability in the laser source and are interpreted by the IS-OCT algorithm as larger diameter particles.
Similarly, in the 8.9 µm sample, the large concentration estimates from small diameters are thought to be caused by the same errors -speckle and phase instability in the laser source -but are amplified by the difference in magnitude of the backscattering coefficients for small and large particles.

Discussion
This study has demonstrated the potential application of IS-OCT to quantify the depth dependent diameter and concentration of scattering particles using backscattering spectra from Mie theory. Mie scattering in OCT [15] and LCS [16] has been used to estimate particle size and concentration in samples without depth dependence, and in S-OCT to estimate the size of individual particles with depth dependence [5,19,20]. The method presented in the current study is both depth dependent and compatible with high concentrations of scatterers, such as organelles in biological cells or highly turbid media, where multiple scattering particles in the imaging volume prohibits an individual assessment of each scattering particle. Noise reduction through incoherent averaging and data filtering approaches are essential because errors from speckle and system instability are interpreted as scattering particles.
The limited variation in the Mie backscattering spectra for different particle diameters when observed across a narrow bandwidth hindered the inverse solution stability and was a major limitation of the method in the current study. Shorter wavelength and broader bandwidth light sources can reduce this instability by capturing more of the oscillations in the Mie backscattering spectra, Fig. 5, and may be achieved using visible wavelength and dual-band (visible and near infrared) OCT systems [42,43]. With some prior knowledge of the sample, another method to improve solution stability is to remove some of the spectra, as was performed in the current study for the microsphere experiments. Prior knowledge of the sample can also be built into the algorithm in much the same way the ( ) was implemented; for example, a reference spectral profile might be generated from modelling or measurements across a population, and then spectral deviations from this reference used to indicate changes in size or concentration of scatterers in the sample, where high frequency deviations would indicate a change in large scatterers and low frequency changes a change in small scatterers.
In Ref. [26], prior information was used to great effect to split the attenuation coefficient into absorption from chlorophyll and scattering from the sample structure.
The frequency of oscillations in the Mie backscattering spectra can also be used to populate the coefficient matrix and is the basis for the interferogram IS-OCT implementation used in Ref. [20]. Here, the authors accurately estimated the diameters of 1.0 and 4.5 µm polystyrene particles in air attached to glass slides though with deteriorating accuracy as the spheres occupied the same imaging volume -highlighting the problem of speckle in quantitative measurements as observed in the current study. The IS-OCT studies that fitted Whittle-Matern functions to the TFDs [24][25][26] did not describe speckle as a problem; possibly because the Whittle-Matern functions characterise the backscattering coefficient using a power law of the form , where is the shape factor and is greater than zero, which limits the number of oscillations in the fitted curve. While this limitation is not a problem for small diameter microspheres [24], and biological tissues in mammal [24,25] and coral [26], it may not extend well to large diameter scatterers or broader bandwidth OCT systems.
Limitations of the current study were the limited number of diameters used in the coefficient matrix because of the long central wavelength and the narrow bandwidth laser source; the use of homogeneous microsphere samples meant that the depth dependent capabilities of the method were not demonstrated; and, multiple scattering events were not considered when extracting the attenuation coefficient. To improve the method in future, models of multiple scattering [15] might be used to improve estimation of the attenuation coefficient, and non-linear inverse problem methods [23] used to find a solution which extracts the backscattering, scattering and attenuation coefficients from the TFDs.
The approach to particle estimation presented in the current study, Eqn. 8 and 9, can also be applied to the attenuation coefficient to estimate the composition of the media in the imaging volume. This application is described briefly here to highlight a potential basis of future studies. Assuming the contribution to the attenuation of light is the coherent sum of individual attenuation events, the spectral modification by sample attenuation can be expressed as the sum of the individual attenuation coefficients, which may be separated into contributions from scattering by particles and absorption by media. While Mie theory can be used to generate particle diameter dependent spectral absorption profiles, the diameter dependence is minor and so this is where ′ is the number density of diameter range , , is the effective scattering cross section for a single particle of diameter range , is the molar concentration of media , and , is the molar extinction coefficient for media . In matrix notation, Eqn. 11 can be expressed as = + , where is a × matrix of scattering cross sections for individual particles of different diameters, is a × 1 vector of number densities for each particle diameter, is a × matrix of absorptivity profiles for different media, and is a × 1 vector of volume concentrations for each of the defined media. The matrix may be populated from existing absorptivity data on different biological tissues and other materials, for which the literature is extensive, while may be populated using Mie theory. Eqn. 11 can be rearranged to solve for and in a similar fashion to Eqn. 9 after grouping the coefficient matrices, and , and the unknown parameter vectors, and , into new matrices: that can be solved using the matrix inversion methods described in section 2.4.

Conclusion
A S-OCT method was presented that estimates depth-dependent particle size and concentration from TFDs and Mie backscattering spectra and is most suited to occurrences where multiple scatterers in the imaging volume prohibit IS-OCT analysis of individual particles. The method was demonstrated on polystyrene microspheres in solution using a coefficient matrix containing a limited set of backscattering spectra to compensate for system errors and narrow source bandwidth.

Disclosures
The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.