Bayesian quantification for coherent anti-Stokes Raman scattering spectroscopy

We model CARS spectral measurements with an additive error model given as

where $y_{k}$ denotes a measurement that has been discretized with spectral sampling resolution $h>0$ at a wavenumber location $\nu_{k}=kh$ with $k\in\mathbb{Z}_{+}$, $f(\nu_{k};p,\bm{\theta})$ is the CARS spectrum model with parameter $p$ controlling the baseline and parameters $\bm{\theta}$ for the Voigt line shape, and with measurement error $\epsilon(\nu_{k})\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$ with known variance. For the spectrum, we use a parameter-wise separable model

where $p$ is the interpolated discrete wavelet transform (DWT) detail level, $\varepsilon_{\rm{m}}(\nu;p)$ is the modulating error function, and $S(\nu;\bm{\theta})$ is the error-corrected CARS signal, similar to the representation used in Ref. 22. The signal $S$ can further be represented as

where the exponential part corresponds to the non-Raman part with $A_{J}$ practically constant (see Ref. 22 for details), $\mathcal{H}$ is the Hilbert transform, and

where $\ast$ denotes convolution. $N$ stands for the number of line shapes, with each line shape having $\bm{\theta}_{n}:=(a_{n},\nu_{n},\sigma_{n},\gamma_{n})^{T}$ 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(\nu;\bm{\theta})$.

Instead of the wavelet prism method Kan et al. 2016, we model the modulating error function as

where $p\in[1,J]$, and $\beta=p-\lfloor p\rfloor$, i.e., as an interpolation between the discrete wavelet reconstruction levels $D_{j}$ to have a continuous model for the background. With the above, we can have an unnormalized posterior formulated as

where $\mathbf{y}:=(y_{1},\dots,y_{K})^{T}\in\mathbb{R}^{K}$ is the vector of observations given via eq (1), $\bm{\theta}\in\mathbb{R}^{4N}_{+}$ is the parameter vector $(\bm{\theta}_{1},\dots,\bm{\theta}_{N})^{T}$ for the $N$ Voigt peaks, $\mathcal{L}(\mathbf{y}\mid p,\bm{\theta})$ represents the likelihood distribution of the forward model, and $\pi_{0}(p,\bm{\theta})$ denotes prior distributions for some or all of the model parameters $p$ and $\bm{\theta}$. 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 (SMC) methods, also known as particle filtering and smoothing, are widely used in statistical signal processing Särkkä 2013. These algorithms provide a general procedure for sampling from Bayesian posterior distributions Chopin 2002; Del Moral et al. 2006. 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 Del Moral et al. 2006 to fit a model of surface-enhanced Raman spectra to measurements.

Assuming additive Gaussian measurement errors $\epsilon(\nu_{k})$, likelihood of the model $f(\nu_{k};p,\bm{\theta})$ fitting measurement data $y_{k}$ can be formulated as

and the posterior distribution for step $t$ of the sequential likelihood tempering is given by

