Modeling of Laser-Irradiated Biological Tissue

Abstract

The laser has been widely used in medical fields. One application of the laser is laser-based photo-thermal therapy, wherein the short-pulsed laser is generally used to destroy the cancerous cells. The efficacy of the laser-based photo-thermal therapy can be improved if we minimize the thermal damage to the surrounding healthy tissue. So, it is essential to understand the laser-tissue interaction and thermal behavior of biological tissue during laser-based photo-thermal therapy. The light propagation through the biological tissue is generally mathematically modeled by the radiative heat transfer equation (RTE). The RTE has been solved using the discrete ordinate method (DOM) to determine the intensity inside the laser-irradiated biological tissue. Consequently, the absorbed photon energy act as the source term in the Fourier/non-Fourier model-based bio-heat transfer equation to determine the temperature distribution inside the biological tissue subjected to short-pulse laser irradiation. The non-Fourier model-based bio-heat transfer equation is numerically solved using the finite volume method (FVM). The numerical results have been compared with the analytical results obtained using the finite integral transform (FIT) technique. A comparative study between the Fourier and non-Fourier heat conduction models has also been carried out.

1. Introduction

Cancer is a group of diseases in which uncontrolled growth and invade of abnormal cells into other parts of the body [1]. The World Health Organization says nearly one in six died from cancer in 2020 [2]. Thus, it is essential to detect the cancerous cells at their early stage and destroy them to minimize possible damage to the surrounding healthy tissue. So, various diagnosis techniques such as computerized tomography (C.T.) scan, magnetic resonance imaging (MRI), ultrasound, X-ray, positron emission tomography (PET) scan, etc., and treatment such as surgery, chemotherapy, laser-based photo-thermal therapy, etc., have been developed in the past.

Laser-based photo-thermal therapy has gained unprecedented growth among all available treatment techniques in the past few decades. In this therapy, the temperature of the tissue is increased above the pre-defined threshold value using the laser, which destroys the cancerous cells [3]. In other words, the photon energy has been absorbed by the laser-irradiated biological tissue, which subsequently increases its internal energy; as a result, the temperature of the tissue increases. Amongst the various types of laser, the short-pulse laser has gained substantial importance because it can transmit high energy in a very short interval of time (order of femtosecond to picosecond) [3]. This advantage helps in increasing the temperature to the desired level in the confined region. Therefore, the short-pulse laser has been usually used to destroy cancerous cells during laser-based photo-thermal therapy. It is worth noting that tissue gets thermally damaged when its temperature crosses 43°C [4]. So, the major challenge in the medical field is to destroy the cancerous cells without damaging the surrounding healthy tissue. Thus, it is essential to understand the thermal characteristics of the laser-irradiated biological tissue to improve the efficacy of laser-based photo-thermal therapy.

The modeling of laser-tissue interaction is challenging because it is absorbing as well as scattering in nature due to its constituents such as water, hemoglobin, melanin RBCs, cell membrane, etc. [3, 5]. So, Various mathematical models such as Beer-Lambert's law, diffusion approximation theory, radiative transfer equation (RTE), etc., have been developed to model the laser-tissue interaction. Among these mathematical models, the Beer-Lambert’s law $\left(I=I_{o}exp\left(-\mathit{\kappa z}\right)\right)$is the simple and relatively straightforward mathematical model, which is solved by the researchers to determine the intensity ($I$) variation inside the laser-irradiated biological tissue [6, 7]. Here, $\kappa$represents the absorption coefficient. However, the limitation of this model is that it assumes that the biological tissue is purely absorbing in nature. Few researchers used the diffusion approximation theory to determine the intensity distribution by considering that the biological tissue is highly scattering in nature [8, 9]. So, assuming the biological tissue is either purely absorbing or highly scattering misleads the results. Therefore, the researchers [4] have suggested using the RTE (Eq. 1), which is the most appropriate mathematical model for modeling the laser-tissue interaction [10].

$$\frac{1}{c}\frac{\partial I}{\partial t}+\frac{\mathit{dI}}{\mathit{ds}}=-\mathit{\beta I}+\kappa I_b+\frac{{\sigma }_s}{4\pi }{\int }_{4\pi }I\varnothing \left(\mathrm{\Omega },{\mathrm{\Omega }}^{\prime }\right)d{\mathrm{\Omega }}^{\prime }$$E1

where, $c,\beta ,{\sigma }_s$ and $\varnothing$ represents the speed of light in the medium, extinction coefficient (sum of scattering and absorption coefficient), scattering coefficient, and scattering phase function, respectively. The first and second term on the left-hand side of Eq. (1) represents the temporal and spatial variation of intensity, respectively. On the other hand, the first, second, and third terms on the right-hand side of Eq. (1) denote the attenuation of intensity due to the absorption and out-scattering, augmentation of intensity due to the emission and in-scattering, respectively.

Eq. (1) is the integrodifferential equation, and getting the analytical solution is quite challenging due to the presence of the in-scattering terms (third term of the right-hand side of Eq. (1)). So, researchers have developed various numerical methods, such as the discrete transfer method (DTM) [11], finite volume method (FVM) [12], discrete ordinate method (DOM) [13], etc., to convert this integrodifferential equation into a simpler form to get the numerical solution. The various researchers performed the comparative studies between these numerical methods and found that the DOM is computationally most efficient, uses comparatively lesser memory and the effective wave speed predicted by DOM comes closer to the actual speed of light [14, 15, 16]. So, the DOM has been employed in the present study to solve the RTE.

To understand the thermal behavior of the laser-irradiated biological tissue, we need to solve the bio-heat transfer equation (Eq. 2) wherein the energy generated by the laser acts as the source terms.

$$\rho c_v\frac{\partial T}{\partial t}=-\nabla \bullet \overrightarrow{q}+{\omega }_b{\rho }_b{c}_b\left(T_b-T\right)+Q_m-\nabla \bullet \overrightarrow{q_r}$$E2

where $\rho$, $c_v$, ${\omega }_b$, ${\rho }_b$, and $c_b$ denote the density of the tissue, specific heat of the tissue, blood perfusion, density of blood, and specific heat of blood, respectively.

