Enhancing eMSOT with nonlinear mixed model approaches for precise blood oxygenation imaging in tissues

Hemoglobin oxygen saturation is clinically significant for patients with tumors.^1,2 Hypoxia, a hallmark of many cancers, is linked with tumor progression and poor prognosis.³ Oxygen saturation (sO₂) imaging plays a crucial role in identifying hypoxic regions within tumors, facilitating targeted treatments, monitoring therapeutic efficacy, and enabling real-time adjustments to treatment strategies, ultimately improving patient outcomes.^4,5 This method represents a key advancement in personalized cancer therapy.

Multispectral optoacoustic tomography (MSOT) leverages the unique absorption characteristics of chromophores across different spectral bands.^6–9 By employing pulsed lasers at specific wavelengths, MSOT selectively excites targeted tissue areas, generating unique optoacoustic signals that reflect the absorption spectra of diverse chromophores.¹⁰ This process allows for the creation of optoacoustic images through linear spectral unmixing algorithms,^11,12 which provide concentration data for different chromophores and support label-free imaging of tissue sO₂. However, the accuracy of this method is limited by spectral corruption—nonlinear changes in optical fluence with tissue depth and wavelength,¹³ which linear unmixing methods fail to account for, leading to significant errors in estimating oxygen saturation in deeper tissues.

To mitigate the effects of spectral corruption, several methodologies^14–21 have been proposed. These approaches primarily involve mathematical modeling of light propagation within biological tissues using the radiative transfer equation and related techniques.¹⁷ By doing so, the distribution of optical fluence within the tissue can be obtained.¹⁵ However, these methods necessitate precise specifications of various conditions, such as illumination and tissue boundaries, to avoid convergence issues during the solution process.^13,22,23 The inability to accurately define these boundary conditions in living tissues can significantly undermine the accuracy of the results.

Vasilis team has proposed an eigenspectra optoacoustic tomography (eMSOT) method aimed at addressing the spectral corruption effect.¹³ This method employs a linear mixing model to characterize the wavelength-dependent optical fluence, thereby transforming the optical fluence correction problem from the spatial to the spectral domain. The eMSOT approach relies solely on eigenfluence parameters within the linear model, rendering it independent of the intrinsic optical properties of the tissue, thus making it suitable for in vivo applications. Furthermore, it has successfully facilitated the quantitative estimation of sO₂ in deep tissue in near-infrared-I (NIR-I). Building upon eMSOT, a series of derived methods have been developed: including the use of spectral reliability maps to estimate the level of noise superimposed onto the recorded OA spectra,^24,25 modeling eMSOT as a Bayesian inverse problem²⁶ to enhance sO₂ estimation accuracy under high noise conditions, and integrating deep-learning techniques to directly solve for the eigenfluence parameters, thereby improving the speed and computational accuracy of sO₂ imaging.²⁷ 

Despite the advancements made by the eMSOT method in the realm of label-free imaging of oxygenation distribution, its reliance on a linear mixing model (LMM) limits its applicability by neglecting the effects of multiple scattering, where photons interact with multiple materials.^28,29 In practical situations, the LMM fails to serve as an adequate approximation. NIR-I can penetrate tissues to a depth on the order of centimeters.^30,31 At these distances, the path taken by the photons greatly exceeds the spatial scale of a single photon, leading to the occurrence of multiple scattering effects.²⁸ The photon banana of NIR-I³² describes a typical example of multiple scattering effects. During this process, the interaction of light with various components in the mixture becomes nonlinear, and the spectral characteristics cannot be simply described as a weighted average of the eigenspectra.

Consequently, more complex mixing models need to be considered to cope with nonlinear interactions. We have employed a nonlinear hybrid model to enhance eMSOT, referred to as NL-eMSOT. The nonlinear mixture model (NLMM) introduces additional degrees of freedom, aiming to portray more accurately the propagation characteristics of photons within tissues by simultaneously accounting for the interactions of both linear and nonlinear components.^{28,29,33–35} Specifically, the NLMM not only incorporates spectrally weighted averages of different substances but also considers the influence of interactions among components on the overall spectral characteristics.^36–38 This approach involves constructing a bilinear matrix to characterize the interactions of photons with different components within the mixture, thereby providing a more precise reflection of spectral characteristics in practical scenarios. To validate the accuracy of this method for quantitative imaging of sO₂, we conducted experiments in four aspects: numerical simulations, numerical simulations of biological tissues, imaging of blood phantoms, and live mouse tumor imaging. The results were then compared with those obtained using the eMSOT method. Experimental findings indicate that the NL-eMSOT method reconstructs quantitative images of sO₂ with greater accuracy and stability, both within uniform tissue interiors and at the boundaries of complex and irregular tissues.

