# Mitigating the unicity problem in quantitative photoacoustic imaging using ultrasound blood flow imaging, a simulation study

Charlotte Constans, Ivana Falco, Bastien Arnal

► To cite this version:

Charlotte Constans, Ivana Falco, Bastien Arnal. Mitigating the unicity problem in quantitative photoacoustic imaging using ultrasound blood flow imaging, a simulation study. Physics in Medicine and Biology, 2025, 70 (15), pp.155003. 10.1088/1361-6560/adebda. hal-05184438

HAL Id: hal-05184438

https://hal.science/hal-05184438v1

Submitted on 29 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER

# Mitigating the unicity problem in quantitative photoacoustic imaging using ultrasound blood flow imaging, a simulation study

To cite this article: Charlotte Constans et al 2025 *Phys. Med. Biol.* 70 155003

View the [article online](#) for updates and enhancements.

You may also like

- [Photoacoustic imaging in percutaneous radiofrequency ablation: device guidance and ablation visualization](#)
Kaloor Joseph Francis and Srirang Manohar

- [Transcranial photoacoustic characterization of neurovascular physiology during early-stage photothrombotic stroke in neonatal piglets in vivo](#)
Jeeun Kang, Xiuyun Liu, Suyi Cao et al.

- [Photoacoustic tomography: applications for atherosclerosis imaging](#)
Gurneet S Sangha and Craig J Goergen

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

CrossMark

RECEIVED

23 January 2025

REVISED

30 June 2025

ACCEPTED FOR PUBLICATION

3 July 2025

PUBLISHED

24 July 2025

PAPER

# Mitigating the unicity problem in quantitative photoacoustic imaging using ultrasound blood flow imaging, a simulation study

Charlotte Constans, Ivana Falco and Bastien Arnal*

<figure><img src="image_4.png" /></figure>

Univ. Grenoble Alpes, CNRS, LIPhy, 38058 Grenoble, France

* Author to whom any correspondence should be addressed.

E-mail: bastien.arnal@univ-grenoble-alpes.fr

**Keywords:** inverse problems, quantitative photoacoustic imaging, blood flow imaging, ultrasound power Doppler, Monte Carlo radiance

## Abstract

**Objective.** Quantitative photoacoustic (PA) imaging aims at recovering optical properties in optically turbid media. Here, we aim to recover uniqueness of the solution in absorption and scattering coefficient in a single illumination and wavelength setting. **Approach.** By considering a numerical tissue model composed of optically contrasted regions surrounded by blood vessels, we combine PA imaging with ultrasound power Doppler (UPD). We show that quantitative absorption coefficients inside the blood vessels can be derived from UPD, serving as a prior for a radiance Monte Carlo inverse model. **Main results.** We show that using this vascular imaging prior mitigates the non-uniqueness issue, providing quantitative reconstructions of absorption and scattering even in the absence of background or blood absorption prior information. We generalize to unresolved capillary beds by degrading the resolution of UPD and PA. **Significance.** Since dual PA-ultrasound systems are already widespread, this approach could become an essential brick of quantitative PA imaging.

## 1. Introduction

Photoacoustic (PA) imaging is an emerging biomedical imaging technique that has proven able to reveal optical contrast at depths in soft tissues (Beard 2011, Wang and Hu 2012, Taruttis and Ntziachristos 2015). This technique has the potential to offer real-time imaging (Deán-Ben et al 2015) at a high spatial resolution and penetration depth of endogenous, exogenous or genetically-encoded absorbers (Jathoul et al 2015). PA imaging theoretically offers quantitative chromophore concentration imaging (Cox et al 2009). This applies to ranges of applications going from oxygenation imaging (Laufer et al 2006, Lutzweiler et al 2014, Tzoumas et al 2016) to molecular imaging (Li et al 2007, Lutzweiler et al 2014). For this task, one needs to use spectral features from images obtained in a PA spectroscopic acquisition and operate spectral unmixing of the different chromophores (Cox et al 2009). However, PA sources are proportional to the absorbed energy density, equal to the product between optical absorption coefficient and optical fluence $\mu_a(\lambda, \mathbf{r}) \times \phi(\lambda, \mathbf{r})$. As optical transport is wavelength-dependent, an effect referred to as spectral coloring (Hochuli et al 2019) requires $\phi(\lambda, \mathbf{r})$ to be computed for each wavelength to apply a fluence correction before any demixing algorithm. When the transport or scattering properties of the background medium can be considered homogeneous, techniques have been proposed to account for the fluence distribution in the images (Ranasinghesagara and Zemp 2010, Bu et al 2012, Brochu et al 2016, Jeng et al 2021). When unknown heterogeneous media are considered, a light propagation model has to be employed (Cox et al 2012) within an inverse algorithm. From a PA image, the basic goal of the inverse model is to recreate a digital twin by finding the optimal optical properties throughout the medium. Ultimately, this allows to accurately recover decoupled absorption and scattering properties, providing further quantitative spectral analysis. Learning techniques based on simulated examples appear as simpler ways than the full inverse problem (Tzoumas et al 2016, Bench et al 2020, Gröhl et al 2021, Kirchner and Frenz 2021, Liang et al 2024, Shi et al 2024) but they
might lack generality since they rely on an ensemble of simulated examples and usually consider uniform background scattering properties.

As PA imaging is both optical and acoustical, complete inverse problems involve two distinct parts (Huang et al 2013, Poudel et al 2019). In this work, we focus on the optical inversion part, considering the acoustic reconstruction ideal. Optical transport is governed by the radiative transport equation (RTE) which can be in some cases simplified to the diffusion approximation (Arridge 1999, Jeng et al 2021). As diffusion approximation can be applied only a few mean free paths away from the illuminated surface, modeling of the RTE is preferable for more accuracy (Cox et al 2012). RTE can be modeled using Finite-element methods (Saratoon et al 2013) but the advent of parallel GPU computing makes Monte Carlo-based methods more attractive for their computational speed (Leino et al 2020, Macdonald et al 2020, Hänninen et al 2022). A radiance Monte Carlo inverse method using the directional information of photons was introduced to inverse the RTE without the diffusion approximation (Hochuli et al 2016, Buchmann et al 2019, Zheng et al 2023).

In quantitative PA imaging, inverse problems aim at recovering the optical scattering $\mu_s(\mathbf{r})$ and absorption $\mu_a(\mathbf{r})$ coefficients, from measured data. It has been demonstrated that the solution to the PA inverse problem is not unique when using a single wavelength and illumination. The reason why is when the calculated absorbed energy $H$ differs from the measured one, it is impossible to know whether it comes from an error in the absorption or scattering distribution, which are both responsible for light attenuation. This means that multiple couples $(\mu_s(\mathbf{r}), \mu_a(\mathbf{r}))$ satisfy the inversion (Cox et al 2009). To solve the non-unicity issue, researchers proposed to use multiple illumination patterns (Zemp 2010, Shao et al 2011, Saratoon et al 2013), but this requires dedicated equipment and illumination calibrations which might not be adequate for medical applications. Another approach would be to constrain the problem by having strong priors. Using multiple wavelengths, priors on the chromophore spectral properties can lead to a unique solution (Cox et al 2009) in both $\mu_s(\mathbf{r})$ and $\mu_a(\mathbf{r})$. There are a few experimental results of multispectral inversions but some assumptions are made on $\mu_s$ and inversions only estimate $\mu_a$ (Buchmann et al 2019, Zheng et al 2023). Clearly, the simultaneous inversion of $\mu_s(\mathbf{r})$ and $\mu_a(\mathbf{r})$ is considered as a bottleneck when the medium is heterogeneous with no prior information. In this issue, 2D simulated data input is considered. Our aim is to propose a method to facilitate the recovery of both scattering and absorption coefficients with single illumination and wavelength. To this end, we will rely on blood vessels.

Blood vessels are ubiquitous in the body. Contained within these vessels, hemoglobin is also the main optical absorber present in tissues for near-infrared light. In this paper, we propose to use a multimodal approach with the coupling of PA imaging and another modality offering blood vessels imaging to gain prior information. This modality is supposed to be of a distinct physical nature and to be insensitive to the optical transport phenomenon. A good candidate as a coupled modality is ultrasound imaging since simultaneous PA and US imaging is common (Niederhauser et al 2005, Kolkman et al 2008, Merčep et al 2017). Ultrasound-informed PA reconstructions were proposed using ultrasound-driven anatomical information (Yang et al 2020, Pattyn et al 2021). While this approach is relevant, it does not exploit quantitative information that can be given by ultrasound. One interesting modality is ultrasound power Doppler (UPD) (Murphy and Rubin 1997) that enables to image specifically the vasculature. Recent developments improved significantly the sensitivity of the technique that is now applied for functional imaging in the brain (Macé et al 2011, Deffieux et al 2021). State-of-the-art UPD was recently combined with PA imaging (Robin et al 2022, Godefroy et al 2023). Different high-resolution instrumentation able to image blood vessels (ex: CT angiography, MR-angiography scanners or optical imaging) may also be coupled to PA imaging but it seems adequate to consider ultrasound-based techniques for their resolution, penetration and as the same equipment can be used for dual PA and ultrasound imaging. Considering blood vessels as a support of valuable information, we propose to use UPD to address the non-unicity problem with a single-illumination and wavelength PA+US system.

UPD consists in evaluating the fluctuations caused by blood flow. It consists in integrating the whole Doppler spectrum. Consequently, UPD measures the agitation of blood and is insensitive to the flow speed. More precisely, it quantifies blood volume (the fraction of blood within one imaging pixel/voxel) under the assumption of uniform hematocrit (Rubin et al 1995, Macé et al 2011). For a given concentration of hemoglobin in blood vessels, it is theoretically possible to convert a UPD picture into a prior optical absorption coefficient due to blood, including in resolved vessels with blood volume of 100% and unresolved vessels with intermediate blood volume down to unresolved capillary beds. The conversion is done through a multiplication by a factor that is a *priori* unknown. We demonstrate it can be optimally determined within our proposed inversion procedure. The absorption due to blood is considered as dominant over other chromophores around 800 nm. In the general case, we then consider UPD gives the effective absorption coefficient locally within vascularized area. Furthermore, while PA imaging and UPD can both image blood, it is essential to account for resolution differences between both techniques. However, resolution degradation will lead to a blurring of the vasculature seen by UPD, which can be problematic if the blurring overlaps with
some heterogeneities. A strategy is further defined to use the prior in the case of an optical heterogeneity surrounding or covering some blood vessels.

We based our inverse algorithm on a previously-reported Monte Carlo radiance inverse problem (Hochuli et al 2016). Thanks to the prior brought by UPD, we show we can simultaneously recover the absorption and scattering coefficients in regions surrounded or covered by blood capillaries with or without prior background information. The media in this paper are composed of blood capillaries (or vessels) among a background medium that includes optical absorption and scattering heterogeneities that are considered as markers of pathologies. For the sake of generality, we assume these pathologies result in local changes of optical absorption, scattering, or both. We consider that these pathologies are not visible from the vascular imaging perspective which implies that they are solely brought to light from their optical properties.

In this issue, we present the inversion model in the general case, where $\mu_s$ and $\mu_a$ are retrieved from a complete PA measurement. We then present the principles of the UPD prior information used in this study and show simulated examples. We demonstrate the UPD prior mitigates the non-uniqueness encountered in single-wavelength and illumination PA imaging. We show the inversion does not need prior information on background, inclusions and hemoglobin concentration and is robust to the presence of noise. Finally, we address the resolution aspects with various decimation of a high resolution medium.

## 2. Theory

### 2.1. Monte Carlo inverse radiance algorithm

The reader can refer to the studies that developed this specific theoretical framework for PA imaging inversion using a Monte Carlo inverse radiance scheme (Saratoon et al 2013, Hochuli et al 2016, Buchmann et al 2019). We briefly expose here the main equations and formulas without repeating their demonstrations. We define the ideal PA reconstruction as $H(\mathbf{r}) = \mu_a(\mathbf{r}) \Phi(\mathbf{r})$ where $\mu_a(\mathbf{r})$ is the optical absorption in the medium and $\Phi(\mathbf{r})$ is the fluence. We distinguish the ‘measured’ reconstruction $H_{\text{meas}}$ and the reconstruction by the model at given iterations $H(\mu_a, \mu_s)$. The radiance $\phi(\mathbf{r}, \hat{s})$ is solution to the RTE, also referred to as forward problem:

$$ (\hat{s} \cdot \nabla + \mu_a(\mathbf{r}) + \mu_s(\mathbf{r})) \phi(\mathbf{r}, \hat{s}) - \mu_s(\mathbf{r}) \int_{S^{N-1}} P_{\theta}(\hat{s}, \hat{s}') \phi(\mathbf{r}, \hat{s}') d\hat{s}' = q(\mathbf{r}, \hat{s}) \quad (1) $$

where $\mathbf{r}$ is the position, $\hat{s}'$ and $\hat{s}$ are the original and scattered propagation directions, $P_{\theta}(\hat{s}, \hat{s}')$ is the scattering phase function, $q(\mathbf{r}, \hat{s})$ is a source term and $S^{N-1}$ indicates integration over angle in $N-1$ dimensions ($N=2$ in 2D). The forward model is run based on an initial guess of the model parameters (priors). The inversion is then performed using a gradient descent algorithm. A classical way is to derive functional gradients from the definition of the adjoint equation:

$$ (-\hat{s} \cdot \nabla + \mu_a(\mathbf{r}) + \mu_s(\mathbf{r})) \phi^*(\mathbf{r}, \hat{s}) - \mu_s(\mathbf{r}) \int_{S^{N-1}} P_{\theta}(\hat{s}, \hat{s}') \phi(\mathbf{r}, \hat{s}')^* d\hat{s}' = Q_{\text{adj}} \quad (2) $$

where $\phi^*$ is the solution to the adjoint equation, $Q_{\text{adj}} = \mu_a(H - H_{\text{meas}})$ is the source term of the adjoint problem. The error function is defined as: $\epsilon(\mu_a(\mathbf{r}), \mu_s(\mathbf{r})) = \frac{1}{2} \int (H(\mathbf{r}, \mu_a(\mathbf{r}), \mu_s(\mathbf{r})) - H_{\text{meas}}(\mathbf{r}))^2 d\mathbf{r}$.

To update the model at each iteration, one has to run the forward and the adjoint problem, allowing to calculate functional gradients as:

$$ \begin{aligned} \frac{d\epsilon}{d\mu_a} &= -\Phi(\mu_a, \mu_s)(H_{\text{meas}} - H(\mu_a, \mu_s)) + \int_{S^{N-1}} \phi(\hat{s}, \mu_a, \mu_s) \phi^*(\hat{s}, \mu_a, \mu_s) d\hat{s} \\ \frac{d\epsilon}{d\mu_s} &= \int_{S^{N-1}} \phi(\hat{s}, \mu_a, \mu_s) \phi^*(\hat{s}, \mu_a, \mu_s) d\hat{s} - \int_{S^{N-1}} \int_{S^1} \phi(\hat{s}, \mu_a, \mu_s) P_{\theta}(\hat{s}, \hat{s}') \phi(\hat{s}', \mu_a, \mu_s)^* d\hat{s}' d\hat{s} \end{aligned} \quad (3) $$

Monte Carlo light simulations usually output the angularly integrated fluence. In order, to compute the radiance scalar products, a decomposition in an angular harmonics basis is performed. In 2D, we use the circular harmonics basis (also called classical Fourier basis). Note that in 3D, one would use the spherical harmonics basis. During the simulation, each photon will impact local Fourier coefficients on their path, depending on their propagation direction resulting in $a_n(\mathbf{r})$ and $b_n(\mathbf{r})$ maps for the forward simulation and $a_n^*(\mathbf{r})$ and $b_n^*(\mathbf{r})$ (Buchmann et al 2019).

During the Monte-Carlo transport simulation, each pixel/voxel accumulates the weights $\delta W$ deposited by the photons along their paths. The Fourier coefficients are updated following, for the 2D case, $a_0 = \frac{\delta W}{\pi}$, $a_n = \frac{\delta W}{\pi} \cos(n\theta')$ and $b_n = \frac{\delta W}{\pi} \sin(n\theta')$, where $\theta'$ is the incoming angle of the photon.
As demonstrated previously (Buchmann et al 2019), using the Fourier decomposition, the gradients (3) become in the 2D case:

$$ 
\begin{aligned}
\frac{d\epsilon}{d\mu_a} &= -\Phi(\mu_a, \mu_s)(H_{\text{meas}} - H(\mu_a, \mu_s)) + \frac{\pi}{2}a_0a_0^* + \pi \sum_n (a_n a_n^* + b_n b_n^*) \\
\frac{d\epsilon}{d\mu_s} &= \pi \sum_{n=1} [a_n a_n^* + b_n b_n^*] (1 - g^n)
\end{aligned}
\quad (4) $$

where g is the anisotropy factor of the medium. The computation of this gradient at each iteration allows updating $\mu_a(\mathbf{r})$ and $\mu_s(\mathbf{r})$ towards the optimal solution.

## 2.2. Quantitative power Doppler

After removing the tissue signals with a filtering step, the power Doppler image $\sigma_{\text{US}}^2(\mathbf{r})$ is calculated as the temporal variance in each pixel. In a previous publication, we have introduced a quantitative power Doppler formula involving the imaging part (Falco et al 2024):

$$ \sigma_{\text{US}}^2(\mathbf{r}) = P_i^2(\mathbf{r}) \times \left(\frac{\Delta Z(\mathbf{r})}{Z_0}\right)^2 \times \eta(\mathbf{r}) W[\eta(\mathbf{r})] V_{\text{RBC}}(\mathbf{r}) \times \int f(\mathbf{r}') |h_{\text{US}}|^2(\mathbf{r}', \mathbf{r}) d\mathbf{r}' \quad (5) $$

where $P_i$ is the incident acoustic pressure, $\Delta Z(\mathbf{r})/Z_0$ is the relative mean impedance fluctuation induced by a single red blood cell, $\eta$ is the volume fraction of blood (hematocrit), $W(\eta)$ is the packing factor, $h_{\text{US}}$ is the US imaging point spread function (PSF) and $f(\mathbf{r})$ defines the blood vessels. This formula is composed of a product between excitation amplitude $P_i$, medium parameter $\Delta Z(\mathbf{r})/Z_0$, statistical medium blood parameters $\eta W(\eta)$ and object-PSF integral (on the right side of the equation). Object-PSF integral involves the squared envelope of the US PSF which, thanks to its positivity, do not allow destructive interferences and enable a faithful representation of the blood vessels at the given resolution. The square factor enables a reduction of the effects of potential side lobes. The object-PSF integral is proportional to the blood volume $BV(\mathbf{r})$, partial volume of vessels in the pixel, or equivalently fractional vessel density.

It is important to note that modern Power Doppler, involving spatio-temporal SVD-filtering (Demené et al 2015), efficiently reconstructs all the blood vessels at the given resolution and is insensitive to the flow speed and direction, unlike other modes such as color Doppler.

In the case the PSF is non-stationary, the blood volume image is derived from a normalization of the object PSF integral by the local squared L2 norm of the PSF: $||h(\mathbf{r})||_2^2 = \int |h_{\text{US}}|^2(\mathbf{r}', \mathbf{r})d\mathbf{r}'$. The normalized object-PSF integral is equal to 1 in vessel larger than the PSF all over the imaging domain. Assuming stationary blood properties (uniform hematocrit and red blood cell impedance and volume) and a correction of the image by a calibrated incident pressure $P_i$, one can assume the whole normalized image $\sigma_{\text{US,n}}^2$ is proportional to the blood volume:

$$ \sigma_{\text{US,n}}^2(\mathbf{r}) = (P_i^2(\mathbf{r}) \cdot ||h(\mathbf{r})||_2^2)^{-1} \sigma_{\text{US}}^2(\mathbf{r}) \propto BV(\mathbf{r}). \quad (6) $$

This holds true in resolved vessels that are delineable on the images but also in unresolved small vessels. Using a wavelength range where hemoglobin absorption dominates over the other chromophores (at 800 nm), Power Doppler reflects the pixel effective absorption coefficient. In a previous study (Falco et al 2024), we demonstrated experimentally that the ratio between PA fluctuation imaging and UPD is proportional to the fluence in blood vessels, proving that Power Doppler detains the blood absorption quantity.

Blood vessels regions are defined, as a domain $M_{\text{BV}}$ satisfying $\mathbf{r} \in M_{\text{BV}}$ if $\tilde{\sigma}_{\text{US}}^2(\mathbf{r}) > \epsilon$, where $\epsilon$ is a threshold set on the UPD image amplitude. Thus, we consider that the absorption coefficient in blood vessels can be derived from the power Doppler image through:

$$ \text{if } \mathbf{r} \in M_{\text{BV}}, \mu_{\text{a,blood}}(\mathbf{r}) = K_1 K_2 \sigma_{\text{US,n}}^2(\mathbf{r}) = \xi \sigma_{\text{US,n}}^2(\mathbf{r}) \quad (7) $$

where $K_1 = 1/\sigma_{\text{US}}^2(BV(\mathbf{r}) = 1)$ is constant dependent on the imaging device, which can be eventually calibrated. It is the inverse of the value taken by the US Doppler picture in a large vessel that has a maximum $BV = 1$ value. Besides, $K_2$ depends on the hemoglobin and blood characteristics of the subject, that can vary among individuals but can be considered uniform within the same individual. Overall, $\xi = K_1 K_2$ is the resulting factor. To summarize, the Power Doppler image only gives a relative $\mu_a$ value where blood vessels are present. $\xi$ is a priori unknown, we further show $\xi$ can be retrieved from a one-dimensional optimization procedure.
When a fluence model can be applied to the measurement, PA imaging can also probe blood volume through the total hemoglobin concentration in the pixel (Wang et al 2006). This does not hold true when spectral coloring affects the measurements or if the medium optical properties are unknown unless inverse problems are used.

In Falco et al (2024), quantitative fluence imaging is achieved experimentally after compensating the resolution differences between the UPD and PA fluctuation images. This was done in the most general case of non-stationary PSF in 3D, using a PSF filter defined to degrade the resolution of UPD to the one of PA fluctuation imaging. For a N-points spatial grid, the PSF filter is a N × N matrix where each column represents a kernel that applies to each position. In the present framework, this PSF filter can be obtained from model-based matrices $H_{PA}$ and $|H_{US}|^2$ containing all the reconstructed PSF for each point of the grid. Consequently, this formalism can be applied to degrade $\sigma_{US}^2$ to the PA resolution. However, in this study, we assume no resolution mismatch between power Doppler and PA imaging. Nonetheless, such a PSF filter could be useful for an experimental validation of our findings.

## 3. Methods

### 3.1. Implementation

To compute the light radiance inside a domain of size $N_z \times N_x$ defined by optical parameters $\mu_a(\mathbf{r})$, $\mu_s(\mathbf{r})$ and $g(\mathbf{r})$, we use a Monte-Carlo implementation of the radiative transfer equation (RTE) called MCXLAB (Yu et al 2018), with GPU acceleration. This open-source software offers light fluence computing in heterogeneous media. We modified the source code in order to project the radiance onto the harmonic basis following the implementation proposed by Hochuli et al (2016). The modified version of MCXLAB was compiled on a UNIX-based operating software (OS) and then on a Microsoft Windows OS to achieve this study. MCXLAB is piloted from MATLAB 2021a (Mathworks, MA, USA) as is the inversion algorithm based on outputs from MCXLAB.

### 3.2. Inverse problem implementation

The forward part simulates a real, physical event where the photons propagate from the source in an absorbing and scattering medium. They lose part of their weight along their path until they reach a threshold (typically $10^{-6}$. $W_0$ where $W_0$ is their initial weight) and die. The adjoint model, however, relates to a virtual situation where the photons propagate backward, due to the minus sign before the convection term (2). With a change of variable $\hat{s} \to -\hat{s}$, we see that the field $\phi(\mathbf{r}, -\hat{s})^*$ is solution to the original RTE equation. Thus, the same model is run for adjoint to get $\phi(\mathbf{r}, -\hat{s})^*$ and the angular variable is taken as its opposite (reversed directions) to obtain $\phi(\mathbf{r}, \hat{s})^*$. Before computing the scalar products, the photon directions are reversed, ie $\theta' \Rightarrow \theta' + \pi$ where $\theta'$ represents the incoming angles of arrival of the photon. Equivalently, the odd Fourier coefficients ($a_1^*, b_1^*, a_3^*, b_3^*$, etc) change sign. To maximize the reach of the photons emitted in the adjoint simulation, we consider a loss less adjoint propagation by conserving the weight of the photons along their paths. The photons die when they leave the domain.

A total variation (TV) regularization (Jensen et al 2012) is used to dampen the effects of noise inherent to MC simulations and potentially to PA measurements. The TV of $\mu_a$ or $\mu_s$ is defined as

$$T(\mu) = \int_{\Omega} \|\nabla \mu(\mathbf{r})\|_{2} d\mathbf{r}. \qquad (8)$$

Note that the Euclidean norm is not taken to the square. A prox-function smoothing version using the Huber function instead of the L2 norm is used in order to make it a derivable function and compute the gradients $\nabla T_{\tau(\mu_s)}$ and $\nabla T_{\tau(\mu_a)}$ where $\tau$ is the proximal smoothing parameter (Jensen et al 2012). Compared to the Tikhonov regularization, the TV presents the advantage of preserving the edges and the high values. The cost function and gradients are changed to:

$$\epsilon(\mu_a(\mathbf{r}), \mu_s(\mathbf{r})) = \frac{1}{2} \int (H(\mathbf{r}) - H_{\text{meas}}(\mathbf{r}))^2 d\mathbf{r} + \alpha_a T(\mu_a(\mathbf{r})) + \alpha_s T(\mu_s(\mathbf{r})) \\ \left[ \frac{d\epsilon}{d\mu_a} \right]_{\text{TV}} (\mathbf{r}) = \frac{d\epsilon}{d\mu_a} (\mathbf{r}) + \alpha_a \nabla T_{\tau}(\mu_a(\mathbf{r})) \\ \left[ \frac{d\epsilon}{d\mu_s} \right]_{\text{TV}} (\mathbf{r}) = \frac{d\epsilon}{d\mu_s} (\mathbf{r}) + \alpha_s \nabla T_{\tau}(\mu_s(\mathbf{r})) \qquad (9)$$

The proximal smoothing $\tau$ was fixed to $5.10^{-4}$. The TV regularization parameters $\alpha$ were adjusted depending on the situation especially if noise is added ($\alpha_s = 1.10^{-5}$ or $1.10^{-4}$ and $\alpha_a = 3.10^{-5}$ or $3.10^{-4}$, resp. without or with noise).
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Model and measurements. (a) Vascular network high resolution binary map. (b) Power Doppler image. (c) and (d) True $\mu_a$ and $\mu_s$ maps. (e) PA imaging $H_{\text{meas}}$. (f) Binary thresholded vascular map. The white color indicates the region where $\mu_a$ is updated through the iterations, $\mu_a$ is locked in the black vessel region. Units are mm$^{-1}$ for absorption and scattering.</figcaption></figure>

An Adam scheme was used to optimize the convergence rate (Kingma and Ba 2014). We consider a 8 × 8 mm 2D medium illuminated by a single homogeneous source on the top of the sample with a 0.1 mm step size. 5.10⁶ photons are launched at each iteration on a standard desktop PC with a NVIDIA GeForce RTX 3090 GPU, each iteration taking about 2.0s.

### 3.3. Tissue model and image generation

Retinal capillary structures are used to represent blood vessels, taken from the DRIVE dataset (https://drive.grand-challenge.org/) (Khanal and Estrada 2019). We consider these images give a binary definition of the vessels at high resolution (see figure 1(a)). Equation (5) states the Power Doppler is obtained by a convolution between the object and the squared envelope of the US PSF $|h_{\text{US}}(\mathbf{r})|^2$. $|h_{\text{US}}(\mathbf{r})|^2$ is a smooth function typically of the form of a 2D Gaussian function. This convolution performs a low-pass filter on the object and the result can be downsampled without loss of information. Effectively, power Doppler represents the blood volume, which is defined as a relative amount of blood in each given pixel. Our Power Doppler image shown in figure 1(b) is obtained by this decimation process, resulting in the blood volume quantity (between 0 and 1). To simplify the test object, resulting values lower than 0.18 are set to zero. Local blood absorption is also governed by the amount of blood in each pixel. Therefore, to obtain the input $\mu_a$ map, we convert the blood volume picture to blood absorption using the $\xi$ coefficient (whole blood: $\mu_a = 4 \text{ cm}^{-1}$ at 800 nm) and we add the inclusions ($\mu_a = 0.8 \text{ cm}^{-1}$ and $\mu_a = 1 \text{ cm}^{-1}$) and background values ($0.25 \text{ cm}^{-1}$) resulting in the figure 1(c). The anisotropy factor and refraction index are supposed homogeneous and taken as $g = 0.9$ and $n = 1.37$. The scattering coefficient of blood is assumed equal to the one of background, as a simplification but it could be set to a more accurate representative value if experiments require this in the future. The $\mu_s$ map is defined as shown in figure 1(d), with one of the scattering inclusions co-aligned with one of the absorbing inclusions ($\mu_s = 30 \text{ cm}^{-1}$, $\mu'_s = 3 \text{ cm}^{-1}$) and another on top left ($\mu_s = 40 \text{ cm}^{-1}$, $\mu'_s = 3 \text{ cm}^{-1}$) and a background value of $\mu_s = 10 \text{ cm}^{-1}$ ($\mu'_s = 1 \text{ cm}^{-1}$). The addition of these scattering inclusions allows demonstrating that the method is efficient with an arbitrary heterogeneous medium in terms of scattering properties. All the optical coefficient values describe a realistic tissue as previously reported in Sandell and Zhu (2011), Jacques (2013).

The incident illumination is defined as a uniform plane wave coming from the top. Following a forward Monte-Carlo simulation giving the fluence $\phi(\mathbf{r})$, the PA image is obtained by the product $\mu_a(\mathbf{r}) \cdot \phi(\mathbf{r})$ (see figure 1(e)).

### 3.4. Priors

To start the inversion, the numerical model needs to be initialized with $\mu_a$ and $\mu_s$ maps, called priors $\mu_{a,0}$ and $\mu_{s,0}$. In the presented examples, we always consider a uniform prior $\mu_{s,0}$ that can be equal to the true
<figure><img src="image_6.png" /><figcaption>Figure 2. Summary of model inversion using a prior coming from ultrasound power Doppler.</figcaption></figure>

background value or not. The same principle is applied to $\mu_{a,0}$ except we use the $\mu_{a,blood}$ values in blood vessels (for $\mathbf{r} \in M_{BV}$). Note that the $\mu_{a,blood}$ and associated conversion factor $\xi$ is initially considered to be the exact one but we further show we can recover $\xi$ if unknown.

### 3.5. Inverse algorithm

The MC forward model is run using the defined medium properties to obtain the ‘measured’ PA image $H_{meas}$ which would be measured in an experiment. In parallel, a numerical model is defined and the previously defined priors allow running a first forward model to obtain the model PA image $H_0$ (see figure 2).

At each iteration, the model PA image $H_n$ is computed from $\mu_{a,n}$ and $\mu_{s,n}$ and compared to the measured value $H_{meas}$, leading to the adjoint source $Q_{adj} = \mu_{a,n} \cdot (H_{meas} - H_n)$. The adjoint propagation is computed from that source, resulting in radiance coefficients which, combined with the forward coefficients, give the gradients expressed in equation (4) (figure 1). Since the source $Q_{adj}$ is signed (positive or negative), 2 simulations are run separately with respectively the positive and the negative part, then the results are subtracted. 4 Fourier harmonics were found sufficient in our model to construct correct gradients.

### 3.6. Optimization of the power Doppler—$\mu_a$ coefficient

Equation (7) expresses the absorption coefficient in blood vessels. As previously stated, the factor $\xi$ presented in the equation depends on the hemoglobin and blood characteristics of the subject, which can vary among individuals but can be considered uniform within the same individual. Quantifying and fixing this factor is crucial for the convergence of the proposed inverse problem, as we will demonstrate in the results section. The absorption coefficient in blood vessels can potentially be quantified by taking a blood sample from the specific patient and measuring it through a spectrophotometer (Thomas et al 2024). However, in an in-vivo application, this approach would make ineffective the non-invasive benefits of PA and ultrasound imaging. Hence, we propose a solution that relies on a one-dimensional optimization where the only unknown is $\xi$. Specifically, we implemented a simple optimization method, which is performed only once before the start of the entire inverse radiance algorithm. Starting from an initial guess, $\xi$ is changed using a fixed step size of $\pm 0.5 \cdot 10^{-3} \text{ mm}^{-1} \text{ a.u.}^{-1}$, in a direction determined by the initial sign of $d\epsilon_{\mu_a}/d\mu_a$. The error function defined in section 2.1 decreases and we let it re-increase for 10 iterations: its minimum is recovered. The routine requires around 60 iterations to find the optimal $\xi$. The validation of this method is discussed in section 4.3 and applied to all the results presented subsequently.

### 3.7. Quantitative aspects

In some of the presented results, some additive noise is added to the input PA image. This noise is generated using a normal distribution of given standard deviation $\sigma_{image}$. We define signal-to-noise ratio (SNR) as the ratio between the maximum of the image and background noise standard deviation. To compare different
situations, a quantitative analysis of results is presented in terms of normalized root mean square error:

$$ NRMSE = \left\langle \left( \frac{\mu(\mathbf{r}) - \mu_{\text{true}}(\mathbf{r})}{\mu_{\text{true}}(\mathbf{r})} \right)^2 \right\rangle \quad (10) $$

where $\mu$ can be either $\mu_s$ or $\mu_a$ and $\langle \rangle$ denotes spatial average over the whole domain.

## 3.8. Resolution aspects

We show the effect of a resolution degradation in results in figure 6. To obtain the UPD picture, the high-resolution binary vessel map is filtered by a 2D Gaussian before getting decimated. In parallel, the true $\mu_a$ map is first defined at high resolution with binary vessels and inclusions. Fluence distribution is simulated in the medium at high resolution and the PA image $H_{0,HR}$ is computed. Resolution is then degraded the same way as UPD to obtain the input $H_0$ of the inversion. The width of the 2D Gaussian is 20, 40, and 60 pixels corresponding to 0.4, 0.8, and 1.2 mm. The decimation factor remained at the same value of five leading to the same image size for the 3 cases.

Resolution degradation blurs the contour of the vascular network which eventually creates pixels containing mixed inclusion and blood. As our algorithm is based on keeping the prior blood absorption values in blood vessels, the inclusion cannot be retrieved this way. To optimize a region containing a supposed inclusion, the prior blood absorption is set to the background $\mu_a$ value in a disk enclosing the inclusion. This means the prior blood absorption is only used around the inclusion. We will show our approach still has a beneficial effect on stabilizing a unique solution in the inversion. The way $\mu_s(\mathbf{r})$ is treated is not affected by this resolution degradation.

# 4. Results

## 4.1. Effect of blood prior

We first consider the case where background $\mu_a$ and $\mu_s$ are known. Figure 3 shows the inversion results with (c-e) and without knowledge (c'-e') of the vascular structure. Without any prior, the reconstruction of both $\mu_a$ and $\mu_s$ is impossible due to the ill-posed problem. With the vascular information, $\mu_a$ is known and fixed in the vessels. This constraint is sufficient to overcome the non-uniquity problem and allows the model to converge to an accurate solution (middle column). We remind that when the prior is considered, $\mu_a$ is not updated in the blood vessels throughout the inversion. However, when no prior is considered, $\mu_a$ is indeed updated everywhere along with $\mu_s$. The NRMSE for $\mu_a$ is 0.76% with prior and 12.9% without prior, and for $\mu_s$ the NRMSE is 8.71% with prior and 24.8% without prior. Interestingly, the cost function minimal value ends up higher with the prior ($1.2 \times 10^{-7}$) than without prior ($0.9 \times 10^{-7}$) which depicts how non-uniqueness can lead to false results. However, using the prior leads to a 16-fold decrease in the initial cost function value.

## 4.2. Effect of inaccurate background priors

We consider an error of +20% for both $\mu_a$ and $\mu_s$ priors in the background. We consider the blood prior as exact. The recovered $\mu_a$ in figure 4(a.2) is very close to the true $\mu_a$ in figure 3(a) with an NRMSE of 4.7%. The effect is more visible for $\mu_s$ with a recovered background slightly perturbated (see figure 4(b.2)) than in figure 3(e). The NRMSE value is increased to 23%. Most importantly, convergence is conserved when inaccurate background priors are considered.

## 4.3. Effect of inaccurate blood absorption value

In section 3.6, we presented a secondary algorithm that can be run to initially optimize the $\mu_a$ blood value through the $\xi$ coefficient 1D optimization. Here, while keeping the situation of figure 4 with wrong background priors, we consider an error in the prior $\xi$ coefficient of 10%. Keeping this false blood absorption still leads to a solution, shown in figures 5(a1)-(b1). The recovered $\mu_a$ is not visually affected (NRMSE = 9.5%) but the $\mu_s$ reconstruction is degraded (NRMSE = 31.9%).

The optimization routine of $\xi$ is performed to lead to improved results (see figures 5(a2)-(b2)), with NRMSE = 5.0% for $\mu_a$ and 26.5% for $\mu_s$. $\xi$ is only optimized at the first iteration of the inversion, based on the homogeneous wrong $\mu_s$ and $\mu_a$ background priors. The evolution of the cost function during this 1D-optimization is shown in figure 5(c). Through this 62 iterations process, the optimal $\xi$ has an absolute error of 1.95% compared to the true value. This result shows considering the blood prior map is effective even if the absolute blood absorption and background priors are erroneous.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Results of the inversion of $\mu_a$ and $\mu_s$ with a prior on $\mu_a$ in the vessels (c)–(e) and without (c')–(e'). True values are displayed in (a) and (b). Units are mm$^{-1}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Inversion of $\mu_a$ and $\mu_s$ with wrong background priors. Units are mm$^{-1}$. (a.1)–(b.1) Priors with 20% error in background values. (a.2)–(b.2) Inversion results. True $\mu_a$ and $\mu_s$ maps are found in figures 3(a) and (b).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Optimization of the UPD to blood absorption coefficient $\xi$ starting with a 10% error and with a 20% error in background properties. (a1)–(b1) Result without optimization of $\xi$. (a2)–(b2) Result with optimization of $\xi$. c Trend of the cost function as the optimization algorithm iterates. The marker indicates the optimal point, representing the value of $\xi$ resulting in the lowest cost.</figcaption></figure>

### 4.4. Effect of resolution

The results presented in this last section account for resolution degradation by considering the effect of the PSF on both PA and ultrasound prior images. The three cases shown in columns 1, 2 and 3 of figure 6 involve, as specified in section 3.8, a degradation of $H_0$ via a 2D Gaussian filtering, with a width of the Gaussian of 20, 40, and 60 pixels, respectively. To simulate a more realistic scenario, these examples include not only the PSF-object convolution but also all the conditions previously mentioned, as a 20% error in the background values of $\mu_a$ and $\mu_s$, and an unknown $\xi$. Additionally, some noise was added to the image before resolution degradation. As an initial step in the inversion, the $\xi$ parameter has been estimated using the optimization algorithm, yielding errors of 2.5%, 5.8%, and 6.6% for the three examples.

Before discussing the findings, we computed the SNRimage on the degraded $H_0$ images using the equation provided in section 3.7. The SNRs were evaluated for the first 15 pixels at the top and the last 15 pixels at the bottom of the image. Across the three degradation scenarios, the SNR at the top was 26.7, 23.9, and 19.9, while at the bottom, it was 9.2, 11.8, and 14.9. The decrease in fluence is responsible for a lower SNR at the deeper region.

From the inversion of $\mu_s$ presented in figures 6(f.1)–(f.3), we can notice that the recovered values are not significantly affected by the resolution degradation (resp. NRMSE = 33.9%, 43.9% and 42.9%). Considering the absorption, figures 6(e.1)–(e.3) shows that the retrieved $\mu_a$ closely matches the PSF-convolved $\mu_a$ in figures 6(c.1)–(c.3). The NRMSE values are equal to 76.5%, 44.1% and 31.0%.

It is worth noting that the inversion is performed using a PSF-convolved $H_0$, and the PSF-convolved $\mu_a$ is not provided to the algorithm. Nevertheless, as fluence is slowly varying, the PSF-convolved $\mu_a$ is a good target to compare with the inverted $\mu_a$.

## 5. Discussion

In this issue, we propose a strategy to add valuable information for quantitative PA imaging. We show how UPD can help to estimate a prior absorption map that reduces the number of unknowns and leads to accurate inversions of $\mu_a$ and $\mu_s$. We show the algorithm is robust to unknown background as well as unknown blood properties. While degrading the imaging resolution, we present an effective strategy when there is no clear separation between optical heterogeneities and blood vessels, in the presence of noise.

The implementation of a radiance Monte Carlo inverse problem was done using MCXLAB. This was made possible by adding features to the source code of MCXLAB to retrieve harmonic coefficients describing the angular information of radiance from the simulation. This preliminary study presents only 2D results. The 3D nature of optical transport and PA imaging requires a 3D framework for working on real data. The exposed implementation and theory are ready to use for 3D data and will be part of a future study. Our algorithm uses an Adam optimizer which speeds up convergence but the computational cost could be
<figure><img src="image_10.png" /></figure>

reduced, for example using Machine learning to adjust the number of needed photons in each Monte Carlo simulation (Macdonald *et al* 2020).

We applied the inverse model to a medium composed of absorption and scattering inclusions surrounded by some blood vessels. We assumed prior $\mu_a$ could be retrieved from a combined vascular imaging modality (based on UPD), allowing the model to start with more knowledge on the medium. The conversion factor $\xi$ can be optimized at different times of the inversion. Different strategies were tried such as optimization every 5 iterations but this turned out to be less efficient compared to a single optimization at the beginning, at least with our medium model. The $\xi$ optimization algorithm was chosen here very simple and there is room to improve its efficiency, which was not our goal here.

Note that the adjoint equation is a modified RTE with a source term $Q_{adj}$ defined as the residual error at each position $\mathbf{r}$, multiplied by $\mu_a(\mathbf{r})$. Hence, $\mu_a(\mathbf{r})$ has a great impact on this source term. This means that without prior, the adjoint source is initially governed by $\Delta H$ which spreads on the whole domain with higher values at shallow depths resulting in parameters optimization from top to bottom, with a high crosstalk between $\mu_a$ and $\mu_s$ and error propagations. However, the $\mu_a$-gradient is dominated by a local term proportional to $\Delta H$ resulting in updates within the right locations. The algorithm converges towards a solution that shows blood vessels and inclusions in $\mu_a$, but with some inaccuracies due to the result in $\mu_s$ that does not show any sign of inclusions. Adversely, the use of the prior $\mu_a$ values in blood generates hot spots in
the adjoint sources $Q_{\text{adj}}$ acting as fluence reference points. Instead of optimizing the coefficients from top to bottom, the inversion happens in a more evenly spatially-distributed manner. The resulting adjoint field propagates from the blood regions and directs towards the regions to be updated. We observed updates in $\mu_s$ can be faster inside blood vessels. Thus, using the prior reduces the crosstalk between $\mu_a$ and $\mu_s$, leading to the convergence of the $\mu_s$ and to a very accurate final $\mu_a$ value. It is then more suited that the prior $\mu_a$ covers the whole region of interest. For future applications, this means the quality of the inversion might depend on the imaged medium.

We show that the sensitivity to noise can be improved using a TV-regularization. This type of regularization uses an empirically-tuned parameter. An automatic algorithm could be useful for general applications. A strong assumption in this study is that the PA reconstructions are ‘perfect’ and so are UPD reconstructions. Behind this assumption, lies the fact that PA imaging and UPD have the same PSF and thus, image in the same manner which is not true in the general case. For example, UPD lateral and axial resolution is different from the one obtained in PA imaging (Godefroy et al 2023) (Falco et al 2024), since ultrasound compounding performs focus in emission, and also since the round-trip nature of ultrasound improves the axial resolution two-fold. If we apply the method proposed here on real data, it means that the prior $\mu_a$ from vascular imaging would not have the same native resolution as $H_{\text{meas}}$. $Q_{\text{adj}}$ being the product of two different quantities with different resolutions, it will have strong spatial fluctuations, for example, on the border of a smooth object, which might prevent converging to an acceptable solution. For application to real data, a resolution-equalization needs to be performed and the easiest manner to proceed is to degrade the best resolution towards the worst of the two techniques as performed in Falco et al (2024). Applying model-based reconstructions for both modalities might still result in different resolutions (Vilov et al 2020).

Blood vessels can have diameters ranging from centimeters down to $10\ \mu\text{m}$. This multiscale nature has been approached in this study with the resolution degradation. In UPD, the evaluated quantity on the image is the blood volume (Macé et al 2011), the same quantity as MRI or CT can evaluate (Zaharchuk 2007). Blood volume images are proportional to the blood volume fraction within each pixel. A maximum value is reached when the pixel is saturated with blood. Particularly, when observing unresolved capillary-bed structures, blood volume is still quantitatively determining the volume fraction of blood. Thus, correct mean effective $\mu_a$ pixel values can still be imposed by this adequate vascular imaging techniques in the case of unresolved vessels.

However, when resolution blurs the contour of the vasculature prior, it is necessary to ignore the prior blood absorption in a region that needs to be updated in $\mu_a$. This strategy supposes that there is a way to identify the presence of a $\mu_a$ heterogeneity not due to flowing blood. This situation can also arise in the presence of additional chromophores (lipid, melanin, contrast agent accumulation) inside or very close to some vessels. The delineation of this region could be done manually or automatically based on cost function values after a first run of the inversion or based on an Ultrasound anatomical image showing a suspicious structure. Herein, we considered vessels in a given range of scales. When applying the method *in vivo*, we could face the case where there is a higher variety of vessel scales, with a capillary bed composing the background. In this case, two possible strategies can be considered. First, if UPD sensitivity allows probing capillary bed structures as an effective vessel density, also responsible for the background optical absorption, this information could lead to a prior optical absorption in the background. Hence, the PA contrast will enable recovering scattering properties and optical absorption contributions other than blood, provided that optical absorption coefficients are left free to some extent during the updates, using regularizations, for example. The other strategy involves considering the UPD prior only in larger vessels, simply by using a threshold on the UPD images. In this case, the same method as proposed here, which locks the prior $\mu_a$ only in larger vessels, will let the inverse problem recover background optical absorption and heterogeneities. The practical use of either strategy depends on the blood vessel structure as well as the PA reconstruction ability to recover given spatial-frequency contents. These considerations will be the focus of further work.

Visibility limited issues were not studied here. The limitations due to the detectors were neglected. It could be interesting to study the effect of the vascular imaging prior (which is not affected by visibility artifacts) in the inversion of PA visibility-limited data. This implies considering a model of detectors and acoustic reconstruction in addition.

Another assumption of our study is that optical absorption in pixels containing blood vessels is only due to hemoglobin. This is not the general case, since other chromophores such as lipids, water and melanin can participate in absorption and spectral unmixing needs to be performed before quantitative blood analysis. However, our approach uses the isobestic point at 800 nm where the effects of other chromophores are minimized (Cox et al 2009). Additionally, oxygenation modulations do not affect blood absorption coefficients at this wavelength. In the advent of a highly-absorbing heterogeneity mixed with blood, the strategy presented in section 3.8 is effective to estimate the total $\mu_a$ in the region. Eventually, from the difference $\mu_a(\mathbf{r}) - \mu_{a,\text{blood}}(\mathbf{r})$, one could estimate the absorption due to the materials other than blood.
Over such a small propagation distance, we consider acoustic attenuation can be neglected. However, depending on the frequency band of the signals, acoustic attenuation effects could be significant and would have to be corrected within the employed reconstruction.

Our tissue model was chosen with a finer resolution in $\mu_a$ than in $\mu_s$. When testing the inversion on $\mu_a$ and $\mu_s$ maps using smaller inclusions (results not shown here), the resolution of the result was significantly higher for $\mu_a$, while such small inclusions were not recovered for $\mu_s$. This happens because the PA signal is significantly more sensitive to $\mu_a$ than to $\mu_s$. Additionally, the gradient analytical formulae contribute to this difference, since $d\epsilon/d\mu_a$ contains a local term proportional to the measurement itself $\Delta H$ whereas $d\epsilon/d\mu_s$ involves forward and adjoint radiance which plays a secondary role.

## 6. Conclusion

As part of an optical inverse problem, we show the combination of two modalities that measure blood vessels in different manners, UPD and PA imaging, can bring the additional information needed for convergence. The fact that the support of information constituted by blood vessels is omnipresent in vivo and holds strong signals is a great opportunity. This methodology brings a promising strategy for quantitative imaging in PAs.

## Data availability statement

The code will be shared after further development and experimental demonstrations. The data that support the findings of this study are available upon reasonable request from the authors.

## Acknowledgments

This work was supported by the Agence National de la Recherche (ANR) under the FULBOX Project (ANR-21-CE19-0003-01). The authors wish to thank Qianqian Fang for technical help with MCXLAB.

## References

Arridge S R 1999 *Inverse Problems* 15 R41

Beard P 2011 Biomedical photoacoustic imaging *Interface Focus* 1 602–31

Bench C, Hauptmann A and Cox B 2020 *J. Biomed. Opt.* 25 085003

Brochu F M, Brunker J, Joseph J, Tomaszewski M R, Morscher S and Bohndiek S E 2016 *IEEE Trans. Med. Imaging* 36 322–31

Bu S, Liu Z, Shiina T, Kondo K, Yamakawa M, Fukutani K, Someda Y and Asao Y 2012 *IEEE Trans. Biomed. Eng.* 59 1354–63

Buchmann J, Kaplan B A, Powell S, Prohaska S and Laufer J 2019 *J. Biomed. Opt.* 24 1

Cox B T, Arridge S R and Beard P C 2009 *J. Opt. Soc. Am. A* 26 443–55

Cox B, Laufer J G, Arridge S R and Beard P C 2012 *J. Biomed. Opt.* 17 061202

Deán-Ben X L, Ford S J and Razansky D 2015 *Sci. Rep.* 5 10133

Deffieux T, Demené C and Tanter M 2021 *Neuroscience* 474 110–21

Demené C et al 2015 *IEEE Trans. Med. Imaging* 34 2271–85

Falco I, Bossy E and Arnal B 2024 *Phys. Med. Biol.* 69 165026

Godefroy G, Arnal B and Bossy E 2023 *Sci. Rep.* 13 2961

Gröhl J, Kirchner T, Adler T J, Hacker L, Holzwarth N, Hernández-Aguilera A, Herrera M A, Santos E, Bohndiek S E and Maier-Hein L 2021 *Sci. Rep.* 11 6565

Hänninen N, Pulkkinen A, Arridge S and Tarvainen T 2022 *Proc. SPIE* 27 083013

Hochuli R, An L, Beard P C and Cox B T 2019 *J. Biomed. Opt.* 24 121914

Hochuli R, Powell S, Arridge S and Cox B 2016 *J. Biomed. Opt.* 21 126004

Huang C, Wang K, Nie L, Wang L V and Anastasio M A 2013 *IEEE Trans. Med. Imaging* 32 1097–110

Jacques S L 2013 *Phys. Med. Biol.* 58 R37

Jathoul A P et al 2015 *Nat. Photon.* 9 239–46

Jeng G S, Li M L, Kim M, Yoon S J, Pitre J J J, Li D S, Pelivanov I and O'Donnell M 2021 *Nat. Commun.* 12 716

Jensen T L, Jørgensen J H, Hansen P C and Jensen S H 2012 *BIT Numer. Math.* 52 329–56

Khanal A and Estrada R 2019 Dynamic deep networks for retinal vessel segmentation (arXiv:1903.07803)

Kingma D P and Ba J 2014 arXiv:1412.6980

Kirchner T and Frenz M 2021 *J. Biomed. Opt.* 26 085001

Kolkman R G, Brands P J, Steenbergen W and van Leeuwen T G 2008 *J. Biomed. Opt.* 13 050510

Laufer J, Delpy D, Elwell C and Beard P 2006 *Phys. Med. Biol.* 52 141

Leino A A, Unttila T, Mozumder M, Pulkkinen A and Tarvainen T 2020 *IEEE Trans. Med. Imaging* 39 2985–95

Li P C, Wei C W, Liao C K, Chen C D, Pao K C, Wang C R C, Wu Y N and Shieh D B 2007 *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* 54 1642–7

Liang Z, Zhang S, Liang Z, Mo Z, Zhang X, Zhong Y, Chen W and Qi L 2024 *Photoacoustics* 37 100601

Lutzweiler C, Meier R, Rummeny E, Ntziachristos V and Razansky D 2014 *Opt. Lett.* 39 4061–4

Macdonald C M, Arridge S and Powell S 2020 *J. Biomed. Opt.* 25 085002

Macé E, Montaldo G, Cohen I, Baulac M, Fink M and Tanter M 2011 *Nat. Methods* 8 662–4

Merčep E, Deán-Ben X L and Razansky D 2017 *IEEE Trans. Med. Imaging* 36 2129–37

Murphy K J and Rubin J M 1997 Power doppler: it's a good thing *Seminars in Ultrasound, CT and MRI* vol 18 (Elsevier) pp 13–21
Niederhauser J J, Jaeger M, Lemor R, Weber P and Frenz M 2005 *IEEE Trans. Med. Imaging* **24** 436–40

Pattyn A, Mumm Z, Alijabbari N, Duric N, Anastasio M A and Mehrmohammadi M 2021 *Photoacoustics* **23** 100275

Poudel J, Lou Y and Anastasio M A 2019 *Phys. Med. Biol.* **64** 14TR01

Ranasinghesagara J C and Zemp R J 2010 *J. Biomed. Opt.* **15** 046016

Robin J, Özbek A, Reiss M, Dean-Ben X L and Razansky D 2022 *IEEE Trans. Med. Imaging* **41** 846–56

Rubin J M, Adler R S, Fowlkes J B, Spratt S, Pallister J E, Chen J F and Carson P L 1995 *Radiology* **197** 183–90

Sandell J L and Zhu T C 2011 *J. Biophoton.* **4** 773–87

Saratoon T, Tarvainen T, Cox B T and Arridge S R 2013 *Inverse Problems* **29** 075006

Shao P, Cox B and Zemp R J 2011 *Appl. Opt.* **50** 3145–54

Shi M, Vercauteren T and Xia W 2024 *Proc. SPIE* **12842** 69–72

Taruttis A and Ntziachristos V 2015 *Nat. Photon.* **9** 219–27

Thomas A, Rietberg M, Akkus M, van Soest G and Francis K J 2024 *Proc. SPIE* **13319** 100673

Tzoumas S, Nunes A, Olefir I, Stangl S, Symvoulidis P, Glasl S, Bayer C, Multhoff G and Ntziachristos V 2016 *Nat. Commun.* **7** 12121

Vilov S, Arnal B, Hojman E, Eldar Y C, Katz O and Bossy E 2020 *Sci. Rep.* **10** 4637

Wang L V and Hu S 2012 *Science* **335** 1458–62

Wang X, Xie X, Ku G, Wang L V and Stoica G 2006 *J. Biomed. Opt.* **11** 024015

Yang H, Jüstel D, Prakash J, Karlas A, Helfen A, Masthoff M, Wildgruber M and Ntziachristos V 2020 *Photoacoustics* **19** 100172

Yu L, Nina-Paravecino F, Kaeli D and Fang Q 2018 *J. Biomed. Opt.* **23** 010504

Zaharchuk G 2007 *Am. J. Neuroradiol.* **28** 1850–8

Zemp R J 2010 *Appl. Opt.* **49** 3566–72

Zheng S, Yingsa H, Meichen S and Qi M 2023 *Phys. Med. Biol.* **68** 065009