The term on the left-hand side of Eq. (2) represents the time rate of change of energy inside the biological tissue. On the other hand, the first, second, third, and fourth terms on the right-hand side of Eq. (2) represent the diffusion term due to the net heat transfer across the boundaries, blood perfusion term due to the energy exchange between the tissue and surrounding capillary blood vessels, metabolic heat generation term, and the divergence of the radiative heat flux which represents net radiative heat loss from the medium.

When Fourier’s law of heat conduction (Eq. 3) is substituted into Eq. (2), the resulting equation is known as the Fourier model-based bio-heat transfer equation, also known as Pennes bio-heat transfer equation.

$$\overrightarrow{q}=-k\nabla T$$E3

The limitation of the Pennes bio-heat transfer equation is that it assumes that arterial blood temperature is almost uniform throughout the tissue domain while the venous blood temperature equals that of the local tissue temperature [17]. However, the assumption of thermal equilibrium between the blood vessel and the surrounding tissue is not valid when the diameter of the blood vessel is more than 500 μm (i.e., large blood vessel). In that case, the energy equation has been separately solved for the tissue and blood vessel domain to get the accurate temperature distribution. The interested reader can find the numerical solution of the energy equation to determine the temperature distribution inside the laser-irradiated biological tissue embedded with the large blood vessel elsewhere [18].

Few researchers experimentally demonstrated that Fourier's law of heat conduction under predicts the temperature value that does not accurately match the experimental data because of the inherent assumption of the infinite speed of thermal wave propagation through the biological tissue [4, 19, 20]. Furthermore, the Fourier model gives inaccurate results in an area where used the short-pulse laser, the temperature is in the range of cryogenic, studying the thermal response of non-homogeneous structures, e.g., biological tissue [21]. To overcome the limitations of Fourier heat conduction models, researchers have modified Eq. (3) by considering the relaxation time associated with the heat flux (${\tau }_q$) and temperature gradient (${\tau }_T$), and it is known as the dual-phase lag (DPL) model (Eq. 4).

$$\overrightarrow{q}+{\tau }_q\frac{\partial \overrightarrow{q}}{\partial t}=-k\left\{\nabla T+{\tau }_T\frac{\partial \left(\nabla T\right)}{\partial t}\right\}$$E4

When ${\tau }_T$ is equal to zero, the DPL model becomes the Cattaneo-Vernotte (C-V) equation. However, the limitation of this model is that it doesn’t consider the effects of microstructural interactions in a non-homogeneous medium, e.g., biological samples [21]. The DPL model becomes the conventional Fourier’s law of heat conduction when both relaxation times are equal to zero. So, the DPL model is the generalized form of the non-Fourier heat conduction model. Thus, the numerical and analytical solutions of the DPL model-based bio-heat transfer equation have been obtained using the FVM and finite integral technique (FIT), respectively, which are discussed in Sections 3 and 4 respectively.

As mentioned earlier, the last term on the right-hand side of Eq. (2) represents the divergence of radiative heat flux (Eq. (5)), which is obtained by solving the Eq. (1).

$$\nabla \bullet \overrightarrow{q_r}=\kappa \left(4\pi I_b-G\right)=\kappa \left(4\sigma T^4-G\right)$$E5

where, $I_b$ and $G={\sum }_{m=1}^M{\omega }^m{I}^m$represent the black body intensity and the incident intensity. Here ${\omega }^m$denotes the angular weight in the particular direction m, and M is the total number of divisions of the $4\pi$ solid angle.

It is to be noted here that the time scale used for solving the RTE is the order of picosecond because the light propagates through the biological tissue generally has the same order of time scale. However, at such time scales, the contribution of the diffusion term, blood perfusion term, and metabolic heat generation term (the first three terms on the right-hand side of Eq. (2)) to the increment of the temperature of the tissue are insignificant as compared to the divergence of the radiative heat flux. However, the time scale of the order of the millisecond is generally used for solving the bio-heat transfer equation to capture the effect of these terms. Therefore, it is a multi-time scale problem because the two different time scales have been used to solve the RTE and bio-heat transfer equation. The solution procedure for solving this multi-time scale problem has been discussed in Section 2.

2. Problem formulation

A two-dimensional squared-shaped biological tissue having the dimensions of 2 mm $\times$ 2 mm was considered in the present study, which is shown in Figure 1. The short-pulse laser is vertically incident at the top surface of the domain. Short laser pulses have been modeled as Gaussian profile in the space domain and Top-Hat profile in the time domain with a pulse width of t_p. When light propagates through the biological tissue, the intensity has two components: collimated and diffuse. The collimated component of the intensity has been mathematically expressed by Eq. (6) [21].

Figure 1.

Schematic diagram of the physical domain under consideration (dimensions in mm).

$$I_c\left(x,y,\Omega ,t\right)=\begin{array}{l}{I}_{\mathit{\text{cmax}}}\times \mathit{exp}\left\{-{\left(x-L/2\right)}^2/d^2\right\}\times \mathit{exp}\left\{-\beta \left(W-y\right)/\left|{\eta }_c\right|\right\}\\ \times \left[H\left\{\mathit{\beta ct}-\beta \left(W-y\right)/\left|{\eta }_c\right|\right\}-H\left\{\left(\mathit{\beta ct}-\beta \left(W-y\right)/\left|{\eta }_c\right|\right)-\mathit{\beta c}{t}_p\right\}\right]\\ \times \delta \left(\mu -{\mu }_c\right)\delta \left(\eta -{\eta }_c\right)\end{array}$$E6

where, $I_{\mathit{\text{cmax}}}$, $L$, $W$, $\mu ,$ and $\eta$ represent the maximum intensity incident on the top surface of the tissue phantom, the length of tissue, the width of the tissue, and direction cosines in x- and y-directions, respectively. Here, $\delta$ denotes the Dirac-delta function, which considers the angle of incidence, and $H$ is the Heaviside function which ensures that the laser energy is available at any location only for the pulse width of the laser. The first and second exponential terms in Eq. (6) represent the Gaussian profile in the space domain and the attenuation of collimated intensity due to absorption and out-scattering, respectively.

The diffuse component of the intensity is obtained by solving the RTE. So, Eq. (1) is simplified for the two-dimensional rectangular domain and also neglects the emission from the tissue phantom because the intensity of the short-pulse laser used for irradiating the sample is much higher than the blackbody radiation intensity:

$$\frac{1}{c}\frac{\partial I_d}{\partial t}+\mu \frac{\partial I_d}{\partial x}+\eta \frac{\partial I_d}{\partial y}=-\beta I_d+\frac{{\sigma }_s}{4\pi }{\int }_{4\pi }I\varnothing \left(\Omega ,{\Omega }^{\prime }\right)d{\Omega }^{\prime }$$E7

The boundary conditions for a diffusely emitting and reflecting wall at a given point r_w on that surface can be expressed as [21]:

$$I\left(r_w,\widehat{s},t\right)=\in \left(r_w\right)I_b\left(r_w,t\right)+\frac{\rho \left(r_w\right)}{\pi }{\int }_{\widehat{n\bullet }{\widehat{s}}^{\prime }<0}I\left(r_w,{\widehat{s}}^{\prime },t\right)\left|\widehat{n}\bullet {\widehat{s}}^{\prime }\right|d{\Omega }^{\prime },\widehat{n}\bullet \widehat{s}>0$$E8

Here $\widehat{n}$ is the inward normal vector to the surface. The first and second terms on the right-hand side of Eq. (8) represent the emission and reflection from the wall, respectively. The expression for boundary conditions given by Eq. (8), in general, holds valid for all the four walls of the rectangular enclosure.

Once the intensity distribution inside the laser-irradiated biological tissue is obtained, then the divergence of the radiative heat flux (Eq. 5) is calculated. As earlier mentioned, finding the temperature distribution inside the biological tissue subjected to the short-pulse laser is a multi-time scale problem. So, we need to divide the Eq. (2) into Eqs. (9) and (10). Eq. (9) is solved for determining local temperature rise due to a single pulse for which the time scale is the order of picoseconds.

$$\rho c_v\frac{\partial T}{\partial t}=-\nabla \bullet \overrightarrow{q_r}$$E9

To determine the temperature distribution inside the biological tissue subjected to a train of pulses, the temperature rise due to the single pulse is added to the temperature distribution at the previous time instant if the time difference between the two consecutive pulses is equal to the time step employed else it is not added. The temperature distribution at the previous time instant is used as the initial condition to determine the temperature distribution at the current time instant by solving the Eq. (10) [21]. It is to be noted that the time scale used for solving Eq. (10) is the order of milliseconds.

$$\rho c_v\frac{\partial T}{\partial t}=-\nabla \bullet \overrightarrow{q}+{\omega }_b{\rho }_b{c}_b\left(T_b-T\right)+Q_m$$E10

The DPL model-based bio-heat transfer equation has been solved in the present study, which is obtained by eliminating temperature from Eq. (10) using the Eq. (4), and the governing equation is written in terms of heat flux which is expressed as:

$$\frac{{\tau }_q}{\alpha }\frac{{\partial }^2\overrightarrow{q}}{\partial t^2}+\frac{1}{\alpha }\frac{\partial \overrightarrow{q}}{\partial t}=\nabla \left(\nabla \bullet \overrightarrow{q}\right)+{\tau }_T\frac{\partial \left[\nabla \left(\nabla \bullet \overrightarrow{q}\right)\right]}{\partial t}$$E11

It is to be noted here that the blood perfusion and metabolic heat generation terms have been neglected while deriving Eq. (11) because their contribution is negligible to temperature rise [21].

The required boundary conditions for solving the Eq. (10) are as follows: the all three walls except the top wall (shown in Figure 1) are maintained at the core body temperature (37°C), while the top wall is subjected to the convective boundary conditions (heat transfer coefficient of 15 W/m²∙K and ambient temperature of 25°C) [21]. The initial condition (t = 0) for temperature is equal to 37°C, and the two initial conditions are required for Eq. (11) because it has a second-order derivative concerning the time, so heat flux and time derivative of heat flux are assumed to be zero.

Eqs. (10) and (11) have been numerically solved using the FVM, and its solution is given in Section 3.

The following laser parameters have been used in this study: beam diameter (d)= 0.025 mm, amplitude of pulse (I_cmax)= 1.6$\times$10⁻³ J/mm²/ps, pulse width (t_p)= 4.6667 ps, repetition rate (fr)= 1 kHz and wavelength= 1100 nm [21]. The absorption coefficient (κ) and scattering coefficient $\left({\sigma }_s\right)$ of tissue are 0.051 mm⁻¹ and 6.14 mm⁻¹, respectively [4, 21]. The thermo-physical properties of the tissue are as follows: density (ρ)= 1000 kg/m³, thermal conductivity (k)= 0.63 W/m∙K and specific heat (c_v)= 4200 J/kg∙K [21].

3. Numerical modeling

This section discusses the methodology for solving the RTE (Eq. 7) using the DOM. Then, the numerical solution of the DPL model-based bio-heat transfer is obtained using the FVM to determine the temperature distribution inside the laser-irradiated biological tissue.

In DOM, the RTE, an integrodifferential equation, has been transformed into a set of partial differential equations [10]. So, the RTE for the diffuse component of intensity (Eq. (12)) with the corresponding boundary conditions (Eq. 13) in any given discrete direction (${\Omega }^{m,n}$) can be written as follows

$$\begin{array}{l}\frac{1}{c}\frac{\partial I_d^{m,n}}{\partial t}+{\mu }^{m,n}\frac{\partial I_d^{m,n}}{\partial x}+{\eta }^{m,n}\frac{\partial I_d^{m,n}}{\partial y}=-\beta I_d^{m,n}+\frac{{\sigma }_s}{4\pi }\sum _{m^{\prime }=1}^{N_{\theta }}\sum _{n^{\prime }=1}^{N_{\phi }}{I}_d^{m^{\prime },n^{\prime }}{\varnothing }^{m^{\prime },n^{\prime };m,n}{\omega }_{\theta }^{m^{\prime }}{\omega }_{\phi }^{n^{\prime }}\\ +\frac{{\sigma }_s}{4\pi }{I}_c^{m_c,n_c}{\varnothing }^{m_c,n_c;m,n}\end{array}$$E12

$$\text{Left wall}\left(\mathit{x}=0,0<\mathit{y}<\mathit{W}\right):I_d^{m,n}=0,{\mu }^{m,n}>0$$E13