The eMSOT model is based on an LMM. During its propagation through tissues, light experiences multiple scattering effects. Figure 1(a) illustrates the interaction of light with various chromophores, which non-linearly influences the OA spectrum. Numerous studies have demonstrated that when multiple scattering effects are present, the LMM is an inadequate modeling approach, necessitating consideration of nonlinear effects. Consequently, employing eMSOT to model light fluence may compromise the accuracy of sO₂ estimation. To enhance eMSOT, we propose a nonlinear model, as detailed below.

Before delving into the specifics, it is essential to clarify that the derivation of both eMSOT and NL-eMSOT methods is predicated on the following assumptions: In the NIR-I, Oxy and deOxy are the primary absorbers within the tissue, while the impact of other absorbers on photons is considered negligible.

The methodology for acquiring the eigenspectra in this model aligns with that of eMSOT, both utilizing PCA. By fitting the light fluence dataset to the eigenspectral model, the eigenspectral fluence parameters, denoted as m and e, are computed and expressed as functions related to tissue background oxygenation levels and tissue depth, as shown in Fig. 1(b). It is evident that, while m₂ exhibits a uniform decreasing trend with tissue depth, the other parameters do not display distinct characteristics. However, the introduction of nonlinear terms in NL-eMSOT increases the number of unknowns by three compared to the eMSOT method, thereby escalating the complexity of the solution.

The nonlinear multivariate function can be minimized using the sequential quadratic programming algorithm from the MATLAB toolbox to obtain eigenfluence parameter values for constraint grid points. By employing interpolation techniques, parameter values for all pixels can be derived. Dividing the initial acoustic pressure image by the light fluence at each position r and wavelength l yields the image corrected by the NL-eMSOT model, denoted as P^NL-eMSOT. Subsequently, linear spectral unmixing can be performed on this image to ultimately determine sO₂, as illustrated in Fig. 1(e).

To validate the precision of modeling light fluence in tissues, we employed a light propagation model to simulate the distribution of light fluence in circular tissue structures with arbitrary optical properties (radius of 1 cm), as depicted in Fig. 2. The first row presents results for highly heterogeneous tissues, while the second row illustrates results for lowly heterogeneous tissues. Figures 2(a) and 2(h) display random images of tissue background sO₂ distributions, with the simulated spatial distribution of tissue absorption coefficients shown in Figs. 2(b) and 2(i) and scattering coefficients in Figs. 2(c) and 2(j). Utilizing the finite element method to solve the diffusion equation for tissue structures with varying absorption and scattering distributions, the resulting light fluence distributions are illustrated in Fig. 2(d) and 2(k). All parameters adhere to Gaussian distributions. Subsequently, the light fluence obtained at each position in the simulated tissue is individually fitted with the eMSOT and NL-eMSOT methods, generating fitting residual distributions illustrated in Figs. 2(e), 2(f), 2(l), and 2(m), respectively. Analysis of these fitting residual distributions reveals that the residual distribution of the NL-eMSOT method is markedly more uniform and smooth and exhibits relatively smaller residual values compared to the eMSOT method. Boxplot statistical analysis of the spectral fitting residuals for NL-eMSOT is displayed in Figs. 2(g) and 2(n), highlighting smaller residual values and affirming the accuracy of the NL-eMSOT method’s forward model in modeling light fluence in tissue.

To assess the accuracy of the NL-eMSOT method for quantitative estimation of tissue sO₂, we conducted numerical simulations on biological tissues using a three-dimensional digital mouse model. This digital mouse model includes various regions such as the heart, liver, pancreas, and others, with distinctively complex and irregular boundaries that mimic real biological tissue structures. We selected a cross-section at the dashed line position shown in Fig. 3(a), as illustrated in Figs. 3(b) and 3(c), encompassing optical properties of regions of the liver, skin, and bone. Within the liver region, a circular hypoxic zone was designated to simulate a tumor’s hypoxic environment.

