# Computational Laser Spectroscopy in a Biological Tissue

Authors: M. Gantri; H. Trabelsi; E. Sediki; R. Ben Salah
Journal: Journal of Biophysics
 (2010)

---

## 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.


---

## Body

## 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 SPN 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]:(1)nc∂Iλ(r,Ω,t)∂t+Ω·∇Iλ(r,Ω,t)+(μa(r)+μs(r))Iλ(r,Ω,t)=μs∫02πp(Ω,Ω′)Iλ(r,Ω′,t)dΩ′+Sλ(r,Ω,t).
In this equation,(i)
Iλ(r,Ω,t)is the monochromatic directional energetic radiance at the spatial position-vector r=(x,y,z). When used in a dimensional form, the radiance Iλ(r,Ω,t)is often expressed in units of W·cm-2sr-1. The photon direction is described using a unitary direction vectorΩ;(ii)
μa(r)andμs(r) are the absorption and the scattering coefficients, respectively. They are both given in units of cm-1. c is the speed of light in vacuum. n is the medium’s refractive index. Also we make use of reduced scattering coefficient:
(2)μs′=(1-g)μs;(iii)
The source termSλ(r,Ω,t)is an injected monochromatic radiance at the medium boundary;(iv)
The phase functionp(Ω,Ω′)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]:(3)p(cos(θ))=1-g24π(1+g2-2gcos(θ))2/3,
where g is the anisotropy factor and θ is the angle between Ω andΩ'. The integral of radiance over all angles Ω at one point r yields the fluence rate:(4)Φλ(r,t)=∫Iλ(r,Ω,t)dΩ.
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(5)I(rb,Ω,t)=S(rb,Ω,t)+R·I(rb,Ωref,t),n·Ω<0,Ωref·n=-Ω·n,
where rb is a vector position on the boundary and n is the unit outer normal vector. The reflectivity R is given for each direction using a Fresnel’s relationship:
(6)R=12(nmcosφtra-n0cosφrefnmcosφtra+n0cosφref)2+12(n0cosφref-nmcosφtran0cosφref+nmcosφtra)2ifφref<φcrandelse,R=1,
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:(7)nmsinφref=n0sinφtra.
The medium refractive index nm depends on the tissue and n0 is taken as unity. The critical angle is determined by the following equation:(8)sinφcr=n0nm.
Thus, the transmitted fluence rate at the boundary can be computed in each detector point using the following integral formula:(9)Φtra(rb,t)=∫Ω·n>0(1-R)I(rb,Ω,t)dΩ.
The fluence rate on the boundary which enters the detector for a given aperture AP is computed as.(10)Φd(rb,t)=∫Ap(1-R)I(rb,Ω,t)dΩ
In this work, Φd will be called detected fluence rate and we will take AP = 45° in all our investigations.


## 2.1. Physical Model

