A Bayesian approach to spectral quantitative photoacoustic tomography

In the photoacoustic effect, the absorption of a short light pulse in a target generates an initial acoustic pressure distribution that is proportional to the absorbed energy density [1]. The connection between the absorbed optical energy density and the initial pressure distribution, the photoacoustic efficiency, is often called the Grüneisen parameter, terminology which is adopted in this paper. Due to the elastic nature of tissue, the initial pressure distribution subsequently propagates as an ultrasonic wave, which can be measured on the boundaries of the target. In photoacoustic imaging, these acoustic measurements are used to deduce properties of the target [2, 3]. As the propagation of light is dependent on the optical properties of the tissue, such as the absorption and the scattering, the ultrasonic wave carries with it information about these parameters, which can in turn provide valuable morphological, pathological and functional data about the underlying tissue. In quantitative photoacoustic tomography (QPAT) the goal is to determine these properties.

The inverse problem associated with QPAT is two-fold. Ultrasonic measurements from the boundary of the target are first used to determine the initial pressure distribution. This is the acoustical inverse problem associated with photoacoustic imaging, which has been studied extensively [4–22]. Secondly, the initial pressure distribution, which is related to the optical properties of the target, is used as the data in an optical inverse problem. In the optical inverse problem, the distribution of the optical parameters is estimated from the initial pressure distribution and knowledge of the illumination scheme. In the optical inverse problem of QPAT two models of light propagation have been used: the radiative transfer equation (RTE) [23–29] and its diffusion approximation (DA) [26, 30–37].

In spectral quantitative photoacoustic tomography (SQPAT), measurements of the target are obtained with illuminations at different optical wavelengths [38–44]. After obtaining the initial pressure distribution, via the acoustical inversion, for each wavelength, spectral parameter models (models of the wavelength dependence of optical absorption and scattering parameters, which couple the distributions at the different wavelengths) can be used to aid the optical inversion.

There are two approaches to using the spectral parameter models. One is to reconstruct the optical absorption and scattering at each different wavelength separately, and then fit the spectral parameter models to these optical parameters [40, 43, 44]. The second approach is to express the optical properties of the forward problem as a function of the spectral parameter model, and estimate the spectral parameters directly [41–44]. The latter approach is used in this work, and is compared to the former. In addition to being able to provide information about the optical properties of the target, SQPAT makes it possible to reconstruct physiological quantities, such as chromophore concentrations or blood oxygenation in biological tissues.

The optical inverse problem in SQPAT has previously been studied in [41–44]. Common assumptions have been to assume that the Grüneisen parameter is a known constant [41–43], or the scattering (or diffusion) to be proportional to a known constant power of the illumination wavelength [42–44]. In this work a Bayesian approach to ill-posed inverse problem is taken [45]. The Grüneisen and scattering parameters are estimated simultaneously with the absorption parameters utilizing quantitative prior information provided by the Bayesian approach. The forward model used here is based on the DA, though the approach could be implemented with the RTE. The spectral parameter model for the absorption is expressed as a spatially dependent chromophore concentration weighted sum of known absorption spectra of different chromophores [41–44]. The spectral model for scattering is treated in a more general fashion to previous studies, as the scattering is expressed as a function of a proportionality factor and a power law exponent, both spatially dependent. The Grüneisen parameter is modeled as spatially dependent, but independent of the illuminating wavelength and independent of the concentrations of the chromophores present. The parameters that are reconstructed are the chromophore concentrations, the exponent and proportionality factor of the scattering parameter model, and the Grüneisen parameter. The estimates obtained by this direct approach are compared to estimates computed by fitting the spectral parameter model to estimates of the absorption and scattering parameters at multiple wavelengths.

The paper is organized as follows. The forward problem associated with SQPAT is described in section 2, along with the spectral model of tissue properties. In section 3, the Bayesian inversion approach in SQPAT is presented. The derived inverse scheme is numerically validated in section 4, and a discussion is provided in section 5.

The measurement obtained in QPAT of a target $\Omega \subset {{\mathbb{R}}^{k}}$ (k = 2 or 3) is the ultrasonic (photoacoustic) pulse as a function of time measured at the boundary $\partial \Omega$. After performing the acoustical inversion, the initial pressure distribution ${{p}_{0}}$ inside $\Omega$ is obtained, which acts as the data in the optical inverse problem. The initial pressure distribution is connected to the absorbed optical energy density through

$$\begin{eqnarray}&&{{p}_{0}}\left( r;\gamma ,{{\mu }_{\text{a}}},\mu _{\text{s}}^{\prime },s \right)=\gamma {{\mu }_{\text{a}}}\Phi (r;{{\mu }_{\text{a}}},\mu _{\text{s}}^{\prime },s),\end{eqnarray} \tag{ 1 }$$

where $r\in \Omega$ is the position, γ is the Grüneisen parameter, ${{\mu }_{a}}\Phi$ is the absorbed optical energy density, ${{\mu }_{\text{a}}}$ and $\mu _{\text{s}}^{\prime }$ are the absorption and (reduced) scattering coefficients, $\Phi$ is the fluence inside the target, and s is the illumination inducing the photoacoustic effect. All parameters are functions of position r and wavelength λ of the illumination, with the exception of γ, being usually treated as wavelength independent.

For highly scattering media the fluence $\Phi$ in (1) is governed by the DA [46, 47]