The results obtained using the NL-eMSOT and eMSOT methods are depicted in Figs. 3(d) and 3(e), respectively. It is evident that the NL-eMSOT method provides a more accurate estimation of sO₂ values in the tumor region compared to the eMSOT method, which significantly underestimates these values. Intensity distribution curves of the reconstructed results obtained by both methods, along with the corresponding real sO₂ image along the dashed lines (Line 1 and Line 2) shown in Fig. 3(c), are presented in Figs. 3(f) and 3(g). Notably, the intensity distribution curve of sO₂ derived from the NL-eMSOT method closely aligns with the real sO₂ distribution curve, exhibiting only minor discrepancies in the tumor region. In contrast, the eMSOT method demonstrates considerable deviations from the actual values at tissue boundaries, particularly where turning points occur in the intensity distribution curve and significantly underestimates sO₂ in the tumor region, with a relative error of up to 30%. The quantitative reconstruction results of sO₂ obtained by the eMSOT and NL-eMSOT methods were compared based on d and PSNR, as summarized in Table I.

To create blood tissue phantoms with varying oxygenation levels, we controlled the addition of sodium hydrosulfite to fresh sodium citrate anti-coagulated bovine blood. This resulted in three different blood phantoms: a blood phantom1 with a high oxygenation background without the addition of sodium hydrosulfite; a blood phantom2 with sodium hydrosulfite added to complete deoxygenation; and a blood phantom3 with a background of high oxygenation without sodium hydrosulfite, with the region of sodium hydrosulfite-added blood inserted. Figure 4 illustrates the performance of the NL-eMSOT method on the phantom experimental data. Figure 4(a) depicts the region of sO₂ imaging, while Figs. 4(b) and 4(c) represent sO₂ estimation imaging on a homogeneous phantom using the eMSOT method; Figs. 4(e) and 4(f) show sO₂ estimation imaging on a homogeneous phantom using NL-eMSOT, revealing the higher uniformity of imaging obtained by NL-eMSOT. Figures 4(d) and 4(g) display sO₂ estimation imaging from both methods on the inserted phantom, with the dashed box indicating the position of the deoxygenated phantom insertion. Experimental results from the non-uniform phantom indicate that the imaging outcomes produced by the eMSOT method exhibit irregular distortions at the boundaries of heterogeneous regions, leading to partially erroneous hypoxic areas in deeper regions of the phantom. In contrast, the sO₂ imaging results obtained through the NL-eMSOT method demonstrate excellent fidelity in both shape and quantitative accuracy across the heterogeneous regions. White arrows in Fig. 4(d) and 4(g) indicate the positions analyzed by the spectroscopy. The spectral fitting curves for the two methods are denoted in Figs. 4(h) and 4(i), respectively. Notably, the spectra generated using the eigenfluence parameters derived from NL-eMSOT are meaningful. Figure 4(j) illustrates the errors in sO₂ estimation for the insertion region and its surrounding area as assessed by both eMSOT and NL-eMSOT.

A BALB/c male nude mouse weighing 20 g and aged 4–5 weeks, bearing a 4T1 mammary tumor on its lower abdomen [as shown in Fig. 5(a)], was selected to validate the capability of NL-eMSOT in assessing sO₂ in live animals. Quantitative images of sO₂ near the tumor region were acquired using the eMSOT and NL-eMSOT methods. A constrained grid was applied to cover the tumor region and adjacent healthy tissues, as shown in Fig. 5(b). The experimental results of the eMSOT and NL-eMSOT methods are illustrated in Figs. 5(c) and 5(d), respectively.

From the images, it is evident that both methods depict a lower distribution of sO₂ at the tumor center, reflecting the hypoxic physiological characteristics typical of actual tumors, which arise from enhanced vascularization and accelerated oxygen metabolism. However, the eMSOT method significantly underestimates sO₂ values in the central tumor region, failing to align with the actual tumor oxygen environment. In healthy tissue regions, NL-eMSOT effectively reconstructs high sO₂ distribution images consistent with normal physiological conditions, while eMSOT continues to present lower sO₂ values. This discrepancy is particularly evident in hypoxic areas deeper within the tissues, mirroring results observed in blood-mimicking experiments [Fig. 4(d)].

