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

The statistical independence between the distributions of different chromophores in tissue has previously been used for linear unmixing with independent component analysis (ICA). In this study, we propose exploiting this statistical property in a nonlinear model-based inversion method. The aim is to reduce the sensitivity of the inversion scheme to errors in the modelling of the fluence, and hence provide more accurate quantification of the concentration of independent chromophores. A gradient-based optimisation algorithm is used to minimise the error functional, which includes a term representing the mutual information between the chromophores in addition to the standard least-squares data error. Both numerical simulations and an experimental phantom study are conducted to demonstrate that, in the presence of experimental errors in the fluence model, the proposed inversion method results in more accurate estimation of the concentrations of independent chromophores compared to the standard model-based inversion.


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][4][5][6][7][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.

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 1 and y 2 , are statistically independent if their joint probability density function (PDF), ρ y 1 ,y 2 (y 1 , y 2 ), is the product of their marginal PDFs, ρ y 1 (y 1 ) and ρ y 2 (y 2 ) [14], such that where y 1 and y 2 denote possible values of y 1 and y 2 . 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 1 and y 2 is given by where H(y k ) and H(y 1 , y 2 ) are the entropy and the joint entropy of y 1 and y 2 respectively. They are defined by where k = 1 or 2, and ρ y 1 ,y 2 (y 1 , y 2 )logρ y 1 ,y 2 (y 1 , y 2 )dy 1 dy 2 .
Similarly, the MI between multiple random variables is defined by H(y k ) − H(y 1 , ..., y K ).
As indicated in Eqs. (2)(3)(4), the MI depends on the probabilities of the variables. We consider the concentrations of K independent chromophores c 1 , ..., c K as continuous random variables with the PDFs ρ c 1 , ..., ρ 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ρ andρ where ξ k are the values at which the PDF is evaluated and κ(x) is a kernel function that satisfies κ(x) ≥ 0 and ∫ ∞ −∞ κ(x)dx = 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: 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 = 4 3M 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 ∆ξ k (11) where the summation symbol denotes a multiple sum of q 1 , ..., q K and ∆ξ k is the spacing between the sampling points for the PDF. Using the approximated entropies, the MI can be estimated by 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] ∂Ȋ The first term in the right hand side of Eq. (13) is where by the product rule, and where κ denotes the derivative of κ, given by for the Gaussian kernel. Similarly, the second term in Eq. (13) is where and

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 = 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 model m,λ n at voxel m and wavelength λ n is calculated based on 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 meas m,λ n , and hence the minimisation problem is given by 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 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.

Generating multiwavelength photoacoustic images of tissue mimicking phantoms
The accuracy of the quantification using ε d+M I 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 2 and NiCl 2 . For both the numerical and the experimental phantom, the distributions of CuCl 2 and NiCl 2 are arranged such that they are statistically independent of each other.

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 −1 NiCl 2 and the second and the fourth tube contain 36gL −1 CuCl 2 . The absorption spectra of the chromophores and scattering spectrum of Intralipid are shown in Fig. 2. The CuCl 2 and NiCl 2 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 2 and NiCl 2 are both also statistically independent from water, even though they can be found at the same locations.) 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 13x13mm 2 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.   [21]. A spectrophotometer (Lambda 750S, Perkin Elmer) was used to measure the transmittance of CuCl 2 , NiCl 2 and India ink in order to determine their absorption spectra. Reprinted from [22].
The iterative time reversal method [27,28] was used to reconstruct 3D images from the photoacoustic time series. A 2D cross-sectional 12x12mm 2 region of interest centred at the tubes was used for the optical inversion, and the dimension was reduced to 72x72 pixels to reduce the computational time and memory requirements. The 2D slices are shown for three wavelengths in Fig. 3.

Numerically simulated images
The numerical phantom has an element spacing of 100µm and represents a 5x5mm 2 area with six insertions arranged in two columns, as illustrated in Fig. 4. The left column of insertions represents solutions of CuCl 2 with concentrations 12, 24 and 36gL −1 in increasing order, where  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 −1 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). 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.

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+M I 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 2 , NiCl 2 , 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+M I . For all inversions using ε d+M I , the MI was calculated between CuCl 2 and NiCl 2 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+M I 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.

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 2 and NiCl 2 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+M I 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.
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+M I 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.

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 2 and NiCl 2 are both known to increase with concentration and are given by where Γ H 2 O is the Grüneisen parameter for water, and β i are 5.80×10 −3 Lg −1 and 2.25×10 −3 Lg −1 for CuCl 2 or NiCl 2 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 2 (top row) and NiCl 2 (bottom row) using ε d and ε d+M I 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 −1 . 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 2 in the second tube and NiCl 2 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 2 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 2 and NiCl 2 are almost completely removed when ε d+M I is used. The absolute concentrations of the CuCl 2 is estimated accurately with an error of 3gL −1 on average for the four tubes. The overestimation of the NiCl 2 concentration in the top tube remains present with ε d+M I . However, this overestimation error is reduced when ε d+M I is used compared to ε d .

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+M I , 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 model m,λ n , 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 model m,λ n . Therefore, the quantification errors were greatly reduced when ε d+M I 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+M I , 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][38][39].

Conclusion
We proposed exploiting the statistical independence between certain chromophores in the modelbased 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.