Computational Laser Spectroscopy in a Biological Tissue

Abstract

We present a numerical spectroscopic study of visible and infrared laser radiation in a biological tissue. We derive a solution of a general two-dimensional time dependent radiative transfer equation in a tissue-like medium. The used model is suitable for many situations especially when the external source is time-dependent or continuous. We use a control volume-discrete ordinate method associated with an implicit three-level second-order time differencing scheme. We consider a very thin rectangular biological-tissue-like medium submitted to a visible or a near infrared light sources. The RTE is solved for a set of different wavelength source. All sources are assumed to be monochromatic and collimated. The energetic fluence rate is computed at a set of detector points on the boundaries. According to the source type, we investigate either the steady-state or transient response of the medium. The used model is validated in the case of a heterogeneous tissue-like medium using referencing experimental results from the literature. Also, the developed model is used to study changes on transmitted light in a rat-liver tissue-like medium. Optical properties depend on the source wavelength and they are taken from the literature. In particular, light-transmission in the medium is studied for continuous wave and for short pulse.

1. Introduction

Diffuse light imaging and spectroscopy aims to investigate tissue physiology in subsurface. There is a spectral window existing within tissues in the near-infrared. The absorption of hemoglobin and water is small in the near-infrared, but elastic scattering from organelles and other microscopic interfaces is large. In tissue optics, any measuring is devoted for either spectroscopy or imaging. Spectroscopy is useful for measurement of time-dependent variations in the absorption and scattering of large tissue volumes. Many spectroscopic methods involving laser interactions in scattering media were employed to study different types of biological tissue [1–3]. Photon propagation in a biological tissue is affected by absorption and scattering. A computational light propagation model which describes the interaction of photons with scattering and absorbing media is essential in biomedical optics. It is useful for setting optical tissue properties and for the development of optical imaging algorithms [4–8].

The purpose of this paper is to produce a computational Laser spectroscopic study in a biological tissue. We aim to examine light-transmission through a biological tissue for different wavelength sources. Our concern is with modeling light propagation using the radiative transfer theory. We derive a model that unifies the basic applications of diffuse optical tomography which differ essentially by the time behavior of the external source. In particular, we will work with two choices for the external source, continuous wave and short-pulse. The simplest and easiest source to use is the continuous wave (CW). In this case, the source amplitude is constant, and the transmitted amplitude is measured as a function of source-detector separation or wavelength. The second type is a short-time-pulse. In this scheme, a short pulse is launched into the medium, and the temporal point spread function of the transmitted pulse is measured. First, we must say a few words about another way to model propagation of visible and near-infrared light in biological media which is the diffusion equation [9–11]. It is a limit behavior of the photon transport equation in a very high scattering media. However, low scattering and highly absorbing regions could not be accurately modeled by diffusion equation. Also, a forward model has to account for a hyperbolic propagation of the collimated laser pulse at the finite velocity of light in the medium, and diffusion approximations family fails to describe accurately this phenomenon [12]. For a general biological tissue, it would be more convenient to use a forward model based on radiative transfer equation (RTE). So, providing solutions to the RTE is a crucial research task especially in this domain. Concerning this problem, many works have used the RTE model with a variety of numerical techniques [13–17]. Most of these works are devoted to the numerical treatment of situations where a unique choice for the external source is considered. When the source is continuous the time-independent RTE is used. The supply of short-pulsed external source requires the use of a time-dependent RTE model. In all cases there, are many difficulties in implementing these models. To overcome some of these difficulties, many techniques of RTE-solution have been used either in time-independent or in time-dependent cases. Klose and Hielsher [18] have developed an upwind-difference discrete ordinate formulation of time-independent RTE. The authors tested their model with experimental data obtained from homogeneous tissue-phantoms and from phantoms that contain void-like regions. In another paper, Klose and Larsen [19] have used a simplified spherical harmonics equation to approximate the equation of radiative transfer for modeling light propagation in biological tissues. In that work, authors tried to demonstrate that the 𝑆𝑃_𝑁 equation provides accurate solutions for a time-independent RTE. Boulanger and Charette [20] have used an upwind advection scheme to treat a propagation wave in homogenous and heterogeneous media submitted to a short-pulsed laser source. They proved that this technique can overcome a difficulty encountered in optical spectroscopy and due to high scattering regions especially in biological tissues. The most complete account of this problem is found in the works of Alfano group which has produced many papers in this field. Among these, we can cite papers [21, 22] where Cai et al. developed an analytical cumulant solution of the time-dependent radiative transfer equation in an infinite uniform medium with an arbitrary phase function. The same group of authors has made improvements in this approach especially by separating the ballistic component from the scattered component of light in [23, 24]. Such improvements allowed the authors to consider complicated situations like a polarized short pulse of light incident on a scattering medium. They even investigated backscattering of circularly polarized light from turbid media. Their theoretical results were validated by experimental data. Recently Pu et al. [25] investigated time-resolved and spectral measurements of the SW emission from cancerous and normal prostate tissues. This contributes to a more comprehension of Laser-tissues interactions and offers new fields to the radiative transfer application.