where the superscript $(t)$ denotes the iteration or “time” step of the algorithm and $\kappa^{(t)}$, $\kappa^{(t-1)}<\kappa^{(t)}<\kappa^{(t+1)}<\dots\leq 1$ with $\kappa^{(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 towards the complete Bayes’ theorem. The tempering parameter $\kappa^{(t)}$ can be defined simply as an strictly increasing sequence so that $\kappa^{(t)}\in[0,1]$ or as done in Ref. 32, the parameter can be determined adaptively according to a given learning rate $\eta$ such that the relative reduction in the ESS between iterations is approximately $\eta$.

Using $Q$ particles, individual weights of each particle at initial step $t=0$ are set as equally important $w_{q}^{(0)}=\frac{1}{Q}$. The weights are then updated, and normalized, at each step $t$ according to

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 when the ESS has fell below a set threshold $Q_{\text{min}}$ according to a chosen resampling algorithm and the particle weights are reset as $w_{q}^{(t)}=\frac{1}{Q}$. Due to resampling, duplicates of the particles are obviously generated. This is undesirable and as such, each particle is additionally updated using Markov chain Monte Carlo (MCMC) with the target distribution given by the tempered posterior distribution defined in eq (8) at the current iteration or “time” step $t$. A pseudo-code of the SMC algorithm used is presented in Algorithm 1.

In line narrowing one aims at spectral line sharpening. We shall use line narrowing for making rough initial estimation of peak locations $\nu_{n}$, amplitudes $a_{n}$, and number of line shapes $N$. This is a preprocessing step for the statistical estimation method described in the previous section. Let us model and approximate a spectrum with Lorentzian line shapes as

where $\widetilde{V}_{N}(\nu_{k},\widetilde{\bm{\theta}})$ denotes a measured spectrum, $\nu_{k}\in\mathbb{R}^{K}$ measurement locations, $\widetilde{\bm{\theta}}:=\linebreak(a_{n},\nu_{n},\gamma_{n})^{T}$, $\gamma$ a single constant width parameter, and $\delta(\nu-\nu_{n})$ is the Dirac delta function.

Our starting point is the LOMEP method Kauppinen et al. 1981; Kauppinen et al. 1991, where the constant $\gamma$ approximation is used. With suitably chosen $\gamma$, we have

where $\mathcal{F}$ denotes the Fourier transform, $t_{k}$ the Fourier domain variable, and $x_{\textrm{LP}}(t_{k};\gamma_{m},N_{\textrm{FIR}})$ is the linearly predicted time signal. In LOMEP, the linear prediction is done finite impulse response filtering with filter length $N_{\textrm{FIR}}-1$Kauppinen et al. 1981; Kauppinen et al. 1991. The major limitation of LOMEP is the heuristic choice of $\gamma$, and additionally, the $q$-curve optimization method fails when number of line shapes $N$ increases. Despite the drawbacks, the potential of the linear prediction scheme is nevertheless attractive for its ability to heavily narrow down line shapes when successful.

We propose an alternative approximative model in eq (11) as a linear combination of $M$ similarly constructed convolutions

using a set of width parameters $\gamma_{m}$ in contrast to fixed $\gamma$. Then, the approximation of the Dirac delta functions is

The squared sum of residuals for a single convolution, denoted here by $d(\gamma_{m},N_{\textrm{FIR}})$, can be given as

We additionally define a constrained squared sum of residuals as

where $\mathbf{1}_{D_{\textrm{A}}>0}=1$, if $D_{\textrm{A}}>0$ and $0$ otherwise, and $c_{\textrm{n}}$ is a normalization constant so that the area under the spectrum is conserved:

With $d_{\textrm{C}}(\gamma_{m},N_{\textrm{FIR}})$ we truncate any negative parts of $D_{\textrm{A}}(\nu_{k},\gamma_{m},N_{\textrm{FIR}})$ and distort the truncated spectrum according to the normalization constant $c_{\textrm{n}}$ the more 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(\nu)$ can be represented as

where $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_{j}$ and, a measure of this concentration of wavelet coefficient energy (we) can be defined as

With the above formulations, we propose Algorithm 2: Define a set of width parameters $\gamma_{m}$, for example, inferred from computational chemistry. Similarly, define an upper bound for the impulse response parameter $N_{\textrm{FIR}}$. Then, compute $D_{\textrm{A}}(\nu_{k},\gamma_{m},N_{\textrm{FIR}})$ using linear prediction for all parameter combinations of $\gamma_{m}$ and $N_{\textrm{FIR}}$ and residuals $d$ and $d_{\textrm{C}}$ along with the wavelet energy concentrations $C_{\textrm{we}}$.

Then using the filtering criterion $f_{\textrm{c}}=d+d_{\textrm{C}}$, narrow down the set of possible solutions by sorting them according to $f_{\textrm{c}}$ and $C_{\textrm{we}}$. Take a percentage $p_{\textrm{we}}$ of the wavelet energy sorted solutions, including the largest energy concentrations. Similarly, take a percentage $p_{f_{\textrm{c}}}$ 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 $\widetilde{M}$ according to $d$. Finally, estimate eq (13) by choosing $M$ so that the sum of residuals $d_{M}$ is minimized:

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

We obtained priors by manually correcting for the experimental artefacts modelled by eq (5) and simultaneously applying phase retrieval Vartiainen 1992; Vartiainen et al. 2006; Liu et al. 2009; Cicerone et al. 2012 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 $\nu_{k}$ as normal distributions

where $\mu_{\nu_{n}}$ and $\sigma_{\nu_{n}}^{2}$ are estimated for each line shape $V(\nu,\bm{\theta}_{n})$ by numerically integrating perceived individual line shapes in the line narrowed spectrum to estimate the means $\mu_{\nu_{n}}$ and variances $\sigma_{\nu_{n}}^{2}$. The line narrowing algorithm utilizes multiple Lorentzian line shapes with differing scale parameters $\gamma_{m}$, thereby giving access to an informative prior for $\gamma_{n}$. As in Ref. 32, we set a prior common for each $\gamma_{n}$ as a log-normal distribution:

where the estimates for the mean and variance, $\mu_{\log(\gamma)}$ and $\sigma_{\log(\gamma)}^{2}$, are obtained from the parameters contained in the intersection set of size $\widetilde{M}$. Priors for the Gaussian shape parameters $\sigma_{n}$ are obtained by scaling $\pi_{0}\big(\log(\gamma_n)\big)$ by $\sqrt{2\log(2)}$. This would correspond to using identical priors for the full-width at half maximums for 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 as

where the mean $\mu_{a_{n}}$ 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 $\sigma_{\epsilon}^{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 (20). This smooth representation of the Raman signal is then transformed to the measurement space by eq (3) and then by eq (2). The resulting residuals between the transformed smooth Raman signal and the measured CARS spectrum were used as an estimate for the noise variance $\sigma_{\epsilon}^{2}$. Detailed descriptions of priors specific for each experimental data set of fructose, glucose, sucrose, and adenosine phosphate can be found in the supplementary material.

The sugar samples used in the multiplex CARS spectroscopy were equimolar aqueous solutions of D-fructose, D-glucose, and their disaccharide combination sucrose. For sample preparation D-fructose, D-glucose, and their disaccharide combination, sucrose ($\alpha$-D-glucopyranosyl-(1$\to$2)-$\beta$-D-fructofuranoside) were dissolved in buffer solutions (50 mM HEPES, pH=7) at equal molar concentrations of 500 mM Müller et al. 2007. The adenosine phosphate sample was an equimolar mixture of AMP, ADP and ATP in water for a total concentration of 500 mM Vartiainen et al. 2006. The adenine ring vibrations Mathlouthi and Luu 1980 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 Rinia et al. 2006. The tri-phosphate 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.

All CARS spectra were recorded using a multiplex CARS spectrometer, the detailed description of which can be found elsewhere Müller and Schins 2002; Rinia et al. 2006. 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 non-resonant 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.

The SMC algorithm was computed using $Q=2000$ particles with the resampling threshold set to $Q_{\text{min}}=1000$ and the learning parameter set as $\eta=0.9$. Resampling was done, as in Ref. 32, via residual resampling Douc and Cappe 2005. 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 seconds to produce the final posterior estimate of the parameters for the fructose, glucose, sucrose, and phosphate samples respectively. For modelling the modulating error function $\varepsilon_{\rm{m}}(\nu;p)$ symlet 34 basis functions were used.

The line narrowing algorithm was run with $\gamma_{m}\in[1,35]$ linearly spaced using 33 points. The maximum number of measurement points $N_{\textrm{FIR}}$ used was 150, meaning that $N_{\textrm{FIR}}=\{1,\dots,150\}$. The length of the extrapolated signal Kauppinen et al. 1981; Kauppinen et al. 1991 was set to equal the number of measurement points in each spectrum. The percentages $p_{\textrm{we}}$ and $p_{f_{\textrm{c}}}$ were set as 50% and 2.5% respectively. To ensure that $\widetilde{M}>0$, $p_{f_{\textrm{c}}}$ was incrementally increased by 2.5% until a minimum intersection set size $\widetilde{M}\geq 50$ was achieved. For computation of the wavelet energy concentration $C_{\textrm{we}}$ symlet 8 basis functions were used.