$$\text{Right wall}\left(\mathit{x}=\mathit{L},0<\mathit{y}<\mathit{W}\right):I_d^{m,n}=0,{\mu }^{m,n}<0$$E14

$$\text{Bottom wall}\left(\mathit{y}=0,0<\mathit{x}<\mathit{L}\right):I_d^{m,n}=0,{\eta }^{m,n}>0$$E15

$$\mathrm{Top}\text{wall}\left(\mathit{y}=\mathit{W},0<\mathit{x}<\mathit{L}\right):I_d^{m,n}=0,{\eta }^{m,n}<0$$E16

The boundary conditions Eqs. (13)–(16) have been obtained using Eq. (8) under the assumption that the boundary is non-reflecting and the magnitude of the light intensity emitted from the surface is negligible compared to the intensity of the short-pulse laser used.

Using a fully implicit backward differencing scheme in time, Eq. (12), after simplification, becomes

$${\mu }^{m,n}\frac{\partial I_d^{m,n}\left(t^{\ast }\right)}{\partial x}+{\eta }^{m,n}\frac{\partial I_d^{m,n}\left(t^{\ast }\right)}{\partial y}+\left[\frac{\beta }{∆t^{\ast }}+\beta \right]I_d^{m,n}\left(t^{\ast }\right)=\frac{\beta }{∆t^{\ast }}{I}_d^{m,n}\left(t^{\ast }-∆t^{\ast }\right)+S_t^{m,n}$$E17

where $t^{\ast }=\mathit{\beta ct},$

$$S_t^{m,n}=\frac{{\sigma }_s}{4\pi }\sum _{m^{\prime }=1}^{N_{\theta }}\sum _{n^{\prime }=1}^{N_{\phi }}{I}_d^{m^{\prime },n^{\prime }}{\varnothing }^{m^{\prime },n^{\prime };m,n}{\omega }_{\theta }^{m^{\prime }}{\omega }_{\phi }^{n^{\prime }}+\frac{{\sigma }_s}{4\pi }{I}_c^{m_c,n_c}{\varnothing }^{m_c,n_c;m,n}$$E18

Integrating Eq. (17) over the control volume ΔV leads to

$$\begin{array}{l}{\mu }^{m,n}∆y\left[I_{d,E}^{m,n}\left(t^{\ast }\right)-I_{d,W}^{m,n}\left(t^{\ast }\right)\right]+{\eta }^{m,n}∆x\left[I_{d,N}^{m,n}\left(t^{\ast }\right)-I_{d,S}^{m,n}\left(t^{\ast }\right)\right]\\ =-\left[\frac{\beta }{∆t^{\ast }}+\beta \right]\mathit{\Delta V}{I}_{d,P}^{m,n}\left(t^{\ast }\right)+\mathit{\Delta V}{S}_{t,P}^{m,n}+\frac{\mathit{\beta \Delta V}}{∆t^{\ast }}{I}_{d,P}^{m,n}\left(t^{\ast }-∆t^{\ast }\right)\end{array}$$E19

where, $I_{d,P}^m$ is the intensity at the cell center P, $S_{t,P}^m$ is the source term at the cell center P and $∆V=∆x\times ∆y$. For relating the cell-surface intensities with cell-center intensity, the step differencing scheme has been used to avoid any possibility of providing unphysical results [10].

Eq. (19) is simplified using the step differencing scheme, and the resulting equation becomes as:

$$I_{d,P}^{m,n}\left(t^{\ast }\right)=\frac{\left[{\mu }^{m,n}∆y{I}_{d,W}^{m,n}\left(t^{\ast }\right)+{\eta }^{m,n}∆x{I}_{d,S}^{m,n}\left(t^{\ast }\right)+∆V{S}_{t,P}^{m,n}+\left(\frac{\beta ∆V}{∆t^{\ast }}\right)I_{d,P}^{m,n}\left(t^{\ast }-∆t^{\ast }\right)\right]}{{\mu }^{m,n}∆y+{\eta }^{m,n}∆x+\left[\frac{\beta }{∆t^{\ast }}+\beta \right]\mathit{\Delta V}}$$E20

As earlier mentioned, the solution of RTE provides the intensity distribution, which is used to determine the divergence of the radiative heat flux (Eq. 5). Subsequently, find the temperature rise due to a single pulse by solving the Eq. (9). Then, Eq. (10) has been solved using the FVM to determine the temperature distribution inside the biological tissue subjected to the train of pulses. The algorithm for solving this multi-time scale problem has been discussed in Section 2.

Under the Cartesian coordinates system, the governing equation (Eq. 11) in x- directions can be expressed as [21]

$$\frac{{\tau }_q}{\alpha }\frac{{\partial }^2{q}_x}{\partial t^2}+\frac{1}{\alpha }\frac{\partial q_x}{\partial t}=\frac{\partial }{\partial x}\left(\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}\right)+{\tau }_T\frac{\partial }{\partial t}\left\{\frac{\partial }{\partial x}\left(\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}\right)\right\}$$E21

Integrating Eq. (21) over a given control volume and over the time step of Δt, i.e., from t to t+Δt followed by discretization based on the backward difference in time leads to the following discretized form of the Eq. (22) [21].

$${a_x}_P{q_x}_P^{t+∆t}=a_W{q_x}_W^{t+∆t}+a_E{q_x}_E^{t+∆t}+b_x$$E22

where$a_W=\frac{∆y}{{\left(\mathit{\delta x}\right)}_w}\left(∆t+{\tau }_T\right)$, $a_E=\frac{∆y}{{\left(\mathit{\delta x}\right)}_e}\left(∆t+{\tau }_T\right)$, ${a_x}_P=\frac{∆x∆y}{\alpha }\left(1+\frac{{\tau }_q}{∆t}\right)+a_W+a_E$ ,

