Finite difference methods for solving the transport equation in the problems of optical biomedical diagnostics

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.


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 30 Apr 2017 © JBPE 010311-3 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][10][11][12][13][14][15][16][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

Problem statement
The propagation of radiation from a pulsed source in a turbid medium is described by the linear integrodifferential transfer equation [12] where The function Q( !r, !Ω, 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 Here S 0 is the source power, 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.

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 where the functions For the non-scattered component Φ 0 ( ! r, !Ω, t) For the singly scattered component Φ 1 ( ! r, !Ω, t) For the component allowing for two or more scattering events Φ( !r, !Ω, t) The solution of the problem (2.2) is found analytically using the technique of generalised functions.Similarly, the right-hand side ŜΦ 0 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).

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][10][11] and Raduga-5.2(P)for the nonstationary 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).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.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 2хAMD 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 10 9 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.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.
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, 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].

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.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 I V 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.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.

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.
photons in the medium.The extinction coefficient Σ t ( ! r ) determines the probability of interaction between the photon and the medium; the scattering coefficient Σ s ( ! r ) specifies the probability of changing the photon direction as a result of collision with atoms of the medium; the scattering indicatrice ρ( !r, !Ω ⋅ !′ Ω ) describes the probability of changing the direction of photon propagation from ! ′ Ω to !Ω .At each spatial point this function is determined by the scalar product !Ω⋅ !′ Ω , i.e., the cosine of photon deflection angle at each scattering event.

δFig. 1
Fig. 1 Geometry of the modelled source of photons.
photons, and the function Φ( !r, !Ω, 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

Fig. 2
Fig.2Scattered 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 = 0-3 mm (а) and r= 0-20 mm (b).

Fig. 5 (
Fig. 5 (а) 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.

Fig. 6 (
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(I 0 / I V ) is presented, where I V and I 0 is the diffuse reflected signal from the phantom with a branching vessel and from a similar phantom without a vessel, respectively.