Quantitative photoacoustic tomography augmented with surface light measurements.

Quantitative photoacoustic tomography is an imaging modality in which distributions of optical parameters inside tissue are estimated from photoacoustic images. This optical parameter estimation is an ill-posed problem and it needs to be approached in the framework of inverse problems. In this work, utilising surface light measurements in quantitative photoacoustic tomography is studied. Estimation of absorption and scattering is formulated as a minimisation problem utilising both internal quantitative photoacoustic data and surface light data. The image reconstruction problem is studied with two-dimensional numerical simulations in various imaging situations using the diffusion approximation as the model for light propagation. The results show that quantitative photoacoustic tomography augmented with surface light data can improve both absorption and scattering estimates when compared to the conventional quantitative photoacoustic tomography.


Introduction
Photoacoustic tomography (PAT) is an imaging modality based on the photoacoustic effect generated through the absorption of an externally introduced light pulse. The method combines optical contrast with high spatial resolution of ultrasound. The optical contrast is provided by distinctive absorption spectra by different chromophores. The chromophores of interest are, for example, haemoglobin, melanin and optionally various contrast agents. The ultrasonic waves carry this optical information directly to the surface with minimal scattering, thus retaining accurate spatial information as well. Nowadays, PAT can be used to provide images of soft biological tissues with high spatial resolution. It has successfully been applied to the visualisation of different structures in biological tissues, such as human blood vessels, microvasculature of tumors, and the cerebral cortex in small animals. For more information about PAT, see e.g. the reviews by [1][2][3][4][5][6][7][8] and the references therein.
Quantitative photoacoustic tomography (QPAT) is a technique which aims at estimating the absolute concentrations of the chromophores from photoacoustic images, i.e. from the reconstructed initial pressure distribution [9]. This is an ill-posed problem which needs to be approached in the framework of inverse problems. The concentrations of the chromophores can be estimated either by directly estimating them from photoacoustic images obtained at various wavelengths [10][11][12][13][14][15][16] or by first recovering the absorption coefficients at different wavelengths and then calculating the concentrations from the absorption spectra [9,10,13,16]. Different approaches for the solution of the optical inverse problem of QPAT have been considered, see e.g. [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35]. As an alternative to the conventional two stage approach, estimation of the optical parameters directly from the photoacoustic time-series has also been considered recently [34,[36][37][38][39][40][41]. In this work, the two stage approach is taken and it is assumed that the acoustic inverse problem, i.e. estimation of the initial pressure, has been solved. Further, one wavelength, i.e. estimation of absorption and scattering, is considered. Extensions of the developed numerical approach to one stage approach and spectral QPAT are straightforward.
It has been shown that, in order to obtain accurate estimates for absorption in QPAT, the scattering effects need to be taken into account [24,25,29,42]. Estimation of scattering is more ill-posed than absorption, and thus more sensitive to errors in data and modelling [25,42]. Furthermore, it has been shown that false scattering values can lead to errors in absorption estimates [29].
In this work, improving QPAT by including information from surface measurements of light is investigated. Surface light measurements have been previously utilised in QPAT using the following approaches. In [43,44], a two-step approach was suggested in which scattering dis-tribution was first solved using diffuse optical tomography (DOT) measurements, and then this information was utilised in the estimation of the absorbed optical energy density in photoacoustic imaging. Furthermore, in [45], also utilising a two-step approach, a DOT experiment was first used to determine constant values for absorption and scattering, and then these were later utilised as background values in QPAT image reconstruction. In [46], a hybrid approach was introduced. In the approach, the vector field method developed in [24] was generalised and utilised in the estimation of the boundary parameters from surface light and the optical parameters inside the domain from absorbed optical energy.
In this work, estimating optical parameters using both absorbed optical energy density and surface light measurements is considered. In the approach, these two data sources are utilised simultaneously in order to solve the inverse problem of QPAT in a Bayesian framework. The approach is investigated in various imaging situations including differently sized domains and various internal structures of the optical parameters. Furthermore, in this paper, simultaneous estimation of the optical parameters and the Grüneisen coefficient is investigated. Simultaneous estimation of the absorption, scattering and Grüneisen parameters is non-unique if only one wavelength of light is used to obtain the QPAT data [13]. In this work, enhancing QPAT with surface light measurements is suggested to overcome this problem.
In QPAT, the most common approach has been to use the diffusion approximation (DA) as the light transport model, although the radiative transfer equation (RTE) has also been utilised [15,21,25,38,47,48]. Recently, utilising Monte Carlo simulation methods have also been considered [35,49]. In this work, the diffusion approximation (DA) is used as the forward model for light propagation. The DA is a valid approximation in a highly scattering medium further than a few scattering lengths from the light source, and thus it can be regarded as a relatively safe approximation for the simulation cases considered in this work.
The rest of the paper is organised as follows. The forward and inverse problem of QPAT are described in Section 2 and the related numerical implementations are given in Section 3. The results of simulations are shown in Section 4 and the conclusions are given in Section 5.

Forward model
In quantitative photoacoustic tomography, a short pulse of laser light is used to illuminate the tissue region of interest. The propagation of light can be described using the radiative transfer equation (RTE) [50]. In biomedical imaging, the RTE is often approximated with the diffusion approximation (DA) where Φ(r , t) [W/mm 2 ] is the photon fluence at time instance t, c is the speed of light in the medium, µ a (r) is a diffuse boundary current at the source position ǫ j ⊂ ∂Ω, ∂Ω is the boundary of the domain Ω, γ d is a dimension-dependent constant which takes values γ 2 = 1/π and γ 3 = 1/4 andn is an outward unit normal. The DA is a valid approximation in situations in which the radiance is almost an uniform distribution, i.e. in a scattering dominated medium further than a few scattering lengths from the light source [50].
As light propagates within the tissue, it is absorbed by chromophores. This generates localised increases in pressure. This pressure increase propagates through the tissue as an acoustic wave and is detected by ultrasound sensors at the surface of the tissue. The propagation of the acoustic wave occurs on a microsecond time scale, about five orders of magnitude slower than the optical propagation, so only the total absorbed optical energy density is of interest and not the rate of the absorption. Thus, in QPAT, light propagation can be modelled using a timeindependent (continuous wave) model of light transport. The time-independent form of the DA is of the form is the time-independent fluence and I s (r) [J/mm 2 ] is a time-independent boundary light source. Furthermore, the absorbed optical energy density and the initial acoustic pressure p 0 (r) [Pa] is [9] p 0 (r) = p(r , t = 0) = G(r)H(r) where G(r) [unitless] is the Grüneisen parameter which is used to identify the photoacoustic efficiency. The time evolution of the resulting photoacoustic wave fields can be modelled using the equations of linear acoustics [51].

Measuring surface light
In this work, utilising surface light measurements in QPAT is investigated. Basically, this corresponds on using both QPAT and diffuse optical tomography data in determining the optical parameters. In a typical DOT measurement setup, the target is illuminated using either a short pulse of light (time-domain measurement systems), intensity modulated light (frequencydomain systems) or continuous light [52]. The measurable quantity is exitance Γ + (r , · ) on the boundary of the target Depending on the measurement system, the exitance can be time-dependent Γ + (r , t) [W/mm 2 ], frequency-dependent Γ + (r , ω) [J/mm 2 ] or intensity only Γ + (r) [J/mm 2 ]. The time-domain and frequency-domain data are related through Fourier-transform. In addition, other moments can be calculated from time-dependent measurements [53]. In this work, we consider frequencydomain data. In practice, these can be obtained either by a frequency-domain measurement device or by measuring the time-response of a time-domain system and Fourier-transforming the data into the frequency domain. Similarly, in order to solve the modelled exitance, the timedomain DA (1)-(2) or its counterpart in frequency domain where ω is the angular modulation frequency of the input light and i is the imaginary unit, is used.

Inverse problem of QPAT
In this work, the optical inverse problem of QPAT, i.e. estimation of distributions of the optical parameters from photoacoustic images, is considered. The problem is approached in the framework of Bayesian inverse problems [28,54]. A discrete observation model for QPAT in the presence of additive noise is where H meas ∈ R m qpat is a measurement vector where m qpat is the number of measurements which in this case is the number of illuminations multiplied with the number of discretised elements to represent the data space, µ a = (µ a 1 , . . . , µ a K ) T ∈ R K and µ ′ s = (µ ′ s 1 , . . . , µ ′ s K ) T ∈ R K are discretised absorption and scattering coefficients, K is the number of discretised parameters, H is the discretised forward operator corresponding to model (3)-(5) and e qpat ∈ R m qpat denotes the noise.
Let us assume that µ a , µ ′ s and H meas are random variables. The solution of the inverse problem is the posterior probability density which is obtained through Bayesian formula and can be written in the form where π(H meas |µ a , µ ′ s ) is the likelihood density and π(µ a , µ ′ s ) is the prior density. Assuming that the unknowns µ a and µ ′ s and noise e are mutually independent and Gaussian distributed, i.e.
, where η µ a , η µ ′ s and η e ,qpat are the means and Γ µ a , Γ µ ′ s and Γ e ,qpat are the covariances for the absorption, scattering and noise, respectively, the posterior density (11) becomes The practical solution for the inverse problem is obtained by calculating point estimates from the posterior density. Since we are interested in computationally efficient inverse problem solvers, we consider here the maximum a posteriori (MAP) estimate. It is obtained as where L e ,qpat is the Cholesky decomposition of the inverse of the noise covariance matrix Γ −1 e ,qpat = L T e ,qpat L e ,qpat , and L µ a and L µ ′ s are the Cholesky decompositions of the inverse of the prior covariance matrices for absorption and scattering, Thus, in the image reconstruction problem of QPAT, we seek to find the discretised distributions of absorption and scattering coefficients (μ a ,μ ′ s ) which solve the minimisation problem (13).

QPAT augmented with surface light data
A discrete observation model for DOT in the presence of additive noise is where Γ + meas ∈ R m dot is a measurement vector of light exitance measured at the detectors on the surface of the target where m dot is the number of measurement, Γ + is the discretised forward operator corresponding to model (7)-(9) and e dot ∈ R m dot denotes the noise. Following a similar derivation as for QPAT, the inverse problem of QPAT augmented with surface light measurements can formulated as a minimisation problem where η e ,dot is the mean of the noise of the surface light measurements and L e ,dot is the Cholesky decomposition of the inverse of the noise covariance matrix of the surface light measurements

Simultaneous estimation of the Grüneisen parameter
Generally, estimation of the absorption, scattering and Grüneisen parameter is non-unique if only one wavelength of light is used [24]. This has been overcome by using multiple optical wavelengths and estimating the optical parameters and the Grüneisen coefficient simultaneously [13,15,16]. In this work, simultaneous estimation of the absorption, scattering and Grüneisen parameter is considered using initial pressure and surface light as data. The minimisation problem is of the form where G = (G 1 , . . . , G K ) T ∈ R K is the discretised distribution of the Grüneisen parameter which was assumed to be Gaussian distributed and independent of the noise, p 0,meas is the initial pressure obtained from measurements, p 0 (µ a , µ ′ s , G) is the modelled initial pressure, and η e ,p is the mean and L e ,p is the Cholesky decomposition of the inverse of the noise covariance matrix Γ −1 e ,p = L T e ,p L e ,p of the initial pressure data. Further, η G is the mean and L G is the Cholesky decomposition of the inverse of the covariance matrix Γ −1 G = L T G L G of the prior for the Grüneisen parameter.

FE-approximation of the DA
In this work, a finite element method (FEM) is used for the numerical solution of the DA. In order to obtain the FE-approximation, first a variational problem is formulated of the DA with its boundary condition, and then the solution of the variational problem is approximated in a piece-wise linear basis. As a result, a matrix equation can be derived. For more details, see e.g. [25,55].
In the case of QPAT, the FE-approximation of the time-independent DA with the diffuse boundary source model, Eqs. (3)-(4), is obtained by solving the matrix equation where Φ qpat is the fluence in the nodes of the FE-mesh and A qpat = K + C + R. In the case of the frequency-domain DA, Eqs. (8)-(9), the FE-approximation of the DA on an angular modulation frequency ω is obtained by solving the matrix equation where Φ dot is the fluence in the nodes of the FE-mesh and A dot = K + C + R + Z. The matrices K, C, R and Z and the source vectors b qpat and b dot are of the form where ϕ q (r) and ϕ p (r) are FE-basis functions, q, p = 1, . . . , N , and N is the number of nodes. The QPAT data, absorbed optical energy density, H can be computed from the fluence Φ qpat obtained with (17) using Eq. (5) and the surface light measurements Γ + can be computed from the fluence Φ dot obtained with (18) using Eq. (7).

Gauss-Newton iteration
In this work, the minimisation problems (13), (15) and (16) are solved using a Gauss-Newton method. In the case of QPAT image reconstruction problem augmented with surface light measurements, Eq. (15), the Gauss-Newton iteration can be written in a form where x = (µ a , µ ′ s ) T = (µ a 1 , . . . , µ a K , µ ′ s 1 , . . . , µ ′ s K ) T ∈ R 2K are the estimated absorption and scattering parameters which in this work are represented in piece-wise constant bases µ a (r) is a characteristic function of the element Ω k . It should be noted that other bases can also be used for the representation of the optical parameters and that different discretisations can be used for different parameters. Furthermore, F meas = H meas , Γ + meas T are the measured absorbed optical energy density and exitance, F = H , Γ + T is the solution of the discretised forward models for QPAT, Eqs. (17) and (5), and DOT, Eqs. (18) and (7), s is the step length of the minimisation algorithm determined using a line-search method, and Further, J is the Jacobian of the form where Jacobians for absorption and scattering are J where the columns k = 1, . . . , K of each matrix are obtained by vectorisation of and M qpat ∈ R m qpat × N is a measurement matrix which contains the discretised measurement operator for QPAT which maps the piece-wise linear solution of the QPAT forward model to piece-wise constant data space and M dot ∈ R m dot × N is a measurement matrix which contains the discretised measurement operator for surface light which maps the piece-wise linear solution of the DOT forward model to detector measurements. The Gauss-Newton algorithm for the solution of the conventional QPAT problem (13), can be formulated from (25) by dropping the DOT related matrices and vectors from the forward model F, measurement vector F meas , noise statistics L e , η e and the Jacobian matrix J.
In the case of estimating the absorption, scattering and Grüneisen parameter simultaneously from the initial pressure and surface light data, i.e. the minimisation problem (16), the Gauss-Newton iteration is of the form where wherex = (µ a , µ ′ s , G) T = (µ a 1 , . . . , µ a K , µ ′ s 1 , . . . , µ ′ s K , G 1 , . . . , G K ) T ∈ R 3K are the estimated absorption, scattering and Grüneisen parameters,F meas = p 0,meas , Γ + meas T are the measured initial pressure and exitance,F = p 0 , Γ + T is the solution of the discretised forward models and Further,J is the Jacobian of the form where Jacobians for the absorption, scattering and Grüneisen parameter areJ qpat µ a = j qpat (1) and J dot µ a and J dot µ ′ s are as in (28) and (30).

Results
The proposed approach was tested with simulations which were implemented in twodimensional rectangular domains. Estimating absorption and scattering was investigated with varying domain size and varying parameter distributions using absorbed optical energy density and exitance as data. Furthermore, simultaneous estimation of the absorption, scattering and Grüneisen parameter was investigated using initial pressure and exitance as data. The results were compared to conventional QPAT reconstructions. The simulations were performed with a Windows workstation with Intel(R) processor with 2 cores, a speed of 3.16 GHz and 8 GB RAM using matlab (R2014a, The MathWorks Inc. Natick, Massachusetts, USA).

Data simulation
In all of the simulations, four planar illuminations, one at each side of the rectangle, with a uniform intensity covering the whole side length and a total energy of 1 J, were used. It should be noted that, since the numerical simulations are performed with a noise amplitude related to the data, the absolute magnitude of the input light energy does not affect the results of these simulations. The absorbed optical energy density, i.e. QPAT data, was simulated by using the FE-solution of the DA (17) to obtain the fluence and then multiplying that with the absorption to get the absorbed optical energy density using Eq. (5). To get the initial pressure distribution, the absorbed optical energy density was multiplied with the Grüneisen parameter using Eq. (6). In all simulations, Gaussian distributed noise with zero mean and standard deviation of 1 % of the maximum absorbed optical energy density or initial pressure was added to the simulated data. Furthermore, the boundary light data was simulated by first solving the FE-approximation of the frequency domain DA (18) for fluence, and then calculating the exitance using Eq. (7). The exitance was solved in 174 equally distributed detector points (58 at each side) on the object boundary excluding the illumination side of the target. In the simulations, the angular modulation frequency of ω = 100 · 10 6 rad/s was used. Gaussian distributed noise with standard deviation of 1 % of the amplitude and phase was added to the simulated data.
The number of nodes and elements of the FE-discretisations used in the simulations are given in Table 1. The optical parameters were represented in piece-wise constant bases using the elements of the FE-discretisation.

Reconstructing the parameters of interest
The absorption and scattering distributions were reconstructed using the suggested augmented QPAT approach by minimising (15). The minimisation problem was solved using the Gauss-Newton method (25) equipped with a line search algorithm for the determination of the step length. The results were compared to the conventional QPAT reconstructions by minimising (13). Furthermore, the absorption, scattering and Grüneisen parameter distributions were reconstructed using the suggested augmented QPAT approach by minimising (16). The minimisation problem was solved using the Gauss-Newton method (33) equipped with a line search algorithm for the determination of the step length. All minimisation problems converged in less than 15 iterations. Typically, the augmented QPAT took few more simulations to converge than the conventional QPAT. The number of nodes and elements of the FE-discretisations used in the reconstructions is given in Table 2. The optical parameters were represented in piece-wise constant bases by the elements of the FE-discretisations. In this work, the prior model for the unknown parameters µ a , µ ′ s and Γ was chosen to be based on the Ornstein-Uhlenbeck process [16,56]. The prior is a Gaussian distribution with mean η and covariance matrix Γ which was defined as being proportional to where σ is the standard deviation of the prior and Ξ is a matrix which has its elements defined as where i and j denote the row and column indices of the matrix, r i and r j denote the coordinates of the centre of elements i and j, and ξ is the characteristic length scale of the prior describing the spatial distance that the parameter is expected to have (significant) spatial correlation for. In this work, the mean and standard deviation of the prior were chosen based on expected mean and variation that the parameters were assumed to have, and the characteristic length scale was chosen based on size of the structures that the medium was assumed to have. The mean, standard deviation and characteristic length scale used in the simulations are given in Table 3.
The validity of the reconstructions was evaluated by computing the mean relative errors for Table 3. The mean and standard deviation of the prior distributions for the absorption η µ a (mm −1 ), σ µ a (mm −1 ), scattering η µ ′ s (mm −1 ), σ µ ′ s (mm −1 ) and Grüneisen parameter η G (unitless), σ G (unitless), and the characteristic length scale ξ(mm) used in the simulations.
where µ TRUE a , µ ′ s TRUE and G TRUE are the simulated parameters interpolated to the reconstruction grid andμ a ,μ ′ s andĜ are the estimated parameters.

Varying domain size
The performance of the augmented QPAT in different size domains was investigated. The sizes of the simulation domains were 20 mm × 20 mm, 40 mm × 40 mm and 60 mm × 60 mm. The absorption and scattering inclusions were stripe-like structures with thickness of 2 mm and length of the side-length of the domain minus 1 mm. The stripes were located with a distance of 3 mm from each other. The simulated parameter distributions and the reconstructions are shown in Fig. 1. The relative errors for the estimated absorption and scattering, Eqs. (40)-(41), are given in Table 4.
As it can been seen, the absorption distributions reconstructed using the augmented QPAT and the conventional QPAT approaches are almost the same quality in the 20 mm × 20 mm size domain. The scattering estimates, on the other hand, obtained with the augmented QPAT resemble the original distribution more than the estimates obtained with the conventional QPAT. The difference between the two approaches becomes more apparent when the domain size increases. Especially in the domain of size 40 mm × 40 mm, the reconstructions are clearly better quality when the surface light data has been utilised in the reconstructions. In the largest domain, 60 mm × 60 mm, the quality of both conventional and augmented QPAT is weaker than in the smaller domains. That is, one is reaching the limit where QPAT reconstructions cannot significantly be improved by including exitance data.
The calculated mean relative errors, Table 4, support these results. The relative errors of the absorption and scattering estimates obtained with the augmented QPAT are smaller than those of the conventional QPAT reconstructions. The most accurate estimates are obtained in the 20 mm × 20 mm domain and as the domain size increases the relative errors of the estimates increase as well. The most significant improvement into the accuracy of the estimates by the augmented QPAT is obtained in 40 mm × 40 mm for absorption estimates.  The convergence of the minimisation problems was visualised by calculating the mean relative errors for absorption and scattering, Eqs. (40)-(41), at each iteration and plotting the results. All the simulations had a similar behaviour, and thus only the results of once case (simulations in 20 mm × 20 mm domain) are shown in Fig. 2. As it can be seen, the augment QPAT required more iterations to converge than the conventional QPAT but it converged to more accurate estimates for absorption and scattering.

Variations in the optical parameter distribution
The capability of the augmented QPAT approach to reconstruct fine structures in internal optical parameter distributions was investigated. The simulation domain size was 20 mm×20 mm. In all simulations, the absorption and scattering were given two different values but their distribution was varied from coarse to fine structures. The total area of absorbing and scattering inclusions  Table 4.
As it can be seen, QPAT augmented with surface light measurements gives better quality reconstructions both for absorption and scattering. In the absorption reconstructions, the shapes of the absorption inclusions are well defined by both methods. However, it seems that the augmented QPAT gives more accurate values for the absorbing inclusions. In the case of scattering, the inclusions are better distinguished if augmented QPAT data is used. The relative errors of the absorption and scattering estimates obtained with the augmented QPAT are smaller than those obtained with the conventional QPAT. That is, the augmented QPAT provides more accurate quantitative estimates for the optical parameters.

Simultaneous estimation of the Grüneisen parameter
The simultaneous estimation of the absorption, scattering and Grüneisen parameter was investigated in 10 mm × 10 mm domain. For comparison, two types of conventional QPAT reconstructions, Eq. (13), were computed. In the first one, the Grüneisen parameter was assumed to have a constant value G = 0.3, which is the mean of the simulated Grüneisen parameter distribution. In the second one, the simulated Grüneisen parameter distribution was mapped to the reconstruction mesh. This can be considered the best possible reference for the reconstructions, which however is unrealistic since, in practice, the Grüneisen parameter distribution is unknown. The simulated and reconstructed parameter distributions are shown in Fig. 4. The mean relative errors, Eqs. (40)-(42), are given in Table 5.
As it can be seen, the approach can produce as good quality reconstructions for absorption and scattering as the reference approach with the correct Grüneisen parameter distribution and better than the reference approach with the constant Grüneisen parameter. The relative errors for the absorption are larger and the relative errors for the scattering are smaller when compared to the the reference approach with the correct Grüneisen parameter distribution. Further, the relative errors are smaller when compared with the reference approach with the constant Grüneisen parameter. Thus, the simulations show that simultaneous estimation of the absorption, scattering and Grüneisen parameter may be possible using the QPAT augmented with surface light measurements.

Conclusions
In this work, improving QPAT by combining information from surface measurements of light was investigated. The QPAT image reconstruction problem was formulated as a minimisation problem in which absorption and scattering distributions were reconstructed utilising both QPAT and surface light data. Furthermore, simultaneous estimation of the Grüneisen parameter was investigated. The approach was tested with simulations.
The results show that including surface light measurements into the QPAT image reconstruction improves both quality of the reconstructed absorption and scattering images as well as their quantitative estimates. The results are in line with other approaches in which surface light data has been utilised in QPAT [46]. Furthermore, the simulations demonstrate some of the situations in which the augmented QPAT is most useful when compared to the conventional QPAT, i.e. in distinguishing scattering, in larger domains, and when the medium consists of fine structures.
The augmented QPAT was also tested in simultaneous estimation of the optical parameters and the Grüneisen parameter. The results show that simultaneous estimation of all parameters results into better estimates for absorption and scattering than the conventional QPAT approach where the Grüneisen parameter was assumed to be constant. This indicates that it may be possible to estimate the Grüneisen parameter simultaneously with the optical parameters if QPAT is augmented with surface light measurements. However, it should be noted that the uniqueness of this minimisation problem remains to be studied.
In this work, the augmented QPAT was tested with numerical simulations in a twodimensional setting using the diffusion approximation as the light transport model at one optical wavelength. The future work includes implementing the methodology using more realistic and numerically more challenging radiative transfer equation as well as extension of the methodology to spectral QPAT.

Funding
Academy of Finland (projects 286247 and 250215 Finnish Centre of Excellence in Inverse Problems Research), Magnus Ehrnrooth foundation, and Jane and Aatos Erkko foundation.

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