CC BY
46
Finite difference methods for solving the transport equation in the problems of optical biomedical diagnostics
Leonid P. Bass1*, Olga V. Nikolaeva1, Alexander V. Bykov2, and Mikhail Yu. Kirillin3
1 M.V. Keldysh Institute of Applied Mathematics RAS, Miusskaya pl., 4, Moscow 125047, Russia
2 Opto-Electronics and Measurement Techniques Lab., University of Oulu, P.O. Box 4500, Oulu 90570, Finland
3 Institute of Applied Physics RAS, Ulyanov St., 46, Nizhny Novgorod 603950, Russia
* e-mail: bass@kiam.ru
Abstract. The present paper is a retrospective review of the development of the approach to numerical modeling of radiation propagation in biological tissues, based on the use of finite difference methods for the solution of the radiative transfer equation (RTE) and its application to the problems of biomedical diagnostics. The advantage of finite difference methods is the possibility to obtain the solution of the direct problem of radiation propagation in turbid media, when the exact analytical solution is impossible. In turn, the possibility to solve the RTE opens wide perspectives for the solution of inverse problem and reconstruction of structural and functional characteristics of objects from the detected scattered radiation. We present a review of the finite difference methods applied to the solution of angiography problems and investigation of blood using optical biomedical diagnostic technologies. © 2017 Journal of Biomedical Photonics & Engineering.
Keywords: radiative transfer equation; finite difference techniques; Monte Carlo simulations; stationary problem; non-stationary problem.
Paper #3176 received 5 Apr 2017; revised manuscript received 14 Apr 2017; accepted for publication 17 Apr 2017; published online 30 Apr 2017. doi: 10.18287/JBPE17.03.010311. [Special Issue. Years in Biophotonics: 70th Anniversary of Prof. A.V. Priezzhev].
References
1. A. Ishimaru, Wave Propagation and Scattering in Random Media, Wiley-IEEE Press (1999).
2. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, "Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction," Commun. Numer. Methods Eng. 25(6), 711-732 (2009).
3. H. Dehghani, B. R. White, B. W. Zeff, A. Tizzard, and J. P. Culver, "Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography," Appl. Optics 48(10), D137-D143 (2009).
4. C. Zhu, and Q. Liu, "Review of Monte Carlo modeling of light transport in tissues," J. Biomed. Opt. 18(5) 050902 (2013).
5. L. Wang, S. L. Jacques, and L. Zheng, "MCML—Monte Carlo modeling of light transport in multilayered tissues," Comput. Methods Prog. Biomed. 47(2), 131-146 (1995).
6. C.-C. Chuang, Y.-T. Lee, C.-M. Chen, Y.-S. Hsieh, T.-C. Liu, and C.-W. Sun, "Patient-oriented simulation based on Monte Carlo algorithm by using MRI data," Biomed. Eng. Online 11, 21 (2012).
7. A. V. Gorshkov, and M. Yu. 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).
8. H. Li, J. Tian, F. Zhu, W. Cong, L. V. Wang, E. A. Hoffman, and G. Wang, "A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo method," Acad. Radiol. 11(9), 1029-1038 (2004).
9. A. V. Bykov, A. V. Priezzhev, A. A. Dergachev, L. P. Bass, O. V. Nikolaeva, and V. S. Kuznetsov, "Simulation of light propagation in strongly scattering media modelling biotissues with the use of different
algorithms," Proceedings of the II Eurasian Congress on Medical Physics and Engineering, Moscow, Russia, 21-24 June 2005, 182 (2005) [in Russian].
10. A. V. Bykov, A. V. Priezzhev, L. P. Bass, V. S. Kuznetsov, and R. A. Myllyla, "Light propagation in highly scattering media mimicking biotissues: comparison of different algorithm," Proc. of Forth International Conference on Photonics and Imaging in Biology and Medicine, Tianjin, China, 3-6 September 2005, 1-8 (2005).
11. L. P. Bass, O. V. Nikolaeva, V. S. Kuznetsov, A. V. Bykov, A. V. Priezzhev, and A. A. Dergachev. "Modelling of propagation of optical radiation in a phantom of biological tissue using the supercomputer MBC1000/M," Matematicheskoye modelirovaniye 18(1), 29-42 (2006) [in Russian].
12. V. S. Kuznetsov, O. V. Nikolaeva, L. P. Bass, A. V. Bykov, and A. V. Priezzhev, "Modelling of ultrashort light pulse propagation through a strongly scattering medium," Matematicheskoye modelirovaniye 21(4), 3-14 (2009) [in Russian].
13. L. P. Bass, V. S. Kuznetsov, O. V. Nikolaeva, A. V. Bykov, and A. V. Priezzhev, "Parallel algorithm of the discrete ordinate method to laser impulse propagation simulation in turbid media," Proc. of 21th International Conference on Transport Theory, Torino, Italy, July 12-17, 9 (2009).
14. L. P. Bass, V. S. Kuznetsov, O. V. Nikolaeva, A. V. Bykov, and A. V. Priezzhev, "Modelling of ultrashort laser pulse propagation in biotissue in application to problems of non-invasive biomedical diagnostics," Proc. of Progress in Electromagnetics Research Symposium, 18-21 August, 562 (2009).
15. L. P. Bass, O. V. Nikolaeva, V. S. Kuznetsov, A. V. Bykov, and A. V. Priezzhev, "Parallel algorithms for simulation of ultrashort pulse propagation in turbid media," Nuovo cimento 33C(1), 39-46 (2010).
16. O. V. Nikolaeva, L. P. Bass, A. V. Bykov, V. S. Kuznetsov, E. M. Loskutov, A. V. Priezzhev, and Yu. V. Yakhno, "Solution of direct and inverse problem in the diagnostics of blood vessels using the computers with parallel architecture," Proceedings of the Conference on Parallel Computing Technologies (PAVT -2011), Moscow, March 28 - April 1, 768 (2011) [in Russian].
17. L. P. Bass, O. V. Nikolaeva, A. V. Priezzhev, A. V. Bykov, and Yu. V. Yakhno, "Using supercomputers to calculate light fields reflected from strongly scattering biotissues for the visualisation of blood vessels and blood sugar assessment," The 1st International Science and Technical Conference "Computer Biology - from Fundamental Science to Biotechnology and Biomedicine." Pushchino, Moscow Region, December 7-9, 2011, 14-16 (2011) [in Russian].
18. A. D. Klose, U. Netz, J. Beuthan, and A. H. Hielscher, "Optical tomography using the time-independent equation of radiative transfer — Part 1: forward model," Journal of Quantitative Spectroscopy & Radiative Transfer 72(5), 691-713(2002).
1 Introduction
The progress of modern optical biomedical diagnostics requires the development of models of probing radiation propagation in a biological object. In the radiative transfer theory, they distinguish between the direct problem, i.e., the calculation of the scattered intensity distribution basing on the object geometry and optical properties, and the inverse problem, i.e., the reconstruction of the optical and geometrical characteristics of the object from the measured scattered radiation. The solution of the inverse problem is the base of optical biomedical diagnostics.
A specific feature of biological objects is their complex geometry, characterised by different optical inhomogeneities. The basic equation of the radiative transfer theory (RTE) generally has no analytical solution for arbitrary geometry. Since biotissues are mostly strong scattering media for the radiation of optical range, the diffusion approximation is commonly used for the RTE, describing the propagation of radiation [1]. In the diffusion approximation, the RTE has analytical solutions for directional and isotropic sources of light and for infinite or semi-infinite homogeneous media. However, the use of diffusion approximation suffers from a number of restrictions; in
particular, it is hardly applicable at the distances from the source smaller than the transport length. For the arbitrary geometry the RTE in the diffusion approximation has also no analytical solution, however, the finite difference technique can be efficiently used for solving it [2]. Such approach has found wide applications in optical diffusion tomography, in particular, for the noninvasive diagnostics of breast cancer and functional diagnostics of human brain [2, 3]. An alternative approach is the numerical simulation using the stochastic Monte Carlo method [4]. This method is based on multiple calculation of random photon trajectories in the medium having the specified geometry and optical characteristics, followed by the statistical processing of the obtained results at the detection points [5, 6]. In particular, the Monte Carlo method was used to simulate the propagation of probing radiation in human head [6, 7] and laboratory animal [8]. The geometry obtained by MRT diagnostics can be used in this case as the initial data. However, both of the above approaches have their drawbacks. The numerical solution of RTE in the diffusion approximation does not provide sufficient accuracy near the source, and the Monte Carlo simulation is often accompanied by the statistical noise, the elimination of which requires high
computational expenses, which is not always acceptable, e.g., in solving inverse problems. In such cases, it is efficient to use finite-difference technique for solving RTE in its original form rather than in the diffusion approximation. The present paper presents a retrospective review of the approaches to the RTE solution by means of finite difference methods, including the cases of complex geometry [9-17]. The validation of approaches was based on the comparison of the results with those of Monte Carlo simulation. The developed approaches were used to solve the problems of angiography and optical diagnostics of blood. The finite difference methods can be also efficiently used to solve the direct problem in optical tomography [18].
2 Finite difference methods for solving the radiative transfer equation
2.1 Problem statement
The propagation of radiation from a pulsed source in a turbid medium is described by the linear integro-differential transfer equation [12]
1 ! ! - ! ! ! ! !
-—+ L¥(r,Q,t) = S¥(r,Q,t) + Q(r,Q,t), r eG, (1.1)
v dt
where
L¥(r,Q, t) = Q ■ grad ¥(r,Q ,t) + Xt (r) ¥(r,Q, t), S ¥(r,Q,t) = Es( r) JQp( r,Q Q')¥(r,Q',t)dQ
Here the function r , Q ,t) describes the radiation flux density at the point r of the spatial domain G along the
direction Qefl . The quantity v is the velocity of photons in the medium. The extinction coefficient Xt(r) determines the probability of interaction between the photon and the medium; the scattering coefficient X (r) specifies the probability of changing the photon direction as a result of collision with atoms of the medium; the scattering indicatrice p(r,Q-Q') describes the probability of changing the direction of photon propagation from Q' to Q. At each spatial point this function is determined by the scalar product Q-Q', i.e., the cosine of photon deflection angle at each scattering event.
The function Q(r,Q,t) determines the density of the source of photons that emits into a spatially limited cone at the initial moment of time (see Fig. 1). This function
Q(r,Q,t) has the form
Q(r,t) =
xn
2n(1 - cosa)
■S(r - rn) x
I r - r | | r - r |
I 0 1 I ax 0 1
- cosa
(1.2)
■S(t - to).
Here S0 is the source power,
S0 = J dt J dQJ dr Q(r, Q, t);
0
S is the Dirac delta function, providing the photon emission from the point rg at the time moment t0; the function 77(a) specifies the radiation cone
77(a) = 0 for a < 0, 77(a) = 1 for a > 0;
a is the opening angle of the source, measured from the cone axis defined by the vector r , see Fig. 1.
Fig. 1 Geometry of the modelled source of photons.
At the boundary, on which the radiation from the source is incident, the boundary condition for Eq. (1.1) is chosen as the total reflection condition. It is also assumed that the photons are free to escape from the calculation domain through the rest boundaries.
2.2 Solution approach
The source singularity makes it impossible to use the finite difference approximations directly for the solution of the considered problem. Therefore, the desired solution is presented as a sum of three functions
Q, t) = O0 (!, Q, t) + (!, Q, t) + O(!, Q, t), (2.1)
where the functions O0(r,Q,t), O^r,Q,t) describe the flux density of non-scattered and singly scattered photons, and the function 0(r, Q,t) describes the flux density of photons that undergo two or more scattering events. Substituting the representation (2.1) into Eq. (1.1) and using the additivity property of the transfer
r-r
r - r
z
equation, we arrive at the boundary problem for each component of the solution.
For the non-scattered component O0(r, Q,t)
1 dO
v + ^ O = Ö'
(rÂt )| ö
= 0.
(2.2)
For the singly scattered component 01(r, Q,t)
1 dO - !! I
-—^ + LOi = SO0,01(r,Q,t)l v dt 1 0 1
= 0.
(2.3)
For the component allowing for two or more scattering events 0(r, Q,t)
1 dO - - ! ! I
-—+ LO = SO + SO,, 0(r,Q,t) = 0. (2.4)
v dt <o
The solution of the problem (2.2) is found analytically using the technique of generalised
functions. Similarly, the right-hand side ,SO0 in Eq.
(2.3) and the solution of this equation are found.
The solution of Eq. (2.4) has already no singularity and, therefore, is found using the integro-interpolational finite difference scheme. The system of finite-difference equations approximating Eq. (2.4) is solved using the iteration process for the collisions. The integrals (the operator S) in the right-hand side of Eq. (2.4) are replaced with quadrature sums.
Note, that since the problem is non-stationary, great computational load is required for time intervals (steps). To reduce the computation time one can use parallel algorithms. The parallelisation of the calculations can be implemented by means of the widely used method of spatial subdomains.
If the problem is stationary, then the first term in the left-hand side of Eqs. (1.1), (2.2)-(2.4) is absent, and the algorithm of solution corresponds to a single time step in the problem (1.1).
3 Comparison of finite difference methods with Monte Carlo method
The algorithm presented above was implemented in the M.V. Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, in the form of the software package Raduga-5.1(P) for the stationary problem [9-11] and Raduga-5.2(P) for the non-stationary one [12-15]. Ref. [11] presents the comparison of the spatial distribution of the scattered intensity calculated using the finite difference method and the Monte Carlo methods for the stationary problem (Fig. 2a, b).
104-
103,
10-
10'
Experiment Monte Carlo
Finite difference approach Diffusion approximation
0.0 0.5
1.0 1.5 2.0 r , mm
2.5 3.0
1041
103]
3 102j
TO
^
to 10'n
0
10°]
10-1]
10-2-
• Experiment
-Monte Carlo
--- Finite difference approach
■ ■ ■ ■ Diffusion approximation
0 2 4 6 8 10 12 14 16 18 20 r , mm
Fig. 2 Scattered intensity versus the distance r from the source, experimentally measured in a biotissue phantom and calculated using the finite difference method (RADUGA-5.1(P) software), the Monte Carlo method, and the diffusion approximation within the ranges r = 03 mm (a) and r= 0-20 mm (b).
Fig. 2 shows that the results of Monte Carlo simulation and the solution of the transfer equation using the finite difference method agree well with the results of experiment in the biotissue model, whereas the solution of diffusion equation is unsatisfactory in the vicinity of the source (Fig. 2a), which gives rise to large errors in attempts to reconstruct the medium properties using the solution of RTE in the diffusion approximation.
Besides that, at large distances from the source the Monte Carlo solution demonstrates oscillations (Fig. 2,b), which are absent in the finite-difference solution of the RTE. These oscillations are due to insufficient number of photon trajectories contributing to the signal.
The comparison of finite-difference and Monte Carlo solutions of a non-stationary problem was presented in Ref. [15]. Fig. 3 shows the results for a delta-pulse, backscattered by the medium, calculated using these methods. One can see good agreement of the results.
a
5.0x10-
4.0x10-
3.0x10"
2.0x10"
1.0x10-
0.0x10°
Monte-C arlo RADUGA-5.2(P)
8 12 time (ps)
16
20
Fig. 3 Comparison of numerical simulations of the medium response to a delta-pulse from a monodirectional point source using the finite difference method (RADUGA-5.2(P) software) and the Monte Carlo method.
To implement the finite difference method, 49 processors 2xAMD Opteron 248 (2.2 GHz) were used; one calculation took 25 minutes. For the validation of the developed algorithm and assessment of its efficiency, the solution of the same problem was
obtained using the Monte Carlo method with 109 trajectories. The Monte Carlo calculation took 10 minutes using 100 Intel Xeon processors (2.8 GHz), so that the both approaches required comparable computational time expenses.
0.005
monodirectional source conic source
Fig. 4 Comparison of numerical simulations of the response of the medium, modelling human skin with the glucose content 500 mg/dl, to the delta-pulse from monodirectional point and conical (focused) radiation sources.
such approach is not applicable to focused beams, since in this case the photons have different initial direction of propagation. In Refs. [12, 13] the finite difference solution of the non-stationary RTE was implemented for the conical-shaped beam. The results for the monodirectional beam and the conical one are compared in Fig. 4 [13].
4 Application of finite difference methods of solving RTE to the problems of optical biomedical diagnostics
The main goal of Refs. [13, 15] was to use numerical simulation for studying the efficiency of time-resolved optical probing of skin aimed at the assessment of glucose content in blood. The typical responses for skin probing with a delta-pulse (the probing geometry is presented in Fig. 5a), the glucose concentration being 0 and 500 mg/dl, are presented in Fig. 5b. It is shown that the difference in the maxima of the recorded pulses amounts to 13%, which is sufficient for in vitro measurements. However, for the in vivo probing this difference can be masked by the typical spread of the skin optical characteristics.
0.005
0.004
0.003
0.002
0.001
0.000
— 0 mg/dl
- - 500 mg/dl
2 4 6 8 10 12 14 16 Time, ps
One more advantage of the considered finite difference approach is the possibility to choose the light beam geometry. In most problems, they consider point monodirectional beams, and the solution for a beam with complex geometry in the case when all rays are similarly directed can be obtained, e.g., by the convolution of the solution for a monodirectional point beam with the beam function (see, e.g., [7]). However,
Fig. 5 (a) Geometry of probing a homogeneous biotissue phantom aimed at the blood sugar measurement by means of the pulsed conical (focused) laser radiation and ring detector. (b) Numerical simulation of the human skin model response to the delta-pulse of conical (focused) light source. The curves correspond to the glucose content of 0 and 500 mg/dl.
0
4
a
0
In Refs. [16, 17] the results of numerical simulation using the Monte Carlo method and finite difference methods were used to develop an approach to the solution of inverse problem of blood vessels diagnostics with the application of a neuron network. Fig. 6a [16] schematically shows the biotissue phantom with vessels, and Fig. 6b illustrates the result of simulating the continuous-wave laser radiation, reflected from this phantom. It is worth nothing that the signal IV from the phantom with a vessel enters the denominator of the function F, shown in Fig. 6b, which provides a clearer result for the vascular diagnostics.
%
E E o 00
Ô /
Source
15mnrh
d=1mm
a=26° 5=0.5 mm
The purpose of Ref. [17] was to develop an approach to the automated determination of the diameter and depth of location of blood vessel in the model biotissue. For this aim using the numerical methods, the signals of spatially resolved reflectometry were simulated for different positions and diameters of the vessel. The results were used as a learning sample for the neuron network with the number of neurons from 2 to 5, depending on the desired accuracy of reconstruction. It was shown that the relative error of determining the diameter and the location depth of the vessel by means of the neuron network does not exceed 2% for the considered model of a cylindrical vessel in homogeneous medium.
5 Conclusion
The finite difference methods of solving the radiative transfer equation (RTE) provide an efficient tool for the solution of the problem of optical-range radiation propagation through biological tissues, which allows the solution of both direct and the inverse problem of optical diagnostics. The comparison of finite difference methods with Monte Carlo simulations (MC) showed good agreement of the obtained results, the finite difference methods providing smoother curves as compared to the MC method. The finite difference methods have been applied to the problems of optical biomedical diagnostics for the solution of both stationary and non-stationary RTE for different beam geometries.
Acknowledgments
The work was carried out within the frameworks of the Governmental project for the Institute of Applied Physics of the Russian Academy of Science (Project No. 0035-2014-0008).
Fig. 6 (a) Schematic top view of biotissue phantom with a branching blood vessel. The phantom thickness is 20 mm, the blood vessel axis is at the depth 1 mm. (b) Surface distribution of the radiation intensity diffusely reflected from the considered phantom. The 2D distribution of F = ln(/0 / IV) is presented, where IV and I0 is the diffuse reflected signal from the phantom with a branching vessel and from a similar phantom without a vessel, respectively.