$$\begin{array}{l}{b}_x=\left(\frac{∆x∆y}{\alpha }\left(1+\frac{2{\tau }_q}{∆t}\right)+\frac{∆y{\tau }_T}{{\left(\mathit{\delta x}\right)}_w}+\frac{∆y{\tau }_T}{{\left(\mathit{\delta x}\right)}_e}\right){q_x}_P^t-\frac{∆y{\tau }_T}{{\left(\mathit{\delta x}\right)}_w}{q_x}_W^t-\frac{∆y{\tau }_T}{{\left(\mathit{\delta x}\right)}_e}{q_x}_E^t\\ -\frac{∆x∆y{\tau }_q}{\alpha ∆t}{q_x}_P^{t-∆t}+\left(∆t+{\tau }_T\right)\left({q_y}_{\mathit{ne}}^{t+∆t}-{q_y}_{\mathit{se}}^{t+∆t}-{q_y}_{\mathit{nw}}^{t+∆t}+{q_y}_{\mathit{sw}}^{t+∆t}\right)\\ -{\tau }_T\left({q_y}_{\mathit{ne}}^t-{q_y}_{\mathit{se}}^t-{q_y}_{\mathit{nw}}^t+{q_y}_{\mathit{sw}}^t\right)\end{array}$$

Here, the symbols t+Δt, t, and t-Δt represent the current time step, previous time step, and previous to previous time step, respectively.

Similarly, the discretized form of the equation for the heat flux component in the y-direction of Eq. (11) can be obtained.

Once the heat flux component in both directions has been obtained, the Eq. (10) has been discretized using the FVM to determine the temperature distribution at the current time instant using Eq. (23) [21].

$$T_P^{t+∆t}=T_P^t+\frac{∆t}{\rho c_v∆x}\left({q_x}_w^{t+∆t}-{q_x}_e^{t+∆t}\right)+\frac{∆t}{\rho c_v∆y}\left({q_y}_s^{t+∆t}-{q_y}_n^{t+∆t}\right)$$E23

The obtained numerical results have been compared with the analytical solution. The analytical solution of the DPL model-based bio-heat transfer equation has been obtained using the FIT technique, which is discussed in Section 4.

4. Analytical modeling

In this section, the analytical solution of the DPL model-based bio-heat transfer equation is obtained using the FIT technique. Then, the algorithm for coupling this analytical solution with a solution of RTE is discussed to determine the laser-irradiated biological tissue.

To determine the temperature distribution inside the laser-irradiated biological tissue, the DPL model-based bio-heat transfer equation in terms of temperature is obtained using Eqs. (4) and (10), which can be expressed as [22]:

$$\rho c_v{\tau }_q\frac{{\partial }^{2}T}{{\partial t}^2}+\left(\rho c_v+{\tau }_q{\omega }_b{\rho }_b{c}_b\right)\frac{\partial T}{\partial t}=\begin{array}{l}k\left[\nabla \bullet \left(\nabla T\right)+{\tau }_T\frac{\partial \left[\nabla \bullet \left(\nabla T\right)\right]}{\partial t}\right]+{\omega }_b{\rho }_b{c}_b\left(T_b-T\right)\\ +{\tau }_q{\omega }_b{\rho }_b{c}_b\frac{\partial T_b}{\partial t}+Q_m+{\tau }_q\frac{\partial Q_m}{\partial t}\end{array}$$E24

Under the assumptions of T_b = constant=37 ^oC, Q_m = constant and$\theta$ defined as $T-T_b$, the governing equation (Eq. 24) in the two-dimensional Cartesian coordinate system can be expressed as [22]:

$$\begin{array}{l}\rho c_v{\tau }_q\frac{{\partial }^2\theta }{{\partial t}^2}+\left(\rho c_v+{\tau }_q{\omega }_b{\rho }_b{c}_b\right)\frac{\partial \theta }{\partial t}=k\left[\frac{{\partial }^2}{{\partial x}^2}\left(\theta +{\tau }_T\frac{\partial \theta }{\partial t}\right)+\frac{{\partial }^2}{{\partial y}^2}\left(\theta +{\tau }_T\frac{\partial \theta }{\partial t}\right)\right]\\ -{\omega }_b{\rho }_b{c}_b\theta +Q_m\end{array}$$E25

Eq. (24) can also be expressed in a cylindrical coordinate system for axisymmetric biological tissue by Eq. (26) [23].

$$\begin{array}{l}\rho c_v{\tau }_q\frac{{\partial }^2\theta }{{\partial t}^2}+\left(\rho c_v+{\tau }_q{\omega }_b{\rho }_b{c}_b\right)\frac{\partial \theta }{\partial t}=k[\frac{{\partial }^2\theta }{\partial r^2}+\frac{1}{r}\frac{\partial \theta }{\partial r}+\frac{{\partial }^2\theta }{\partial z^2}\\ +{\tau }_T\frac{\partial }{\partial t}\left(\frac{{\partial }^2\theta }{\partial r^2}+\frac{1}{r}\frac{\partial \theta }{\partial r}+\frac{{\partial }^2\theta }{\partial z^2}\right)]-{\omega }_b{\rho }_b{c}_b\theta +Q_m\end{array}$$E26

The FIT technique [24] converts the partial differential equation into the ordinary differential equation by removing the space derivative using the integral transformation pair. Thus, the ordinary differential equation is only the function of time. This equation with the transformed initial conditions is solved. The transformed temperature can be inverted back into the temperature using the inversion formula.

The boundary conditions Eqs. (27)–(30) and initial conditions Eqs. (31) and (32) required for solving Eq. (25) are given below:

$$\theta \left(0,y,t\right)=0$$E27

$$\theta \left(L,y,t\right)=0$$E28

$$\theta \left(x,0,t\right)=0$$E29

$$\frac{\partial \theta \left(x,L_y,t\right)}{\partial y}+\mathit{H\theta }\left(x,L_y,t\right)=H{\theta }_{\infty }$$E30

$$\theta \left(x,y,0\right)=0$$E31

$$\frac{\partial \theta \left(x,y,0\right)}{\partial t}=0$$E32

$$\text{Here}H=\frac{h}{k}.$$

The integral transform pair for a given function $\theta \left(x,y,t\right)$concerning the x variable can be expressed as below

Inversion formula:

$$\theta \left(x,y,t\right)=\sum _{p=1}^{\infty }\frac{X\left({\nu }_p,x\right)}{N\left({\nu }_p\right)}\overline{\theta }\left({\nu }_p,y,t\right)$$E33

Integral transform:

$$\overline{\theta }\left({\nu }_p,y,t\right)={\int }_{x^{\prime }=0}^{\mathrm{L}}X\left({\nu }_p,x^{\prime }\right)\theta \left(x^{\prime },y,t\right)d{x}^{\prime }$$E34