In this work, we derive a control volume discrete-ordinate method associated with an implicit three-level second-order time differencing scheme to solve a two-dimensional time-dependent radiative transfer equation in a general biological tissue medium. We investigate many cases concerning the properties of the medium and the characteristics of the source. We will prove that the used model suits for all these applications. First, the proposed model is tested in the case of a heterogeneous tissue-like phantom illuminated by a continuous wave source. Experimental results are taken from [18]. Numerical results concerning the steady-state transmitted light on the boundary are compared with measurements. Next, the proposed model is used to establish a spectroscopic study with a continuous wave source. Results concerning transmitted light through a rat-liver tissue-like medium are presented for different wavelengths. Also, we present time-dependent spectroscopy for the same medium submitted to a one short-pulse. Results concerning transmitted light and internal light distribution are presented for some source wavelengths. Before all, theoretical and numerical background of our simulations is presented below.

2. General Time-Dependent RTE

2.1. Physical Model

Light propagation in a human biological tissue is modeled by using a general time-dependent RTE [26]:

𝑛𝑐𝜕𝐼_𝜆(𝐫,Ω,𝑡)𝜕𝑡+Ω⋅∇𝐼_𝜆(𝐫,Ω,𝑡)+(𝜇_𝑎(𝐫)+𝜇_𝑠(𝐫))𝐼_𝜆(𝐫,Ω,𝑡)=𝜇_𝑠0^2𝜋𝑝(Ω,Ω^′)𝐼_𝜆(𝐫,Ω^′,𝑡)𝑑Ω^′+𝑆_𝜆(𝐫,Ω,𝑡). (1) In this equation,

(i)𝐼_𝜆(𝐫,Ω,𝑡)is the monochromatic directional energetic radiance at the spatial position-vector 𝐫=(𝑥,𝑦,𝑧). When used in a dimensional form, the radiance 𝐼_𝜆(𝐫,Ω,𝑡)is often expressed in units of W·cm^-2sr^-1. The photon direction is described using a unitary direction vectorΩ;(ii)𝜇_a(𝐫)and𝜇_s(𝐫) are the absorption and the scattering coefficients, respectively. They are both given in units of cm^-1. 𝑐 is the speed of light in vacuum. 𝑛 is the medium’s refractive index. Also we make use of reduced scattering coefficient: 𝜇^′_s=(1−g)𝜇_s; (2)(iii)The source term 𝑆_𝜆(𝐫,Ω,𝑡)is an injected monochromatic radiance at the medium boundary; (iv)The phase function 𝑝(Ω,Ω^′)describes the probability that during a scattering event a photon with directionΩ^′ is scattered into the directionΩ. In this paper, we use a commonly-applied phase function in tissue optics. It is the Henyey-Greenstein phase function [27, 28]:

