Bayesian Quantification for Coherent Anti-Stokes Raman Scattering Spectroscopy.

We propose a Bayesian statistical model for analyzing coherent anti-Stokes Raman scattering (CARS) spectra. Our quantitative analysis includes statistical estimation of constituent line-shape parameters, the underlying Raman signal, the error-corrected CARS spectrum, and the measured CARS spectrum. As such, this work enables extensive uncertainty quantification in the context of CARS spectroscopy. Furthermore, we present an unsupervised method for improving spectral resolution of Raman-like spectra requiring little to no a priori information. Finally, the recently proposed wavelet prism method for correcting the experimental artifacts in CARS is enhanced by using interpolation techniques for wavelets. The method is validated using CARS spectra of adenosine mono-, di-, and triphosphate in water, as well as equimolar aqueous solutions of d-fructose, d-glucose, and their disaccharide combination sucrose.

Introduction

Coherent anti-Stokes Raman scattering (CARS) spectroscopy offers a unique microscopic tool in biophysics, biology, and materials research.[1−14] In addition to being ideally suited for qualitatively label-free microscopy,[2,3,6] the multiplex approach of CARS can also provide complete (position-dependent) vibrational spectra. In principle, this would allow a quantitative, local analysis of chemical composition.[1,11,13,15−18] However, a CARS measurement does not directly provide any quantitative information. Sophisticated analytical methods are therefore required in order to extract this information from the spectroscopic measurements.

An observed CARS spectrum arises from a coherent addition of both resonant contributions from different vibrational modes and a constant, nonresonant (NR) background contribution. This results in a complex line shape, where the positions, amplitudes, and line widths of each vibrational mode are generally hidden. This is particularly true for condensed-phase samples, where the vibrational spectra are highly congested with strongly overlapping vibrational modes.[19] At a minimum, quantitative analysis requires extracting the Raman line shapes from CARS spectra. This can be done by using a suitable phase retrieval method[19,20] on the normalized CARS spectrum. However, the technology is still limited in terms of comparable and quantitative analysis methods, which remain active and ongoing topics of research.[18,19,21,22] Moreover, the analysis is complicated by experimental errors encountered in obtaining a normalized CARS line-shape spectrum, which leads to an erroneous, nonadditive, and nonconstant background component to the NR background due to the reference CARS intensity not arising from a purely nonresonant third-order susceptibility or due to broadband laser behavior inside the sample.[21,22] If it remains uncorrected, this artifact in the NR background can prevent any quantitative information from being obtained from a CARS measurement. Recently, a procedure based on the wavelet prism decomposition algorithm was proposed to address this issue.[22]

Sequential Monte Carlo (SMC) methods have been successfully applied in a wide variety of contexts, including motion tracking,[23,24] satellite image analysis,[25] medical applications,[26] and geophysics.[27] In spectroscopy, Bayesian methods such as SMC have recently been gaining significant attention from the research community. Bayesian statistical inference has been applied to electrochemical impedance,[28] double electron–electron resonance,[29] time-resolved analysis of gamma-ray bursts,[30] and the estimation of elastic and crystallographic features by resonance ultrasound spectroscopy[31] to name a few. In particular, a hierarchical Bayesian approach combining modeling of individual line shapes with a continuous background model, with estimation done via SMC methods, has been introduced for Raman spectroscopy.[32]

The contributions of this study are 3-fold. We introduce a method for correcting experimental artifacts in raw CARS measurements, extending further the existing method based on wavelet prism decomposition.[22] Second, we propose a line-narrowing method with improved properties compared to the line shape optimized maximum entropy linear prediction (LOMEP) method.[33,34] Our method utilizes linear prediction, as in LOMEP, but in contrast circumvents the need of assuming a single a priori common line shape for all spectral lines. This constitutes a major improvement over the LOMEP method. Third and foremost, Bayesian inference is introduced to CARS spectrum analysis, extending previously available analysis methods. We formulate a Bayesian inference model that is capable of estimating predictive distributions of the underlying Raman signal, the error-corrected CARS spectrum, and the measurement CARS spectrum. This is enabled by parametric modeling of Voigt line shapes, along with a continuous, wavelet-based model for experimental artifacts.