Similarly, the integral transform pair for the function $\overline{\theta }\left({\nu }_p,y,t\right)$ concerning the y variable can be expressed as

Inversion formula:

$$\overline{\theta }\left({\nu }_p,y,t\right)=\sum _{m=1}^{\infty }\frac{Y\left({\eta }_m,y\right)}{N\left({\eta }_m\right)}\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t\right)$$E35

Integral transform:

$$\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t\right)={\int }_{y^{\prime }=0}^{\mathrm{W}}Y\left({\eta }_m,y^{\prime }\right)\overline{\theta }\left({\nu }_p,y^{\prime },t\right)d{y}^{\prime }$$E36

Here, ${\nu }_p$,${\eta }_m$ are eigenvalues, $X\left({\nu }_p,x\right)$ , $Y\left({\eta }_m,y\right)$ are Eigen function and $N\left({\nu }_p\right)$, $N\left({\eta }_m\right)$ are the norm. The eigenvalues, Eigen function, and norm for the givens set boundary conditions Eqs. (27)–(30) can be found elsewhere [24].

Multiplying both sides of Eq. (25) by the operator${\int }_{x^{\prime }=0}^{L}X\left({\nu }_p,x^{\prime }\right)d{x}^{\prime }$ and using the Eq. (34) and Eqs. (27) and (28), the derivative concerning x variables is removed from Eq. (25), and the resulting equation is in terms of $\overline{\theta }\left({\nu }_p,y,t\right)$. Similarly, the derivative concerning y variables is removed from this transformed governing equation by multiplying both sides by ${\int }_{y^{\prime }=0}^{W}Y\left({\eta }_m,y^{\prime }\right)d{y}^{\prime }$ and using the Eq. (36) and the transformed boundary conditions of Eqs. (29) and (30), we get [22]

$$\frac{d^2\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t\right)}{{\mathit{dt}}^2}+a\frac{d\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t\right)}{\mathit{dt}}+b\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t\right)=\frac{1}{\rho c_v{\tau }_q}\left(\mathit{kH}{\theta }_{\infty }sin{\eta }_m{L}_y\frac{\left[1-{\left(-1\right)}^p\right]}{{\nu }_p}\right)$$E37

where $a=\frac{k{\tau }_T\left({\nu }_p^2+{\eta }_m^2\right)+\left(\rho c_v+{\tau }_q{\omega }_b{\rho }_b{c}_b\right)}{\rho c_v{\tau }_q}$ , $b=\frac{k\left({\nu }_p^2+{\eta }_m^2\right)+{\omega }_b{\rho }_b{c}_b}{\rho c_v{\tau }_q}$

Eq. (37) is the second-order non-homogenous ordinary differential equation, and its solution using the transformed initial conditions has been obtained, which is given below [22]:

$$\begin{array}{l}\theta \left(x,y,t_{n+1}\right)=\sum _{p=1}^{\infty }H{\theta }_{\infty }\frac{\left[1-{\left(-1\right)}^p\right]}{{\nu }_p}\frac{\mathit{sinh}{b}_1\left(p\right)y}{b_1\left(p\right)\mathit{cosh}{b}_1\left(p\right)W+\mathit{\text{Hsinhh}}{b}_1\left(p\right)W}\frac{\mathit{sin}{\nu }_{p}x}{N\left({\nu }_p\right)}\\ +\sum _{p=1}^{\infty }\sum _{m=1}^{\infty }\left[\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t_n\right)-\left(\frac{\mathit{kH}{\theta }_{\infty }sin{\eta }_{m}W}{\rho c_v{\tau }_{q}b}\right)\left[\frac{1-{\left(-1\right)}^p}{{\nu }_p}\right]\right]\left[\left(\frac{\gamma +{\omega }_1}{2{\omega }_1}\right)\mathit{exp}\left\{-\left(\gamma -{\omega }_1\right)∆t\right\}\right.\\ \left.-\left(\frac{\gamma -{\omega }_1}{2{\omega }_1}\right)\mathit{exp}\left\{-\left(\gamma +{\omega }_1\right)∆t\right\}\right]\frac{\mathit{sin}{\nu }_{p}x}{N\left({\nu }_p\right)}\frac{\mathit{sin}{\eta }_{m}y}{N\left({\eta }_m\right)}\\ +\sum _{p=1}^{\infty }\sum _{m=1}^{\infty }\frac{d\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t_n\right)}{\mathit{dt}}\left[\frac{\mathit{exp}\left\{-\left(\gamma -{\omega }_1\right)∆t\right\}-\mathit{exp}\left\{-\left(\gamma +{\omega }_1\right)∆t\right\}}{2{\omega }_1}\right]\frac{\mathit{sin}{\nu }_{p}x}{N\left({\nu }_p\right)}\frac{\mathit{sin}{\eta }_{m}y}{N\left({\eta }_m\right)},\\ {\left(\frac{a}{2}\right)}^2>b\end{array}$$E38

$$\begin{array}{l}\theta \left(x,y,t_{n+1}\right)=\sum _{p=1}^{\infty }H{\theta }_{\infty }\frac{\left[1-{\left(-1\right)}^p\right]}{{\nu }_p}\frac{\mathit{sinh}{b}_1\left(p\right)y}{b_1\left(p\right)\mathit{cosh}{b}_1\left(p\right)W+\mathit{\text{Hsinhh}}{b}_1\left(p\right)W}\frac{\mathit{sin}{\nu }_{p}x}{N\left({\nu }_p\right)}\\ +\sum _{p=1}^{\infty }\sum _{m=1}^{\infty }\left[\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t_n\right)-\left(\frac{\mathit{kH}{\theta }_{\infty }sin{\eta }_{m}W}{\rho c_v{\tau }_{q}b}\right)\left[\frac{1-{\left(-1\right)}^p}{{\nu }_p}\right]\right]\left[\frac{\mathit{exp}\left(-\gamma ∆t\right)}{{\omega }_2}\right]\left[{\omega }_{2}cos{\omega }_2∆t\right.\\ \left.+\gamma sin{\omega }_2∆t\right]\frac{\mathit{sin}{\nu }_{p}x}{N\left({\nu }_p\right)}\frac{\mathit{sin}{\eta }_{m}y}{N\left({\eta }_m\right)}+\sum _{p=1}^{\infty }\sum _{m=1}^{\infty }\frac{d\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t_n\right)}{\mathit{dt}}\left[\frac{sin{\omega }_2∆\mathit{\text{texp}}\left(-\gamma ∆t\right)}{{\omega }_2}\right]\frac{\mathit{sin}{\nu }_{p}x}{N\left({\nu }_p\right)}\frac{\mathit{sin}{\eta }_{m}y}{N\left({\eta }_m\right)},\\ {\left(\frac{a}{2}\right)}^2<b\end{array}$$E39

