Boundary integral method for simulating laser short-pulse penetration into biological tissues

Mohammad Ali Ansari; Reza Massudi

doi:10.1117/1.3526351

1 November 2010 Boundary integral method for simulating laser short-pulse penetration into biological tissues

Mohammad Ali Ansari, Reza Massudi

Author Affiliations +

Journal of Biomedical Optics, Vol. 15, Issue 6, 065009 (November 2010). https://doi.org/10.1117/1.3526351

Abstract

The study of short-pulse propagation through biological tissues is important due to the medical applications of short-pulse lasers. Techniques used for numerical study of short pulses through human tissues include the Monte Carlo (MC) method, the finite-element method, and the finite-difference time-domain (FDTD), but these are often time consuming. Recently, the boundary integral method (BIM) was applied to overcome this problem. The literature shows that the BIM is faster than the other mentioned methods. We first investigate the precision of results obtained by the BIM by comparison with those results obtained by the MC and FDTD methods. Then we use the BIM to investigate the short-pulse penetration into biological tissues. We also study the effects of optical properties of tissues such as scattering, the absorption coefficient, the anisotropic factor on the penetrating pulse. We also, consider the propagation of pulses emitted from extended sources with different temporal evolutions.

1. Introduction

Recently, short-pulse lasers have been widely used in noninvasive optical tomography.^{1, 2, 3, 4, 5, 6} As compared with continuous wave imaging, more information can be achieved in short-pulse imaging by studying the temporal distribution of reflectance/transmittance.^{7, 8, 9, 10, 11} The temporal distribution of a reflected/transmitted signal is broadened due to multiple scattering of photons into tissues.^{2, 12} The level of broadening depends on optical properties of the tissue, i.e., absorption and scattering; consequently, variations of those parameters can be investigated by studying the temporal distribution of the scattered light. Thus, the temporal distribution of the reflectance/transmittance can be exploited to detect normal or malignant biological tissue, because the optical properties of normal tissues vary during malignant progress.¹³

Photon transport through biological tissues can be described by a diffuse equation. Sometimes, this equation can be solved analytically.¹⁴ But in most cases, numerical methods such as the Monte Carlo (MC) method, the finite-element method (FEM), and the finite-difference time-domain (FDTD) method are most often used to solve that equation.^{15, 16, 17, 18, 19, 20, 21, 22} In those methods, the whole sample should be discritized and the calculation is time consuming.¹⁵

The boundary integral method (BIM) can also be used to solve a diffuse equation and to simulate a reflected pulse on the surface of biological tissues. Since this method requires surface tessellation, the computation time is reduced and the accuracy of the results increases as compared to other numerical methods.^{23, 24, 25, 26} The effects of optical properties of biological tissues on the intensity of a reflected pulse were investigated and reported in Ref. 26, but the behavior of a pulse penetrating through biological tissues has not yet been studied by this method.

In this paper, pulse penetration into biological tissue is studied using the BIM. First, the appropriate Green's function is obtained to convert the diffuse equation to integral form using Green's second theorem. Then, the surface integral is discretized using the boundary element method (BEM). In this method, the boundary of the sample is first discretized to elements. Then, the observation point is located on the surface of tissue and an equation containing fluence at that point is achieved. Locating the observation point on different nodes, a system of equations is obtained that gives the fluence at those points.^{23, 24, 25} By using this technique, the photon intensity inside the sample and the intensity of the diffusely penetrated pulse inside tissue are calculated. To investigate the accuracy and precision of results, they are compared with those obtained analytically and by other numerical methods. Furthermore, the effects of the absorption and scattering coefficients and the anisotropic factor on the diffusely penetrating pulse are also studied. Finally, the penetration of the short pulse with a arbitrary time shape is studied.

2. Review of Theory

Short-pulse laser propagation into biological tissue Ω with boundary [TeX:] $\Gamma _s$ $Γ_{s}$ can be studied by a diffuse equation and Robin's boundary condition, which are, respectively, given by¹⁶

Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{\partial }{{c\partial t}}\varphi \left({{\bm r}\!,t} \right) - D\nabla ^2 \varphi \left({{\bm r}\!,t} \right) + a\varphi \left({{\bm r}\!,t} \right) = S\left({{\bm r}\!,t} \right)\quad{\bm r} \in \Omega, \end{equation}\end{document}

\frac{\partial}{c \partial t} ϕ (r, t) - D \nabla^{2} ϕ (r, t) + a ϕ (r, t) = S (r, t) r \in Ω,

Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \hspace*{25pt}\varphi \left({{\bm r}\!,t} \right) - 2C_R D\frac{{\partial \varphi \left({{\bm r}\!,t} \right)}}{{\partial z}} = 0\quad {\bm r} \in \Gamma _s, \end{equation}\end{document}