Overall, the eMSOT method consistently underestimates sO₂ values in both tumor and deep healthy tissue regions, while NL-eMSOT excels in faithfully reconstructing sO₂ values at various locations, aligning with numerical simulations and phantom experimental outcomes. Figures 5(e) and 5(f) present the spectral fitting curves for both methods at the location indicated by the white arrow in Fig. 5(b), demonstrating the superior fitting performance of the NL-eMSOT method over eMSOT. The spectral fitting residual value is 2.85% for NL-eMSOT and 5.57% for eMSOT. In vivo mouse experiments highlight the potential of NL-eMSOT in achieving high-precision quantitative imaging in complex biological tissues, showcasing its promise for biomedical research applications.

During the experiments, we observed that the computational time required by NL-eMSOT was longer than that for eMSOT. To compare the computational resources demanded by the two methods, we utilized an i7-9700k CPU and conducted tests on digital mouse model data. The specific results of the time comparison are presented in Table II:

From the table, it is evident that as the number of computation points increases, the time required by both methods exhibits a linear growth. Although NL-eMSOT incurs a greater computational time compared to eMSOT, it demonstrates a lower computational error.

The introduction of the NL-eMSOT method provides a certain improvement in blood oxygenation imaging within heterogeneous tissues. By incorporating a nonlinear mixing model, NL-eMSOT addresses the limitations of the linear mixture model employed in eMSOT, particularly concerning multiple scattering effects and complex tissue interactions. The Hadamard product, integral to NL-eMSOT, enhances the representation of interactions among different chromophores, thereby increasing sensitivity to variations in tissue oxygenation.

Experimental results from numerical simulations, phantom studies, and in vivo mouse models show that NL-eMSOT obviously outperforms eMSOT in terms of accuracy and stability. Notably, in vivo studies, the estimation of sO₂ in healthy tissues by NL-eMSOT closely aligns with actual values reported in the literature,^39–41 highlighting its remarkable capability in addressing tissue boundaries and complex structures.

The increased accuracy in sO₂ estimation is vital for clinical applications, especially in oncology, as it can significantly impact treatment planning and prognosis. By providing more reliable data, NL-eMSOT has the potential to enhance the effectiveness of targeted therapies and personalized medicine.

However, NL-eMSOT necessitates more sophisticated modeling and greater computational resources to effectively manage its complexity and ensure efficient data processing. Consequently, future investigations should focus on developing more efficient computational methodologies. In addition, while the NL-eMSOT method currently considers interactions between two chromophores, it is essential to recognize that melanin also serves as a significant absorber within the NIR-I. Future research should incorporate this aspect to fully capture the complex dynamics of tissue interactions.

In summary, NL-eMSOT represents a substantial advancement over eMSOT methods in accurately quantifying tissue oxygenation in deep and heterogeneous tissues. Its ability to handle nonlinear interactions within tissues offers precise imaging, pushing the boundaries of optoacoustic technology and showing promise for improving clinical outcomes in oncology and beyond.

See the supplementary material for the light fluence fitting results of NL-eMSOT and eMSOT; the mathematical expression of the evaluation metrics..

This work was supported by the National Natural Science Foundation of China (Grant No. 82171989 and 62075156) and the 2022 Shenzhen College and University Stability Support Plan.

The authors have no conflicts to disclose.

The animal use protocol listed below has been reviewed and approved by the Tianjin University Animal Ethical and Welfare Committee (AEWC), Hereby certify.

P.H. and Y.L. contributed equally to this work.

Pengwei Han: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Yiping Lv: Data curation (equal); Formal analysis (equal); Resources (equal); Validation (equal); Writing – review & editing (equal). Binxue Zhang: Data curation (equal); Validation (equal); Writing – original draft (equal). Shuang Wu: Formal analysis (equal); Validation (equal); Writing – original draft (equal). Jiayue Wang: Writing – review & editing (equal). Hao Zhang: Writing – review & editing (equal). Feng Gao: Funding acquisition (equal); Supervision (equal); Writing – original draft (equal). Jiao Li: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.