where ${\omega }_1=\sqrt{{\left(\frac{a}{2}\right)}^2-b}$, and ${\omega }_2=\sqrt{b-{\left(\frac{a}{2}\right)}^2}$

It is to be noted here that the subscript “n” denotes the previous time instant while subscript “n+1” represents the current time instant. Here, $∆t=t_{n+1}-t_n$ and $\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t_n\right)$is the temperature distribution at the previous time instant $t_n$. If $t_n=0$, Eqs. (31) and (32) is used as the initial conditions.

When multiplying both sides of Eq. (37) by ${\tau }_q$ and then substituting ${\tau }_q={\tau }_T=0$ into Eq. (37), the DPL model-based bio-heat transfer equation becomes the Fourier model-based bio-heat transfer equation, and its solution is expressed by Eq. (40).

$$\begin{array}{l}\theta \left(x,y,t_{n+1}\right)=\sum _{p=1}^{\infty }H{\theta }_{\infty }\frac{\left[1-{\left(-1\right)}^p\right]}{{\nu }_p}\frac{\mathit{sinh}{b}_1\left(p\right)y}{b_1\left(p\right)\mathit{cosh}{b}_1\left(p\right)W+\mathit{\text{Hsinhh}}{b}_1\left(p\right)W}\frac{\mathit{sin}{\nu }_{p}x}{N\left({\nu }_p\right)}\\ +\sum _{p=1}^{\infty }\sum _{m=1}^{\infty }\left[\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t_n\right)\right.\\ \left.-\left(\frac{\mathit{kH}{\theta }_{\infty }sin{\eta }_{m}W}{\rho c_{v}b}\right)\left[\frac{1-{\left(-1\right)}^p}{{\nu }_p}\right]\right]\mathit{exp}\left(-b∆t\right)\frac{\mathit{sin}{\nu }_{p}x}{N\left({\nu }_p\right)}\frac{\mathit{sin}{\eta }_{m}y}{N\left({\eta }_m\right)}\end{array}$$E40

Following the approach discussed in Section 2, the temperature rise due to a single pulse ($∆T$) was first obtained by solving Eq. (9). After that, this temperature rise has been added to the temperature distribution already obtained for the previous time instant to determine the temperature distribution at the current time instant:

$$\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t_n\right)=\stackrel{\sim }{\overline{\theta }}\left({\nu }_p,{\eta }_m,t_n\right)+∆\theta \left({\nu }_p,{\eta }_m\right)$$E41

where

$$∆\theta \left({\nu }_p,{\eta }_m\right)={\int }_{x^{\prime }=0}^L{\int }_{y^{\prime }=0}^W∆T\left(x^{\prime },y^{\prime }\right)sin{\nu }_p{x}^{\prime }sin{\eta }_m{y}^{\prime }{\mathit{dx}}^{\prime }{\mathit{dy}}^{\prime }$$E42

The trapezoidal rule has been employed for evaluating the double integral appearing on the right-hand side of Eq. (42).

Similarly, the analytical solution of the DPL model-based bio-heat transfer equation in the cylindrical coordinate system Eq. (26) can be obtained using FIT, and the details can be given in elsewhere [23].

5. Results and discussion

The results obtained using the solutions presented in Sections 3 and 4 have been discussed in this section. First, the spatial variation of the two components of the total intensity (collimate and diffused) was calculated at the location of the point of laser irradiation (Location 1 having coordinates of x = 0.98 mm and y = 2 mm) and the center of the domain (Location 2 having coordinates of x = 0.98 mm and y = 0.98 mm) to understand the laser-tissue interaction. After that, the temperature rise due to single-pulse and train of pulses was calculated at the same two locations. Then, the comparative study between the Fourier and non-Fourier (C-V and DPL) heat conduction model were carried out. Furthermore, numerical results were compared with the analytical results obtained using the FIT.

When light propagates through the biological tissue, the intensity has two components: collimated and diffused. The collimated component of the intensity represents the variation of the ballistic photons, which don't undergo any multiple scattering events, and its spatial variation at different time instants is shown in Figure 2a. The figure shows that the collimated intensity is maximum at y = 2 mm because it is the point of laser irradiation (Location 1). As expected, the collimated intensity decreases exponentially following the Beer-Lambert law. The collimated intensity is only available when time is less than or equal to the pulse width; otherwise, the value of the collimated intensity is zero.

Figure 2.

Spatial variation of collimated (a) and diffuse (b) components of light along a section passing through x = 0.98 mm [3].

Variation of the diffuse component of light intensity with respect to the depth (y) of the tissue phantom at various time instants is shown in Figure 2b. At the initial time instant (t = 4.6667 ps), the magnitude of diffuse light is significantly much smaller than the collimated intensity. Furthermore, in contrast to the profile seen for the collimated component wherein the maxima were observed at the top surface of the tissue phantom (the point of laser irradiation), the maxima of the diffuse component are observed at points below the top surface inside the body of the tissue phantom. This trend can be attributed to the scattering properties of the tissue domain, which affect the propagation of the diffuse component of the light intensity.

Figure 3a shows the temporal variation of the temperature at the same two locations considered in the biological tissue subjected to single pulse laser irradiation. As expected, the change in temperature at Location 1 (the point of laser irradiation) is more significant in comparison to Location 2. It is also to be seen from the figure that the rate of temperature rise is quite rapid during the first 5 ps, which is almost equal to the width of the laser pulse employed for irradiation since the laser power is available for this complete duration. Furthermore, the temperature rise achieves a constant value after nearly 20 ps. Because of this observation, the local temperature profile at t = 180 ps was used as the initial temperature field for solving the bio-heat transfer equation (Eq. 10) for determining the temperature distribution within biological tissue subjected to a train of laser pulses.

Figure 3.