𝑝(cos(𝜃))=1−𝑔²4𝜋(1+𝑔²−2𝑔cos(𝜃))^2/3, (3) where 𝑔 is the anisotropy factor and 𝜃 is the angle between Ω andΩ^′. The integral of radiance over all angles Ω at one point 𝐫 yields the fluence rate:

Φ_𝜆(𝐫,𝑡)=𝐼_𝜆(𝐫,Ω,𝑡)𝑑Ω. (4) In order to lighten the equations, the subscript λ will be omitted in all the rest of this paper.

2.2. Boundary Condition

The radiance in the boundary is the sum of the external source contribution and the partly reflected radiance due to the refractive index mismatch. This gives

𝐼(𝐫_𝐛,Ω,𝑡)=𝑆(𝐫_𝐛,Ω,𝑡)+𝑅⋅𝐼(𝐫_𝐛,Ω_ref,𝑡),𝑛⋅Ω<0,Ω_ref⋅𝑛=−Ω⋅𝑛, (5) where 𝐫_𝐛 is a vector position on the boundary and 𝑛 is the unit outer normal vector. The reflectivity 𝑅 is given for each direction using a Fresnel’s relationship: 1𝑅=2(𝑛_𝑚cos𝜑_tra−𝑛₀cos𝜑_ref𝑛_𝑚cos𝜑_tra+𝑛₀cos𝜑_ref)²+12(𝑛₀cos𝜑_ref−𝑛_𝑚cos𝜑_tra𝑛₀cos𝜑_ref+𝑛_𝑚cos𝜑_tra)²if𝜑_ref<𝜑_crandelse,𝑅=1, (6) where 𝜑_refand 𝜑_tra are the reflection and the transmission angle respectively. In the air, the transmission angles, is calculated through the following Snell’s law:

𝑛_𝑚sin𝜑_ref=𝑛₀sin𝜑_tra. (7) The medium refractive index 𝑛_𝑚 depends on the tissue and 𝑛₀ is taken as unity. The critical angle is determined by the following equation:

sin𝜑_cr=𝑛₀𝑛_𝑚. (8) Thus, the transmitted fluence rate at the boundary can be computed in each detector point using the following integral formula:

Φ_tra(𝐫_𝐛,𝑡)=_Ω⋅𝑛>0(1−𝑅)𝐼(𝐫_𝐛,Ω,𝑡)𝑑Ω. (9) The fluence rate on the boundary which enters the detector for a given aperture AP is computed as.

Φ_𝑑(𝐫_𝐛,𝑡)=_Ap(1−𝑅)𝐼(𝐫_𝐛,Ω,𝑡)𝑑Ω (10) In this work, Φ_𝑑 will be called detected fluence rate and we will take AP = 45° in all our investigations.

3. Solution Method and Numerical Implementation

3.1. Control Volume Discrete Ordinates Method: CVDOM

In this paper, we use a control volume discrete-ordinate method [29, 30] associated with an implicit three-level second-order time differencing scheme to solve (1) in a general biological tissue medium. At first, this requires the integration of (1) on an elementary control volume noted ▵𝑉. The spatial calculus domain is divided into a set of 𝐼×𝐽 elementary uniform volumes ▵𝑉 with unitary depth. This yields a calculus domain subdivided in 𝐼∗𝐽 rectangular cells. Each (𝑖,𝑗)-cell is Δ𝑥▵𝑦-sized. So, we can write.