In what follows, we introduce the Bayesian statistical model for CARS. The Raman signal of the CARS spectrum is modeled using a linear combination of Voigt line shapes. Using the Hilbert transform, we construct the modulus of the resonant part of the CARS spectrum. A nonresonant part, estimated from the data,[22] is added to obtain an error-free CARS spectrum, which is finally modulated with a slowly varying error function. Next, we describe the numerical algorithms used for statistical inference and line narrowing, along with our Bayesian prior distributions. We then present experimental details along with obtained results for the means of the constituent line shapes and the predictive intervals for the resonant Raman signal, modulating error function, error-corrected CARS spectrum, and the measurement CARS spectrum. Lastly, the key aspects of the study are briefly remarked upon, thereby concluding the paper.

Methods

Statistical Model

We model CARS spectral measurements with an additive error model given aswhere y denotes a measurement that has been discretized with spectral sampling resolution h > 0 at a wavenumber location ν = kh with , f(ν; p, θ) is the CARS spectrum model with parameter p controlling the baseline and parameters θ for the Voigt line shape, and with measurement error with known variance. For the spectrum, we use a parameter-wise separable modelwhere p is the interpolated discrete wavelet transform (DWT) detail level, εm(ν; p) is the modulating error function, and S(ν; θ) is the error-corrected CARS signal, similar to the representation used in ref (22). The signal S can further be represented aswhere the exponential part corresponds to the non-Raman part with A practically constant (see ref (22) for details), is the Hilbert transform, andwhere * denotes convolution. N stands for the number of line shapes, with each line shape having θ ≔ (a, ν, σ, γ) parameters standing for the amplitude, location, scale of the Gaussian shape, and scale of the Lorentzian shape, respectively. Thus, we have 4N parameters in total for our model of S(ν; θ).

Instead of the wavelet prism method,[22] we model the modulating error function aswhere p ∈ [1, J], and β = p – ⌊p⌋, i.e., as an interpolation between the discrete wavelet reconstruction levels D to have a continuous model for the background in contrast to the wavelet prism method where only the discrete wavelet reconstruction levels are used. With the above, we can have an unnormalized posterior formulated aswhere is the vector of observations given via eq , is the parameter vector (θ1, ..., θ) for the N Voigt peaks, represents the likelihood distribution of the forward model, and π0(p, θ) denotes prior distributions for some or all of the model parameters p and θ. As such, the total number of parameters in the model is 4N + 1. The solution of (6) is unavailable in closed form, but following ref (32), we can use Monte Carlo methods to obtain samples from this distribution, as described in the following section.

Sequential Monte Carlo

Sequential Monte Carlo (SMC) methods, also known as particle filtering and smoothing, are widely used in statistical signal processing.[35] These algorithms provide a general procedure for sampling from Bayesian posterior distributions.[36,37] SMC methods utilize a collection of weighted particles, initialized from a prior distribution, which are ultimately transformed to represent a posterior distribution under investigation. The methodology used in this study is similar to the one used in ref (22), where they use sequential likelihood tempering[37] to fit a model of surface-enhanced Raman spectra to measurements.