(a) Temporal profile of temperature rise in tissue subjected to single pulse laser irradiation at x = 0.98 mm; (b) Temporal profile of temperature subjected to a train of pulse at x = 0.98 mm [3].

The Pennes bio-heat transfer equation was solved to determine the temperature distribution inside the biological tissue subjected to the train of pulses (repetition rate: 1 kHz). Figure 3b shows the temporal variation of temperature at the same two locations. The figure shows that the temperatures at these locations increase during the first one second, and thereafter a decay in the temperature values can be seen. It is expected since the total duration of laser irradiation is only one second. The temperature rise is comparatively much higher at Location 1, where the laser pulse strikes the sample, than in the rest of the region within the biological tissue. The temperature drop after one second is also steeper at the top wall of the tissue phantom as it is directly exposed to the ambient conditions.

The efficacy of laser-based photo-thermal therapy for selective destruction of cancerous cells depends on accurately predicting the temperature distribution inside the laser-irradiated biological tissue. Because of this, various heat conduction models, e.g., Fourier, C-V, DPL, etc., have been developed by researchers in the past. The need for the non-Fourier heat conduction model is realized because the photons travel with a finite speed through the biological tissues, contrary to the infinite speed of thermal wave propagation assumed in the Fourier model. So, we compared the temperature distribution obtained using the Fourier and non-Fourier heat conduction models to study the thermal response of laser-irradiated biological tissue phantom. Figure 4 shows the two-dimensional temperature distribution inside the laser-irradiated biological tissue at different time instants for Fourier and non-Fourier (C-V and DPL) models. While a nearly uniform diffusion of heat is to be seen in the contour plots corresponding to the Fourier model, a distinct wave nature in the thermal profiles is associated with the predictions of non-Fourier (C-V and DPL) heat conduction models. The oscillations in temperature values penetrate within the body of the tissue phantom as time progresses. The C-V heat conduction model predicts the maximum temperature values at any given time instant, followed by the DPL-based model. The C-V model also results in maximum amplitude of oscillations seen in the thermal wavefront compared to the DPL model.

Figure 4.

Two-dimensional temperature distribution at different time instants (i) Fourier (ii) C-V (iii) DPL model [21].

The numerical results obtained using the FVM were compared with analytical results obtained using the FIT. Figure 5 shows the analytical and numerical results for different heat conduction models (Fourier and non-Fourier). The temporal profiles of temperature without markers represent the results based on the FIT-based analytical solution, and those with markers represent the numerical results. The figure shows that the results obtained using numerical simulations predict relatively lower values of temperatures at Locations 1 and 2 in comparison with those obtained based on the analytical approach. It is also to be seen from the figure that the Fourier heat conduction model predicts relatively lower temperature values than those calculated using non-Fourier conduction models, which is to be attributed to the infinite speed of propagation of thermal waves considered in the conventional Fourier model.

Figure 5.

Comparison of temporal profiles of temperature predicted using Fourier and non-Fourier heat conduction models at Location 1 (a) and 2 (b) (Without marker: Analytical results; With marker: numerical results) [22].

Furthermore, the C-V heat conduction model predicts the highest temperature values (shown in the inset of Figure 5a and b for better clarity), while the DPL model predicts the temperatures that lie in between Fourier and C-V heat conduction models. The observed trend may be explained based on the C-V heat conduction model only considers the effect of the thermal relaxation time associated with heat flux (τ_q). In contrast, the dual-phase lag model considers the coupled effects of non-zero values of relaxation times related to the temperature gradient (τ_T) and heat flux (τ_q). In physical terms, the phase lag associated with the temperature gradients, i.e., τ_T tends to suppress the amplitude of the thermal wavefront. Because the value of τ_T is zero in the C-V heat conduction model, the resultant temperature profiles are relatively free of such dampening effects. Hence, the absolute values of temperatures predicted based on this form of non-Fourier heat conduction model are expected to be higher than predicted using the DPL model, wherein these two relaxation times are considered non-zero. The profiles shown in Figure 5 support this observation wherein the magnitude of temperatures as predicted using the DPL model lies between those obtained using the C-V heat conduction model (maximum values) and the Fourier-based heat conduction model (minimum values).

Figure 5b shows the temporal variation of temperature at Location 2. The profiles corresponding to the Fourier heat conduction model show a sudden drop in temperature values immediately after the laser power is switched off (t > 1.0 s). On the other hand, because of the phase lag terms associated with temperature gradients and/or heat flux, the temperature profiles predicted using the non-Fourier heat conduction models (hyperbolic and DPL) show nearly constant values of temperature for a relatively long period of time (t$\le$ 11 s) before the drop in temperature values starts. It indicates the prolonged (sustained) effects of thermal energy deposited at a given spatial location within the body of the tissue phantom, according to non-Fourier heat conduction models.

6. Conclusion

The present chapter discusses the modeling of the laser-irradiated biological tissue to understand its thermal behavior, which may help improve the efficacy of the laser-based photo-thermal therapy to destroy the cancerous cells with minimal damage to the surrounding healthy tissue. The light propagation through the biological tissue was mathematically modeled using the RTE. The RTE was solved using the DOM to determine the intensity distribution inside the biological tissue subjected to short-pulse laser irradiation. Once the intensity distribution was obtained, the divergence of radiative heat flux was calculated, which acts as the source term in the Fourier/non-Fourier model-based bio-heat transfer equation to determine the temperature distribution inside the laser-irradiated biological tissue. However, it is a multi-time scale problem because of the two different time scales used to solve the RTE and bio-heat transfer equation. So, the algorithm for solving this multi-time scale was presented. The non-Fourier model-based bio-heat transfer equation was numerically solved using the FVM to determine the temperature distribution. The numerical results were compared with the analytical results obtained using the FIT technique and found that the numerical solution predicted relatively lower temperature values than the analytical solution. A comparative study between the Fourier and non-Fourier (C-V and DPL) models was conducted and found that the temperature predicted using the DPL model lies between the C-V and Fourier models.

Acknowledgments

The author sincerely acknowledges the technical input provided by Prof. Atul Srivastava, IIT Bombay, India, during the preparation of this book chapter. The constant support and encouragement received from him over the last several years are gratefully acknowledged. The author also acknowledges the continuous support received from IIT Bombay and NIT Rourkela during the preparation of this work.