𝑛𝑐𝜕𝐼_𝑃,𝑚𝜕𝑡+Δ𝑦⋅𝜉_𝑚(𝐼_𝐸,𝑚−𝐼_𝑊,𝑚)+Δ𝑥⋅𝜂_𝑚(𝐼_𝑁,𝑚−𝐼_𝑆,𝑚)+Δ𝑥⋅Δ𝑦(𝜇_𝑎+𝜇_𝑠)𝐼_𝑃,𝑚=Δ𝑥⋅Δ𝑦(𝑆_𝑃,𝑚+𝜇_𝑠^𝑀_𝑚′=1𝑤_𝑚′𝑝_𝑚,𝑚′𝐼_𝑃,𝑚′), (11) where _𝜉𝑚=cos(Ω_𝑚⋅𝑛_𝑥)and𝜂_𝑚=cos(Ω_𝑚⋅𝑛_𝑦)with 𝑛_𝑥and𝑛_𝑦being two normal vectors to the (𝑥,𝑧)-plane and the (𝑦,𝑧)-plane, respectively. 𝑤_𝑚 is a weighting factor depending on the chosen quadrature formula. In this work, we use a constant weight quadrature. After angular discretisation, the phase function term becomes

𝑝_𝑚𝑚′=1−𝑔²4𝜋(1+𝑔²−2𝑔(𝜉_𝑚𝜉_𝑚′+𝜂_𝑚𝜂_𝑚′))^3/2. (12) If the direction cosines are positive, the directional energetic radiances are known on the faces 𝑊 and 𝑆 and they are unknown on the faces 𝐸 and 𝑁 of the (𝑖,𝑗)-cell and also in the centre 𝑃. Therefore, we need two complementary relations to eliminate 𝐼_𝐸,𝑚 and 𝐼_𝑁,𝑚.These relations can be obtained by using interpolation formula:

𝐼_𝑃,𝑚=𝛼𝐼_𝐸,𝑚+(1−𝛼)𝐼_𝑊,𝑚,𝐼_𝑃,𝑚=𝛼𝐼_𝑁,𝑚+(1−𝛼)𝐼_𝑆,𝑚, (13) where 𝛼 is an interpolation parameter. Using the above formulae, (11) can be written as

𝑛𝑐𝜕𝐼_𝑃,𝑚+𝜕𝑡Δ𝑦⋅𝜉_𝑚𝛼(𝐼_𝑃,𝑚−𝐼_𝑊,𝑚)+Δ𝑥⋅𝜂_𝑚𝛼(𝐼_𝑃,𝑚−𝐼_𝑆,𝑚)+Δ𝑥⋅Δ𝑦(𝜇_𝑎+𝜇_𝑠)𝐼_𝑃,𝑚=Δ𝑥⋅Δ𝑦(𝑆_𝑃,𝑚+𝜇_𝑠^𝑀_𝑚′=1𝑤_𝑚′𝑝_𝑚,𝑚′𝐼_𝑃,𝑚′). (14) For time discretization, an implicit three-level second-order time differencing scheme [31] is used:

𝜕𝐼_𝑃,𝑚=𝜕𝑡3𝐼^𝑛+1_𝑃,𝑚−4𝐼^𝑛_𝑃,𝑚+𝐼^𝑛−1_𝑃,𝑚2▵𝑡,𝑛=1,…𝑛_max, (15) whith Δ𝑡 being the discrete time step. So we can rewrite (14) as follows:

𝐼^𝑛+1_𝑃,𝑚(3𝑛+2𝑐▵𝑡Δ𝑦⋅𝜉_𝑚𝛼+Δ𝑥⋅𝜂_𝑚𝛼+Δ𝑥⋅Δ𝑦(𝜇_𝑎+𝜇_𝑠=))4𝑛𝐼2𝑐▵𝑡^{𝑛𝑃,𝑚}−𝑛𝐼2𝑐▵𝑡^{𝑛−1𝑃,𝑚}+Δ𝑦⋅𝜉_𝑚𝛼𝐼^𝑛+1_𝑊,𝑚+Δ𝑥⋅𝜂_𝑚𝛼𝐼^𝑛+1_𝑆,𝑚+Δ𝑥⋅Δ𝑦(𝑆_𝑃,𝑚+𝜇_𝑠^𝑀_𝑚′=1𝑤_𝑚′𝑝_𝑚,𝑚′𝐼^𝑛+1_𝑃,𝑚′). (16) Theoretically, if we know 𝐼_𝑃,𝑚 on the (𝑖,𝑗)-cell, we can obtain results over the cells (𝑖+1,𝑗) and (𝑖,𝑗+1) by using the boundary conditions and the following relations:

𝐼_𝑊,𝑚(𝑖+1,𝑗)=𝐼_𝐸,𝑚𝐼(𝑖,𝑗),𝑖=1,…,𝐼−1,_𝑆,𝑚(𝑖,𝑗+1)=𝐼_𝑁,𝑚(𝑖,𝑗),𝑗=1,…,𝐽−1. (17) If (𝜉_𝑚>0and𝜂_𝑚>0), then we get the following equation:

𝐼^𝑛+1_{𝑖,𝑗,𝑚}=[3𝑛+2𝑐▵𝑡Δ𝑦⋅𝜉_𝑚𝛼+Δ𝑥⋅𝜂_𝑚𝛼+Δ𝑥⋅Δ𝑦(𝜇_𝑎+𝜇_𝑠)]⁻¹×[4𝑛𝐼2𝑐▵𝑡^{𝑛𝑖,𝑗,𝑚}−𝑛𝐼2𝑐▵𝑡^{𝑛−1𝑖,𝑗,𝑚}+Δ𝑦⋅𝜉_𝑚𝛼𝐼^𝑛+1_{𝑖−1,𝑗,𝑚}+Δ𝑥⋅𝜂_𝑚𝛼𝐼^𝑛+1_{𝑖,𝑗−1,𝑚}+Δ𝑥⋅Δ𝑦(𝑆_{𝑖,𝑗,𝑚}+𝜇_𝑠^𝑀_𝑚′=1𝑤_𝑚′𝑝_𝑚,𝑚′𝐼^𝑛+1_{𝑖,𝑗,𝑚′})]. (18) In (18), we have used these notations:

(i)𝐼^𝑛_{𝑖,𝑗,𝑚}=𝐼(𝑖Δ𝑥,𝑗Δ𝑦,𝜉_𝑚,𝜂_𝑚,𝑛Δ𝑡)(ii)𝑆^𝑛_{𝑖,𝑗,𝑚}=𝑆(𝑖Δ𝑥,𝑗Δ𝑦,𝜉_𝑚,𝜂_𝑚,𝑛Δ𝑡).

To solve (18), we use successive iterations to actualise the implicit internal source term in the right member. So, we obtain for𝑛≥1,