Assuming additive Gaussian measurement errors ϵ(ν) as in eq , the likelihood of the model f(ν; p, θ) fitting measurement data y can be formulated asand the posterior distribution for step t of the sequential likelihood tempering is given bywhere the superscript (t) denotes the iteration or “time” step of the algorithm and κ(, κ( < κ( < κ( < ··· ≤ 1 with κ(0) = 0, being a parameter controlling the degree of tempering of the likelihood, with the initial state being equal to the prior distribution while increasingly tempering the total likelihood toward the complete Bayes’ theorem. The tempering parameter κ( can be defined simply as an strictly increasing sequence so that κ( ∈ [0, 1] or as done in ref (32), the parameter can be determined adaptively according to a given learning rate η such that the relative reduction in the ESS between iterations is approximately η.

Using Q particles, with Q being the number of parameter values used to approximate the posterior distribution, individual weights of each particle at initial step t = 0 are set as equally important, so w(0) = 1/Q. The weights are then updated at each step t according toand then normalized so that ∑w( = 1. However, updating the particle weights gradually impoverishes the sample distribution. This degradation is measured by the effective sample size (ESS)To counteract this, the particles are resampled according to a chosen resampling algorithm when the ESS has fallen below a set threshold Qmin. The particle weights are then reset as w( = 1/Q. Some duplication of the particles is inevitably introduced due to the resampling procedure. To remove these duplicates, each particle is additionally updated using Markov chain Monte Carlo (MCMC) targeting the invariant distribution given by the tempered posterior defined in eq at the current iteration or “time” step t. The pseudocode of our SMC algorithm is presented in Algorithm 1.

Line Narrowing

We employ a line-narrowing method to obtain an initial estimation of peak locations ν, amplitudes a, and number of line shapes N. This is a preprocessing step for the statistical estimation method described in the previous section. A spectrum with Lorentzian line shapes can be approximately modeled aswhere Ṽ(ν, θ̃) denotes a spectrum measured at location ν with parameters θ̃ ≔ (a, ν, γ), γ is a single, constant parameter for the line width, and δ(ν – ν) is the Dirac delta function.

Our starting point is the LOMEP method,[33,34] where the constant γ approximation is used. With suitably chosen γ, we havewhere denotes the Fourier transform, t the Fourier domain variable, and xLP(t; γ, NFIR) is the linearly predicted time signal. In LOMEP, the linear prediction is done using finite impulse response filtering with filter length NFIR – 1.[33,34] The major limitation of LOMEP is the heuristic choice of γ. Additionally, the q-curve optimization method fails when the number of line shapes N increases. Despite these drawbacks, the potential of the linear prediction scheme is nevertheless attractive for its ability to substantially sharpen the line shapes when it is successful.

As an alternative to the approximation model in eq , we propose a linear combination of M similarly constructed convolutionsusing a set of width parameters γ in contrast to fixed γ. Then, the approximation of the Dirac delta functions isThe squared sum of residuals for a single convolution, denoted here by d(γ, NFIR), can be given asWe additionally define a constrained squared sum of residuals aswhere 1 = 1, if DA > 0 and 0 otherwise, and cn is a normalization constant so that the area under the spectrum is conserved:With dC(γ, NFIR), we truncate any negative parts of DA(ν, γ, NFIR) and distort the truncated spectrum according to the normalization constant cn depending on how much signal energy is present on the negative parts. By Parseval’s theorem, and by using an orthonormal wavelet basis, the energy of a signal g(ν) can be represented aswhere a and b are the scaling function and wavelet coefficients obtained using DWT. Given a signal with sharp features, the energy of the signal should be concentrated on the wavelet coefficients b and, a measure of this concentration of wavelet coefficient energy (we) can be defined asWith the above formulations, we propose Algorithm 2: Define a set of width parameters γ, for example, inferred from computational chemistry. Similarly, define an upper bound for the impulse response parameter NFIR. Then, compute DA(ν, γ, NFIR) using linear prediction for all parameter combinations of γ and NFIR and residuals d and dC along with the wavelet energy concentrations Cwe.

Using the filtering criterion fc = d + dC, narrow down the set of possible solutions by sorting them according to fc and Cwe. Take a percentage pwe of the wavelet energy sorted solutions, including the largest energy concentrations. Similarly, take a percentage p of the filtering criterion sorted solutions, including the smallest filtering criteria. Thus, an intersection of these sets should include solutions with mostly positive and sharp line shapes. Sort this intersection set of size M̃ according to d. Finally, estimate eq by choosing M so that the sum of residuals d is minimized:

As needed, smooth the obtained line-narrowed spectrum with a smoothing function.

Priors

We obtained priors by manually correcting for the experimental artifacts modeled by eq and simultaneously applying phase retrieval[19,20,38,39] and computation of the resonant imaginary component of the CARS spectrum until a reasonable Raman signal was observed. The line-narrowing algorithm was applied on the manually estimated Raman signal, producing a line-narrowed spectrum from which individual line shapes could be identified. We follow ref (32) in setting informative priors for the line shape locations ν as normal distributionswhere μν and σν2 are estimated for each line shape V(ν, θ) by numerically integrating perceived individual line shapes in the line-narrowed spectrum to estimate the means μν and variances σν2. The line-narrowing algorithm utilizes multiple Lorentzian line shapes with differing scale parameters γ, thereby giving access to an informative prior for γ. As in ref (32), we set a common prior for each γ as a log-normal distribution:where the estimates for the mean and variance, μlog(γ) and σlog(γ)2, are obtained from the parameters contained in the intersection set of size M̃. Priors for the Gaussian shape parameters σ are obtained by scaling π0(log(γ)) by . This would correspond to using identical priors for the full-width-at-half-maximum of both the Gaussian and Lorentzian line shapes. For the amplitudes, we can obtain an estimate for the areas straight-forwardly by the same numerical integration used to estimate the priors for the locations, as described above. We set a fairly wide prior for the amplitude by setting them aswhere the mean μ is the numerically integrated area of each line shape. A prior for the background parameter p is set as a uniform prior:

An estimate for the noise level σϵ2 was also obtained using the line-narrowing algorithm. The algorithm fits a smooth representation of the Raman spectrum to the manually corrected data according to eq . This smooth representation of the Raman signal is then transformed to the measurement space by eq and then by eq . The resulting residuals between the transformed smooth Raman signal and the measured CARS spectrum were used as an estimate for the noise variance σϵ2. Detailed descriptions of priors specific for each experimental data set of fructose, glucose, sucrose, and adenosine phosphate can be found in the Supporting Information.

Experimental Details

Samples

The sugar samples used in the multiplex CARS spectroscopy were equimolar aqueous solutions of d-fructose, d-glucose, and their disaccharide combination, sucrose (α-d-glucopyranosyl-(1→2)-β-d-fructofuranoside). For sample preparation, the sugar samples were dissolved in buffer solutions (50 mM HEPES, pH = 7) at equal molar concentrations of 500 mM.[16] The adenosine phosphate sample was an equimolar mixture of AMP, ADP and ATP in water for a total concentration of 500 mM.[19] The adenine ring vibrations[40] are found at identical frequencies for either for AMP, ADP, or ATP around 1350 cm–1. The phosphate vibrations between 900 and 1100 cm–1 can be used to discriminate between the different nucleotides.[15] The triphosphate group of ATP shows a strong resonance at 1123 cm–1, whereas the monophosphate resonance of AMP is found at 979 cm–1. For ADP a broadened resonance is found in between at 1100 cm–1.

Multiplex CARS Spectroscopy

All CARS spectra used to validate our methodology were recorded using a multiplex CARS spectrometer, the detailed description of which can be found elsewhere.[1,15] In brief, a 10 ps and an 80 fs mode-locked Ti:sapphire lasers were electronically synchronized and used to provide the narrowband pump/probe and broadband Stokes laser pulses in the multiplex CARS process. The center wavelengths of the pump/probe and Stokes pulses were 710 nm. The Stokes laser was tunable between 750 and 950 nm. The sugar spectra were probed within a wavenumber range from 700 to 1250 cm–1, and the AMP/ADP/ATP spectrum within a range from 900 to 1700 cm–1. The linear and parallel polarized pump/probe and Stokes beams were made collinear and focused with an achromatic lens into a tandem cuvette. The latter could be translated perpendicular to the optical axis to perform measurements in either of its two compartments, providing a multiplex CARS spectrum of the sample and of a nonresonant reference under near-identical experimental conditions. Typical average powers used at the sample were 95 mW (75 mW, in case of AMT/ADP/ATP) and 25 mW (105 mW) for the pump/probe and Stokes laser, respectively. The anti-Stokes signal was collected and collimated by a second achromatic lens in the forward-scattering geometry, spectrally filtered by short-pass and notch filters, and focused into a spectrometer equipped with a CCD camera. The acquisition time per CARS spectrum was 200 ms for sugar spectra and 800 ms for the AMP/ADP/ATP spectrum.

Computational Details

The SMC algorithm was computed using Q = 2000 particles with the resampling threshold set to Qmin = 1000 and the learning parameter set as η = 0.9. Resampling was done, as in ref (32), via residual resampling.[41] Target MCMC acceptance rate was set to 0.23 and the number of MCMC updates at each iteration was 200. An AMD Ryzen 3950X processor was used with 27 CPU threads utilized, with the SMC estimation taking 580, 522, 413, and 688 s to produce the final posterior estimate of the parameters for the fructose, glucose, sucrose, and phosphate samples, respectively. For modeling the modulating error function εm(ν; p), symlet 34 basis functions were used.

The line-narrowing algorithm was run with γ ∈ [1, 35] linearly spaced using 33 points. The maximum number of measurement points NFIR used was 150, meaning that NFIR = {1, ..., 150}. The length of the extrapolated signal[33,34] was set to equal the number of measurement points in each spectrum. The percentages pwe and p were set as 50% and 2.5% respectively. To ensure that M̃ > 0, p was incrementally increased by 2.5% until a minimum intersection set size M̃ ≥ 50 was achieved. For computation of the wavelet energy concentration Cwe symlet, eight basis functions were used.

Results and Discussion

In what follows, 95% predictive intervals for the forward model f(ν; p, θ), the modulating error function εm(ν; p), and the error-corrected spectra S(ν; θ) are presented in Figures a, 2a, 3a, and 4a for the experimental spectra of fructose, glucose, sucrose, and phosphate, respectively. Similarly, in Figures b, 2b, 3b, and 4b the 95% predictive interval for the Raman signal represented by V(ν, θ) is presented along with the means for each constituent line shape V(ν, θ). To illustrate how the priors were estimated, the manually corrected Raman signal and the result obtained via the proposed line narrowing method are shown in Figure . Additionally, the obtained posterior distributions for θ, alongside their respective prior distributions, are presented in the Supporting Information.

Figure 1

(a) Obtained 95% predictive intervals for y, f, S, and εm shown in blue, red, yellow, and purple respectively for a CARS measurement of a fructose sample. (b) Obtained 95% predictive intervals for V(ν; θ) and means of each individual line shape V(ν; θ) for the fructose sample.

Figure 2

(a) Obtained 95% predictive intervals for y, f, S, and εm shown in blue, red, yellow, and purple respectively for a CARS measurement of a glucose sample. (b) Obtained 95% predictive intervals for V(ν; θ) and means of each individual line shape V(ν; θ) for the glucose sample.

Figure 3

(a) Obtained 95% predictive intervals for y, f, S, and εm shown in blue, red, yellow, and purple respectively for a CARS measurement of a sucrose sample. Some discrepancies between y and f can be seen around the boundaries. These areas of the data should be ignored in the optimization. (b) Obtained 95% predictive intervals for V(ν; θ) and means of each individual line shape V(ν; θ) for the sucrose sample.

Figure 4

(a) Obtained 95% predictive intervals for y, f, S, and εm shown in blue, red, yellow, and purple respectively for a CARS measurement of a adenosine phosphate sample. (b) Obtained 95% predictive intervals for V(ν; θ) and means of each individual line shape V(ν; θ) for the adenosine phosphate sample.

Figure 5

Manually estimated Raman signal, according to the procedure described in the section Priors, of the fructose sample and the line-narrowed Raman signal are shown in blue and red, respectively. The perceivable individual line shapes were numerically integrated to yield informative prior estimates for Voigt line shape parameters.

The inference model proposed here was found to adequately model the CARS measurements along with perceived noise levels in the spectra. For future work, it would be interesting to include heteroscedasticity in the model instead of assuming a constant measurement error variance. Comparing the estimated predictive intervals of the obtained Raman signal showed clear correspondence to measured Raman intensities of aqueous solutions for fructose and glucose.[42] To validate the potential of the line narrowing method, we considered the number of line shapes identified for the aqueous solution of sucrose to resemble the 18 line shapes reported for solid sucrose.[43] The estimated priors were considered not to restrict the parameter posterior which can be observed in the posterior distributions when seen alongside the respective priors, especially so for the cases of fructose, sucrose, and adenosine phosphate. The obtained Raman signals for fructose, glucose, and sucrose are similar to results obtained ref (22), which further supports the applicability of the methodology presented in this study. As our method is immediately applicable to more complex samples, such as solutions with multiple solutes, applying the method to such samples provides further interesting future work.

Obtaining informative priors can be approached in different ways for chemically known samples, as was done in ref (32), where the authors use results obtained by density functional theory (DFT) software to derive estimates for the location priors and existing studies on structural properties of a known sample such as observed in refs (42 and 43). Naturally, any other forms of information on the underlying line shapes could just as well be used for the prior distributions. Here we have considered estimating the priors purely from the data using a line-narrowing algorithm, requiring minimal a priori information on the sample under study. Although this information would clearly be available in this case,[42] there are many potential applications of our method where much less is known about the molecules in question. Additionally, the use of maximum-entropy methods in improving spectral resolution can cause individual line shapes to split.[44,45] In our proposed line-narrowing method, the averaging together of multiple, resolution-enhanced spectra is postulated to possibly lessen the effect of this undesired spectral line splitting.

Conclusion

A Bayesian inference model applicable to coherent anti-Stokes Raman spectroscopy is proposed and numerically implemented. This work extends the current methodology of analyzing CARS spectra by introducing Bayesian inference in the field, enabling uncertainty quantification of spectral features. The statistical inference model is able to produce posterior distributions for physically informative parameters, line shape amplitudes, widths, and locations, for each constituent line shape along with predictive distributions for the estimated resonant Raman signal contained in the CARS measurement spectrum, the error-corrected CARS measurements, and the CARS measurement spectrum, as well as extending currently existing methodology for modeling experimental artifacts present in CARS measurements. Additionally, we have developed a line-narrowing algorithm requiring minimal a priori information on the underlying line shapes, which is readily applicable to various spectral measurements. We have successfully used this algorithm to obtain informative priors purely from the measurement data for the Bayesian inference model. The applicability of the methods is demonstrated with experimental CARS spectra of sucrose, fructose, glucose, and adenosine phosphate.