ϕ (r, t) - 2 C_{R} D \frac{\partial ϕ (r, t)}{\partial z} = 0 r \in Γ_{s},

where [TeX:] $\varphi \left({{\bm r}\!,t} \right)$ $ϕ (r, t)$ and [TeX:] $S\left({{\bm r}\!,t} \right)$ $S (r, t)$ are, respectively, the fluence and the isotropic source term at position [TeX:] ${\bm r}$ $r$ and at moment t. The velocity of light is shown by c. The parameter [TeX:] $D = {1 / {3({a + \sigma^ \prime })}}$ $D = 1 / 3 (a + σ^{'})$ is the diffusion coefficient, where a and [TeX:] $\sigma^ \prime = \sigma ({1 - g})$ $σ^{'} = σ (1 - g)$ are the absorption and reduced scattering coefficients, respectively; and σ and g are also the scattering coefficient and the anisotropic factor, respectively. In Eq. 2, [TeX:] $C_R = {{({1 + R})} / {({1 - R})}}$ $C_{R} = (1 + R) / (1 - R)$ , where R is the Fresnel reflection coefficient.

The boundary integral method is based on using of Green's function.²³ The Green's function of Eq. 1 in domain Ω is the solution of

Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\frac{\partial }{{c\partial t}}g({{\bm r}\!,{\bm r'};t,t'}) - D\nabla ^2 g({{\bm r}\!,{\bm r'};t,t'}) + ag({{\bm r}\!,{\bm r'};t,t'})\nonumber \\ && \quad= - \delta ({{\bm r} - {\bm r}'})\delta ({t - t'}). \end{eqnarray}\end{document}

\begin{matrix} \frac{\partial}{c \partial t} g (r, r^{'}; t, t^{'}) - D \nabla^{2} g (r, r^{'}; t, t^{'}) + a g (r, r^{'}; t, t^{'}) \\ = - δ (r - r^{'}) δ (t - t^{'}) . \end{matrix}

Applying a Laplace transform on t in Eq. 3 and rearranging the resulting equation, gives

Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\hspace*{-6pt} \nabla ^2 G({{\bm r}\!,{\bm r'};s,t'})-\kappa ^2 G({{\bm r}\!,{\bm r'};s,t'}) = \frac{{\delta ({{\bm r} - {\bm r'}})}}{D}\,\exp\, ({ - st'}),\phantom{00} \end{eqnarray}\end{document}

\begin{matrix} \nabla^{2} G (r, r^{'}; s, t^{'}) - κ^{2} G (r, r^{'}; s, t^{'}) = \frac{δ (r - r^{'})}{D} \exp (- s t^{'}), \end{matrix}

where [TeX:] $G({{\bm r}\!,{\bm r^\prime}; s,t^\prime }) = L[ {g\,({{\bm r}\!,{\bm r^\prime}\!;s,t^\prime })} ]$

G (r, r^{'}; s, t^{'}) = L [g (r, r^{'}; s, t^{'})]

, [TeX:] $\kappa ^2 = (s/cD) + a/D,$

κ^{2} = (s / c D) + a / D,

and [TeX:] $\nabla ^2$

\nabla^{2}

operates on [TeX:] ${\bm r}$

r

. Next, applying Fourier transform on [TeX:] ${\bm r} $

r

in Eq. 4 and doing some mathematics results in

Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \tilde G({k,{\bm r}';s,t'}) = \frac{{\exp\, ({ - ik.{\bm r}'})}}{{D({k^2 + \kappa ^2 })}}\exp\, ({ - st'}), \end{equation}\end{document}

\tilde{G} (k, r^{'}; s, t^{'}) = \frac{\exp (- i k . r^{'})}{D (k^{2} + κ^{2})} \exp (- s t^{'}),

where [TeX:] $\tilde G$

\tilde{G}

is Fourier transform of G.

The inverse Fourier and Laplace transforms of Green's function stated in Eq. 5 gives the Green's function as

Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \hspace*{-16pt}g({{\bm r},{\bm r}';t,t'}) &=& H({t - t'})\frac{c}{{{\{ {[ {4\pi {\rm }Dc{\rm }({t - t'})} ]} ^3\}}^{1/2}}} \nonumber\\ && \times\, \exp[ { - ca({t - t'})} ]\exp \left[ { - \frac{{\left| {{\bm r} - {\bm r}'} \right|^2 }}{{4Dc{\rm }\left({t - t'} \right)}}} \right], \end{eqnarray}\end{document}

\begin{matrix} g (r, r^{'}; t, t^{'}) & = & H (t - t^{'}) \frac{c}{{{[4 π D c (t - t^{'})]}^{3}}^{1 / 2}} \\ \times \exp [- c a (t - t^{'})] \exp [- \frac{{|r - r^{'}|}^{2}}{4 D c (t - t^{'})}], \end{matrix}

where [TeX:] $H({t - t^\prime })$

H (t - t^{'})

is the Heaviside function.

For a short-point pulse where [TeX:] $S\left({{\bm r}\!,t} \right) = \delta \left({\bm r} \right)\delta \left(t \right)$ $S (r, t) = δ (r) δ (t)$ , the fluence at any arbitrary point inside an infinite medium for t > 0 is equal to²⁶

Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm }\varphi ({{\bm r}\!,t}) = \frac{c}{{ {[ {({4\pi {\rm }Dct})} ^3]^{1/2} } }}\exp ({ - cat})\exp \left({ - \frac{{| {\bm r} |^2 }}{{4Dct}}} \right). \end{equation}\end{document}

ϕ (r, t) = \frac{c}{{[{(4 π D c t)}^{3}]}^{1 / 2}} \exp (- c a t) \exp (- \frac{{| r |}^{2}}{4 D c t}) .

Assuming the sample as a slab and by using Eq. 7, the transmited intensity from it is given by²⁷

Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} T\left({d,t} \right) &=& \left({4\pi Dct} \right)^{{{ - 1}/2}} t^{{{ - 3}/2}} \,{\rm exp}\,\left({ - act} \right) \nonumber \\ && \times\, \left\{ {\left({d - z_0 } \right)\,\exp \,\left[ { - \displaystyle\frac{{\left({d - z_0 } \right)^2 }}{{4Dct}}} \right]} \right. \nonumber \\ && -\, \left({d + z_0 } \right)\,\exp \,\left[ { - \displaystyle\frac{{\left({d + z_0 } \right)^2 }}{{4Dct}}} \right] \nonumber\\ && +\, \left({3d - z_0 } \right)\,\exp \,\left[ { - \displaystyle\frac{{\left({3d - z_0 } \right)^2 }}{{4Dct}}} \right] \nonumber \\ && -\, \left.\left({3d + z_0 } \right)\,\exp \,\left[ { - \displaystyle\frac{{\left({3d + z_0 } \right)^2 }}{{4Dct}}} \right] \right\}, \end{eqnarray}\end{document}

\begin{matrix} T (d, t) & = & {(4 π D c t)}^{- 1 / 2} t^{- 3 / 2} \exp (- a c t) \\ \times \{(d - z_{0}) \exp [- \frac{{(d - z_{0})}^{2}}{4 D c t}] \\ - (d + z_{0}) \exp [- \frac{{(d + z_{0})}^{2}}{4 D c t}] \\ + (3 d - z_{0}) \exp [- \frac{{(3 d - z_{0})}^{2}}{4 D c t}] \\ - (3 d + z_{0}) \exp [- \frac{{(3 d + z_{0})}^{2}}{4 D c t}]\}, \end{matrix}

where d is the thickness of the slab, and [TeX:] $z_0 = 1/\sigma ^\prime $

z_{0} = 1 / σ^{'}

. For an extended source, the fluence is obtained by using Green's second theorem:

Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \varphi ({{\bm r}\!,t}) &=& \displaystyle\int\nolimits^t_\Omega\!\! {\int\nolimits_0 {S({{\bm r}',t'})} } g ({{\bm r}\!,t;{\bm r}',t'})\,{\rm d}{\bm r}'dt' \nonumber\\ && - \displaystyle\int\nolimits^t_{\Gamma _s }\! \int\nolimits_0 \Bigg\{\left[ \frac{c}{C_R D}\varphi ({\bm r}_s, t')g({\bm r}\!,t;{\bm r}_s, t') \right]\nonumber \\ && - \left[ \varphi ({\bm r}_s, t')\frac{\partial }{\partial n}g ({\bm r}\!,t;{\bm r}_s, t') \right]\Bigg\} {\rm d}{\bm r}_s\, dt', \end{eqnarray}\end{document}

\begin{matrix} ϕ (r, t) & = & \int_{Ω}^{t} \int_{0} S (r^{'}, t^{'}) g (r, t; r^{'}, t^{'}) d r^{'} d t^{'} \\ - \int_{Γ_{s}}^{t} \int_{0} {[\frac{c}{C_{R} D} ϕ (r_{s}, t^{'}) g (r, t; r_{s}, t^{'})] \\ - [ϕ (r_{s}, t^{'}) \frac{\partial}{\partial n} g (r, t; r_{s}, t^{'})]} d r_{s} d t^{'}, \end{matrix}

where [TeX:] ${\bm r}_s$

r_{s}

is the observation point vector on the boundary [TeX:] $\Gamma _s$

Γ_{s}

.

The intensity of a diffused short pulse can be calculated from Eq. 9, which is the integral form of diffuse equation. The BEM can be used to solve Eq. 9. In this method, the boundary of the sample is first discretized to elements. Then, observation point [TeX:] ${\bm r}_s$ $r_{s}$ is located on the surface of tissue and an equation containing fluence at that point is achieved. Locating observation point on different nodes, a system of equations is obtained which gives the fluence at those points. Thus, the boundary [TeX:] $\Gamma _s^{}$ $Γ_{s}^{}$ is discretized to square elements and the fluence φ is approximated as

Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \varphi \left({{\bm r}\!,t} \right) = \sum\limits_{k = 1}^n {\sum\limits_{l = 1}^m {N_k \left({\bm r} \right)\varphi _{kl} \left({{\bm r}\!,t} \right)} }, \end{equation}\end{document}

ϕ (r, t) = \sum_{k = 1}^{n} \sum_{l = 1}^{m} N_{k} (r) ϕ_{k l} (r, t),

where k and l refer to node k and time step l, and [TeX:] $N_i ({{\bm r}_j }) = \delta _{ij}$

N_{i} (r_{j}) = δ_{i j}

is the Kronecker symbol. If [TeX:] ${\bm r}_j$

r_{j}

spans all the nodes on the surface of the boundary, Eqs. 1, 2 will give, respectively, the following set of algebraic equations:

Eq. 11

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{l} H{\bm u} + \Gamma {\bm v} = \bar {\bm S}, \\[2pt] {\bm v} = - R{\bm u} + P, \\ \end{array} \end{equation}\end{document}

\begin{matrix} H u + Γ v = \bar{S}, \\ v = - R u + P, \end{matrix}

where [TeX:] $R = C_R^{ - 1}$

R = C_{R}^{- 1}

, and [TeX:] ${\bm u}$

u

, [TeX:] ${\bm v}$

v

, [TeX:] ${\bm P}$

P

, and [TeX:] $\bar {\bm S}$

\bar{S}

are column vectors of the nodal values of the fluence φ, normal derivative q, the prescribed boundary flux p, and the volume source S, respectively, which are given by

Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} {\bm u}_{n \times 1} &=& \left[ \begin{array}{l} \phi _1 \\ {\cdots} \\ {\cdots} \\ \phi _n \\ \end{array} \right], \quad {\bm v}_{n \times 1} = \left[ \begin{array}{l} q_1 \\ {\cdots} \\ {\cdots} \\ q_n \\ \end{array} \right] , \nonumber\\ {\bm P}_{n \times 1} & =&\left[ \begin{array}{l} p_1 \\ {\cdots} \\ {\cdots} \\ p_n \\ \end{array} \right] , \quad\bar {\bm S}_{n \times 1} = \left[ \begin{array}{l} s_1 \\ {\cdots} \\ {\cdots} \\ s_n \\ \end{array} \right]{\rm, } \end{eqnarray}\end{document}

\begin{matrix} u_{n \times 1} & = & [\begin{matrix} φ_{1} \\ \dots \\ \dots \\ φ_{n} \end{matrix}], v_{n \times 1} = [\begin{matrix} q_{1} \\ \dots \\ \dots \\ q_{n} \end{matrix}], \\ P_{n \times 1} & = & [\begin{matrix} p_{1} \\ \dots \\ \dots \\ p_{n} \end{matrix}], {\bar{S}}_{n \times 1} = [\begin{matrix} s_{1} \\ \dots \\ \dots \\ s_{n} \end{matrix}], \end{matrix}

where

Eq. 13

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} s_j = \int\nolimits_{\Gamma _s }\!\! {\int\limits_{} {g({\rho; t,t'})S({\bm r},{t})\ {\rm d}rdt',} } {\rm } \end{equation}\end{document}

s_{j} = \int_{Γ_{s}} \int_{} g (ρ; t, t^{'}) S (r, t) d r d t^{'},

and [TeX:] $\rho = | {{\bm r} - {r}_j } |$

ρ = | r - r_{j} |

. The spatial and temporal shape of [TeX:] $S\left({{\bm r}\!,t} \right)$

S (r, t)

can be arbitrarily chosen. The elements of matrices [TeX:] ${{\mbox{\sffamily\bfseries H}}}_{n \times n \times m}\break = \{ {h_{i,j,l} } \}$

H_{n \times n \times m} = {h_{i, j, l}}

and [TeX:] ${\bm \Gamma} _{n \times n \times m} = \{ {\xi _{i,j,l} }\}$

Γ_{n \times n \times m} = {ξ_{i, j, l}}

are as

Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} h_{i,j,l} &=& \delta _{ij} I + \displaystyle\int\nolimits^t_{\Gamma _s }\! {\int\nolimits_0 {\frac{{\partial g\left({\rho; \tau } \right)}}{{\partial n}}N_{kl} \left({{\bm r}\!,t} \right)\,{\rm d}{\bm r}\,{\rm d}t'} }, \nonumber \\[-9pt] \\[-9pt]\nonumber \xi _{i,j,l} &=& -\displaystyle \int\nolimits^t_{\Gamma _s }\! {\int\nolimits_0 {g\left({\rho; \tau } \right)N_{kl} \left({{\bm r}\!,t} \right){\rm d}{\bm r}\,{\rm d}t'}, } \end{eqnarray}\end{document}

\begin{matrix} h_{i, j, l} & = & δ_{i j} I + \int_{Γ_{s}}^{t} \int_{0} \frac{\partial g (ρ; τ)}{\partial n} N_{k l} (r, t) d r d t^{'}, \\ ξ_{i, j, l} & = & - \int_{Γ_{s}}^{t} \int_{0} g (ρ; τ) N_{k l} (r, t) d r d t^{'}, \end{matrix}

where [TeX:] $\tau = t - t^\prime $

τ = t - t^{'}

. Locating observation point on different nodes and using Eq. 11, a system of equations are obtained that give the fluence at those surface points. Next, one can calculate the fluence at any arbitrary point inside and outside the sample. According to Eq. 9, the fluence [TeX:] $\varphi \left({{\bm r}\!,t} \right)$

ϕ (r, t)

inside the sample can be calculated by using

Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} %\begin{array}{l} \varphi ^* ({{\bm r}\!,t}) &=&\displaystyle \int\nolimits^t_\Omega \!\!{\int\nolimits_0 {S ({{\bm r}',t'})} } g ({{\bm r}\!,t;{\bm r}',t'})\,{\rm d}{\bm r}'\,{\rm d}t' \nonumber \\ && - \displaystyle\int\nolimits^t_{\Gamma _s } \!\int\nolimits_0 \Bigg\{\left[ \frac{c}{{C_R D}}\bar \varphi ({{\bm r}_s, t'})g ({{\bm r}\!,t;{\bm r}_s, t'}) \right] \nonumber \\ &&- \left[ \bar \varphi ({{\bm r}_s, t'})\frac{\partial }{{\partial n}}g ({{\bm r}\!,t;{\bm r}_s, t'}) \right]\Bigg\} \,{\rm d}{\bm r}_s \,{\rm d}t', % \end{array} \end{eqnarray}\end{document}

\begin{matrix} ϕ^{*} (r, t) & = & \int_{Ω}^{t} \int_{0} S (r^{'}, t^{'}) g (r, t; r^{'}, t^{'}) d r^{'} d t^{'} \\ - \int_{Γ_{s}}^{t} \int_{0} {[\frac{c}{C_{R} D} \bar{ϕ} (r_{s}, t^{'}) g (r, t; r_{s}, t^{'})] \\ - [\bar{ϕ} (r_{s}, t^{'}) \frac{\partial}{\partial n} g (r, t; r_{s}, t^{'})]} d r_{s} d t^{'}, \end{matrix}

where [TeX:] $\varphi ^* \left({{\bm r}\!,t} \right)$

ϕ^{*} (r, t)

and [TeX:] $\bar \varphi \left({{\bm r}_s, t} \right)$

\bar{ϕ} (r_{s}, t)

, respectively, represent the value of the fluence at point [TeX:] ${\bm r}$

r

inside the sample and at point [TeX:] ${\bm r}_s$

r_{s}

on the surface at specified moment t. Note that to calculate of the fluence at each internal points, we only need to use the solution of the system of equations that were obtained on the boundary, the values of [TeX:] $\bar \varphi _{}$

{\bar{ϕ}}_{}

, which we now know, so Eq. 15 can be used as many times as desired, with only the surface integration required to be calculated for every new internal point.

3. Numerical Results

The diffuse equation is used to study propagation of a short-pulse laser in biological tissues. We use the BIM to solve diffuse equation and to study behavior of light inside the biological tissues.

First, to confirm precision of this method, the results obtained by the BIM are compared with those obtained analytically using Eq. 8 and the MC and FDTD methods. An interesting parameter in vivo time domain tomography is 〈ct〉, i.e., the mean distance traversed by photons before exiting the tissues. Using Eq. 8, 〈ct〉 can be analytically found to be equal to 77.6 mm for a 10-mm slab tissue with [TeX:] $a = 0.00434\ {\rm mm}^{{\rm - 1}}$ $a = 0.00434 {mm}^{- 1}$ , [TeX:] $\sigma = 6\ {\rm mm}^{{\rm - 1}}$ $σ = 6 {mm}^{- 1}$ , and g = 0.72. Patterson calculated this quantity as 80.6 mm, an error of 3.8% by the MC method.²⁷ The same quantity is calculated as method 81.6 mm, an error of only 5.1%, and 89.3 mm, an error of 15%, using the BIM and FDTD method, respectively. However, the computational time of the FDTD method is more 4 time longer than the BIM, whereas the mesh resolution in the FDTD method and the BIM are 2.00 pS × 0.45 mm and 12.19 pS × 1.36 mm, respectively.

Note that Eqs. 7, 8 can only be used for point sources. However, it is not appropriate to study propagation of a pulse originating from sources that are temporally extended. We assume a semi-infinite sample with [TeX:] $a = 0.02\ {\rm mm}^{{\rm - 1}}$ $a = 0.02 {mm}^{- 1}$ , [TeX:] $\sigma = 10\ {\rm mm}^{{\rm - 1}}$ $σ = 10 {mm}^{- 1}$ , and g = 0.9, which is illuminated by a Gaussian pulse that is presented by

Eq. 16

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} I\left({{\bm r}\!,t} \right) = I_0 \,{\rm exp}\,\left[ { - \left({\frac{{t - t_0 }}{\tau }} \right)^2 } \right]{\rm }\delta \left({\bm r} \right), \end{equation}\end{document}

I (r, t) = I_{0} \exp [- {(\frac{t - t_{0}}{τ})}^{2}] δ (r),

with duration of τ = 10.0 ps. In this paper, we assume resolution of the sample is 12.19 pS × 1.36 mm, which gives the number of 51 × 41 nodes (Fig. 1).

Fig. 1

Temporal evolution of diffuse photon intensity calculated using the BIM at different depths.

Next, we study the effect of different parameters on the temporal evolution of a diffusely penetrating pulse into a phantom like breast tissue with optical properties similar to that represented in Ref. 28. Temporal broadening of the penetrating pulse is due to scattering and absorption of diffused photons in tissue. Therefore, it is important to study such an effect in more detail. To study the effect of multiple scattering on the diffused pulse, we assume a semi-infinite sample with [TeX:] $a = 0.02\ {\rm mm}^{{\rm - 1}}$ $a = 0.02 {mm}^{- 1}$ to be illuminated by a Gaussian pulse with duration of τ = 10.0 ps. Figure 2 shows that the peak value of the pulse decreases for larger scattering coefficients. This is because by increasing the scattering coefficient, more photons are scattered from the illumination direction, and consequently, the intensity of penetrating pulse decreases. The results are in agreement with those reported in Ref. 29.

Fig. 2

Temporal distribution of diffuse photon intensity for different scattering coefficient values.

Figure 3 illustrates the temporal evolution of a penetrating pulse for two different values of the anisotropic factors g. One can see that the peak value of the pulse increases for larger anisotropic factor values. This is because for a larger anisotropic factor, the majority of launched photons are scattered along the illumination direction, and therefore, its intensity increases.

Fig. 3

Temporal distribution of diffuse photon intensity for two different anisotropic factor values.

Furthermore, to study effect of the absorption coefficient on the penetrating pulse, the temporal distribution of penetrating pulse is calculated and the results for three different absorption coefficients values are calculated and depicted in Fig. 4. The results show that by increasing the absorption coefficient, more diffused photons are absorbed in tissue and the density of the remaining photons decreases.

Fig. 4

Temporal distribution of diffuse photon intensity for three different absorption coefficient values.

As mentioned earlier, using this technique, the temporal evolution of the source can take arbitrary shape, and we also consider a source emitting pulses with the following temporal evolution:

Eq. 17

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} I ({{\bm r}\!,t}) = I_0 \Bigg\{ {\rm exp}\!\left[ {-} \!\left(\frac{t - t_0 }{\tau } \right)^2 \right] + {\rm exp}\!\left[ - \left (\frac{t - 2t_0 }{2\tau } \right) \!\right]\!\Bigg\} \delta ({\bm r}),\nonumber\\ \end{eqnarray}\end{document}

\begin{matrix} I (r, t) = I_{0} \{\exp [- {(\frac{t - t_{0}}{τ})}^{2}] + \exp [- (\frac{t - 2 t_{0}}{2 τ})]\} δ (r), \end{matrix}

and the penetration of arbitrary temporal source is calculated by the BIM (Fig. 5). The shown results in Fig. 5 have been calculated for three different depths. This result shows the BIM can simulate the diffused photons inside biological sample with arbitrary temporal evaluation.

Fig. 5

Variation of the diffuse photon intensity with z for an arbitrary input source.

4. Conclusion

Recently, the diffuse equation was solved by the BIM (Refs. 24 to 26), and note that another paper³⁰ solved the transient diffuse equation for turbine component by means of the BEM.³⁰ In Ref. 26, the BIM solution of transient diffuse equation for biological tissues was presented and the accuracy of this method was inspected and the showed resulting emphasizing that the BIM method is faster than the MC and FDTD methods and more precise than the FDTD method. In Ref. 26, only the reflected pulse on the surface of the tissue was calculated.

In this paper, the propagation of short pulses inside biological tissue such as breast was studied. The BIM was used to study the diffusion of short pulses penetrating inside biological samples and to calculate the photon density at any arbitrary point. The effects of the optical properties of tissues on diffusely penetrating pulses inside it were studied and it was observed those parameters can alter the intensity and temporal distribution of the penetrating pulse. Since the optical properties of malignant tissues gradually changes, this offers an effective technique to trace the progress of cancer.

Using this method, the propagation of pulse emitted from extended sources with an arbitrary temporal evolution can be studied, where in this paper, sources with Gaussian shape were studied.

This method is suitable to study the interaction of a short-pulse laser with biological tissue, especially for regions where shock waves are produced.

References

1.

F. Liu, K. M. Yoo, and R. R. Alfano, “Ultrafast laser-pulse transmission and imaging through biological tissues,” Appl. Opt., 32 554 –558 (1993). https://doi.org/10.1364/AO.32.000554 Google Scholar

2.

M. Sakami, K. Mitra, and T. Vo-Dinh, “Analysis of short-pulse laser photon transport through tissues for optical tomography,” Opt. Lett., 27 336 –338 (2002). https://doi.org/10.1364/OL.27.000336 Google Scholar

3.

M. Braun, P. Gilch, and W. Zinth, Ultrashort Laser Pulse in Biology and Medicine, Springer, New York (2007). Google Scholar

4.

E. Alerstam, S. Andersoon-Engels, and T. Svensson, “``Improved accucuracy in timedomain diffuse reflectance spectroscopy,” Opt. Express., 16 10434 –10448 (2008). https://doi.org/10.1364/OE.16.010440 Google Scholar

5.

F. Martelli, A. Sassarolli, S. Del Bianco, and G. Zaccanti, “Solution of time-dependent diffuse equation for a three-layer medium: application to study photon migration through a simplified adult head model,” Phys. Med. Biol., 52 2827 –2843 (2007). https://doi.org/10.1088/0031-9155/52/10/013 Google Scholar

6.

Z. Guo and S. Kumar, “Equitant isotropic scattering formulation for transient short-pulse radiative transfer in anisotropic scattering planar media,” Appl. Opt., 39 4411 –4417 (2000). https://doi.org/10.1364/AO.39.004411 Google Scholar

7.

A. Trivedi, S. Basu, and K. Mitra, “Time-resolved optical tomography using short pulse laser for tumor detection,” J. Quantum Spectrosc. Radiat. Transf., 93 337 –348 (2005). https://doi.org/10.1016/j.jqsrt.2004.08.033 Google Scholar

8.

J. C. Chai, “Time-resolved optical tomography using short pulse laser for tumor detection,” Opt. Rev., 10 609 –610 (2003). https://doi.org/10.1007/s10043-003-0609-3 Google Scholar

9.

C. Das, A. Trivedi, K. Mitra, and T. Vo-Dinhg, “Short pulse laser propagation through tissues for biomedical imagin,” J. Phys. D Appl. Phys., 38 1714 –1721 (2003). https://doi.org/10.1088/0022-3727/36/14/310 Google Scholar

10.

K. Mitra and S. Kumar, “Development and comparison of models for light pulse transport through scattering absorbing media,” Appl. Opt., 38 188 –196 (1999). https://doi.org/10.1364/AO.38.000188 Google Scholar

11.

D. J. Hall, J. C. Hebden, and D. T. Delpy, “Imaging very-low-contrast objects in breastlike scattering media with a time resolved method,” Appl. Opt., 36 7270 –7276 (1997). https://doi.org/10.1364/AO.36.007270 Google Scholar

12.

K. M. Yoo and R. R. Alfano, “Time resolved coherent and incoherent component of forward light scattering in random media,” Opt. Lett., 15 320 –322 (1990). https://doi.org/10.1364/OL.15.000320 Google Scholar

13.

G. Pal, S. Basu, K. Mitra, and T. Vo-Dinh, “Time-resolved optical tomography using short-pulse laser for tumor detection,” Appl. Opt., 45 6270 –6282 (2006). Google Scholar

14.

C. Y. Hopen, “Optical pulse propagation through a slab of random medium,” Waves Rand. Med., 9 551 –560 (1999). https://doi.org/10.1088/0959-7174/9/4/307 Google Scholar

15.

D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express., 10 159 –170 (2002). Google Scholar

16.

L. V. Wang and H. I. Wu, Biomedical Optics: Principles and Imaging, Wiley-Interscience, Hoboken, NJ (2007). Google Scholar

17.

Y. Hasegwa, Y. Yamada, M. Tamura, and Y. Nomura, “Monte Carlo simulation of light transmission through living tissues,” Appl. Opt., 30 4515 –4520 (1991). https://doi.org/10.1364/AO.30.004515 Google Scholar

18.

Z. Guo, S. Kumar, and K. C. San, “``Multidimensional Monte Carlo simulation of short-pulse laser transport in scattering media,” J. Thermophys. Heat Transf., 14 504 –511 (2000). https://doi.org/10.2514/2.6573 Google Scholar

19.

A. H. Hielscher, R. L. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol., 43 1285 –1302 (1998). Google Scholar

20.

Z. Guo, J. Aber, B. A. Garetz, and S. Kumar, “Monte Carlo simulation and experiments of pulsed radiative transfer,” J. Quant. Spectrosc. RA, 73 159 –168 (2002). https://doi.org/10.1016/S0022-4073(01)00203-5 Google Scholar

21.

T. Tanifuji and M. Hijikata, “Finite difference time domain (FDTD) analysis of optical pulse responses in biological tissues for spectroscopic diffused optical tomography,” IEEE. Trans. Med. Imaging, 21 181 –184 (2002). https://doi.org/10.1109/42.993136 Google Scholar

22.

A. Sassaroli, F. Martelli, D. Imai, and Y. Yamada, “Study on the propagation of ultra-short pulse light in cylindrical optical phantoms,” Phys. Med. Biol, 44 2747 –2763 (1999). https://doi.org/10.1088/0031-9155/44/11/305 Google Scholar

23.

F. Paris and J. Canas, Boundary Element Method: Fundamentals and Applications, Oxford University Press, Oxford (1997). Google Scholar

24.

S. Srinivasan, B. W. Pogue, C. Carpenter, P. K. Yalavarthy, and K. Paulsen, “Estimation,” Med. Phys., 34 4545 –4557 (2007). https://doi.org/10.1118/1.2795832 Google Scholar

25.

M. A. Ansari and R. Massudi, “Study of light propagation in Asian and Caucasian skins by means of the boundary element method,” Opt. Laser. Eng., 47 965 –970 (2009). https://doi.org/10.1016/j.optlaseng.2009.04.006 Google Scholar

26.

M. A. Ansari and R. Massudi, “Study of short pulse laser in biological tissue by means of boundary element method,” Google Scholar

27.

M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt., 28 2331 –2336 (1989). https://doi.org/10.1364/AO.28.002331 Google Scholar

28.

N. Ghosh, S. K. Mohanty, S. K. Majumder, and P. K. Gupta, “Measurement of optical transport properties of normal and malignant human breast tissue,” Appl. Opt., 40 176 –184 (2001). https://doi.org/10.1364/AO.40.000176 Google Scholar

29.

P. Rath, S. C. Mishra, P. Mahanta, U. K. Saha, and K. Mitra, “Discrete transfer method applied to transient radiative transfer problems in participating medium,” Numer. Heat Transfer A, 44 183 –197 (2003). Google Scholar

30.

D. A. S. Curran, B. A. Lewis, and M. Cross, “A boundary element method for the solution of the transient diffusion equation in two dimensions,” Appl. Math. Model., 10 107 –113 (1986). https://doi.org/10.1016/0307-904X(86)90080-6 Google Scholar

Citation Download Citation

Mohammad Ali Ansari and Reza Massudi "Boundary integral method for simulating laser short-pulse penetration into biological tissues," Journal of Biomedical Optics 15(6), 065009 (1 November 2010). https://doi.org/10.1117/1.3526351

Published: 1 November 2010

Access the abstract

JOURNAL ARTICLE
5 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 7 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Tissues

Photons

Finite-difference time-domain method

Absorption

Scattering

Monte Carlo methods

Optical properties

1.

Introduction

2.

Review of Theory