(𝐼^𝑛+1_{𝑖,𝑗,𝑚})^𝑘+1=[3𝑛+2𝑐▵𝑡Δ𝑦⋅𝜉_𝑚𝛼+Δ𝑥⋅𝜂_𝑚𝛼+Δ𝑥⋅Δ𝑦(𝜇_𝑎+𝜇_𝑠)]⁻¹×[4𝑛2𝑐▵𝑡𝐼^𝑛_{𝑖,𝑗,𝑚}−𝑛2𝑐▵𝑡𝐼^𝑛−1_{𝑖,𝑗,𝑚}+Δ𝑦⋅𝜉_𝑚𝛼(𝐼^𝑛+1_{𝑖−1,𝑗,𝑚})^𝑘+1+Δ𝑥⋅𝜂_𝑚𝛼(𝐼^𝑛+1_{𝑖,𝑗−1,𝑚})^𝑘+1+Δ𝑥⋅Δ𝑦(𝑆^𝑛+1_{𝑖,𝑗,𝑚}+𝜇_𝑠^𝑀_𝑚′=1𝑤_𝑚′𝑝_𝑚,𝑚′(𝐼^𝑛+1_{𝑖,𝑗,𝑚′})^𝑘)]. (19) The iteration process is repeated until a convergence criterion is attempted and to improve convergence speed, we use a successive over relaxation technique [32]. So the updated value (𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘+1_updated becomes a linear combination of the iterated value (𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘+1 and the previously computed value (𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘 as

(𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘+1_updated=(1−𝜌)(𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘+𝜌(𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘+1, (20) where 𝜌 being a relaxation parameter whose value is usually between 1 and 2. In most of our calculus, we have used: 𝜌=1.1. The solution is obtained when the relative discrepancy value, 𝜖=|(𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘+1−(𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘/(𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘|, is smaller than a tolerance value tol. In that case, we take (𝐼^𝑛_{𝑖,𝑗,𝑚})^𝑘 as a solution. It will be noted 𝐼^𝑛_{𝑖,𝑗,𝑚}.

It is worth noting that, in all our simulations, we take a tolerance value: tol=10⁻⁸ and for initial condition, we take 𝐼⁰_{𝑖,𝑗,𝑚}=0. We use Δ𝑥=Δ𝑦 and a uniform time step Δ𝑡=10⁻¹³s. Also it is important to mention that when a continuous source is injected, another iterating process is implemented to attempt the steady-state solution. This is obtained by iterating the solution 𝐼^𝑛_{𝑖,𝑗,𝑚} until|𝐼^𝑛_{𝑖,𝑗,𝑚}−𝐼^𝑛−1_{𝑖,𝑗,𝑚}/𝐼^𝑛−1_{𝑖,𝑗,𝑚}|≤tol. At the equilibrium, all time derivatives will be vanished and we will get a solution of the time-independent radiative transfer equation. This procedure could be a validation criterion for the implemented algorithm convergence. Now, if the direction cosines are not positive, the precedent equations are valid provided that the orientation WESN of cells is done according to the direction of propagation. To sweep the calculus domain, we use negative increments. So, the set of all angular directions (Ω_𝑚)is divided into four subsets depending on the sign of 𝜉_𝑚 and 𝜂_𝑚. By the use of direction cosine sign, all cases can be given through the following general equation:

(𝐼^𝑛+1_{𝑖,𝑗,𝑚})^𝑘+1=[3𝑛+2𝑐▵𝑡Δ𝑦⋅sign(𝜉_𝑚)𝜉_𝑚𝛼+Δ𝑥⋅sign(𝜂_𝑚)𝜂_𝑚𝛼+Δ𝑥⋅Δ𝑦(𝜇_𝑎+𝜇_𝑠)]⁻¹×[4𝑛2𝑐▵𝑡𝐼^𝑛_{𝑖,𝑗,𝑚}−𝑛2𝑐▵𝑡𝐼^𝑛−1_{𝑖,𝑗,𝑚}+Δ𝑦⋅sign(𝜉_𝑚)𝜉_𝑚𝛼⋅(𝐼^𝑛+1_{𝑖−1×sign(𝜉𝑚),𝑗,𝑚})^𝑘+1+Δ𝑥⋅sign(𝜂_𝑚)𝜂_𝑚𝛼(𝐼^𝑛+1_{𝑖,𝑗−1×sign(𝜂𝑚),𝑚})^𝑘+1+Δ𝑥⋅Δ𝑦(𝑆^𝑛+1_{𝑖,𝑗,𝑚}+𝜇_𝑠^𝑀_𝑚′=1𝑤_𝑚′𝑝_𝑚,𝑚′(𝐼^𝑛+1_{𝑖,𝑗,𝑚′})^𝑘)]. (21)

3.2. External Source and Computed Quantities

In all our investigations, the power laser-source is assumed to be equivalent to a forward collimated radiance injected at a pinpoint on the bottom side of the boundary. We use two types of source:

(i)a continuous wave source with a uniform equivalent intensity value of 20 mWcm^-2sr^-1.(ii)a short-Laser-pulse source, applied only at the first instant of calculus. It is a 100fs-pulse having an equivalent radiance value of 200 mWcm^-2sr^-1,

In our presented results on the boundary, we use either calculated detected fluence rate or normalized detected fluence rate:

Φ^𝑁_𝑑=𝜙_𝑑1/𝐷∑^𝐷_𝑑𝑤_𝑑𝜙_𝑑, (22) where𝜙_𝑑is the computed detected fluence rate at a detector point 𝑑,𝐷 is the number of the detector points on the concerned side of the boundary, and 𝑤_𝑑 is a weighting factor from the generalized trapezoidal integration rule [34]. Also, distribution of light in the medium is presented through the calculated fluence rate matrix which is given at the instant, 𝑛Δ𝑡, as:

𝜙^𝑛_𝑖𝑗=^𝑀_𝑚=1𝑤_𝑚𝐼^𝑛_{𝑖,𝑗,𝑚}. (23)

4. Results and Discussion

4.1. Model Validation

To test our model, we have considered a case of 678 nm-light punctual continuous wave source. Our Model is tested by comparing obtained results with phantom experimental data reported in [18]. Figure 1 shows the phantom, it is assumed to be 4 cm×4 cm sized heterogeneous medium which contains a water-filled, void-like ring with an inner diameter of 2.8 cm and an outer diameter of 3.0 cm. A continuous source is placed on the middle of the bottom side of the boundary. The numerically predicted results are compared with corresponding experimentally measured data extracted from in [18, Figure  7]. The convergence of our algorithm is conditioned by some considerations in grid-mesh characteristics. In particular, the choice of the cell dimension should be taken too small. More precisely, to move from a cell to another neighboring cell, we should travel less than the total mean free path of photons. So, it should be more convenient to do calculations under the condition: Δ𝑥<1/(𝜇_𝑎+𝜇^′_𝑠). In these calculations, we have used Δ𝑥=Δ𝑦=0.05cm.

Figure 1 

Geometry and optical properties of the referencing medium [18].

Figure 2(a) represents the calculated steady state normalized transmitted fluence rate on the opposite side to the source for different values of the interpolation parameter. The case (𝛼=0.5) shows the closest results to the experimental measurements. In the rest of investigations, the interpolation parameter value will be putted 0.5. Also Figure 2(b) shows an almost angular grid-independency of our numerical results. We will take a set of 16 uniformly distributed discrete angles which is sufficient to do the rest of our calculus. These results and other results shown in [14] constitute both numerical tests and experimental validations of our model.

(a)

(b)

Figure 2 

Comparison of calculated normalized steady state transmitted fluence rate and experimental data of a heterogeneous tissue-like medium. Calculus is done in the case of a continuous wave source for different values of the interpolation parameter (a) and for different angular grids (b).

4.2. Computational Spectroscopy with a Continuous Wave Source

In this investigation, we consider a rat liver tissue-like medium as it is shown in Figure 3. Optical properties of the tissue are taken from [33] from a set of several wavelengths. Figures 4(a) and 4(b) illustrate transmission through the considered medium on the top side and on the left side for different laser-source wavelengths. Calculus is carried out as in the previous paragraph for a continuous wave source. Transmitted directional radiances are calculated in 18 detector points on the top side and on the left side of the boundary. Figures 4(a) and 4(b) show that transmission for 488 nm and 2100 nm is very weak. The maximal transmission from the tissue is observed in the red and the near infrared part of radiation spectrum. Also, Figure 5 shows the time behavior of the detected fluence rate in a point on the top side of the studied medium. During the transient period, Figure 4 shows no response of the medium for 488 nm and 2100 nm and a very weak transmission in 1320 nm and a relative larger global transmission in red and near infrared zones.

(a)

(b)

Figure 3 

Geometry and optical properties of the rat-liver-tissue at different wavelengths [33].

(a)

(b)

Figure 4 

Detected fluence rate for different source-wavelengths. Calculus is done in the case of a continuous wave source. Results are shown for top side (a) and left side (b).

Figure 5 

Time-dependent detected fluence rate at a detector point on top side of the rat-liver tissue-like medium. Calculus is done in the case of a continuous wave source.

4.3. Spectroscopy with a Short-Laser Pulse

In this investigation, we study the same medium as in the just precedent paragraph but we use a short-pulsed source. In all cases, it is a 100fs-laser pulse injected in the medium through a point on the middle of the bottom side of the boundary. Figure 6 shows the transmitted light at a detector point on the top side of the boundary. Results are presented for a set of different wavelengths. Only near infrared light displays a relative appreciable transmission. A maximum of detected fluence rate is noted at 800 nm-light. A little delay of the maximum detection of radiation is observed when wavelength increases. Figure 7 shows fluence rate distribution into the medium at different moments after the pulse. In all moments, symmetric propagation of light is observed into the medium. Also, it can be observed that scattering prevails in 488 nm while more and more absorption is observed in 2100 nm. At 488 nm, the phenomenon of multiple scattering is dominant into the medium. The dispersion of the pulse is accentuated in all sides. Light disappears after 250ps. At 2100 nm, more light is absorbed into the medium so a fraction of energy persists into the medium for more time.

Figure 6 

Time-dependent detected fluence rate on the top side of the rat-liver tissue-like medium. Calculus is done in the case of a 100fs-laser pulse injected in the bottom side with different wavelengths.

Figure 7 

Internal light distribution in the rat-liver tissue-like medium for different instants after the pulse. Calculus is done in the case of a 100fs-laser pulse for different wavelengths.

Figure 8 shows typical movement of different photons in the medium at three different wavelengths. Calculations in our studied medium are carried out by using Monte Carlo simulations as it is described in [35]. It can be observed that a typical 800 nm-photon is almost snake. For 488 nm-photon, the mean scattering free path is very small. Multiple scattering is dominant. The typical trajectory in this case is not ballistic and the photons are almost diffusive. For 2100 nm, the mean scattering free path is longer than the 488 nm-photon but the absorption mean free path is narrower. A typical photon in this wavelength is almost absorbed into the medium.

Figure 8 

Typical trajectory of photons for different values of the wavelength source as obtained through Monte Carlo simulations.

5. Conclusion

This paper investigated situations concerning biological tissue interaction with Laser radiation by using radiative transfer theory. More precisely, a computational radiative transfer model is used to establish a numerical Laser spectroscopic study in a rat-liver-tissue-like medium. The response of the medium to both time-dependent and time-independent sources is explored. Obtained Results are tested by using experimental measurements from the literature and by using Monte Carlo simulations. The maximal transmission from the tissue is observed in the red and the near infrared part of radiation spectrum. Also multiple scattering is prevailing in 488 nm while more and more absorption is observed in 2100 nm.

Nomenclature

𝑐:	Speed of light, (in vacuum), cm/s
𝑔:	Anisotropy factor
I:	Intensity, W⋅cm−2sr−1
𝑛:	Refractive index
𝑝(Ω,Ω′):	Phase function
𝐫:	Position vector, cm
𝑅:	Reflectivity coefficient
𝐫𝐛:	Position vector in the boundary
𝑆:	Source term, W⋅cm−3sr−1
𝑤𝑑:	Weighting factor from trapezoidal integration
𝑤𝑚:	Weight associated with discrete direction
𝑥,𝑦:	Cartesian coordinates, cm.

Greek Symbols

𝛼∶	Interpolation parameter
Φ:	Fluence rate, W⋅cm−2
Φ𝑑:	Detected fluence rate, W⋅cm−2
𝜆:	Wavelength, nm
Ω:	Direction solid angle, sr
𝜇𝑎:	Absorption coefficient, cm-1
𝜇𝑠:	Scattering coefficient, cm-1
𝜇′𝑠:	Reduced scattering coefficient, cm-1
𝜌:	Relaxation parameter
𝜃:	Scattering angle, rad
𝜉,𝜂:	Direction cosines.