Light propagation in a human biological tissue is modeled by using a general time-dependent RTE [26]:(1)nc∂Iλ(r,Ω,t)∂t+Ω·∇Iλ(r,Ω,t)+(μa(r)+μs(r))Iλ(r,Ω,t)=μs∫02πp(Ω,Ω′)Iλ(r,Ω′,t)dΩ′+Sλ(r,Ω,t).
In this equation,(i)
Iλ(r,Ω,t)is the monochromatic directional energetic radiance at the spatial position-vector r=(x,y,z). When used in a dimensional form, the radiance Iλ(r,Ω,t)is often expressed in units of W·cm-2sr-1. The photon direction is described using a unitary direction vectorΩ;(ii)
μa(r)andμs(r) are the absorption and the scattering coefficients, respectively. They are both given in units of cm-1. c is the speed of light in vacuum. n is the medium’s refractive index. Also we make use of reduced scattering coefficient:
(2)μs′=(1-g)μs;(iii)
The source termSλ(r,Ω,t)is an injected monochromatic radiance at the medium boundary;(iv)
The phase functionp(Ω,Ω′)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]:(3)p(cos(θ))=1-g24π(1+g2-2gcos(θ))2/3,
where g is the anisotropy factor and θ is the angle between Ω andΩ'. The integral of radiance over all angles Ω at one point r yields the fluence rate:(4)Φλ(r,t)=∫Iλ(r,Ω,t)dΩ.
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(5)I(rb,Ω,t)=S(rb,Ω,t)+R·I(rb,Ωref,t),n·Ω<0,Ωref·n=-Ω·n,
where rb is a vector position on the boundary and n is the unit outer normal vector. The reflectivity R is given for each direction using a Fresnel’s relationship:
(6)R=12(nmcosφtra-n0cosφrefnmcosφtra+n0cosφref)2+12(n0cosφref-nmcosφtran0cosφref+nmcosφtra)2ifφref<φcrandelse,R=1,
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:(7)nmsinφref=n0sinφtra.
The medium refractive index nm depends on the tissue and n0 is taken as unity. The critical angle is determined by the following equation:(8)sinφcr=n0nm.
Thus, the transmitted fluence rate at the boundary can be computed in each detector point using the following integral formula:(9)Φtra(rb,t)=∫Ω·n>0(1-R)I(rb,Ω,t)dΩ.
The fluence rate on the boundary which enters the detector for a given aperture AP is computed as.(10)Φd(rb,t)=∫Ap(1-R)I(rb,Ω,t)dΩ
In this work, Φd 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 ∆V. The spatial calculus domain is divided into a set of I×J elementary uniform volumes ∆V with unitary depth. This yields a calculus domain subdivided in I*J rectangular cells. Each (i,j)-cell is Δx∆y-sized. So, we can write.(11)nc∂IP,m∂t+Δy·ξm(IE,m-IW,m)+Δx·ηm(IN,m-IS,m)+Δx·Δy(μa+μs)IP,m=Δx·Δy(SP,m+μs∑m′=1Mwm′pm,m′IP,m′),
where ξm=cos(Ωm·nx)andηm=cos(Ωm·ny)with nxandnybeing two normal vectors to the (x,z)-plane and the (y,z)-plane, respectively. wm 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(12)pmm′=1-g²4π(1+g²-2g(ξmξm′+ηmηm′))3/2.
If the direction cosines are positive, the directional energetic radiances are known on the faces W and S and they are unknown on the faces E and N of the (i,j)-cell and also in the centre P. Therefore, we need two complementary relations to eliminate IE,m and IN,m.These relations can be obtained by using interpolation formula:(13)IP,m=αIE,m+(1-α)IW,m,IP,m=αIN,m+(1-α)IS,m,
where α is an interpolation parameter. Using the above formulae, (11) can be written as(14)nc∂IP,m∂t+Δy·ξmα(IP,m-IW,m)+Δx·ηmα(IP,m-IS,m)+Δx·Δy(μa+μs)IP,m=Δx·Δy(SP,m+μs∑m′=1Mwm′pm,m′IP,m′).
For time discretization, an implicit three-level second-order time differencing scheme [31] is used:(15)∂IP,m∂t=3IP,mn+1-4IP,mn+IP,mn-12∆t,n=1,…nmax,
whith Δt being the discrete time step. So we can rewrite (14) as follows:(16)IP,mn+1(3n2c∆t+Δy·ξmα+Δx·ηmα+Δx·Δy(μa+μs))=4n2c∆tIP,mn-n2c∆tIP,mn-1+Δy·ξmαIW,mn+1+Δx·ηmαIS,mn+1+Δx·Δy(SP,m+μs∑m′=1Mwm′pm,m′IP,m′n+1).
Theoretically, if we know IP,m on the (i,j)-cell, we can obtain results over the cells (i+1,j) and (i,j+1) by using the boundary conditions and the following relations:(17)IW,m(i+1,j)=IE,m(i,j),i=1,…,I-1,IS,m(i,j+1)=IN,m(i,j),j=1,…,J-1.
If (ξm>0andηm>0), then we get the following equation:(18)Ii,j,mn+1=[3n2c∆t+Δy·ξmα+Δx·ηmα+Δx·Δy(μa+μs)]-1×[4n2c∆tIi,j,mn-n2c∆tIi,j,mn-1+Δy·ξmαIi-1,j,mn+1+Δx·ηmαIi,j-1,mn+1+Δx·Δy(Si,j,m+μs∑m′=1Mwm′pm,m′Ii,j,m′n+1)].
In (18), we have used these notations:(i)
Ii,j,mn=I(iΔx,jΔy,ξm,ηm,nΔt)(ii)
Si,j,mn=S(iΔx,jΔy,ξm,ηm,nΔt).To solve (18), we use successive iterations to actualise the implicit internal source term in the right member. So, we obtain forn≥1,(19)(Ii,j,mn+1)k+1=[3n2c∆t+Δy·ξmα+Δx·ηmα+Δx·Δy(μa+μs)]-1×[4n2c∆tIi,j,mn̅-n2c∆tIi,j,mn-1̅+Δy·ξmα(Ii-1,j,mn+1)k+1+Δx·ηmα(Ii,j-1,mn+1)k+1+Δx·Δy(Si,j,mn+1+μs∑m′=1Mwm′pm,m′(Ii,j,m′n+1)k)].
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 (Ii,j,mn)updatedk+1 becomes a linear combination of the iterated value (Ii,j,mn)k+1 and the previously computed value (Ii,j,mn)k as(20)(Ii,j,mn)updatedk+1=(1-ρ)(Ii,j,mn)k+ρ(Ii,j,mn)k+1,
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, ϵ=|(Ii,j,mn)k+1-(Ii,j,mn)k/(Ii,j,mn)k|, is smaller than a tolerance value tol. In that case, we take (Ii,j,mn)k as a solution. It will be noted I̅i,j,mn.It is worth noting that, in all our simulations, we take a tolerance value:tol=10-8 and for initial condition, we take I̅i,j,m0=0. We use Δx=Δy and a uniform time step Δt=10-13s. 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 I̅i,j,mn until|I̅i,j,mn-I̅i,j,mn-1/I̅i,j,mn-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 (Ωm)is divided into four subsets depending on the sign of ξm and ηm. By the use of direction cosine sign, all cases can be given through the following general equation:(21)(Ii,j,mn+1)k+1=[3n2c∆t+Δy·sign(ξm)ξmα+Δx·sign(ηm)ηmα+Δx·Δy(μa+μs)3n2c∆t]-1×[∑m′=1Mwm′pm,m′4n2c∆tIi,j,mn̅-n2c∆tIi,j,mn-1̅+Δy·sign(ξm)ξmα·(Ii-1×sign(ξm),j,mn+1)k+1+Δx·sign(ηm)ηmα(Ii,j-1×sign(ηm),mn+1)k+1+Δx·Δy(Si,j,mn+1+μs∑m′=1Mwm′pm,m′(Ii,j,m′n+1)k)].