$$\begin{eqnarray}&&\left\{ \begin{array}{} -\nabla \cdot \kappa \nabla \Phi +{{\mu }_{\text{a}}}\Phi =0,\qquad r\in \Omega \\ {{\zeta }_{k}}\Phi +\frac{1}{2}A\kappa \nabla \Phi \cdot \nu =s,\qquad r\in \partial \Omega , \\ \end{array} \right.\end{eqnarray} \tag{ 2 }$$

where $\kappa ={{\left( k\left( {{\mu }_{\text{a}}}+\mu _{\text{s}}^{\prime } \right) \right)}^{-1}}$ is the optical diffusion, ${{\zeta }_{k}}$ is a parameter related to the number of spatial dimensions (${{\zeta }_{2}}={{\pi }^{-1}}$, ${{\zeta }_{3}}={{4}^{-1}}$), A describes the reflectivity at the boundary $\partial \Omega$ and ν is the outward normal on $\partial \Omega$. Throughout this work value A = 1 is used. s in (2) is the spatially dependent illumination pattern (inward light current) on the boundary $\partial \Omega$.

The optical properties of the target can be expressed using their spectral representation. Commonly used spectral models for the optical properties can be written as [39–44]

$$\begin{eqnarray}\left\{ \begin{array}{} {{\mu }_{\text{a}}}\left( r;{{c}_{1}},\ldots ,{{c}_{N}},\lambda \right) & = & \underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,{{c}_{n}}\left( r \right){{\mu }_{\text{a},n}}\left( \lambda \right) \\ \mu _{\text{s}}^{\prime }\left( r;\mu _{\text{s},\text{REF}}^{\prime },b,\lambda \right) & = & \mu _{\text{s},\text{REF}}^{\prime }\left( r \right){{\left( \frac{\lambda }{{{\lambda }_{\text{REF}}}} \right)}^{-b\left( r \right)}}, \\ \end{array} \right.\end{eqnarray} \tag{ 3 }$$

where the absorption coefficient ${{\mu }_{\text{a}}}$ is expressed as chromophore concentration ${{c}_{n}}$ weighted sum of N known chromophore absorption spectra ${{\mu }_{\text{a},n}}$. The scattering coefficient $\mu _{\text{s}}^{\prime }$ is expressed, according to Mie scattering theory, as being proportional to $\mu _{\text{s},\text{REF}}^{\prime }\left( r \right)=\mu _{\text{s}}^{\prime }\left( r;{{\lambda }_{\text{REF}}} \right)$ and power $-b\left( r \right)$ of a relative wavelength $\lambda /{{\lambda }_{\text{REF}}}$ [48–51]. $\mu _{\text{s},\text{REF}}^{\prime }$ is the reduced scattering coefficient at a reference wavelength ${{\lambda }_{\text{REF}}}$.

Substituting the spectral parameter model (3) into (1) results in the SQPAT forward model. In this work the SQPAT forward model is discretized using the finite element method. The model is discretized using linear basis functions ${{\phi }_{q}}$ defined at points ${{r}_{q}}\in \Omega$ into $q=1,\ldots ,Q$ discretization points. The parameters are expressed as

$$\begin{eqnarray}&&h\left( r \right)=\underset{q=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,h\left( {{r}_{q}} \right){{\phi }_{q}}\left( r \right),\qquad h={{p}_{0}},\;{{c}_{1}},\ldots ,{{c}_{N}},\;\gamma ,\;\mu _{\text{s},\text{REF}}^{\prime }\;\text{or}\;b,\end{eqnarray} \tag{ 4 }$$

where $h\left( {{r}_{q}} \right)$ is the discrete representation of the parameter. The resulting SQPAT forward model for measurements $m=1,\ldots ,M$ performed with illuminations ${{s}_{m}}$ at wavelengths ${{\lambda }_{m}}$ can be expressed as

$$\begin{eqnarray}&&{{z}_{m}}={{f}_{m}}\left( x \right),\qquad m=1,\;\ldots ,\;M,\end{eqnarray} \tag{ 5 }$$

where ${{z}_{m}}={{\left( {{p}_{0,m}}\left( {{r}_{1}} \right),\ldots ,{{p}_{0,m}}\left( {{r}_{Q}} \right) \right)}^{\top }}\in {{\mathbb{R}}^{Q}}$ is a discrete representation of ${{p}_{0}}$ associated with measurement $m,x=({{c}_{1}}\left( {{r}_{1}} \right),\ldots ,{{c}_{1}}\left( {{r}_{Q}} \right),\ldots ,{{c}_{N}}\left( {{r}_{1}} \right),\ldots ,{{c}_{N}}\left( {{r}_{Q}} \right),\gamma \left( {{r}_{1}} \right),\ldots ,\gamma \left( {{r}_{Q}} \right)$, $\mu _{\text{s},\text{REF}}^{\prime }\left( {{r}_{1}} \right),$ $\ldots ,\mu _{\text{s},\text{REF}}^{\prime }\left( {{r}_{Q}} \right)$, $b\left( {{r}_{1}} \right),\ldots$, $b\left( {{r}_{Q}} \right){{)}^{\top }}\in {{\mathbb{R}}^{Q\left( N+3 \right)}}$ are the parameters of the target, and ${{f}_{m}}:{{\mathbb{R}}^{Q\left( N+3 \right)}}\to {{\mathbb{R}}^{Q}}$ is the discretized model of (1)–(3). Forward model (5) can be further simplified to

$$\begin{eqnarray}&&z=f\left( x \right),\end{eqnarray} \tag{ 6 }$$

where $z={{\left( z_{1}^{\top },\ldots ,z_{M}^{\top } \right)}^{\top }}\in {{\mathbb{R}}^{MQ}}$, and $f\left( x \right)={{\left( {{f}_{1}}{{\left( x \right)}^{\top }},\ldots ,{{f}_{M}}{{\left( x \right)}^{\top }} \right)}^{\top }}:{{\mathbb{R}}^{Q\left( N+3 \right)}}\to {{\mathbb{R}}^{MQ}}$. Details of the finite element approximation of the forward model are presented in appendix A.

The discrete measurement data z in (6) is polluted with some additive noise

$$\begin{eqnarray}&&z=f\left( x \right)+e,\end{eqnarray} \tag{ 7 }$$

with e being the noise vector of each measurement m such that $e={{\left( {{e}_{1}},\ldots ,{{e}_{M}} \right)}^{\top }}\in {{\mathbb{R}}^{MQ}}$ and ${{e}_{m}}\in {{\mathbb{R}}^{Q}}$. In the case of SQPAT, the noise e is the residual noise remaining after the application of the acoustical inversion to the noisy acoustical measurement on $\partial \Omega$ [37].

In the Bayesian approach, the inverse problem is treated as a problem of statistical inference [45]. All variables are modeled as random variables and the measurements are used to determine the posterior density of the parameters of primary interest. Assuming that z and x are random variables, the likelihood distribution associated with the observation model (7) is

$$\begin{eqnarray}&&\pi \left( z\;|\;x,e \right)=\delta \left( z-f\left( x \right)-e \right),\end{eqnarray} \tag{ 8 }$$

where $\delta \left( \cdot \right)$ is the Dirac δ-distribution. The solution to the inverse problem is the posteriori distribution $\pi \left( x\;|\;z \right)$, which can be obtained using Bayes' theorem

$$\begin{eqnarray} \pi \left( x\;|\;z \right) & = & \pi {{\left( z \right)}^{-1}}{{\mathop{\int }^{}}_{{{\mathbb{R}}^{MQ}}}}\;\pi \left( z,x,e \right)\;\text{d}e \\ {} & = & \pi {{\left( z \right)}^{-1}}{{\mathop{\int }^{}}_{{{\mathbb{R}}^{MQ}}}}\pi \left( z\;|\;x,e \right)\;\pi \left( x \right)\;\pi \left( e \right)\;\text{d}e \\ {} & = & \pi {{\left( z \right)}^{-1}}{{\mathop{\int }^{}}_{{{\mathbb{R}}^{MQ}}}}\delta \left( z-f\left( x \right)-e \right)\;\pi \left( x \right)\;\pi \left( e \right)\;\text{d}e \\ {} & = & \pi {{\left( z \right)}^{-1}}\;{{\pi }_{e}}\left( z-f\left( x \right) \right)\;\pi \left( x \right), \\ \end{eqnarray} \tag{ 9 }$$

where ${{\pi }_{e}}\left( \cdot \right)$ is used to denote the distribution of e, $\pi \left( x \right)$ is the prior distribution of x, and it has been assumed that e and x are uncorrelated. Since, for fixed measurement z, $\pi \left( z \right)$ is a fixed scalar, (9) can be expressed as a proportionality

$$\begin{eqnarray}&&\pi \left( x\;|\;z \right)\propto {{\pi }_{e}}\left( z-f\left( x \right) \right)\;\pi \left( x \right).\end{eqnarray} \tag{ 10 }$$

For normal distributions $e\sim \mathcal{N}\left( {{\eta }_{e}},{{\Gamma }_{e}} \right)$, $x\sim \mathcal{N}\left( {{\eta }_{x}},{{\Gamma }_{x}} \right)$ the posteriori distribution (10) can be written as

$$\begin{eqnarray}&&\pi \left( x\;|\;z \right)\propto \exp \left( -\frac{1}{2}||{{L}_{e}}\left( z-f\left( x \right)-{{\eta }_{e}} \right)|{{|}^{2}}-\frac{1}{2}||{{L}_{x}}\left( x-{{\eta }_{x}} \right)|{{|}^{2}} \right),\end{eqnarray} \tag{ 11 }$$

where $L_{e}^{\top }{{L}_{e}}=\Gamma _{e}^{-1}$ and $L_{x}^{\top }{{L}_{x}}=\Gamma _{x}^{-1}$ are the Cholesky decompositions of the inverse covariance matrices, and $||\cdot ||$ is the Euclidean norm.

A pointwise estimate of x is the maximum a posteriori (MAP) estimate, which can be defined for (11) as

$$\begin{eqnarray}&&{{x}_{\text{MAP}}}=\underset{x}{\mathop{\arg \;\min }}\,\;||{{L}_{e}}\left( z-f\left( x \right)-{{\eta }_{e}} \right)|{{|}^{2}}+||{{L}_{x}}\left( x-{{\eta }_{x}} \right)|{{|}^{2}}.\end{eqnarray} \tag{ 12 }$$

The MAP estimate is used in this work due to its computational efficiency.

Priors for the parameters ${{c}_{1}},\ldots ,{{c}_{N}}$, γ, $\mu _{\text{s},\text{REF}}^{\prime }$ and b are based on the Ornstein–Uhlenbeck process [52]. The covariance matrices for the estimated parameters are defined in relation to covariance matrix $\Xi \in {{\mathbb{R}}^{Q\times Q}}$, with sth row and pth column defined with the Ornstein–Uhlenbeck covariance function as

$$\begin{eqnarray}&&{\Xi}_{sp}=\exp \left( -||{{r}_{s}}-{{r}_{p}}||\;/\;l \right),\end{eqnarray} \tag{ 13 }$$

where l describes the corralation length (i.e. the smoothness) of the prior. The specific choice of the prior mean and the scaling of the prior covariance matrix is described later in section 4.2, where the combined prior of the estimated parameters and its relation to ${\Xi}$ is defined. Other possible choices for the prior could be, for example, a white noise prior with diagonal covariance matrix or a squared exponential prior with the matrix elements defined as $\exp \left( -\frac{1}{2}||{{r}_{s}}-{{r}_{p}}|{{|}^{2}}\;/\;{{l}^{2}} \right)$. A white noise prior is not a smooth prior, and is more suited for estimation of parameters which are independent or which have no spatial correlation. A squared exponential prior, on the other hand, is infinitely differentiable. The choice of the prior distribution, and the parameters of the distribution, will affect the MAP estimate. White noise prior supports nonsmooth or very noise sensitive MAP estimates whereas the squared exponential prior supports estimates which can be unrealistically smooth considering the true parameters. Both of these prior choices may be too harsh a presumption for applications where it is possible that the target is composed of heterogeneities separated by sharp edges (an example of such could be blood vessels in photoacoustic biomedical imaging). Ornstein–Uhlenbeck provides a compromise between the two: while the prior is not differentiable (and as such not a smooth prior), it still has finite cross-correlation between spatially different points promoting distributions which can be close to homogeneous locally with sharp changes between different areas.

In this work the minimization problem (12) is solved using Gauss–Newton method augmented with an approximative line search algorithm. The minimization algorithm is started from the prior mean. Details of the computation of the partial derivatives of the forward model required by the Gauss–Newton algorithm are presented in appendix B. The line search is performed on fixed interval $\left[ 0,{{\alpha }_{\text{MAX}}}/2 \right]$. ${{\alpha }_{\text{MAX}}}$ was chosen as the largest $\alpha \prime$ such that ${{\mu }_{\text{a}}}\left( r;{{c}_{1}},\ldots ,{{c}_{N}},\lambda \right)+\mu _{\text{s}}^{\prime }\left( r;\mu _{\text{s},\text{REF}}^{\prime },b,\lambda \right)>0$ for $\forall \alpha \in [0,\alpha \prime ],\forall r\in \Omega ,$ $\forall \lambda \in \left\{ {{\lambda }_{1}},\ldots ,{{\lambda }_{M}} \right\}$ with ${{\mu }_{\text{a}}}$ and $\mu _{\text{s}}^{\prime }$ corresponding to optical parameters at $x+\alpha {{\Delta }_{\text{GN}}}$, where ${{\Delta }_{\text{GN}}}$ is the Gauss–Newton update direction. This choice of step length avoids issues that could arise with the singularity of the optical diffusion κ in (2) when ${{\mu }_{\text{a}}}+\mu _{\text{s}}^{\prime }$ approaches zero.

4.1. Simulation of the measurement data

The data was simulated using the DA in a two dimensional square domain $\Omega$ with sidelength of $10\;\text{mm}$. The data $z={{\left( z_{1}^{\top },\ldots ,z_{M}^{\top } \right)}^{\top }}$ was formed based on six measurements (M = 6) in total, such that two different illumination patterns were used at three different wavelengths. The wavelengths were ${{\lambda }_{1}}={{\lambda }_{2}}=700\;\text{nm}$, ${{\lambda }_{3}}={{\lambda }_{4}}=800\;\text{nm}$, ${{\lambda }_{5}}={{\lambda }_{6}}=900\;\text{nm}$. The illumination patterns were defined as

$$\begin{eqnarray} {{s}_{1}}\left( r \right)={{s}_{3}}\left( r \right)={{s}_{5}}\left( r \right)=\left\{ \begin{array}{} 1, & r\in \partial {{\Omega }_{\text{L}}}\partial {{\Omega }_{\text{R}}} \\ 0, & r\in \partial {{\Omega }_{\text{T}}}\partial {{\Omega }_{\text{B}}} \\ \end{array} \right. \\ {{s}_{2}}\left( r \right)={{s}_{4}}\left( r \right)={{s}_{6}}\left( r \right)=\left\{ \begin{array}{} 0, & r\in \partial {{\Omega }_{\text{L}}}\partial {{\Omega }_{\text{R}}} \\ 1, & r\in \partial {{\Omega }_{\text{T}}}\partial {{\Omega }_{\text{B}}} \\ \end{array} \right. \\ &&,\end{eqnarray} \tag{ 14 }$$

where $\partial {{\Omega }_{\text{L}}}$, $\partial {{\Omega }_{\text{R}}}$, $\partial {{\Omega }_{\text{T}}}$, $\partial {{\Omega }_{\text{B}}}$ are the left, right, top and bottom segments, respectively, forming the square boundary such that $\partial \Omega =\partial {{\Omega }_{\text{L}}}\partial {{\Omega }_{\text{R}}}\partial {{\Omega }_{\text{T}}}\partial {{\Omega }_{\text{B}}}$. The optical absorption coefficient was composed of three different chromophores (N = 3): fat, deoxygenated and oxygenated blood corresponding to chromophore concentrations ${{c}_{1}}$, ${{c}_{2}}$ and ${{c}_{3}}$ and known absorption spectra ${{\mu }_{\text{a},1}}$, ${{\mu }_{\text{a},2}}$, ${{\mu }_{\text{a},3}}$ respectively. Values for the absorption spectra are shown in table 1. Two problems were considered: a case with smoothly varying spectral parameters and a case with sharply varying spectral parameters. The values for the spectral parameters ${{c}_{1}}$, ${{c}_{2}}$, ${{c}_{3}}$, γ, $\mu _{\text{s},\text{REF}}^{\prime }$ and b used to simulate the data in the two cases are presented in the first columns of figures 1 and 2 for the smoothly and sharply varying parameter cases respectively. $\mu _{\text{s},\text{REF}}^{\prime }$ was defined at a reference wavelength ${{\lambda }_{\text{REF}}}=700\;\text{nm}$. The spectral parameters correspond to ${{\mu }_{\text{a}}}$ varying within range 0.10–0.92 $\text{m}{{\text{m}}^{-1}}$ and $\mu _{\text{s}}^{\prime }$ varying within range $0.36\text{--}0.90\;\text{m}{{\text{m}}^{-1}}$ at the investigated wavelengths for both of the cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Absorption spectra used in this paper: ${{\mu }_{\text{a},1}}$ corresponds to fat [53], ${{\mu }_{\text{a},2}}$ and ${{\mu }_{\text{a},3}}$ to oxygenated and deoxygenated blood [54].

$\lambda$ (nm)	${{\mu }_{\text{a},1}}$ ($\lambda$) $(\text{m}{{\text{m}}^{-1}})$	${{\mu }_{\text{a},2}}\;(\lambda )\;(\text{m}{{\text{m}}^{-1}})$	${{\mu }_{\text{a},3}}(\lambda )\;(\text{m}{{\text{m}}^{-1}})$
700	0.0700	0.9781	0.1713
800	0.0750	0.4496	0.4632
900	0.0800	0.4754	0.7155

The computational domain $\Omega$ was composed of 9016 triangular elements and 4629 nodes. After simulating the noiseless data $z={{\left( {{z}_{1}},\ldots ,{{z}_{M}} \right)}^{\top }}$, it was interpolated onto an inversion grid composed of 5022 triangular elements and Q = 2598 nodes. The interpolation was performed in order to avoid the inverse crime. Additive error with a Gaussian white noise distribution was added to the noiseless data such that ${{z}_{m}}\to {{z}_{m}}+{{e}_{m}}$. It was drawn from a distribution with ${{e}_{m}}\sim \mathcal{N}\left( 0,{{(0.01/3)}^{2}}\cdot \Delta _{m}^{2}\cdot I \right)$, with ${{\Delta }_{m}}=\max \{{{z}_{m}}\}-\min \{{{z}_{m}}\}$, $0\in {{\mathbb{R}}^{Q}}$ being a vector of zeros and $I\in {{\mathbb{R}}^{Q\times Q}}$ being the identity matrix.

4.2. Reconstructions

The MAP estimates were obtained by minimizing (12) as described in section 3. The parameters of the combined prior $x\sim \mathcal{N}\left( {{\eta }_{x}},{{\Gamma }_{x}} \right)$ were

$$\begin{eqnarray} {{\eta }_{x}} & = & {{\left( {{\eta }_{c}}\cdot {{1}^{\top }},{{\eta }_{c}}\cdot {{1}^{\top }},{{\eta }_{c}}\cdot {{1}^{\top }},{{\eta }_{\gamma }}\cdot {{1}^{\top }},{{\eta }_{\mu _{\text{s},\text{REF}}^{\prime }}}\cdot {{1}^{\top }},{{\eta }_{b}}\cdot {{1}^{\top }} \right)}^{\top }}, \\ {{\Gamma }_{x}} & = & \text{diag}\left\{ {{\Gamma }_{c}},{{\Gamma }_{c}},{{\Gamma }_{c}},{{\Gamma }_{\gamma }},{{\Gamma }_{\mu _{\text{s},\text{REF}}^{\prime }}},{{\Gamma }_{b}} \right\}, \\ \end{eqnarray} \tag{ 15 }$$

where $\text{diag}\left\{ \cdot \right\}$ indicates a block diagonal matrix formed with the matrices in the braces, and the parameters correspond to priors of the estimated parameters such that ${{c}_{1}},{{c}_{2}},{{c}_{3}}\sim$ $\mathcal{N}\left( {{\eta }_{c}}\cdot 1,{{\Gamma }_{c}} \right)$, $\gamma \sim \mathcal{N}\left( {{\eta }_{\gamma }}\cdot 1,{{\Gamma }_{\gamma }} \right)$, $\mu _{\text{s},\text{REF}}^{\prime }\sim \mathcal{N}\left( {{\eta }_{\mu _{\text{s},\text{REF}}^{\prime }}}\cdot 1,{{\Gamma }_{\mu _{\text{s},\text{REF}}^{\prime }}} \right)$, $b\sim \mathcal{N}\left( {{\eta }_{b}}\cdot 1,{{\Gamma }_{b}} \right)$, $1\in {{\mathbb{R}}^{Q}}$ is a vector of ones, and each covariance matrix of the priors has the form

$$\begin{eqnarray}&&{{\Gamma }_{w}}={{\alpha }_{w}}{\Xi},\qquad w=c,\;\gamma ,\;\mu _{\text{s},\text{REF}}^{\prime },\;\text{or}\;b,\end{eqnarray} \tag{ 16 }$$

where ${\Xi}$ is given by (13), and ${{\alpha }_{w}}$ are scaling parameters for the marginal priors. Values used for the prior means ${{\eta }_{w}}$ and the scaling parameters ${{\alpha }_{w}}$ are given in table 2. The prior parameters were chosen by presuming good quantitative prior information of the form ${{\eta }_{w}}=\left( \max \left\{ w \right\}+\min \left\{ w \right\} \right)/2$, and ${{\alpha }_{w}}={{\left( \left( \max \left\{ w \right\}-\min \left\{ w \right\} \right)/2 \right)}^{2}}$, where w is the true γ, $\mu _{\text{s},\text{REF}}^{\prime }$, or b. Parameters ${{\eta }_{c}}$ and ${{\alpha }_{c}}$ were chosen similarly, but assuming that $\max \left\{ c \right\}=1$ and $\min \left\{ c \right\}=0$. Correlation length of $l=1\;\text{mm}$ was used in computing ${\Xi}$.

Table 2.  Prior parameters used for the MAP estimate.

w	${{\eta }_{w}}$	${{\alpha }_{w}}$
c	0.5000	0.2500
γ	0.1100	0.0001
$\mu _{\text{s},\text{REF}}^{\prime }$	$0.6750\;\text{m}{{\text{m}}^{-1}}$	$0.0506\;\text{m}{{\text{m}}^{-2}}$
b	0.7500	0.2500

For the noise statistics, a distribution of $e\sim \mathcal{N}\left( {{\eta }_{e}},{{\Gamma }_{e}} \right)$ was used in the inversion, with accurately defined mean and covariance

$$\begin{eqnarray} {{\eta }_{e}} & = & {{\left( {{0}^{\top }},{{0}^{\top }},{{0}^{\top }},{{0}^{\top }},{{0}^{\top }},{{0}^{\top }} \right)}^{\top }}, \\ {{\Gamma }_{e}} & = & {{\left( \frac{0.01}{3} \right)}^{2}}\text{diag}\left\{ \Delta _{1}^{2}I,\Delta _{2}^{2}I,\Delta _{3}^{2}I,\Delta _{4}^{2}I,\Delta _{5}^{2}I,\Delta _{6}^{2}I \right\}. \\ \end{eqnarray} \tag{ 17 }$$

The second column of figure 1 shows MAP estimates obtained with (12) for the case of smoothly varying spectral parameters. The estimate of ${{c}_{1}}$ is clearly worse than the estimates of the other chromophore concentrations ${{c}_{2}}$ and ${{c}_{3}}$. This is due to the low absorption spectrum of the first chromophore in comparison to the other two, as evident from the values listed in table 1. The main contributors to the optical absorption in the domain are the oxygenated and deoxygenated components of blood, with concentrations ${{c}_{2}}$ and ${{c}_{3}}$. Their reconstructions are comparable in visual quality to each other, and their shape resembles that of the true parameters. The Grüneisen parameter γ and the reduced scattering parameters $\mu _{\text{s},\text{REF}}^{\prime }$ and b are in their overall shape similar to the true parameters, but show more distortion in comparison to ${{c}_{2}}$ and ${{c}_{3}}$.

The second column of figure 2 shows MAP estimates in the case of sharply varying spectral parameters. The reconstructions are similar to the case of smoothly varying spectral parameters. Estimates of ${{c}_{2}}$ and ${{c}_{3}}$ show the sharp changes in the true parameters. However, the estimates of ${{c}_{1}}$, γ, $\mu _{\text{s},\text{REF}}^{\prime }$ and b do not show the sharp changes, while the overall shape of the estimates are similar to the true parameters.

4.3. Reference reconstructions

For a comparison with the direct estimation of the spectral parameters, the same parameters were estimated indirectly, and more conventionally, by first estimating the optical coefficients at each wavelength. A Bayesian approach [37] was used to obtain the optical coefficient estimates, then the spectral parameters were fitted to the estimated absorption and scattering parameters based on (3). The approach was performed under two conditions: with a known and unknown γ. These are referred to as LS and LSγ approaches for the known and unknown γ respectively. The approaches were performed for the smoothly and sharply varying material parameter cases.

Both in the case of known and unknown γ, the prior parameters for the absorption and scattering in the LS and LSγ approaches were defined similarly as in section 4.2 as ${{\mu }_{\text{a}}}\sim \mathcal{N}\left( {{\eta }_{{{\mu }_{\text{a}}}}},{{\alpha }_{{{\mu }_{\text{a}}}}}{\Xi} \right)$ and $\mu _{\text{s}}^{\prime }\sim \mathcal{N}\left( {{\eta }_{\mu _{\text{s}}^{\prime }}},{{\alpha }_{\mu _{\text{s}}^{\prime }}}{\Xi} \right)$ with ${{\eta }_{w}}=\left( \max \left\{ w \right\}+\min \left\{ w \right\} \right)/2$ and ${{\alpha }_{w}}={{\left( \left( \max \left\{ w \right\}-\min \left\{ w \right\} \right)/2 \right)}^{2}}$ with w being either ${{\mu }_{\text{a}}}$ or $\mu _{\text{s}}^{\prime }$. The prior parameters of ${{\mu }_{\text{a}}}$ and $\mu _{\text{s}}^{\prime }$ were different for each of the investigated wavelengths. In the LSγ approach, where γ was assumed to be unknown, its prior was defined as $\gamma \sim \mathcal{N}\left( {{\eta }_{\gamma }},{{\alpha }_{\gamma }}{\Xi} \right)$ where ${{\eta }_{\gamma }}=\left( \max \left\{ \gamma \right\}+\min \left\{ \gamma \right\} \right)/2$ and ${{\alpha }_{\gamma }}={{\left( \left( \max \left\{ \gamma \right\}-\min \left\{ \gamma \right\} \right)/2 \right)}^{2}}$. Values of the prior parameters are shown in table 3. As in section 4.2, accurate mean and covariance were used for the additive error in both of the approaches.

Table 3.  Prior parameters used for LS and LSγ approaches.

w	${{\eta }_{w}}$	${{\alpha }_{w}}$
$\lambda =700\;\text{nm}$
${{\mu }_{\text{a}}}$	$0.512\;\text{m}{{\text{m}}^{-1}}$	$0.1584\;\text{m}{{\text{m}}^{-2}}$
$\mu _{\text{s}}^{\prime }$	$0.675\;\text{m}{{\text{m}}^{-1}}$	$0.0506\;\text{m}{{\text{m}}^{-2}}$
$\lambda =800\;\text{nm}$
${{\mu }_{\text{a}}}$	$0.3696\;\text{m}{{\text{m}}^{-1}}$	$0.0727\;\text{m}{{\text{m}}^{-2}}$
$\mu _{\text{s}}^{\prime }$	$0.6107\;\text{m}{{\text{m}}^{-1}}$	$0.0506\;\text{m}{{\text{m}}^{-2}}$
$\lambda =900\;\text{nm}$
${{\mu }_{\text{a}}}$	$0.4882\;\text{m}{{\text{m}}^{-1}}$	$0.1294\;\text{m}{{\text{m}}^{-2}}$
$\mu _{\text{s}}^{\prime }$	$0.5590\;\text{m}{{\text{m}}^{-1}}$	$0.0347\;\text{m}{{\text{m}}^{-2}}$
γ	0.1100	0.0001

The estimates of ${{\mu }_{\text{a}}}$ and $\mu _{\text{s}}^{\prime }$ were computed at the three wavelengths using the true value of γ in the LS approach. For LS γ approach γ was estimated as well. Following the estimation of the absorption and scattering using the LS and LSγ approaches, the spectral parameter model (3) was fitted to the estimates using a least squares method.

The third column of figure 1 shows the estimates of the spectral properties of the target for the LS approach (known γ) for the smoothly varying parameter case. It is observed that the estimates of ${{c}_{1}}$ and b do not resemble the true parameter distributions, while the estimates of ${{c}_{2}}$, ${{c}_{3}}$ and $\mu _{\text{s},\text{REF}}^{\prime }$ do. Similar behaviour is observed in the case of sharply varying material parameters shown in the third column of figure 2.

The fourth column of figure 1 shows the estimated spectral parameters of the target for the LSγ approach (unknown γ) for the smoothly varying parameter case. The estimates of ${{c}_{2}}$ and ${{c}_{3}}$ seem to resemble the true parameters. It is oberserved that the estimate of γ is not very good and that the estimate of $\mu _{\text{s},\text{REF}}^{\prime }$ is worse than in the LS approach. As in the LS approach, the estimate of b using the LSγ does not resemble the true parameter distribution. The fourth column of figure 2 shows the estimates for the sharply varying material parameter case. In this case $\mu _{\text{s},\text{REF}}^{\prime }$ resembles the true parameter distribution even less than in the case of smoothly varying material parameters. Otherwise the estimates behave as in the case of the smooth material parameters.

4.4. Comparison of the reconstructions

The relative error of the estimated spectral parameters was computed using

$$\begin{eqnarray}&&{{E}_{h}}=100\%\cdot \frac{||h-h\prime ||}{||h||},\qquad h,h\prime ={{c}_{1}},\;{{c}_{2}},\;{{c}_{3}},\;\gamma ,\;\mu _{\text{s},\text{REF}}^{\prime },b,{{\mu }_{\text{a}}}\;\text{or}\;\mu _{\text{s}}^{\prime },\end{eqnarray} \tag{ 18 }$$

where h is the true parameter and $h\prime$ is the estimate. Table 4 shows the relative errors of the MAP estimates and the estimates obtained with the LS and LSγ approaches. The relative errors are shown for the smoothly and sharply varying parameter cases. Based on the relative errors it is observed that the direct estimation of the spectral parameters results in better estimates than by first estimating the optical parameters at each wavelength separately and then fitting the spectral parameters based on the spectral parameter model. Estimates of $\mu _{\text{s},\text{REF}}^{\prime }$ with the LS approach (known γ assumption) were better than when using the MAP estimate. However, this assumption is not usually valid as the Grüneisen parameter is typically unknown.

Table 4.  Relative errors of the estimates in percentage for the smoothly and sharply varying parameter cases. Relative error shown for the MAP estimate, as well as for the LS and LSγ approaches.

	Smooth parameter case			Sharp parameter case
	MAP	LS	LSγ	MAP	LS	LSγ
${{c}_{1}}$	91.8	331.0	442.1	86.8	345.0	386.5
${{c}_{2}}$	6.5	15.1	17.1	8.3	13.7	19.7
${{c}_{3}}$	7.9	17.4	21.4	11.8	21.4	26.3
γ	2.9		4.9	4.4		5.0
$\mu _{\text{s},\text{REF}}^{\prime }$	9.4	7.7	12.8	12.6	8.7	13.7
b	17.3	61.0	69.2	25.2	58.4	56.8

In this work a Bayesian approach to SQPAT was described. It employed generic and widely-used spectral tissue models, and was tested in two problems, each composed of three different chromophores with biologically relevant tissue parameters. The first problem had spatially smoothly varying spectral parameters, while the second had sharp changes. In both cases all the parameters (three chromophore concentrations, Grüneisen parameter, and two scattering parameters from Mie scattering theory) were reconstructed to good accuracy from the simulated measurement data formed with six measurement composed of two illuminations at three different wavelengths.

The superiority of the direct estimation of the spectral parameters was demonstrated numerically, by comparing the estimates to ones obtained by first estimating the absorption, scattering and the Grüneisen parameters separately at the three wavelengths and then fitting the spectral parameter model to these estimates. The latter approach failed to provide as accurate estimates as the direct estimation of the spectral parameters.

The approach presented allows the inclusions of quantitative prior information of the estimated parameters into the reconstruction. In addition, although not investigated in this paper, this approach also makes it possible to include information regarding the noise statistics in the reconstruction, as in [37]. As the inverse problem in SQPAT is two-fold, any error in the initial pressure distribution estimate (formed from the noisy measured transient ultrasonic waves) is certain to affect the estimates obtained in the optical inverse problem. The error in the initial pressure distribution estimate might have a nondiagonal and spatially nonuniform distribution in terms of its covariance structure, when it is modeled as a random normally distributed noise process. This could be a consequence of nonideal measurements of the acoustic signal, for example, limited view measurement on the boundary, large area transducers causing averaging of the received signal, or band-limited response of the transducer. The error could also be a consequence of neglecting acoustic inhomogeneities in the investigated target, or using approximative reconstruction algorithms to obtain the initial pressure distribution estimate. The error might be of numerical source by nature as well. An example of such could arise from errors caused by insufficient discretization in the numerical estimation of the intial pressure distribution. Utilizing the Bayesian approach should allow compensation of the said types of imperfections.

For this investigation, the DA was used as the model for light propagation in tissues. A more accurate model for the situation, considering the size of the domain and the values of absorption and scattering, is the RTE [46, 47]. However, practice has shown that the DA works as well in QPAT, at least when the main interest is the optical absorption parameter [35, 38, 40, 42, 55]. This is probably due to the nature of the data in QPAT as it is more sensitive to the absorption, than it is to the scattering [26]. Therefore, it was chosen as the light transport model for this work. In the future, this same approach will be tested using the RTE. Furthermore, compensating the modelling error caused by using the DA with respect to RTE through Bayesian approximation error modelling could be considered, as in [56].

Extension of the presented approach from two dimensions to three is not without issues. As the number of spatial dimensions increases, so does the total number of discretization points involved in the computation. This means that the size of the covariance matrices in the computation and the simulation time of the forward model increase as well. These factors will make the reconstructions slower and more memory demanding in three dimensions. In this work, the time to compute the estimates of the spectral parameters using the approach presented was around 31 min using MATLAB (R2011b, The Mathworks) code on a workstation with two Intel Xeon E5649 processors.

This work has been supported by the Academy of Finland (projects 136220, 140984, 272803, and 250215 Finnish Centre of Excellence in Inverse Problems Research) and by the strategic funding of the University of Eastern Finland.

Variational form of (2) is

$$\begin{eqnarray}&&{{\mathop{\int }^{}}_{\Omega }}\kappa \nabla \Phi \cdot \nabla \psi \;\text{d}r+{{\mathop{\int }^{}}_{\Omega }}\kappa {{\mu }_{\text{a}}}\Phi \psi \;\text{d}r+{{\mathop{\int }^{}}_{\partial \Omega }}\frac{2{{\zeta }_{k}}}{A}\Phi \psi \;\text{d}r={{\mathop{\int }^{}}_{\partial \Omega }}\frac{2}{A}s\psi \;\text{dr},\end{eqnarray} \tag{ A.1 }$$

where $\psi =\psi \left( r \right)$ is a test function. Approximating κ, ${{\mu }_{\text{a}}}$, s and $\Phi$ as sum of linear basis functions ${{\phi }_{n}}$ (n = i or q)

$$\begin{eqnarray} \kappa \left( r \right) & {} & \approx \underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\kappa \left( {{r}_{i}} \right){{\phi }_{i}}\left( r \right),\qquad {{\mu }_{a}}\left( r \right)\approx \underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,{{\mu }_{\text{a}}}\left( {{r}_{i}} \right){{\phi }_{i}}\left( r \right), \\ s\left( r \right) & {} & \approx \underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,s\left( {{r}_{i}} \right){{\phi }_{i}}\left( r \right),\qquad \Phi \left( r \right)\approx \underset{q=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\Phi \left( {{r}_{q}} \right){{\phi }_{q}}\left( r \right), \\ \end{eqnarray} \tag{ A.2 }$$

and inserting (A.2) into (A.1) with choice $\psi \left( r \right)={{\phi }_{p}}\left( r \right)$ for the test function results in the discretized form

$$\begin{eqnarray} &&\underset{q=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\left[ \kappa \left( {{r}_{i}} \right){{\mathop{\int }^{}}_{\Omega }}{{\phi }_{i}}\nabla {{\phi }_{q}}\cdot \nabla {{\phi }_{p}}\;\text{d}r \right]\Phi \left( {{r}_{q}} \right) \\ &&\quad +\underset{q=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\left[ {{\mu }_{\text{a}}}\left( {{r}_{i}} \right){{\mathop{\int }^{}}_{\Omega }}{{\phi }_{i}}{{\phi }_{q}}{{\phi }_{p}}\;\text{d}r \right]\Phi \left( {{r}_{q}} \right) \\ &&\quad +\underset{q=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\left[ \frac{2{{\zeta }_{k}}}{A}{{\mathop{\int }^{}}_{\partial \Omega }}{{\phi }_{q}}{{\phi }_{p}}\;\text{d}r \right]\Phi \left( {{r}_{q}} \right)=\underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\left[ \frac{2}{A}{{\mathop{\int }^{}}_{\partial \Omega }}{{\phi }_{i}}{{\phi }_{p}}\;\text{d}r \right]s\left( {{r}_{i}} \right). \\ \end{eqnarray} \tag{ A.3 }$$

Note that $\kappa \left( {{r}_{i}} \right)={{\left( k\left( {{\mu }_{\text{a}}}\left( {{r}_{i}} \right)+\mu _{\text{s}}^{\prime }\left( {{r}_{i}} \right) \right) \right)}^{-1}}$ in (A.3). Though s only needs to be defined on the boundary $\partial \Omega$, in the above it is defined in the entire domain $\Omega$ for convenience. Defining matrices ${{S}_{i}},\;{{M}_{i}},\;B\in {{\mathbb{R}}^{Q\times Q}}$ by their value in pth row and qth column, and vertical vector ${{R}_{i}}\in {{\mathbb{R}}^{Q}}$ by its value in the pth row as

$$\begin{eqnarray} &&{{\left\{ {{S}_{i}} \right\}}_{pq}}={{\mathop{\int }^{}}_{\Omega }}{{\phi }_{i}}\nabla {{\phi }_{q}}\cdot \nabla {{\phi }_{p}}\;\text{d}r,\qquad {{\left\{ {{M}_{i}} \right\}}_{pq}}={{\mathop{\int }^{}}_{\Omega }}{{\phi }_{i}}{{\phi }_{q}}{{\phi }_{p}}\;\text{d}r, \\ &&{{\left\{ B \right\}}_{pq}}=\frac{2{{\zeta }_{k}}}{A}{{\mathop{\int }^{}}_{\partial \Omega }}{{\phi }_{q}}{{\phi }_{p}}\;\text{d}r,\qquad {{\left\{ {{R}_{i}} \right\}}_{p}}=\frac{2}{A}{{\mathop{\int }^{}}_{\partial \Omega }}{{\phi }_{i}}{{\phi }_{p}}\;\text{d}r, \\ \end{eqnarray} \tag{ A.4 }$$

equation (A.3) can be expressed in matrix notation

$$\begin{eqnarray}&&K\left( {{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime } \right)\hat{\Phi}\,=v\left( \hat{s}\, \right),\end{eqnarray} \tag{ A.5 }$$

where

$$\begin{eqnarray} &&{{\hat{\mu}\,}_{\text{a}}}={{\left( {{\mu }_{\text{a}}}\left( {{r}_{1}} \right),\ldots ,{{\mu }_{\text{a}}}\left( {{r}_{Q}} \right) \right)}^{\top }} \\ &&\hat{\mu}\,_{\text{s}}^{\prime }={{\left( \mu _{\text{s}}^{\prime }\left( {{r}_{1}} \right),\ldots ,\mu _{\text{s}}^{\prime }\left( {{r}_{Q}} \right) \right)}^{\top }} \\ &&\hat{s}\,={{\left( s\left( {{r}_{1}} \right),\ldots ,s\left( {{r}_{Q}} \right) \right)}^{\top }} \\ &&K\left( {{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime } \right)=\underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,\frac{1}{k\left( {{\mu }_{\text{a}}}\left( {{r}_{i}} \right)+\mu _{\text{s}}^{\prime }\left( {{r}_{i}} \right) \right)}{{S}_{i}}+\underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,{{\mu }_{\text{a}}}\left( {{r}_{i}} \right){{M}_{i}}+B, \\ &&\hat{\Phi}\,={{\left( \Phi \left( {{r}_{1}} \right),\ldots ,\Phi \left( {{r}_{Q}} \right) \right)}^{\top }}\in {{\mathbb{R}}^{Q}}, \\ &&v\left( \hat{s}\, \right)=\underset{i=1}{\overset{Q}{\mathop{\mathop{\sum }^{}}}}\,s\left( {{r}_{i}} \right){{R}_{i}}. \\ \end{eqnarray} \tag{ A.6 }$$

A finite element approximation of fluence $\Phi$ can readily be obtained from (A.5) with matrix inversion, from which an approximation of the initial pressure distribution (1) is obtained as

$$\begin{eqnarray}&&{{p}_{0}}\left( {{r}_{q}};\hat{\gamma}\,,{{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right)=\gamma \left( {{r}_{q}} \right){{\mu }_{\text{a}}}\left( {{r}_{q}} \right)\Phi \left( {{r}_{q}};{{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right).\end{eqnarray} \tag{ A.7 }$$

Defining ${{\hat{c}\,}_{n}}={{\left( {{c}_{n}}\left( {{r}_{1}} \right),\ldots ,{{c}_{n}}\left( {{r}_{Q}} \right) \right)}^{\top }}\in {{\mathbb{R}}^{Q}}$ ($n=1,\ldots ,N$), $\hat{\gamma}\,={{\left( \gamma \left( {{r}_{1}} \right),\ldots ,\gamma \left( {{r}_{Q}} \right) \right)}^{\top }}$ $\in {{\mathbb{R}}^{Q}}$, $\hat{\mu}\,_{\text{s},\text{REF}}^{\prime }={{\left( \mu _{\text{s},\text{REF}}^{\prime }\left( {{r}_{1}} \right),\ldots ,\mu _{\text{s},\text{REF}}^{\prime }\left( {{r}_{Q}} \right) \right)}^{\top }}\in {{\mathbb{R}}^{Q}}$ and $\hat{b}\,={{\left( b\left( {{r}_{1}} \right),\ldots ,b\left( {{r}_{Q}} \right) \right)}^{\top }}\in {{\mathbb{R}}^{Q}}$, and substituting ${{\hat{\mu}\,}_{\text{a}}}$ and $\hat{\mu}\,_{\text{s}}^{\prime }$ with their discrete spectral representatations at wavelength λ

$$\begin{eqnarray} &&{{\mu }_{\text{a}}}\left( {{r}_{q}} \right)\to {{\mu }_{\text{a}}}\left( {{r}_{q}};{{\hat{c}\,}_{1}},\ldots ,{{\hat{c}\,}_{N}},\lambda \right)=\underset{n=1}{\overset{N}{\mathop{\mathop{\sum }^{}}}}\,{{c}_{n}}\left( {{r}_{q}} \right){{\mu }_{\text{a},n}}\left( \lambda \right) \\ &&\mu _{\text{s}}^{\prime }\left( {{r}_{q}} \right)\to \mu _{\text{s}}^{\prime }\left( {{r}_{q}};\hat{\mu}\,_{\text{s},\text{REF}}^{\prime },\hat{b}\,,\lambda \right)=\mu _{\text{s},\text{REF}}^{\prime }\left( {{r}_{q}} \right){{\left( \frac{\lambda }{{{\lambda }_{\text{REF}}}} \right)}^{-b\left( {{r}_{q}} \right)}} \\ \end{eqnarray} \tag{ A.8 }$$

in (A.7) results in

$$\begin{eqnarray}&&{{p}_{0}}\left( {{r}_{q}};\hat{\gamma}\,,{{\hat{c}\,}_{1}},\ldots ,{{\hat{c}\,}_{N}},\hat{\mu}\,_{\text{s},\text{REF}}^{\prime },\hat{b}\,,\lambda ,\hat{s}\, \right)=\gamma \left( {{r}_{q}} \right){{\mu }_{\text{a}}}\left( {{r}_{q}} \right)\Phi \left( {{r}_{q}};{{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right),\end{eqnarray} \tag{ A.9 }$$

where

$$\begin{eqnarray} &&{{\mu }_{\text{a}}}\left( {{r}_{q}} \right)={{\mu }_{\text{a}}}\left( {{r}_{q}};{{\hat{c}\,}_{1}},\ldots ,{{\hat{c}\,}_{N}},\lambda \right) \\ &&\mu _{\text{s}}^{\prime }\left( {{r}_{q}} \right)=\mu _{\text{s}}^{\prime }\left( {{r}_{q}};\hat{\mu}\,_{\text{s},\text{REF}}^{\prime },\hat{b}\,,\lambda \right). \\ \end{eqnarray} \tag{ A.10 }$$

Forward models (5) and (6) follow from (A.9)–(A.10) by replacing λ and s with ${{\lambda }_{m}}$, ${{s}_{m}}$, and defining $x={{\left( \hat{c}\,_{1}^{\top },\ldots ,\hat{c}\,_{N}^{\top },{{\hat{\gamma}\,}^{\top }},\hat{\mu}\,_{\text{s},\text{RE}{{\text{F}}^{\top }}}^{\prime },{{\hat{b}\,}^{\top }} \right)}^{\top }}$.

Partial derivatives of (A.9) forming the Jacobians with respect to ${{c}_{n}}\left( {{r}_{p}} \right)$, $\gamma \left( {{r}_{p}} \right)$, $\mu _{\text{s},\text{REF}}^{\prime }\left( {{r}_{p}} \right)$ and $b\left( {{r}_{p}} \right)$ are

$$\begin{eqnarray} &&\frac{\partial {{p}_{0}}\left( {{r}_{q}} \right)}{\partial {{c}_{n}}\left( {{r}_{p}} \right)}=\gamma \left( {{r}_{q}} \right){{\mu }_{\text{a},n}}\left( \lambda \right)\Phi \left( {{r}_{q}};{{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right)+\gamma \left( {{r}_{q}} \right){{\mu }_{\text{a}}}\left( {{r}_{q}} \right){{\mu }_{\text{a},n}}\left( \lambda \right){{D}_{{{\mu }_{\text{a}}}}}\left( {{r}_{q}};{{r}_{p}} \right) \\ && \frac{\partial {{p}_{0}}\left( {{r}_{q}} \right)}{\partial \gamma \left( {{r}_{p}} \right)}=\left\{ \begin{array}{} {{\mu }_{\text{a}}}\left( {{r}_{q}} \right)\Phi \left( {{r}_{q}};{{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right), & q=p \\ 0, & q\ne p \\ \end{array} \right. \\&& \frac{\partial {{p}_{0}}\left( {{r}_{q}} \right)}{\partial \mu _{\text{s},\text{REF}}^{\prime }\left( {{r}_{p}} \right)}=\gamma \left( {{r}_{q}} \right){{\mu }_{\text{a}}}\left( {{r}_{q}} \right){{\left( \frac{\lambda }{{{\lambda }_{\text{REF}}}} \right)}^{-b\left( {{r}_{p}} \right)}}{{D}_{\mu _{\text{s}}^{\prime }}}\left( {{r}_{q}};{{r}_{p}} \right) \\ && \frac{\partial {{p}_{0}}\left( {{r}_{q}} \right)}{\partial b\left( {{r}_{p}} \right)}=-\gamma \left( {{r}_{q}} \right){{\mu }_{\text{a}}}\left( {{r}_{q}} \right)\mu _{\text{s},\text{REF}}^{\prime }\left( {{r}_{p}} \right)\ln \left( \frac{\lambda }{{{\lambda }_{\text{REF}}}} \right){{\left( \frac{\lambda }{{{\lambda }_{\text{REF}}}} \right)}^{-b\left( {{r}_{p}} \right)}}{{D}_{\mu _{\text{s}}^{\prime }}}\left( {{r}_{q}};{{r}_{p}} \right), \\ \end{eqnarray} \tag{ B.1 }$$

where ${{p}_{0}}\left( {{r}_{q}} \right)={{p}_{0}}\left( {{r}_{q}};\hat{\gamma}\,,{{\hat{c}\,}_{1}},\ldots ,{{\hat{c}\,}_{N}},\hat{\mu}\,_{\text{s},\text{REF}}^{\prime },\hat{b}\,,\lambda ,\hat{s}\, \right)$, and

$$\begin{eqnarray} &&{{D}_{{{\mu }_{\text{a}}}}}\left( {{r}_{q}};{{r}_{p}} \right)=\frac{\partial \Phi \left( {{r}_{q}};{{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right)}{\partial {{\mu }_{\text{a}}}\left( {{r}_{p}} \right)} \\ &&{{D}_{\mu _{\text{s}}^{\prime }}}\left( {{r}_{q}};{{r}_{p}} \right)=\frac{\partial \Phi \left( {{r}_{q}};{{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right)}{\partial \mu _{\text{s}}^{\prime }\left( {{r}_{p}} \right)} \\ &&{{\mu }_{\text{a}}}\left( {{r}_{q}} \right)={{\mu }_{\text{a}}}\left( {{r}_{q}};{{\hat{c}\,}_{1}},\ldots ,{{\hat{c}\,}_{N}},\lambda \right) \\ &&\mu _{\text{s}}^{\prime }\left( {{r}_{q}} \right)=\mu _{\text{s}}^{\prime }\left( {{r}_{q}};\hat{\mu}\,_{\text{s},\text{REF}}^{\prime },\hat{b}\,,\lambda \right). \\ \end{eqnarray} \tag{ B.2 }$$

Partial derivatives ${{D}_{{{\mu }_{\text{a}}}}}\left( {{r}_{q}};{{r}_{p}} \right)$ and ${{D}_{\mu _{\text{s}}^{\prime }}}\left( {{r}_{q}};{{r}_{p}} \right)$ are computed by differentiating (A.5) with respect to ${{\mu }_{\text{a}}}\left( {{r}_{p}} \right)$ and $\mu _{\text{s}}^{\prime }\left( {{r}_{p}} \right)$ and solving for the partial derivatives of the fluence

$$\begin{eqnarray} \frac{\partial \hat{\Phi}\,\left( {{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right)}{\partial {{\mu }_{\text{a}}}\left( {{r}_{p}} \right)} & = & -K{{\left( {{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime } \right)}^{-1}}\left[ \frac{1}{k{{\left( {{\mu }_{\text{a}}}\left( {{r}_{p}} \right)+\mu _{\text{s}}^{\prime }\left( {{r}_{p}} \right) \right)}^{2}}}{{S}_{p}}+{{M}_{p}} \right]\hat{\Phi}\,, \\ \frac{\partial \hat{\Phi}\,\left( {{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right)}{\partial \mu _{\text{s}}^{\prime }\left( {{r}_{p}} \right)} & = & -K{{\left( {{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime } \right)}^{-1}}\left[ \frac{1}{k{{\left( {{\mu }_{\text{a}}}\left( {{r}_{p}} \right)+\mu _{\text{s}}^{\prime }\left( {{r}_{p}} \right) \right)}^{2}}}{{S}_{p}} \right]\hat{\Phi}\,, \\ \end{eqnarray} \tag{ B.3 }$$

where the left-hand sides are

$$\begin{eqnarray} \frac{\partial \hat{\Phi}\,\left( {{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right)}{\partial {{\mu }_{\text{a}}}\left( {{r}_{p}} \right)} & = & {{\left( {{D}_{{{\mu }_{\text{a}}}}}\left( {{r}_{1}};{{r}_{p}} \right),\ldots ,{{D}_{{{\mu }_{\text{a}}}}}\left( {{r}_{Q}};{{r}_{p}} \right) \right)}^{\top }} \\ \frac{\partial \hat{\Phi}\,\left( {{\hat{\mu}\,}_{\text{a}}},\hat{\mu}\,_{\text{s}}^{\prime },\hat{s}\, \right)}{\partial \mu _{\text{s}}^{\prime }\left( {{r}_{p}} \right)} & = & {{\left( {{D}_{\mu _{\text{s}}^{\prime }}}\left( {{r}_{1}};{{r}_{p}} \right),\ldots ,{{D}_{\mu _{\text{s}}^{\prime }}}\left( {{r}_{Q}};{{r}_{p}} \right) \right)}^{\top }}. \\ \end{eqnarray} \tag{ B.4 }$$