Statistical independence in nonlinear model-based inversion for quantitative photoacoustic tomography

Lu An; Teedah Saratoon; Martina Fonseca; Robert Ellwood; Ben Cox

doi:10.1364/BOE.8.005297

1. Introduction

Photoacoustic tomography is a non-invasive biomedical imaging technique [1] relying on the absorption of optical energy and the generation of ultrasound waves. It has a relatively low cost of implementation and has the advantage of combining the relatively large penetration depth and high resolution of ultrasound imaging with optical absorption based contrast which provides high specificity. In quantitative photoacoustic tomography, the unique optical absorption spectra of the chromophores are exploited in spectroscopic inversions of the multiwavelength photoacoustic images in order to estimate the concentration of each chromophore. The key endogenous chromophores of interest for quantitative photoacoustic tomography are oxy- and deoxyhaemoglobin, because the ratio of their concentrations is related to the blood oxygenation, which is an important physiological parameter. Spectroscopic techniques are also valuable tools for contrast-enhanced photoacoustic molecular imaging applications, where the detection and quantification of the local accumulation of genetically encoded probes and extrinsically administered contrast agents [2] can provide information on biological processes, drug delivery, disease development and treatment response.

Quantifying the chromophore concentrations is a challenging task because the absorbed optical energy density is not linearly related to the chromophore concentrations, but a product of the absorption coefficient and light fluence, which varies both spatially and spectrally and depends on the unknown chromophore concentrations. The model-based inversion scheme [3–8] is a quantification method that accounts for the effect of the fluence by including a numerical model of the fluence distribution in an iterative scheme to solve the inverse problem. Thus, it has the potential to provide accurate estimation of the absolute chromophore concentrations in complex tissue structures. However, in experimental settings, both the approximate nature of the model, and the difficulty in determining all the “known” model parameters, such as the absorption and scattering spectra and the intensity profile of the excitation beam, to high accuracy, lead to model-mismatch. This poses challenges on the practical implementation of model-based inversion schemes for in vivo imaging. Previous studies [5,8] have relied on reducing the number of unknown variables by segmenting the images into regions with piece-wise constant optical properties to increase the robustness of the inversion scheme. In this study, we exploit instead the statistical independence between certain chromophores to improve the accuracy and robustness of model-based inversion.

Statistical independence has been utilised to spectrally decompose the chromophores via independent component analysis (ICA) [9]. ICA is a fast and simple unmixing algorithm based on the assumption that the photoacoustic images are linear mixtures of the independent source components representing the chromophores. Glatz et al [10] showed that ICA can result in more accurate unmixing than a linear spectroscopic inversion. ICA has subsequently been used to unmix genetic reporters [11] and contrast agents targeted to apoptotic cells [12] or cells expressing selectin [13]. However, ICA is unable to estimate the absolute concentrations of the chromophores and it does not account for the nonlinear fluence distortion.

In this paper, we propose incorporating statistical independence as additional information in the nonlinear model-based inversion method. This involves including a measure of the independence in the error functional in addition to the least-squares data error. The aim is to reduce the quantification errors caused by inaccurate forward modelling of the fluence. The effect of using the statistical independence in the inversion is investigated using both numerical and experimental tissue mimicking phantoms and the quantification results are compared to the model-based inversion using only the data error.

2. Mutual information error functional

Detailed descriptions of statistical independence and its applications in imaging can be found in Refs [14,16,17]. However, as the concept has not yet been widely used in photoacoustic imaging, a brief summary of the essential aspects will be given here. The spatial distributions of certain biomarkers or molecular probes are statistically independent from other tissue chromophores in some cases of practical interest. Examples include fluorescent probes that can be targeted to disease-specific receptors whose spatial distribution is unrelated to that of the blood and the background tissue, and cells that have been genetically modified to express optical absorbers which can be found at locations independent from other absorbers. The statistical independence of two chromophores is related to the histogram of the values of the individual chromophore concentrations, and the joint histogram of the concentrations of both chromophores. The mathematical definition of statistical independence states that two random variables, y₁ and y₂, are statistically independent if their joint probability density function (PDF), ρ_y₁,y₂(y₁, y₂), is the product of their marginal PDFs, ρ_y₁(y₁) and ρ_y₂(y₂) [14], such that

ρ_{y_{1}, y_{2}} (y_{1}, y_{2}) = ρ_{y_{1}} (y_{1}) ρ_{y_{2}} (y_{2}),

where y₁ and y₂ denote possible values of y₁ and y₂. The degree of independence between variables can be measured using the mutual information (MI). MI is an estimate of the amount of information one variable provides on another variable. Variables with higher independence have lower MI, which means that they contain less information about each other. In quantitative photoacoustic tomography, two chromophores are considered independent if the knowledge of the concentration of one chromophore at a location does not affect the estimate of the other chromophore’s concentration at the same location. For example, if a contrast agent is independent of the blood, then the estimation of the contrast agent concentration at a voxel does not in any way predict the blood concentration there. As a counter example, oxy- and deoxyhaemoglobin are not independent chromophores: if a voxel is found to contain a high concentration of deoxyhaemoglobin, then the likelihood that a high concentration of oxyhaemoglobin will be found at the same pixel increases, as the voxel is likely to represent a blood vessel. The MI, I, between y₁ and y₂ is given by