### 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:(22)ΦdN=ϕd1/D∑dDwdϕd,
whereϕdis the computed detected fluence rate at a detector point d,D is the number of the detector points on the concerned side of the boundary, and wd 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, nΔt, as:(23)ϕijn=∑m=1MwmI̅i,j,mn.


## 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 ∆V. The spatial calculus domain is divided into a set of I×J elementary uniform volumes ∆V with unitary depth. This yields a calculus domain subdivided in I*J rectangular cells. Each (i,j)-cell is Δx∆y-sized. So, we can write.(11)nc∂IP,m∂t+Δy·ξm(IE,m-IW,m)+Δx·ηm(IN,m-IS,m)+Δx·Δy(μa+μs)IP,m=Δx·Δy(SP,m+μs∑m′=1Mwm′pm,m′IP,m′),
where ξm=cos(Ωm·nx)andηm=cos(Ωm·ny)with nxandnybeing two normal vectors to the (x,z)-plane and the (y,z)-plane, respectively. wm 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(12)pmm′=1-g²4π(1+g²-2g(ξmξm′+ηmηm′))3/2.
If the direction cosines are positive, the directional energetic radiances are known on the faces W and S and they are unknown on the faces E and N of the (i,j)-cell and also in the centre P. Therefore, we need two complementary relations to eliminate IE,m and IN,m.These relations can be obtained by using interpolation formula:(13)IP,m=αIE,m+(1-α)IW,m,IP,m=αIN,m+(1-α)IS,m,
where α is an interpolation parameter. Using the above formulae, (11) can be written as(14)nc∂IP,m∂t+Δy·ξmα(IP,m-IW,m)+Δx·ηmα(IP,m-IS,m)+Δx·Δy(μa+μs)IP,m=Δx·Δy(SP,m+μs∑m′=1Mwm′pm,m′IP,m′).
For time discretization, an implicit three-level second-order time differencing scheme [31] is used:(15)∂IP,m∂t=3IP,mn+1-4IP,mn+IP,mn-12∆t,n=1,…nmax,
whith Δt being the discrete time step. So we can rewrite (14) as follows:(16)IP,mn+1(3n2c∆t+Δy·ξmα+Δx·ηmα+Δx·Δy(μa+μs))=4n2c∆tIP,mn-n2c∆tIP,mn-1+Δy·ξmαIW,mn+1+Δx·ηmαIS,mn+1+Δx·Δy(SP,m+μs∑m′=1Mwm′pm,m′IP,m′n+1).
Theoretically, if we know IP,m on the (i,j)-cell, we can obtain results over the cells (i+1,j) and (i,j+1) by using the boundary conditions and the following relations:(17)IW,m(i+1,j)=IE,m(i,j),i=1,…,I-1,IS,m(i,j+1)=IN,m(i,j),j=1,…,J-1.
If (ξm>0andηm>0), then we get the following equation:(18)Ii,j,mn+1=[3n2c∆t+Δy·ξmα+Δx·ηmα+Δx·Δy(μa+μs)]-1×[4n2c∆tIi,j,mn-n2c∆tIi,j,mn-1+Δy·ξmαIi-1,j,mn+1+Δx·ηmαIi,j-1,mn+1+Δx·Δy(Si,j,m+μs∑m′=1Mwm′pm,m′Ii,j,m′n+1)].
In (18), we have used these notations:(i)
Ii,j,mn=I(iΔx,jΔy,ξm,ηm,nΔt)(ii)
Si,j,mn=S(iΔx,jΔy,ξm,ηm,nΔt).To solve (18), we use successive iterations to actualise the implicit internal source term in the right member. So, we obtain forn≥1,(19)(Ii,j,mn+1)k+1=[3n2c∆t+Δy·ξmα+Δx·ηmα+Δx·Δy(μa+μs)]-1×[4n2c∆tIi,j,mn̅-n2c∆tIi,j,mn-1̅+Δy·ξmα(Ii-1,j,mn+1)k+1+Δx·ηmα(Ii,j-1,mn+1)k+1+Δx·Δy(Si,j,mn+1+μs∑m′=1Mwm′pm,m′(Ii,j,m′n+1)k)].
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 (Ii,j,mn)updatedk+1 becomes a linear combination of the iterated value (Ii,j,mn)k+1 and the previously computed value (Ii,j,mn)k as(20)(Ii,j,mn)updatedk+1=(1-ρ)(Ii,j,mn)k+ρ(Ii,j,mn)k+1,
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, ϵ=|(Ii,j,mn)k+1-(Ii,j,mn)k/(Ii,j,mn)k|, is smaller than a tolerance value tol. In that case, we take (Ii,j,mn)k as a solution. It will be noted I̅i,j,mn.It is worth noting that, in all our simulations, we take a tolerance value:tol=10-8 and for initial condition, we take I̅i,j,m0=0. We use Δx=Δy and a uniform time step Δt=10-13s. 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 I̅i,j,mn until|I̅i,j,mn-I̅i,j,mn-1/I̅i,j,mn-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 (Ωm)is divided into four subsets depending on the sign of ξm and ηm. By the use of direction cosine sign, all cases can be given through the following general equation:(21)(Ii,j,mn+1)k+1=[3n2c∆t+Δy·sign(ξm)ξmα+Δx·sign(ηm)ηmα+Δx·Δy(μa+μs)3n2c∆t]-1×[∑m′=1Mwm′pm,m′4n2c∆tIi,j,mn̅-n2c∆tIi,j,mn-1̅+Δy·sign(ξm)ξmα·(Ii-1×sign(ξm),j,mn+1)k+1+Δx·sign(ηm)ηmα(Ii,j-1×sign(ηm),mn+1)k+1+Δx·Δy(Si,j,mn+1+μs∑m′=1Mwm′pm,m′(Ii,j,m′n+1)k)].


## 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:(22)ΦdN=ϕd1/D∑dDwdϕd,
whereϕdis the computed detected fluence rate at a detector point d,D is the number of the detector points on the concerned side of the boundary, and wd 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, nΔt, as:(23)ϕijn=∑m=1MwmI̅i,j,mn.


## 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: Δx<1/(μa+μs′). In these calculations, we have used Δx=Δy=0.05cm.Figure 1
Geometry and optical properties of the referencing medium [18].Figure2(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.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).
(a)(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 Figure3. 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.Geometry and optical properties of the rat-liver-tissue at different wavelengths [33].
(a)(b)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).
(a)(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. Figure6 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.Figure8 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.


## 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: Δx<1/(μa+μs′). In these calculations, we have used Δx=Δy=0.05cm.Figure 1
Geometry and optical properties of the referencing medium [18].Figure2(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.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).
(a)(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 Figure3. 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.Geometry and optical properties of the rat-liver-tissue at different wavelengths [33].
(a)(b)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).
(a)(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. Figure6 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.Figure8 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.



---

*Source: 253763-2010-04-07.xml*