CC BY
16
Combined Monte Carlo and k-Wave Simulations for Reconstruction of Blood Oxygen Saturation in Optoacoustics: A Pilot Study
Valeriya Perekatova*, Daria Kurakina, Aleksandr Khilov, and Mikhail Kirillin
Institute of Applied Physics RAS, 46 Ul'yanov str., Nizhny Novgorod 603950, Russia
* e-mail: valeriyaperekatova@gmail.com
Abstract. Optoacoustic (OA) imaging of biological tissues is a modern technique allowing for three-dimensional blood oxygen saturation mapping based on OA spectroscopy data. Since biological tissues are optically inhomogeneous and the spatial distribution of optical parameters within a biological tissue is a priori unknown, Monte Carlo simulation technique is traditionally used to estimate the distribution of probing illumination within tissues in quantitative OA reconstruction. Currently, machine learning techniques are actively employed for reconstructing 3D distribution of blood oxygen saturation or estimating optical properties of biological tissues based on training datasets. In this paper, systemic calculations of synthetic OA images of a medium with embedded vessel-like structures were performed to create a training dataset for machine learning employing combined application of the Monte Carlo technique for direct solution of optical problem and difference-space pseudo-spectral approach implemented through k-Wave Toolbox calculations for the acoustical part. The calculations were performed for probing wavelengths of 532 nm, 658 nm and 1064 nm, which are commonly employed in spectral OA imaging. Simulated OA data for different orientation, diameter and embedding depth of blood vessels allows analyzing the effect of these parameters on the formation of OA image and the reconstruction of blood oxygen saturation. The ratio of OA signals corresponding to probing wavelengths of 658 nm and 1064 nm was employed for simple reconstruction of blood oxygen saturation in silico for different vessel geometries with the precision of < 3-15% for the most of blood vessels diameters and embedding depths and the range of blood oxygen saturation values > 0.8. The obtained set of synthetic OA data has high potential as a training set for employment in machine learning techniques aiming at mapping blood oxygenation based on spectral OA data. © 2022 Journal of Biomedical Photonics & Engineering.
Keywords: optoacoustic imaging; Monte Carlo modeling; k-Wave modeling; blood oxygen saturation mapping.
Paper #3562 received 10 Nov 2022; revised manuscript received 7 Dec 2022; accepted for publication 12 Dec 2022; published online 24 Dec 2022. doi: 10.18287/JBPE22.08.040511.
1 Introduction
Blood oxygen saturation (blood oxygenation level, StO) is an important physiological indicator actively
used in monitoring of functional parameters in biomedical studies, such as the evaluation of chemo-and radiation therapy impact on tumor tissue [1, 2], monitoring of wounds healing [3], or estimation of
tissue reaction to photodynamic therapy [4]. In this regard, the development of a non-invasive visualization technique capable of providing fast in vivo three-dimensional mapping of blood oxygenation level is of high importance. Optical imaging techniques [5, 6] benefit from both non-invasiveness and high sensitivity to blood oxygenation owing to significant difference in absorption spectra of oxy- and deoxyhemoglobin in optical spectral range. Traditionally, blood oxygenation level is estimated with diffuse optical spectroscopy (DOS) [2, 7, 8] based on the registration of the spectra of the probing radiation passed through the biological tissue and the subsequent reconstruction of the medium absorption spectra from the measurement data [9, 10]. This approach, however, provides with the estimation of blood oxygenation value averaged over a particular measurement volume within biological tissue. Blood saturation mapping provides a more detailed information regarding tissue oxygenation and could be performed using diffuse optical tomography (DOT) [11, 12], however, at the expense of spatial resolution of order of 0.5 cm. Nevertheless, the estimation of local oxygenation at tissue layer level in certain blood vessels may be required, for example, for the studies of hemodynamics [13, 14]. Optoacoustic (OA) imaging [15-19] is a modern hybrid imaging technique with high potential in saturation mapping [20] owing to high spatial resolution combined with increased probing depth compared to purely optical modalities.
This technique is often employed for superficial vasculature mapping providing resolution down to tens of microns and imaging depth of several centimeters, which allows obtaining in vivo morphological and functional information on the vascular bed of biological tissue. The contrast of OA images is determined by the difference in the absorption coefficient of blood and surrounding components of the biological tissue revealing high potential for in vivo angiography. Since the optical absorption spectra of oxy- and deoxyhemoglobin differ significantly, spectral separation of these chromophores can be performed based on OA spectroscopy data with consequent determination of blood oxygen saturation [20, 21].
However, the local pressure increment that occurs in the medium due to absorption of pulsed laser radiation by optical inhomogeneities is proportional to the local optical absorption coefficient and the fuence of the probing optical radiation at a given point in the medium. Since biological tissues are optically inhomogeneous and the distribution of optical parameters of a biological tissue is a priori unknown, Monte Carlo simulation technique is traditionally used to estimate the distribution of probing illumination within tissues in quantitative OA imaging [22]. Reconstruction of 3D distribution of functional tissue parameters from OA data benefits from current development of machine learning (ML) techniques. Currently, ML methods (for example, artificial neural networks) are actively employed for reconstructing maps of blood oxygen saturation distribution or estimating
optical properties of biological tissues based on training OA datasets [23-25]. Some of approaches using neural networks estimate oxygenation from single-pixel pressure spectra [26, 27]. Employment of convolutional neural networks allows to process the spectra of entire 2D (3D) images utilizing spatial and spectral information at the same time [28-31]. Due to the lack of experimental training data with a priori known saturation distribution within tissue, synthetic data are commonly used for training the neural network and proof-of-concept verification [23, 24]. In order to accumulate synthetic datasets for training algorithms, it is essential to realistically model the physical processes of optoacoustic imaging.
In this paper, to create a training sample dataset, systemic calculations of OA images were performed employing sequential application of the Monte Carlo method for direct solution of optical problem and k-space pseudo-spectral approach for acoustics. A simple approach to blood oxygen saturation reconstruction based on the revealed monotonous dependence of registered OA signal ratio for probing wavelengths of 658 nm and 1064 nm was proposed and tested at the synthetic data. Obtained results indicated high potential of the developed approach to OA images simulations for further development of ML-based algorithms of the estimation of blood oxygen saturation.
2 Materials and Methods
2.1 Monte Carlo Modeling
Traditional Monte Carlo (MC) technique for light transport studies is based on modeling of a large number of random photon trajectories in turbid media with following statistical analysis of the collected data [32]. In this work previously developed platform for three-dimensional MC modeling of light propagation in biological tissues [22, 33, 34] was customized by implementing elongated absorbers mimicking different vessels filled with blood within biotissue. It was employed for the generation of the maps of the absorbed light dose distribution at probing wavelength lex H(x,y,z,Xex) in a flat tissue-like homogenous medium containing cylinders mimicking blood vessels of different diameters, embedding depths and orientations corresponding to typical morphological parameters.
The full size of a tissue sample considered in simulations was 20 x 20 x 10 mm3. The light source was considered to be a plane wave illuminating the sample from the top. A total of 107 photons were launched into the medium perpendicular to the tissue surface and used to calculate H(x,y,z,Xex) for each of three probing wavelengths: Xex = 532, 658, and 1064 nm. The wavelengths of 532 nm and 1064 nm are the wavelengths commonly employed for OA imaging [35], while the pair of 658 nm and 1064 nm is optimal for more precise estimation of blood oxygen saturation [36]. The refractive index was set equal to 1.38 for blood vessels as well as for background tissue.
Table 1 Ranges of the optical properties used in Monte Carlo simulations.
Human skin in vivo Blood (StO2 varies from 0.1 to 1)
Wavelength, nm
1064
mm-1 mm-1
g
mm
mm
g
532 0.39 13.29 0.67 19.04-22.48 92.07-98.65 0.989-0.988
658 0.06 9.96 0.67 1.51-0.15 82.94-92.74 0.986-0.985
0.04
5.9
0.67 0.23-0.45 42.19-47.75 0.973-0.972
Table 2 Parameters of the simulated medium and ultrasonic detector for acoustic simulations.
Biotissue
Detector
Acoustic attenuation, mm"1MHz"1
Time scale, ns
Grid element size, mm
Volume size, mm3
Numerical aperture
Maximum supported frequency, MHz
Radius, mm
0.01
10
0.05
17.8 x 7.8 x 12.05
0.71
15
• « . a>
Configl 1 mm
b)
• • • •
Config2 1 mm
• • • C) •
Config3 1 mm
Fig. 1 Three examples of side views of blood vessels configuration: a) 0.1 mm < zo < 0.7 mm and 0.1 mm < d < 0.25 mm, b) 0.9 mm < z0 < 1.5 mm and 0.3 mm < d < 0.45 mm, c) 0.5 mm < z0 < 1.1 mm and 0.75 mm < d < 0.9 mm.
Blood oxygen saturation in blood vessels was varied from 0.1 to 1 with the step of 0.1, resulted in different values of the blood optical properties calculated based on absorption and scattering spectra for oxygenated and deoxygenated whole blood taken from [37]. Optical properties for the background tissue were taken from paper [38] for human skin in vivo. All the values of optical properties at corresponding wavelengths are shown in Table 1.
2.2 k-Wave Modeling of OA Images
Modeling of OA images of biotissue based on MC-calculated absorption map H(x,y,z,Xex) was performed using the k-Wave Toolbox [39], which is a standard for acoustic calculations. The toolbox implements a k-spatial pseudo-spectral method for solving first-order acoustic equations for homogeneous and inhomogeneous media. The employed toolbox can take into account arbitrary distribution of inhomogeneities and acoustic absorption. Monte Carlo simulated maps of the absorbed light dose H(x,y, z,Xex) were employed as distributed sources of ultrasonic waves with initial pressure distribution
p0(x,y,z,Xex) = G • H(x,y,z,Xex),
where G is the Gruneisen parameter.
The k-Wave modeling of OA microscopy images was implemented in three-dimensional geometry by the solution of the forward acoustic problem for each detector position in XZ-plane. The ultrasonic detector was supposed to be a spherically focused antenna with the parameters close to those for real detector employed in an OA microscope [40-42]. The acoustic properties, grid size and parameters of the biotissue and the ultrasonic detector are shown in Table 2. After k-Wave modeling, each OA microscopy B-scan was processed with the reconstruction algorithm described in paper [43].
Both MC modeling and k-Wave simulations were performed with the employment of workstation equipped with 64-core CPU AMD Ryzen Threadripper 3990X with frequency of 2.9 GHz with 256 Gb RAM onboard.
2.3 Building Synthetic Dataset for Saturation Reconstruction
Numerical simulations were performed for two different geometry configurations of blood vessels distribution. The first one was a simplified model of skin (Fig. 1), containing four blood vessels with parallel axes with different diameters d and embedding depths zo.
1
1
6
Fig. 2 Complex configuration of vessel net for MC modeling with randomly distributed blood vessels: (a) enface schematic (top view, XY-plane) of vessels orientation (color encodes number of pixels belonging to vessels in corresponding projection); (b) 2D projection (side view, XZ-plane) of blood vessels corresponding to the projection marked with dashed line in Fig. (a).
Fig. 1 shows side views of three examples of such skin models: 0.1 mm < z o < 0.7 mm and
0.1 mm < d < 0.25 mm (Config 1), 0.9 mm < zo < 1.5 mm and 0.3 mm < d < 0.45 mm (Config 2), 0.5 mm < z0 < 1.1 mm and 0.75 mm < d < 0.9 mm (Config 3). The total number of combinations of parameters considered in simulations amounted 200, which covered ten different embedding depths z0 from 0.1 to 1.9 mm and twenty different diameters d from 0.05 to 1 mm. One configuration contained 4 vessels with different parameters resulting in 50 different configurations. Employment of this geometry aims at creating a training dataset that cover the entire range of possible morphological parameters of superficial blood vessels. For all the combinations, ten different values of blood oxygen saturation varying from 0.1 to 1 were considered. Thus, the total number of calculated maps of
H(x,y,z) for different 50 geometries, 3 probing wavelengths and 10 oxygenation values amounted 1500.
Fig. 3 shows the schematic of spherically focused ultrasound detector and initial pressure distribution obtained by MC modeling in projection to XZ-plane for vessel orientation configurations from Fig. 1c and Fig. 2b.
Based on simulated optoacoustic images, a training dataset was created for the implementation of reconstruction methods aiming at the estimation of local blood oxygenation level with a voxel-by-voxel resolution.
2.4 Estimation of Blood Oxygen Saturation
For considered probing wavelengths, the ratios of OA signals p(x,y, z,Xex) in central cross-sections (y = 0) at
i . . p(x,y=0,z,XPV=532 nm)
wavelengths ( , ' ' -- ,
p(X'y=0'Z'Xex=658 nm) p(X'V=0'Z'XPY=658 nm)
^ ' ' ex ' ) were
different
p(x,y=0,z,Xex=532 nm)
p(x,y=0,z,Xex=1064 nm) ' p(x,y=0,z,Xex=1064 nm)
calculated for simple geometries and all considered blood oxygen saturation values. The trends of the OA signal ratios vs blood oxygen saturation were monotonous, and the ratio
R(x,y,z) =
p(x,y,z,Aex=658 nm) p(x,y,z,Aex=1064 nm)'
(2)
revealed the smallest variation within all cases (vessel depth and location) for particular values of StO2. In this connection, the average value
ß _ (p(x,y=0,z,Xex=6S8 nm))\(Xy=oz)Eves (p(x,y=0,z,Xex=1064 nm))\(x,y=o,z)eves ,
where (p(x,y = 0,z,Xex))li (x,y=o,z)eves means the OA signal in central cross-section of blood vessel corresponding to probing wavelength lex averaged over all available cases, was employed to construct the mean dependence R(StO2). The inverse function StO2(R) was further employed for voxel-by-voxel estimation of blood oxygen saturation based on the calculated ratio R(x,y = 0,z):
StO2(x,y = 0,z) = StO2(R(x,y = 0,z)). (4)
1 mm • • • • . 1 mm a) — -r-- ~ b)
Fig. 3 Schematic of three-dimensional detector and initial pressure distribution obtained from MC modeling in projection to XZ-plane for vessel configuration from: a) Fig. 1c; b) Fig. 2b.
St02=0.9 Configl
•
1 mm
St02=0.9 Config2
g)
St02=0.9 1 mm
Config3
j)
St02=0.9
CG 1 mm
• • C) 1 mm
• • • • f) 1 mm
—
• • • • ¡1 1 mm
i)
— ...... 1 mm
10.8 0.6 0.4
I 0.2
'o
Fig. 4 Three-dimensional maps (side view, XZ-plane) ofthe distribution of absorbed light dose in XZ-plane H(x,y = 0,z,Xex) for probing wavelengths of Xex = 532 nm (a, d, g, j), Xex = 658 nm (b, e, h, k), and Xex = 1064 nm (c, f, i, l) and for configurations shown in Fig. 1 (Config1 - a, b, c; Config2 - d, e, f; Config3 - g, h, i) and 2D projection of complex geometry (CG) from Fig. 2b (j, k, l). Notations v1-v4 indicate vessels chosen for further detailed analysis.
Fig. 5 Reconstructed OA B-scans p(x,y = 0,z,Xex) obtained from k-Wave modeling for probing wavelengths of Xex = 532 nm (a, d, g, j), Xex = 658 nm (b, e, h, k), and Xex = 1064 nm (c, f, i, l) and for configurations shown in Fig. 1 (Config1 - a, b, c; Config2 - d, e, f; Config3 - g, h, i) and 2D projection of complex geometry (CG) from Fig. 2b (j, k l).
3 Results and Discussion
3.1 Monte Carlo Modeling
Fig. 4 shows central cross-sections of three-dimensional maps of absorbed light dose H(x,y = 0,z,Xex) in blood vessels for their various diameters and embedding depths for probing wavelengths of 532 nm (Fig 4a), 658 nm (Fig. 4b) and 1064 nm (Fig. 4b) and three different configurations from Fig. 1 and 2D projection of complex geometry (Fig. 2b).
Absorption coefficient of blood for 532 nm is significantly higher than those for 658 nm or 1064 nm, which significantly limits the penetration of radiation of this wavelength into blood vessels, which is manifested by primary absorption of light on the boundaries of blood vessels (Fig. 4a), while for 658 nm and 1064 nm probing light is absorbed in the entire volume of vessels (Figs. 4b, c).
3.2 k-Wave Modeling of OA Images
At the second stage, based on the obtained absorption maps (Fig. 4), the acoustic response of the medium was calculated using the k-Wave Toolbox [39] and reconstructed by the earlier proposed algorithm [43] for all probing wavelengths. The results corresponding to configurations from Fig. 1 and Fig. 2b are shown in Fig. 5.
Due to limitations of RAM capacity of the employed workstation (although it is highest available for non-server solutions) the grid element size was chosen equal to 50 ^m, which exceeds the scanning step in a real OA imaging system [42]; in addition, the central frequency of ultrasound detector is lower for numerical modelling as compared to the real system, which leads to the presence of artifacts in the reconstructed OA images.
z (mm)
Fig. 6 Dependencies of average OA signals vs depth obtained in vivo [41] (Experimental data) and from k-Wave modeling (Simulated data).
Moreover, real OA systems operate with complicated illumination geometry (circular or through the center of the ultrasound detector [42, 44]), while this study is limited by a plane wave illumination case. These limitations lead to the presence of pronounced artifacts around the vessels location areas. In order to analyze the adequacy of the employed simulation approach we have compared the averaged A-scan (OA signal vs Z coordinate) for the experimental OA image of a human palm in vivo data [41] with the same dependence for a simulated image (Fig. 6). This comparison reveals good
agreement of OA signal attenuation, and the only discrepancy is observed in the top tissue layer, where simulated signal exceeds the experimental one. This discrepancy is explained by a simplified skin model employed in simulations, in which morphological skin layers are not distinguished, while in real skin upper stratum corneum and epidermis layers contain no blood thus producing weak OA signal. Nevertheless, good correspondence of the general trend allows concluding on the adequacy of the employed model given that proper normalization is done.
3.3 Reconstruction of Blood Oxygen
Saturation from Ratio of Signal Intensities
Figs. 7a-c show the dependencies of OA signal intensities averaged over each of four blood vessels (denoted as v1-v4 in Fig. 4h) on blood oxygen saturation, which were obtained in numerical simulation. Figs. 7d-f shows the ratios of the signals registered at different probing wavelengths. It is worth mentioning that the signals show monotonous dependence on StO2 owing to monotonous dependence of the medium optical properties on this parameter. However, they show different trends owing to difference in absorption spectra of oxy- and deoxyhemoglobin. The revealed monotonous dependence of signal ratio in this respect is more important, since it demonstrates the feasibility of employing the ratiometric approach in saturation reconstruction.
100-|
90-
13 80-
03
A 70-
t-
1_
CM 60-
Ol
IX)
CL V 50-
40-
30-
160-,
v1 140-
v2
v3 13 120-
v4 03
A, 100-
EE
£=
CO 80-
m
CO
CL V 60-
40-
20-
0.4 0.6 StO,
110-1 C 100
=3
® 90
/T
E 80 £=
g 70
If 60
50
40 0.0
StO,
2.0-1 1.8-
A
¡1.6-
CO • 10 1 A
(O 1.4-CL • £1.2-H 1.0-
<N • 8.0.8-CL •
v 0.60.4-
0.4 0.6
StO,
0.4 0.6 StO,
2.41
A
IT20
c
<X> £ 1.6
.v
2^1.2-1
0.4 0.6 StO,
- v1
v2 v3 v4
Fig. 7 Dependencies of OA signal intensity on blood oxygen saturation for probing depths of a) 532 nm, b) 658 nm, and c) 1064
1 1. 1. p(x,y=0,z,XPY=532 nm) . p(x,y=0,z,XPY=532 nm) ~ p(x,y=0,z,XPY=658 nm)
nm and corresponding ratios: d) ,-i-^ro__, , e) i; , f
p(x,y=0,z,Xex=65B nm)
Notations v1-v4 indicate vessels in Fig. 4h from left to right.
p(x,y=0,z,Xex=1064 nm)
p(x,y=0,z,Xex=1064 nm)
b
a
StO
d
e
f
io:
3.02.52.01.51.00.50.0-
a
1.00.8-
O 0.6-
0.4-
0.2-
0.0
b
0.0 0.2
0.4 0.6 StO2
0.8
1.0
Calculated data Linear fit
0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 R
Fig. 8 (a) Mean dependence of OA signals ratio R vs blood oxygen saturation StO2 calculated for all 50 simple geometry configurations; (b) inverse function StO2 (R) and its approximation with linear fit.
The ratio of OA signals in central cross-section of blood vessels corresponding to probing wavelengths of 658 nm and 1064 nm revealed the smallest variance within all three ratios (Fig. 6d-f), thus, it was employed for the generation of synthetic curve
^^(fi - (P(x,y=°^ex=6S8 nm))\{X}y=o,Z)eves p. 8a
2\ (p(x,y=0,z,Xex=1064 nm))\(x,y=o,z)eves)' .
shows mean dependence of ratio R vs blood oxygen saturation averaged over all fifty considered geometries of blood vessels, which could be considered as a generalized dependence that could be used in saturation map construction.
Fig. 8a shows that the variance of the signal ratio R decreases with the increase in blood oxygen saturation, hence, it is smaller for the blood oxygen saturation values > 0.8, which allows considering the obtained dependence to be used in a simple algorithm of blood oxygenation map reconstruction. The linear fit of inverse function StO2 (R) depicted in Fig. 8b was employed for voxel-by-voxel reconstruction of blood oxygen saturation map based on the map of OA signal ratios ( R (x,y = 0,z) =
p(x,y=0,z,Aex=658 nm) . . . _ . . ,
—----) derived from simulations results.
p(xfy=0fzfAex=1064 nm)
Fitting was performed within Matlab Curve Fitting Toolbox and revealed that dependence StO2(R) is described with the following Eq.:
StO2(R) = -0.59R + 1.37,
(5)
with high precision (coefficient of determination r2 = 0.9991). Reconstructed saturation value StO
2,rec
was
calculated voxel-by-voxel from corresponding values of R. For the cases of reconstructing StO2,rec value outside the physiological range, below 0 or above 1, they were assumed to be equal to 0 or 1, respectively.
1 mm • ■ . a>
St02=0.8 sto2,tnjc
1 mm • d)
St02=0.8 R
1 mm • . »
St02=0.8 St02 rec in vessels
1 mm b)
• • • •
St02=0.9 sto2,true
1 mm k)
• • * •
St02=0.9 sto2.,.c in vessels
1 mm f)
^ m
StO,=l R
1 mm 1)
m
sto2=i St02,rec
1 mm o)
m
sto2=i 6St02,rK
Fig. 9 Blood vessel configurations with StO2 true(x,y = 0,z) value (a, b, c), calculated ratio R(x,y = 0,z) (d, e, f), reconstructed blood oxygen saturation StO2rec(x,y = 0,z) (g, h, i), reconstructed blood oxygen saturation StO2rec(x,y = 0,z) in vessels only (j, k, l) and the discrepancy of reconstructed blood oxygen saturation SStO2 rec(x,y = 0,z) (m, n, o) for initial blood oxygen saturation 0.8 (a, d, g, j, m), 0.9 (b, e, h, k, n), and 1 (c, f, i, l, o). Vessel configurations correspond to Fig. 1.
Fig. 10 Blood vessel configurations with StO2true(x,y = 0,z) value (a, b, c), calculated ratio R(x,y = 0,z) (d, e, f), reconstructed blood oxygen saturation StO2rec(x,y = 0,z) (g, h, i), reconstructed blood oxygen saturation StO2rec(x,y = 0,z) in vessels only (j, k, l) and the discrepancy of reconstructed blood oxygen saturation SStO2 rec (x, y = 0, z) (m, n, o) for initial blood oxygen saturation 0.8 (a, d, g, j, m), 0.9 (b, e, h, k, n) and 1 (c, f, i, l, o). Vessel configurations correspond to Fig. 2b.
The results of reconstruction for several simple geometry cases (shown in Fig. 1) are shown in Fig. 9, while the results for the complex geometry (Fig. 2b) are shown in Fig. 10. Each figure contains vessel location schematic, calculated XZ projections maps of OA signals ratio R (x, y = 0, z), projections of full reconstructed maps of StO2 and reconstructed StO2 values in vessels only, and the discrepancy between reconstructed (StO2,rec) and true (StO2,true) (blood oxygen saturation values, calculated as:
discrepancy 5StO2,rec. For thick vessels SStO2,rec value increases with the increase in their embedding depth due to attenuation of probing light and corresponding decrease in OA signal. This effect is demonstrated in details in Fig. 11, which shows the dependencies of mean error in vessels on their diameter d and embedding depth zo calculated for all vessels in all considered simple geometry OA images (Fig. 1) calculated as mean discrepancy over all voxels belonging to a certain vessel:
SStO2,rec(x, y = 0, z) = = StO2,rec(x, >" = 0, Z) - StO2,tme(X y = 0, z).
SStO
(6)
2,rec
= <8StO2,rec(*.y.z)>|
(7)
For the calculation of SStO2,rec the true blood oxygen saturation value in surrounding tissue outside the vessels was assumed equal to that in vessels, since the distribution of the light field generating OA signal is governed by absorption in vessels, although the considered signals originate from the nearby voxels.
As can be seen from Figs. 9-10, the proposed simple reconstruction algorithm provides with adequate estimation of blood oxygen saturation in vessels with no a priori knowledge of vessel embedding depth required. The discrepancy 5StO2,rec between reconstructed and true blood oxygen saturation, StO2,rec and StO2,true, is the highest for thin and subcutaneous vessels (Figs. 9-10). The ratio R is higher for subcutaneous vessels (Fig. 9d-f, Fig. 10d-f), which leads to smaller values of StO2,rec and, thus, to larger
From Fig. 11 one can see that described effects lead to the presence of the minimal 5StO2,rec (precision < 3%) area, which corresponds to the combinations of parameters that provide the mean values of the calibration curve (Fig. 8a) employed for the reconstruction algorithm. For the combination of vessel diameter d and embedding depth zo above that area the precision is <15% for StO2,true > 0.8. The reconstruction algorithm employs mean values of R, which correspond to the intermediate values of vessel embedding depth thus providing better accuracy for these depths (Fig. 11). For blood vessels located below those depths, the ratio R is smaller due to stronger attenuation of probing light at 658 nm compared to 1064 nm. This leads to overestimation of StO2 and positive values of SStO2,rec (Fig. 11).
Fig. 11 Dependencies of discrepancy of blood oxygen saturation reconstruction SStO2rec in vessels on their diameter d and embedding depth zo for blood oxygen saturation value of (a) 0.6, (b) 0.7, (c) 0.8, and (d) 0.9 calculated for all considered vessels in simple geometry.
On the contrary, the algorithm reveals underestimation of StO2 for subcutaneous vessels. For estimating saturation at minimal and maximal depths among all considered values zo it is reasonable to account for the value of zo in the reconstruction algorithm for a more precise estimation of StO2.
It should also be noted that the values of ratio R are comparable for blood vessels and reconstruction artifacts (Figs. 9-10), which may lead to uncertainties in interpreting the constructed maps. The employment of illumination geometries of real OA systems and three-dimensional reconstruction instead of simplified planar wave case considered in this study has potential to eliminate the reconstruction artifacts and provide with more precise estimation of oxygen saturation in surrounding tissues as well as in blood vessels, however, it requires larger computational power.
4 Conclusion
In this study, a set of OA images of biotissue with embedded vessels was numerically simulated for different orientations, diameters, embedding depths and blood oxygen saturation of vessels. Numerical simulations are based on sequential employment of developed MC algorithm and k-Wave Toolbox. The ratio of averaged OA signals corresponding to different probing wavelengths revealed monotonous dependence on blood oxygen saturation for different pairs of the wavelengths of 532, 680, and 1064 nm. The dependence of ratio R of OA signals corresponding to probing wavelengths of 658 nm and 1064 nm on blood oxygen saturation demonstrated the smallest variance for all
diameters and depths of vessels and was employed for a simple estimation of StO2 value. The algorithm was based on the construction of inverse function StO2 (R) and voxel-by-voxel estimation of blood oxygen saturation value based on the calculated ratio R. Proposed algorithm was tested on simulated OA images for simple and complex vessel geometry. In the space of vessel diameters and embedding depths for a given true value of blood oxygen saturation there is an area of vessel diameters d and embedding depths zo that are closed for the values contributing to the average ratio dependence, which provides the reconstruction accuracy better than 3%. For the combination of vessel diameter d and embedding depth zo above that area the precision is below 15% for true StO2 values > 0.8. Three-dimensional OA reconstruction employed instead of 2-dimensional modality considered in this study has potential to improve the algorithm performance, since it allows eliminating artifacts from OA images. Moreover, depth-dependent StO2 reconstruction is expected to provide even higher accuracy. In this connection, employing a more advanced, however, time-consuming approach to OA image simulations in combination with machine learning techniques instead of the considered simple algorithm has high potential in accurate mapping of blood oxygenation based on spectral OA data.
Disclosures
The authors declare no conflict of interest.
Acknowledgements
The study is supported by Russian Science Foundation (project #22-29-00074).
References
1. C. Menon, D. L. Fraker, "Tumor oxygenation status as a prognostic marker," Cancer Letters 221(2), 225-235 (2005).
2. C. Zhou, R. Choe, N. S. Shah, T. Durduran, G. Yu, A. Durkin, D. Hsiang, R. Mehta, J. A. Butler, and A. E. Cerussi, "Diffuse optical monitoring of blood flow and oxygenation in human breast cancer during early stages of neoadjuvant chemotherapy," Journal of Biomedical Optics 12(5), 051903 (2007).
3. A. A. Tandara, T. A. Mustoe, "Oxygen in wound healing—more than a nutrient," World Journal of Surgery 28, 294300 (2004).
4. A. Orlova, Y. Perevalova, K. Pavlova, N. Orlinskaya, A. Khilov, D. Kurakina, M. Shakhova, M. Kleshnin, E. Sergeeva, and I. Turchin, "In Diffuse optical spectroscopy monitoring of experimental tumor oxygenation after red and blue light photodynamic therapy," Photonics 9(1), 19 (2022).
5. C. Balas, "Review of biomedical optical imaging—a powerful, non-invasive, non-ionizing technology for improving in vivo diagnosis," Measurement Science and Technology 20(10), 104020 (2009).
6. A. P. Dhawan, B. D'Alessandro, and X. Fu, "Optical imaging modalities for biomedical applications," IEEE Reviews in Biomedical Engineering 3, 69-92 (2010).
7. T. Durduran, G. Yu, M. G. Burnett, J. A. Detre, J. H. Greenberg, J. Wang, C. Zhou, and A. G. Yodh, "Diffuse optical measurement of blood flow, blood oxygenation, and metabolism in a human brain during sensorimotor cortex activation," Optics Letters 29(15), 1766-1768 (2004).
8. S. E. Boebinger, R. O. Brothers, S. Bong, B. Sanders, C. McCracken, L. H. Ting, and E. M. Buckley, "Diffuse optical spectroscopy assessment of resting oxygen metabolism in the leg musculature," Metabolites 11(8), 496 (2021).
9. D. Kurakina, V. Perekatova, E. Sergeeva, A. Kostyuk, I. Turchin, and M. Kirillin, "Probing depth in diffuse reflectance spectroscopy of biotissues: A monte carlo study," Laser Physics Letters 19, 035602 (2022).
10. A. Orlova, M. Y. Kirillin, A. Volovetsky, N. Y. Shilyagina, and E. Sergeeva, G. Y. Golubiatnikov, and I. Turchin, "Diffuse optical spectroscopy monitoring of oxygen state and hemoglobin concentration during skbr-3 tumor model growth," Laser Physics Letters 14, 015601 (2016).
11. A. F. Khan, F. Zhang, H. Yuan, and L. Ding, "Brain-wide functional diffuse optical tomography of resting state networks," Journal of Neural Engineering 18, 046069 (2021).
12. J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. Yodh, "Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia," Journal of Cerebral Blood Flow & Metabolism 23(8), 911-924 (2003).
13. L. D. Liao, M. L. Li, H. Y. Lai, Y. Y. I. Shih, Y. C. Lo, S. Tsang, P. C. P. Chao, C. T. Lin, F. S. Jaw, and Y. Y. Chen, "Imaging brain hemodynamic changes during rat forepaw electrical stimulation using functional photoacoustic microscopy," Neuroimage 52(2), 562-570 (2010).
14. E. Liapis, A. Karlas, U. Klemm, and V. Ntziachristos, "Chemotherapeutic effects on breast tumor hemodynamics revealed by eigenspectra multispectral optoacoustic tomography (emsot)," Theranostics 11(16), 7813 (2021).
15. A. Orlova, M. Sirotkina, E. Smolina, V. Elagin, A. Kovalchuk, I. Turchin, and P. Subochev, "Raster-scan optoacoustic angiography of blood vessel development in colon cancer models," Photoacoustics 13, 25-32 (2019).
16. P. Subochev, E. Smolina, E. Sergeeva, M. Kirillin, A. Orlova, D. Kurakina, D. Emyanov, and D. Razansky, "Toward whole-brain in vivo optoacoustic angiography of rodents: Modeling and experimental observations," Biomedical Optics Express 11(3), 1477-1488 (2020).
17. L. V. Wang, "Multiscale photoacoustic microscopy and computed tomography," Nature Photonics 3, 503-509 (2009).
18. L. V. Wang, S. Hu, "Photoacoustic tomography: In vivo imaging from organelles to organs," Science 335, 14581462 (2012).
19. X. L. Dean-Ben, D. Razansky, "Optoacoustic imaging of the skin," Experimental Dermatology 30, 1598-1609 (2021).
20. V. Perekatova, P. Subochev, M. Y. Kirillin, E. Sergeeva, D. Kurakina, A. Orlova, A. Postnikova, and I. Turchin, "Quantitative techniques for extraction of blood oxygenation from multispectral optoacoustic measurements," Laser Physics Letters 16, 116201 (2019).
21. M. Schwarz, A. Buehler, J. Aguirre, and V. Ntziachristos, "Three-dimensional multispectral optoacoustic mesoscopy reveals melanin and blood oxygenation in human skin in vivo," Journal of Biophotonics 9, 55-60 (2016).
22. M. Kirillin, V. Perekatova, I. Turchin, and P. Subochev, "Fluence compensation in raster-scan optoacoustic angiography," Photoacoustics 8, 59-67 (2017).
23. J. Gröhl, M. Schellenberg, K. Dreher, and L. Maier-Hein, "Deep learning for biomedical photoacoustic imaging: A review," Photoacoustics 22, 100241 (2021).
24. H. Deng, H. Qiao, Q. Dai, and C. Ma, "Deep learning in photoacoustic imaging: A review," Journal of Biomedical Optics 26(4), 040901 (2021).
25. C. Yang, H. Lan, F. Gao, and F. Gao, "Review of deep learning for photoacoustic imaging," Photoacoustics 21, 100215 (2021).
26. J. Gröhl, T. Kirchner, T. J. Adler, L. Hacker, N. Holzwarth, A. Hernández-Aguilera, M. A. Herrera, E. Santos, S. E. Bohndiek, and L. Maier-Hein, "Learned spectral decoloring enables photoacoustic oximetry," Scientific Reports 11, 6565 (2021).
27. J. H. Nölke, T. Adler, J. Gröhl, T. Kirchner, L. Ardizzone, C. Rother, U. Köthe, and L. Maier-Hein, "Invertible neural networks for uncertainty quantification in photoacoustic imaging," Bildverarbeitung Für Die Medizin 2021, 330-335 (2021).
28. C. Yang, F. Gao, "EDA-Net: Dense aggregation of deep and shallow information achieves quantitative photoacoustic blood oxygenation imaging deep in human breast," in International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, Cham, 246-254 (2019).
29. K. Hoffer-Hawlik, G. P. Luke, "Abso2luteu-net: Tissue oxygenation calculation using photoacoustic imaging and convolutional neural networks," ENGS 88 Honors Thesis (2019).
30. G. P. Luke, K. Hoffer-Hawlik, A. C. Van Namen, and R. Shang, "O-net: A convolutional neural network for quantitative photoacoustic image segmentation and oximetry," arXiv Preprint, 1911.01935 (2019).
31. C. Bench, A. Hauptmann, and B. T. Cox, "Toward accurate quantitative photoacoustic imaging: Learning vascular blood oxygen saturation in three dimensions," Journal of Biomedical Optics 25(8), 085003 (2020).
32. L. Wang, S. L. Jacques, and L. Zheng, "Mcml—monte carlo modeling of light transport in multi-layered tissues," Computer Methods and Programs in Biomedicine 47(2), 131-146 (1995).
33. M. Kirillin, A. Khilov, D. Kurakina, A. Orlova, V. Perekatova, V. Shishkova, A. Malygina, A. Mironycheva, I. Shlivko, S. Gamayunov, I. Turchin, and E. Sergeeva, "Dual-wavelength fluorescence monitoring of photodynamic therapy: From analytical models to clinical studies," Cancers 13(22), 5807 (2021).
34. A. V. Gorshkov, M. Y. Kirillin, "Acceleration of monte carlo simulation of photon migration in complex heterogeneous media using intel many-integrated core architecture," Journal of Biomedical Optics 20(8), 085002 (2015).
35. V. Perekatova, S. Nemirova, A. Orlova, M. Kirillin, A. Kurnikov, K. Pavlova, A. Khilov, A. Kovalchuk, and P. Subochev, "Three-dimensional dual-wavelength optoacoustic angiography reveals arteriovenous anastomoses," Laser Physics Letters 18, 045601 (2021).
36. V. Perekatova, P. Subochev, M. Kleshnin, and I. Turchin, "Optimal wavelengths for optoacoustic measurements of blood oxygen saturation in biological tissues," Biomedical Optics Express 7(10), 3979-3995 (2016).
37. N. Bosschaart, G. J. Edelman, M. C. Aalders, T. G. van Leeuwen, and D. J. Faber, "A literature review and novel theoretical approach on the optical properties of whole blood," Lasers in Medical Science 29, 453-479 (2014).
38. T. Kono, J. Yamada, "In vivo measurement of optical properties of human skin for 450-800 nm and 950-1600 nm wavelengths," International Journal of Thermophysics 40, 51 (2019).
39. B. E. Treeby, B. T. Cox, "K-wave: Matlab toolbox for the simulation and reconstruction of photoacoustic wave fields," Journal of Biomedical Optics 15(2), 021314 (2010).
40. I. Turchin, S. Bano, M. Kirillin, A. Orlova, V. Perekatova, V. Plekhanov, E. Sergeeva, D. Kurakina, A. Khilov, A. Kurnikov, P. Subochev, M. Shirmanova, A. Komarova, D. Yuzhakova, A. Gavrina, S. Mallidi, and T. Hasan, "Combined fluorescence and optoacoustic imaging for monitoring treatments against ct26 tumors with photoactivatable liposomes," Cancers 14(1), 197 (2021).
41. V. Perekatova, M. Kirillin, S. Nemirova, A. Orlova, A. Kurnikov, A. Khilov, K. Pavlova, V. Kazakov, V. Vildanov, I. Turchin, and P. Subochev, "Quantitative characterization of age-related changes in peripheral vessels of a human palm using raster-scan optoacoustic angiography," Photonics 9(7), 482 (2022).
42. P. Subochev, "Cost-effective imaging of optoacoustic pressure, ultrasonic scattering, and optical diffuse reflectance with improved resolution and speed," Optics Letters 41(5), 1006-1009 (2016).
43. V. V. Perekatova, M. Y. Kirillin, I. V. Turchin, and P. V. Subochev, "Combination of virtual point detector concept and fluence compensation in acoustic resolution photoacoustic microscopy," Journal of Biomedical Optics 23(9), 091414 (2018).
44. V. Perekatova, M. Kirillin, P. Subochev, A. Kurnikov, A. Khilov, A. Orlova, D. Yuzhakova, and I. Turchin, "Quantification of microvasculature parameters based on optoacoustic angiography data," Laser Physics Letters 18, 035602 (2021).