I (y_{1}, y_{2}) = ℋ (y_{1}) + ℋ (y_{2}) - ℋ (y_{1}, y_{2}),

where ℋ(y_k) and ℋ(y₁, y₂) are the entropy and the joint entropy of y₁ and y₂ respectively. They are defined by

ℋ (y_{k}) = - \int_{y_{k}} ρ_{y_{k}} (y_{k}) log ρ_{y_{k}} (y_{k}) d y_{k},

where k = 1 or 2, and

ℋ (y_{1}, y_{2}) = - \int_{y_{1}} \int_{y_{2}} ρ_{y_{1}, y_{2}} (y_{1}, y_{2}) log ρ_{y_{1}, y_{2}} (y_{1}, y_{2}) d y_{1} d y_{2} .

Similarly, the MI between multiple random variables is defined by

I (y_{1}, \dots, y_{K}) = \sum_{k = 1}^{K} ℋ (y_{k}) - ℋ (y_{1}, \dots, y_{K}) .

As indicated in Eqs. (2–4), the MI depends on the probabilities of the variables. We consider the concentrations of K independent chromophores c₁, ..., c_K as continuous random variables with the PDFs ρ_c₁, ..., ρ_{c_K}. However, the underlying probabilities of the chromophore concentrations are unknown and therefore ρ_{c_k} needs to be estimated based on the instantiations of the random variable, which are the concentrations of the chromophore at different voxels [c_k,1, ..., c_k,M]. The total number of voxels, which is equal to the total number of instantiations, is denoted with M. Using these data points, the PDF can be approximated with a kernel density estimator. Conceptually, the kernel density estimation method involves placing a smoothly varying spread of values, or kernel, on the value of each data point. The probability is then approximated by the sum of all kernels, such that

{\overset{⌣}{ρ}}_{c_{k}} (ξ_{k}) = \frac{1}{M} \sum_{m = 1}^{M} κ (ξ_{k} - c_{k, m})

and

{\overset{⌣}{ρ}}_{c_{1}, \dots, c_{K}} (ξ_{1}, \dots, ξ_{K}) = \frac{1}{M} \sum_{m = 1}^{M} \prod_{k = 1}^{K} κ (ξ_{k} - c_{k, m}),

where ξ_k are the values at which the PDF is evaluated and κ(x) is a kernel function that satisfies κ(x) ≥ 0 and

\int_{- \infty}^{\infty} κ (x) d x = 1

. The breve symbol (˘) is used to denote the estimation of a quantity. The motivation for using the kernel density estimator instead of a simple histogram approximation for the probabilities is that the former ensures continuity and differentiability of the estimated PDFs, provided that a continuous and differentiable kernel function is used. To satisfy those criteria, the Gaussian kernel is chosen for this study:

κ (x) = \frac{1}{h \sqrt{2 π}} exp (- x^{2} / 2 h^{2}),

where h is the kernel width, which determines the amount of smoothing in the estimations. It is calculated using the expression for the optimal kernel width normally distributed data, given by [15]

h = {(\frac{4}{3 M})}^{\frac{1}{5}} σ \approx 1.06 σ M^{- \frac{1}{5}},

where σ is the standard deviation of the data.

To estimate the entropies of the chromophore concentrations, the PDF is evaluated at a set of equally spaced points denoted by [ξ_k,1, ..., ξ_k,Q], where Q is the total number of points. The integral in the definition of the entropy and the joint entropy for continuous variables in Eq. (3) and (4) can be approximated as a sum using the trapezium rule, such that

\bar{ℋ} (c_{k}) = - \sum_{q_{k} = 1}^{Q} {\overset{⌣}{ρ}}_{c_{k}} (ξ_{k, q_{k}}) log {\overset{⌣}{ρ}}_{c_{k}} (ξ_{k, q_{k}}) Δ ξ_{k}

and

\overset{⌣}{ℋ} (c_{1}, \dots, c_{K}) = - \sum_{q_{1}, \dots, q_{K} = 1}^{Q} {\overset{⌣}{ρ}}_{c_{1}, \dots, c_{K}} (ξ_{1, q_{1}}, \dots, ξ_{K, q_{K}}) log {\overset{⌣}{ρ}}_{c_{1}, \dots, c_{K}} (ξ_{1, q_{1}}, \dots, ξ_{K, q_{K}}) \prod_{k = 1}^{K} Δ ξ_{k}

where the summation symbol denotes a multiple sum of q₁, ..., q_K and Δξ_k is the spacing between the sampling points for the PDF. Using the approximated entropies, the MI can be estimated by

\overset{⌣}{I} (c_{1}, \dots, c_{K}) = \sum_{k = 1}^{K} \overset{⌣}{ℋ} (c_{k}) - \overset{⌣}{ℋ} (c_{1}, \dots, c_{K}) .

To find the most independent chromophores requires minimising the MI between the chromophores, which can be done efficiently using a gradient-based optimisation approach. The partial derivative of the MI with respect to the chromophore concentration at each voxel is given by [18]

\frac{\partial \overset{⌣}{I} (c_{1}, \dots, c_{K})}{\partial c_{k, m}} = \frac{\overset{⌣}{ℋ} (c_{k})}{\partial c_{k, m}} - \frac{\partial \overset{⌣}{ℋ} (c_{1}, \dots, c_{K})}{\partial c_{k, m}} .

The first term in the right hand side of Eq. (13) is

\frac{\partial \overset{⌣}{ℋ} (c_{k})}{\partial c_{k, m}} = \frac{\partial ℋ (c_{k})}{\partial {\overset{⌣}{ρ}}_{c_{k}}} \frac{\partial {\overset{⌣}{ρ}}_{c_{k}}}{\partial c_{k, m}},

where

\frac{\partial {\overset{⌣}{ℋ}}_{c_{k}}}{\partial {\overset{⌣}{ρ}}_{c_{k}}} = - \sum_{q_{k} = 1}^{Q} [1 + log {\overset{⌣}{ρ}}_{c_{k}} (ξ_{k, q_{k}})] Δ ξ_{k}

by the product rule, and

\frac{\partial {\overset{⌣}{ρ}}_{c_{k}}}{\partial c_{k, m}} = κ^{'} (ξ_{k, q_{k}} - c_{k, m}),

where κ′ denotes the derivative of κ, given by

κ^{'} (ξ_{k, q_{k}} - c_{k, m}) = \frac{ξ_{k, q_{k}} - c_{k, m}}{h^{2}} κ (ξ_{k, q_{k}} - c_{k, m})

for the Gaussian kernel. Similarly, the second term in Eq. (13) is

\frac{\partial \overset{⌣}{ℋ} (c_{1}, \dots, c_{K})}{\partial c_{k, m}} = \frac{\partial \overset{⌣}{ℋ} (c_{1}, \dots, c_{K})}{\partial {\overset{⌣}{ρ}}_{c_{1}, \dots, c_{K}}} \frac{\partial {\overset{⌣}{ρ}}_{c_{1}, \dots, c_{K}}}{\partial c_{k, m}},

where

\frac{\partial \overset{⌣}{ℋ} (c_{1}, \dots, c_{K})}{\partial {\overset{⌣}{ρ}}_{c_{1}, \dots, c_{K}}} = - \sum_{q_{1}, \dots q_{K} = 1}^{Q} [1 + log {\overset{⌣}{ρ}}_{c_{1}, \dots, c_{K}} (ξ_{1, q_{1}}, \dots, ξ_{K, q_{K}})] \prod_{i = 1}^{K} Δ ξ_{i}

and

\frac{\partial {\overset{⌣}{ρ}}_{c_{1}, \dots, c_{K}}}{\partial c_{k, m}} = κ^{'} (ξ_{k} - c_{k, m}) \prod_{i = 1, i \neq k}^{K} κ (ξ_{i} - c_{i, m}) .

3. Model-based inversion with statistical independence

The first step in the model-based inversion scheme used in this study is to make an initial guess for the unknown chromophore concentrations. This initial guess and the known specific absorption spectra α are used to calculate the absorption coefficient, $μ_{a} = \sum_{k} α_{k} c_{k}$ . The fluence, Φ, is then modelled using the diffusion approximation to the radiative transfer equation for the distribution of the absorption coefficient and the scattering coefficient. The scattering coefficient is assumed to be known in this study. Using the modelled fluence, the absorption coefficient and the Grüneisen parameter Γ, the modelled initial pressure $p_{m, λ_{n}}^{model}$ at voxel m and wavelength λ_n is calculated based on

p_{m, λ_{n}}^{model} = S Γ_{m} Φ_{m, λ_{n}} μ_{a_m, λ_{n}},

where S is the system calibration factor that depends on the acoustic response of the sensors, which can be assumed to be spatially homogeneous and independent of the optical wavelength. The data error, ε_d, is defined as the sum of squared differences between the modelled and the measured initial pressure,

p_{m, λ_{n}}^{meas}

, and hence the minimisation problem is given by

\underset{c_{1}, \dots, c_{K_{t}}}{argmin} (c_{1}, \dots, c_{K_{t}}) = \sum_{n = 1}^{N} \sum_{m = 1}^{M} {[p_{m, λ_{n}}^{model} (c_{1}, \dots, c_{K_{t}}) - p_{m, λ_{n}}^{meas}]}^{2},

where the total number of chromophores, wavelengths and voxels are denoted by K_t, N, and M respectively, and c_i denotes the i^th chromophore. By iteratively updating the values of the chromophore concentrations until ε_d is minimised, accurate quantification can be achieved in an idealised scenario. However, in practical applications, model-mismatch arises due to the approximations in the model and uncertainty in the model parameters (which are typically determined experimentally). This leads to the minimum of the data error occurring at erroneous chromophore concentrations, and results in inaccurate quantification. Unlike the data error, the statistical independence is a property of the distribution of the chromophores alone, rather than a function of the forward modelling. Therefore, the errors in the fluence model do not affect the MI between the chromophore concentrations, which will always have a minimum at the true solution, provided that the chromophores are statistically independent. Hence, by including a term representing the MI between the independent chromophores in the error functional, the quantification errors can be reduced. The new minimisation problem with both the data error and the MI is given by

\underset{c_{1}, \dots, c_{K_{t}}}{argmin} ε_{d + M I} (c_{1}, \dots, c_{K_{t}}) = \sum_{n = 1}^{N} \sum_{m = 1}^{M} {[p_{m, λ_{n}}^{model} (c_{1}, \dots, c_{K_{t}}) - p_{m, λ_{n}}^{meas}]}^{2} + γ \overset{⌣}{I} (c_{1}, \dots, c_{K}),

where γ is the weight parameter for the MI, and total number of chromophores may be larger than or equal to the number of independent chromophores, K_t ≥ K.

4. Generating multiwavelength photoacoustic images of tissue mimicking phantoms

The accuracy of the quantification using ε_d+MI compared to ε_d was investigated using experimental and numerical tissue mimicking phantoms containing aqueous solutions of copper(II) chloride dihydrate, nickel(II) chloride hexahydrate and black India ink (Winsor & Newton, London, UK), which represent different absorbers in the tissue, and Intralipid, which provides scattering in the medium. The copper(II) chloride dihydrate and nickel(II) chloride hexahydrate will be referred to as CuCl₂ and NiCl₂. For both the numerical and the experimental phantom, the distributions of CuCl₂ and NiCl₂ are arranged such that they are statistically independent of each other.

4.1. Experimentally acquired images

A schematic of the experimental set-up is shown in Fig. 1. The tissue mimicking phantom consists of four polythene tubes with 0.58mm inner diameter and 0.19mm wall thickness (Scientific Laboratory Supplies Ltd, Nottingham, UK) submerged in a background solution of highly diluted India ink and 1% (w/v) Intralipid, which give rise to an absorption and scattering amplitude comparable to that of typical biological tissue [19]. The tubes are arranged in a line at depths of approximately 3.6, 6.1, 8.1 and 9.8mm from the top surface of the phantom, which are all within the diffusive regime. The first and third tube from the top contain 399gL⁻¹ NiCl₂ and the second and the fourth tube contain 36gL⁻¹ CuCl₂. The absorption spectra of the chromophores and scattering spectrum of Intralipid are shown in Fig. 2. The CuCl₂ and NiCl₂ are statistically independent of each other in this phantom, which is clear from the fact that they are contained in distinct regions that are spatially separated. (However, the spatial separation is not a necessary criterion for statistical independence. For example, CuCl₂ and NiCl₂ are both also statistically independent from water, even though they can be found at the same locations.)

Fig. 1 Experimental set-up and phantom structure. The four tubes containing CuCl₂ or NiCl₂ are fixed in a vertical column and submerged in the India ink and Intralipid solution. Two orthogonal Fabry-Perot interferometer sensors are used for increased detection aperture. The fibre tip at the top of the phantom delivers the pulsed excitation beam.

Download Full Size | PDF

Fig. 2 (a) The absorption coefficients of 36gL⁻¹ CuCl₂ (squares) and 399gL⁻¹ NiCl₂ (circles). (b) The absorption coefficient of the background solution (crosses), which is a sum of the absorption of water [20] (dotted) and the India ink (dashed). (c) The scattering amplitude of 1% Intralipid as a function of wavelength [21]. A spectrophotometer (Lambda 750S, Perkin Elmer) was used to measure the transmittance of CuCl₂, NiCl₂ and India ink in order to determine their absorption spectra. Reprinted from [22].

Download Full Size | PDF

The phantom is imaged in a V-shaped photoacoustic imaging system [24,25] consisting of two orthogonal Fabry-Perot interferometer sensors [26]. This sensor geometry increases the detection aperture of the photoacoustic signals compared to a single planar detector array and hence reduces the limited-view artefacts [25]. The fibre tip was positioned vertically above the phantom to deliver the pulsed excitation light from a Nd:YAG-pumped optical parametric oscillator (GWU, Spectra-Physics, Santa Clara, USA) with 10Hz repetition rate and a pulse energy of 15–19mJ depending on wavelength. The phantom was imaged at 8 wavelengths with equal spacing between 750nm and 890nm by recording the photoacoustic time series at a 13×13mm² area with 100μm step size for both sensors. A small fraction of the light was directed to an integrating sphere to measure the pulse energy, which was used to normalise the measured signals. The beam position and the spatial distribution of the beam intensity at the surface of the phantom was found by acquiring an additional photoacoustic image at 780nm of a transparent sheet with a printed grid of highly absorbing dots (with 1mm grid spacing) which was resting on the surface. Based on the reconstruction of this image, the illumination source was approximated as a Gaussian beam with a 1/e diameter of 6.6mm.

The iterative time reversal method [27, 28] was used to reconstruct 3D images from the photoacoustic time series. A 2D cross-sectional 12×12mm² region of interest centred at the tubes was used for the optical inversion, and the dimension was reduced to 72×72 pixels to reduce the computational time and memory requirements. The 2D slices are shown for three wavelengths in Fig. 3.

Fig. 3 The 2D cross-sectional slices of the 3D reconstructed photoacoustic images which are used for the optical inversion at wavelengths 750, 830 and 890nm. The size of this region of interest is 12×12mm² and the element spacing is 166μm. As expected, the intensity of the tubes decreases with depth for all wavelengths due to the decay of the fluence.

Download Full Size | PDF

4.2. Numerically simulated images

The numerical phantom has an element spacing of 100μm and represents a 5×5mm² area with six insertions arranged in two columns, as illustrated in Fig. 4. The left column of insertions represents solutions of CuCl₂ with concentrations 12, 24 and 36gL⁻¹ in increasing order, where the top insertion has the lowest concentration. The right column of insertions contain solutions of NiCl₂ with concentrations 133, 266 and 399gL⁻¹, also in increasing order from the top insertion. The phantom is designed with increasing concentrations at the insertions with depth to improve the signal to noise ratio at the deeper insertions. The absorption of CuCl₂ and NiCl₂ are based on the measured spectra shown in Fig. 2(a) and assumed to follow linear dependence on concentration. The concentrations are chosen such that the average absorption of both columns is 0.52mm⁻¹ over the wavelength range between 750nm to 890nm, which is similar to the absorption of blood over the same range of wavelengths. Water is present in the whole phantom and the background region outside the insertions represents a solution of India ink and Intralipid, which gives rise to the same absorption and scattering amplitude as shown in Fig. 2(b–c).

Fig. 4 Diagram of the 2D numerical phantom. The phantom contains regions with different concentrations of CuCl₂ and NiCl₂. The background region contains India ink and Intralipid.

Download Full Size | PDF

The top surface of the domain was illuminated with a radially-symmetric light source with a Gaussian intensity profile with a 1/e width of 3mm. The light fluence distributions for the same 8 wavelengths as the experimental measurement in Sec. 4.1 were simulated using the diffusion approximation with the MATLAB software Toast++ [23]. The system calibration factor and the Grüneisen parameter are assumed to be known and equal to one, and the acoustic reconstruction is assumed to be perfect, such that the simulated photoacoustic images were equal to the product of the fluence and the absorption coefficient. Gaussian noise with an amplitude equal to 10% of the mean of the data, which was comparable to the magnitude of the noise in the experimental images (Sec. 4.1), was added to the simulated images.

5. Inverting for the chromophore concentrations

To investigate the effect of incorporating the MI term, the model-based inversion scheme was applied to the simulated and the experimentally acquired multiwavelength 2D photoacoustic images using both ε_d and ε_d+MI error functionals. The error functionals were minimised with the limited-memory BFGS [29] algorithm which searches for the optimal chromophore concentrations using the functional gradients of the data error term [3] and the MI term [18]. The unknown parameters in the inversion are the concentrations of CuCl₂, NiCl₂, India ink and water, and the known model parameters are the absorption spectra, the scattering distribution, the system calibration factor, the Grüneisen parameter and the light source position and width. The unknown chromophore concentrations were initialised with a spatially homogeneous value equal to the true concentration at the background. The iterative update was run for 300 iterations for the inversion of both the simulated and the experimental images using ε_d or ε_d+MI. For all inversions using ε_d+MI, the MI was calculated between CuCl₂ and NiCl₂ and the weight parameter γ of the MI term was set to zero for the first 200 iterations to avoid the algorithm being trapped in the local minima of the MI term. The difference in computation time for ε_d+MI and ε_d was negligible. The computationally efficient calculation of the MI term and its gradient was achieved using fast Fourier transforms [30], which takes advantage of the fact that Eqs. (6) and (7) have convolution structures.

5.1. Effect of model-mismatch

Two case studies were conducted using the simulated images to investigate the effect of the uncertainty in different model parameters on the quantification accuracy. In the first case, the beam diameter was set to be up to 75% smaller or larger than the true value in the inversion. In the second case, an error of up to ±75% was included in the scattering amplitude. The average errors of the estimated concentrations of CuCl₂ and NiCl₂ at the insertions (ROI) using the erroneous beam diameter are shown in Fig. 5(a), where the circles and asterisks correspond to the inversions using ε_d or ε_d+MI respectively. The results show that the increase in the percentage error in the beam diameter leads to larger quantification errors, as expected. The quantification errors for the simulated images are relatively small because only one model parameter contains error at the time. In an experimental setting, there is likely to be a combination of modelling errors, resulting in larger quantification errors. Including the MI term results in a reduction in error compared to using only the standard data error for all data points.

Fig. 5 The average errors of the estimated concentrations of CuCl₂ and NiCl₂ at the insertions as a function of errors in (a) the beam diameter or (b) the scattering amplitude in the inversion. As expected, the quantification errors increase for larger errors in the beam diameter or scattering amplitude. However, the inversions using ε_d+MI (asterisks) result in smaller errors compared to using ε_d (circles) for all data points. The individual errors for CuCl₂ and NiCl₂ show similar general trends as the average of the two. The average errors outside the ROI are <2% for inversions using both ε_d and ε_d+MI. Reprinted from [22].

Download Full Size | PDF

The inaccuracies in the scattering amplitude used in the inversion resulted in similar trends for the quantification error, as shown in Fig. 5(b). The errors are generally larger in Fig. 5(b) than 5(a), which suggests that the changes in scattering amplitude have a larger impact on the fluence distribution than changes in the beam diameter for this numerical phantom. Nonetheless, using ε_d+MI is shown to provide more accurate estimations compared to using ε_d also for this case. The relative improvement in accuracy varied between 37% and 8% with an average of 22% over all data points for both cases.

5.2. Experimental results

The multiwavelength experimental images were divided by the calibration factor and the spatially varying Grüneisen parameter before the inversion. The calibration factor was determined using a forward simulation with the true concentrations. The Grüneisen parameter of aqueous solutions of CuCl₂ and NiCl₂ are both known to increase with concentration and are given by

Γ_{i} = Γ_{H_{2} O} (1 + β_{i} c_{i}), i = {CuCl}_{2} or {NiCl}_{2}

where Γ_H₂O is the Grüneisen parameter for water, and β_i are 5.80×10⁻³Lg⁻¹ and 2.25×10⁻³Lg⁻¹ for CuCl₂ or NiCl₂ respectively [31]. In practical applications, where the true concentrations are unknown, the calibration factor can be obtained by measuring the acoustic sensitivity of the sensors [32], and the Grüneisen parameter can be included in the model as a parameter that is linearly dependent on the estimated chromophore concentrations [5,32].

The results from the model-based inversion of the experimental data are presented in Fig. 6. The estimated concentrations of CuCl₂ (top row) and NiCl₂ (bottom row) using ε_d and ε_d+MI are shown in the left and centre columns respectively in Fig. 6(a), while the true concentrations are shown in the right column. The colour scale indicates the concentrations in units of gL⁻¹. Figure 6(b) compares the estimated with the true concentrations along a line profile across the tubes. The average estimated and expected concentrations for each tube are presented in Table 1. The inversions using ε_d resulted in high overestimation of CuCl₂ in the second tube and NiCl₂ in the first tube, where the estimated concentrations are 94% and 149% larger than the true values respectively. There are also large cross-talk errors in the estimation of both contrast agents. This is most clearly seen for the estimated CuCl₂ concentration, where the third tube shows a high false-positive concentration with comparable magnitude to the concentration in the fourth tube. The accuracy of the quantification is significantly improved when the MI term is included in the error functional. The cross-talk errors for both CuCl₂ and NiCl₂ are almost completely removed when ε_d+MI is used. The absolute concentrations of the CuCl₂ is estimated accurately with an error of 3gL⁻¹ on average for the four tubes. The overestimation of the NiCl₂ concentration in the top tube remains present with ε_d+MI. However, this overestimation error is reduced when ε_d+MI is used compared to ε_d.

Fig. 6 (a) The estimated concentrations of CuCl₂ (top row) and NiCl₂ (bottom row) in units of gL⁻¹. The results from the inversions using ε_d (left column) show overestimation of the the upper tubes and large cross-talk errors, while using ε_d+MI results in more accurate quantification without cross-talk errors. The true concentrations are shown in the column to the right for comparison. (b) The estimated concentration of the CuCl₂ (top) and NiCl₂ (bottom) using ε_d (crosses) and ε_d+MI (circles) along a line across the tubes. The solid curves represent the true concentrations.

Download Full Size | PDF

Table 1. The average estimated and true concentrations of CuCl₂ (left) and NiCl₂ (right) in gL⁻¹ for each tube. The largest improvements using ε_d+MI are mostly seen for the tubes that are not expected to contain the relevant chromophore, as they suffer from significant cross-talk errors when ε_d is used. The average expected concentrations are lower than the true concentrations in the solutions due to the interpolation from the original images.

View Table

6. Discussion

Minimising the MI as well as the data error led to improved quantification accuracy for both simulated and experimental multiwavelength photoacoustic images, compared to using the data error alone. In the numerical study, despite the significant fractional decrease of the quantification errors using ε_d+MI, in absolute terms, the error was decreased only by a few per cent. This was mainly due to the fact that only one model parameter was erroneous in each inversion, while all other assumptions in the model were accurate, which resulted in relatively small quantification errors, even when only ε_d was used. In the experimental study, on the other hand, a combination of different types of modelling errors was likely to have been present simultaneously, leading to poorer quantification results in the absence of the MI term. The main causes of model-mismatch may be due to the limited size of the modelled domain, which does not account for the backscattered light from outside of the domain, and the 2D modelling of the light fluence, which assumes that the light source is constant in the direction along the tubes, while in the experiment, the beam was of circular cross-section. Other possible errors may include uncertainty in the scattering spectra and amplitude, as different values have been reported in the literature [21,33]. These errors in the model affect the calculation of $p_{m, λ_{n}}^{model}$ , which consequently also affect ε_d. The MI is not affected by the fluence modelling errors, because MI is calculated based on only the distribution of the estimated chromophore concentrations in each iteration, and does not require the forward modelling of $p_{m, λ_{n}}^{model}$ . Therefore, the quantification errors were greatly reduced when ε_d+MI was used in the experimental study. These results suggest that incorporating the statistical independence can improve the robustness of model-based inversion schemes for independent chromophores and thus potentially enhance their applicability to pre-clinical or clinical imaging studies.

In order to obtain accurate results with the inversion using ε_d+MI, it is necessary to use an appropriate weight parameter γ for the MI term. The weight was determined through manual trial and error. The solution was continuously dependent on the weight parameter, in the sense that small changes in the weight parameter resulted in small changes in the solution. The same weight was used for all inversions of the simulated and the experimental data, despite the differences in the data and/or the model parameters. This suggests that the concentration estimates were not highly sensitive to small variations of the weighting of the MI term around this value and, although non-trivial [34, 35], it may be possible to develop a general method for finding the optimal weight parameter for different types of applications.

The 2D fluence model based on the diffuse approximation assumes that the features are constant in the third dimension and located at depths within the diffusive regime. These assumptions are appropriate for the phantom geometry used in this study. However, full 3D fluence modelling will be required for applications of the model-based inversion in biological tissue with complex structures. More accurate modelling of the fluence for the superficial layer can be achieved by incorporating the delta-Eddington approximation [5] or using the radiative transfer equation [36]. The calculation of the MI can be straightforwardly extended to 3D without causing significant increase in computation time. Furthermore, using the MI term in the error functional is also compatible with other regularisation methods such as the total-variation regulariser [37–39].

7. Conclusion

We proposed exploiting the statistical independence between certain chromophores in the model-based inversion method by minimising the MI between the independent chromophores in addition to the data error. The improvement in the accuracy of the estimated chromophore concentrations was demonstrated using both numerical simulations and an experimental phantom. The results suggest that the sensitivity of the model-based inversion to model-mismatch can be reduced by incorporating the additional information of statistical independence. Thus, the robustness and hence usefulness of the inversion scheme can potentially be improved for in vivo imaging experiments.

Funding

This work was funded by the UCL EPSRC Centre for Doctoral Training in Medical Imaging.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. P. Beard, “Biomedical photoacoustic imaging,” Interface Focus 1, 602–631 (2011). [CrossRef]

2. J. Weber, P. C. Beard, and S. E. Bohndiek, “Contrast agents for molecular photoacoustic imaging,” Nat. Methods 13, 639–650 (2016). [CrossRef] [PubMed]

3. B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A 26(2), 443–455 (2009). [CrossRef]

4. Y Sun, E Sobel, and H Jiang, “Quantitative three-dimensional photoacoustic tomography of the finger joints: an in vivo study,” J. Biomed. Opt. 14(6), 064002 (2009). [CrossRef]

5. J. Laufer, B. Cox, E. Zhang, and P. Beard, “Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt. 49(8), 1219–1233 (2010). [CrossRef] [PubMed]

6. G. Bal and K. Ren, “On multi-spectral quantitative photoacoustic tomography in diffusive regime,” Inverse Probl. 28(2), 025010 (2012). [CrossRef]

7. T. Tarvainen, B. T. Cox, J. P. Kaipio, and S. R. Arridge, “Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography,” Invserse Prob. 28(8), 084009 (2012). [CrossRef]

8. F. M. Brochu, J. Brunker, J. Joseph, M. R. Tomaszewski, S. Morscher, and S. E. Bohndiek, “Towards quantitative evaluation of tissue absorption coefficients using light fluence correction in optoacoustic tomography,” IEEE Trans. Med. Imag. 36(1), 322–331 (2017). [CrossRef]

9. A. Hyvarinen and E. Oja, “Independent component analysis: algorithms and applications,” Neural Netw. 13, 411–430 (2000). [CrossRef] [PubMed]

10. J. Glatz, N. C. Deliolanis, A. Buehler, D. Razansky, and V. Ntziachristos, “Blind source unmixing in multi-spectral optoacoustic tomography,” Opt. Express 19(4), 3175–3184 (2011). [CrossRef] [PubMed]

11. N. C. Deliolanis, A. Ale, S. Morscher, N. C. Burton, K. Schaefer, K. Radrich, D. Razansky, and V. Ntziachristos, “Deep-tissue reporter-gene imaging with fluorescence and optoacoustic tomography: A performance overview,” Mol. Imaging Biol. 16(5), 652–660 (2014). [CrossRef] [PubMed]

12. A. Buehler, E. Herzog, A. Ale, B. D. Smith, V. Ntziachristos, and D. Razansky, “High resolution tumor targeting in living mice by means of multispectral optoacoustic tomography,” EJNMMI Research 2(1), 1–6 (2012). [CrossRef]

13. A. Taruttis, M. Wildgruber, K. Kosanke, N. Beziere, K. Licha, R. Haag, M. Aichler, A. Walch, E. Rummeny, and V. Ntziachristos, “Multispectral optoacoustic tomography of myocardial infarction,” Photoacoustics 1(1), 3–8 (2013). [CrossRef] [PubMed]

14. A. Hyvarinen, J. Karhunen, and E. Oja, Independent Component Analysis (John Wiley & Sons, 2004), Chap. 2.

15. B. W. Silverman, Density Estimation for Statistics and Data Analysis (Chapman and Hall1986). [CrossRef]

16. J. Wang and C. I. Chang, “Applications of independent component analysis in endmember extraction and abundance quantification for hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens. 44(9), 2601–2616 (2006). [CrossRef]

17. V. D. Calhoun, J. Liu, and T. Adal, “A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and ERP data,” NeuroImage 45(1, Supplement 1), S163–S172 (2009). [CrossRef]

18. C. Panagiotou, S. Somayajula, A. P. Gibson, M. Schweiger, R. M. Leahy, and S. R. Arridge, “Information theoretic regularization in diffuse optical tomography,” J. Opt. Soc. Am. 26(5), 1277–1290 (2009). [CrossRef]

19. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58, R37–R61 (2013). [CrossRef] [PubMed]

20. L. Kou, D. Labrie, and P. Chylek, “Refractive indices of water and ice in the 0.65-to 2.5-μm spectral range,” Appl. Opt. 32(19), 3531–3540 (1993). [CrossRef] [PubMed]

21. H. J. Van Staveren, C. J. Moes, J. van Marie, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in lntralipid-10% in the wavelength range of 400–1100nm,” Appl. Opt. 30(31), 4507–4514 (1991). [CrossRef] [PubMed]

22. L. An, T. Saratoon, M. Fonseca, R. Ellwood, and B. Cox, “Exploiting statistical independence for quantitative photoacoustic tomography,” Proc. SPIE 10064, 1006419 (2017). [CrossRef]

23. M. Schweiger and S. Arridge, “The Toast++ software suite for forward and inverse modeling in optical tomography,” J. Biomed. Opt. 19(4), 040801 (2014). [CrossRef] [PubMed]

24. R. Ellwood, O. Ogunlade, E. Z. Zhang, P. C. Beard, and B. T. Cox, “Orthogonal Fabry-Pérot sensors for photoacoustic tomography,” Proc. SPIE 9708, 97082N (2016).

25. R. Ellwood, O. Ogunlade, E. Zhang, P. Beard, and B. Cox, “Photoacoustic tomography using orthogonal Fabry-Pérot sensors,” J. Biomed. Opt. 22(4), 041009 (2017). [CrossRef]

26. E. Zhang, J. Laufer, and P. Beard, “Backward-mode multiwavelength photoacoustic scanner using a planar Fabry-Perot polymer film ultrasound sensor for high-resolution threedimensional imaging of biological tissues,” Appl. Opt. 47, 561–577 (2008). [CrossRef] [PubMed]

27. G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am. 112(4), 1536–1544 (2002). [CrossRef] [PubMed]

28. P. Stefanov and G. Uhlmann, “Thermoacoustic tomography with variable sound speed,” Inverse Probl. 25(7), 075011 (2009). [CrossRef]

29. J. Nocedal and S. Wright, Numerical Optimization (Springer Science & Business Media, 2006).

30. S. Shwartz, M. Zibulevsky, and Y. Y. Schechner, “ICA using kernel entropy estimation with NlogN complexity,” in Independent Component Analysis and Blind Signal Separation: Fifth International Conference, ICA 2004, Granada, Spain, September 22–24, 2004. Proceedings, C. G. Puntonet and A. Prieto, eds. (Springer, 2004), pp. 422–429. [CrossRef]

31. J. Laufer, E. Zhang, and P. Beard, “Evaluation of absorbing chromophores used in tissue phantoms for quantitative photoacoustic spectroscopy and imaging,” IEEE J. Sel. Top. Quantum Electron. 16(3), 600–607 (2010). [CrossRef]

32. M. Fonseca, E. Malone, F. Lucka, R. Ellwood, L. An, R. Arridge, P. Beard, and B. Cox, “Three-dimensional photoacoustic imaging and inversion for accurate quantification of chromophore distributions,” Proc. SPIE 10064, Photons Plus Ultrasound: Imaging and Sensing 2017, 1006415 (2017). [CrossRef]

33. S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star, and M. J. van Gement, “Optical properties of intralipid: A phantom medium for light propagation studies,” Lasers. Surg. Med. 12(5), 510–519 (1992). [CrossRef] [PubMed]

34. C. R. Vogel, Computational Methods for Inverse Problems (Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA: 2002). Chap. 7. [CrossRef]

35. A. M. Thompson, J. C. Brown, J. W. Kay, and D. M. Titterington, “A study of methods of choosing the smoothing parameter in image restoration by regularization,” IEEE Trans. Pattern Anal. Mach. Intell. 13(13),326–339 (1991). [CrossRef]

36. T. Saratoon, T. Tarvainen, B. T. Cox, and S. R. Arridge, “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Invserse Prob. 29(7), 075006 (2013). [CrossRef]

37. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D. 60(1–4), 259–268 (1992). [CrossRef]

38. C. Huang, K. Wang, L. Nie, L. Wang, and M. Anastasio, “Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media,” IEEE Trans. Med. Imaging. 32(6), 1097–1110 (2013). [CrossRef] [PubMed]

39. S. Arridge, P. Beard, M. Betcke, B. Cox, N. Huynh, F. Lucka, O. Ogunlade, and E. Zhang, “Accelerated high-resolution photoacoustic tomography via compressed sensing,” Phys. Med. Biol. 61(24), 8908 (2016). [CrossRef] [PubMed]

NiCl₂	ε_d	ε_d+MI	Expected
Tube 1	860	698	345
Tube 2	185	−5	0
Tube 3	367	258	352
Tube 4	167	14	0

Statistical independence in nonlinear model-based inversion for quantitative photoacoustic tomography

Abstract

1. Introduction

2. Mutual information error functional

3. Model-based inversion with statistical independence

4. Generating multiwavelength photoacoustic images of tissue mimicking phantoms

4.1. Experimentally acquired images

4.2. Numerically simulated images

5. Inverting for the chromophore concentrations

5.1. Effect of model-mismatch

5.2. Experimental results

6. Discussion

7. Conclusion

Funding

Disclosures

References and links

Cited By

Figures (6)

Tables (1)

Equations (24)

Biomedical Optics Express