Modeling hemoglobin at optical frequency using the unconditionally stable fundamental ADI-FDTD method

doi:10.1364/BOE.2.001169

« Show journal navigation

Modeling hemoglobin at optical frequency using the unconditionally stable fundamental ADI-FDTD method

Ding Yu Heh and Eng Leong Tan »View Author Affiliations

Biomedical Optics Express, Vol. 2, Issue 5, pp. 1169-1183 (2011)
http://dx.doi.org/10.1364/BOE.2.001169

View Full Text Article

Acrobat PDF (1497 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Complex relative permittivity model for hemoglobin
3. Implementation using fundamental ADI-FDTD (FADI-FDTD) method
4. Numerical results
5. Conclusion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (29)
Cited By
Figures (7)
Metrics

Abstract

This paper presents the modeling of hemoglobin at optical frequency (250 nm – 1000 nm) using the unconditionally stable fundamental alternating-direction-implicit finite-difference time-domain (FADI-FDTD) method. An accurate model based on complex conjugate pole-residue pairs is proposed to model the complex permittivity of hemoglobin at optical frequency. Two hemoglobin concentrations at 15 g/dL and 33 g/dL are considered. The model is then incorporated into the FADI-FDTD method for solving electromagnetic problems involving interaction of light with hemoglobin. The computation of transmission and reflection coefficients of a half space hemoglobin medium using the FADI-FDTD validates the accuracy of our model and method. The specific absorption rate (SAR) distribution of human capillary at optical frequency is also shown. While maintaining accuracy, the unconditionally stable FADI-FDTD method exhibits high efficiency in modeling hemoglobin.

1. Introduction

The information on optical properties of human blood is crucial for many medical applications. In particular, the hemoglobin in red blood cells plays an important role by transporting oxygen throughout the human body. Past literatures have successfully determined the molar extinction coefficient [1

W. G. Zijlstra, A. Buursma, and W. P. Meeuwsen-van der Roest, “Absorption spectra of human fetal and adult oxyhemoglobin, deoxyhemoglobin, carboxyhemoglobin, and methemoglobin,” Clin. Chem. 37, 1633–1668 (1991).

–7

J. G. Kim and H. Liu, “Variation of haemoglobin extinction coefficients can cause errors in the determination of haemoglobin concentration measured by near-infrared spectroscopy,” Phys. Med. Biol. 52, 6295–6332 (2007). [CrossRef] [PubMed]

] of hemoglobin at optical frequency. The real part of its complex refractive index has also been found using transmittance and reflectance measurements [8

M. Meinke and M. Friebel, “Complex refractive index of hemoglobin in the wavelength range from 250 to 1100 nm,” Proc. SPIE 5862, 1–7 (2005).

–10

M. Meinke and M. Friebel, “Model function to calculate the refractive index of native hemoglobin in the wavelength range of 250 nm to 1100 nm dependent on concentration,” Appl. Opt. 45, 2838–2842 (2006). [CrossRef] [PubMed]

]. With such knowledge on the complex refractive index, a better understanding of the interaction of light with the biological media is greatly desired. Owing to this, the finite-difference time-domain (FDTD) method [11

K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]

, 12

A. Taflove and S. C. Hagness, Computational Electrodynamics: the Finite-Difference Time-Domain Method , (Artech House, 2005).

] has been widely used for solving electromagnetic problems and modeling electromagnetic wave propagation in various media. Unfortunately, the usage of such explicit FDTD method is limited by the fact that its time step is restricted by the Courant-Friedrich-Lewy (CFL) condition, given below as

Δ t \leq \sqrt{\frac{μ ɛ}{\frac{1}{Δ x^{2}} + \frac{1}{Δ y^{2}} + \frac{1}{Δ z^{2}}}}

(1)

where ε and μ are the permittivity and permeability of the media, respectively. Δt, Δx, Δy and Δz are the time step and spatial steps in x, y and z-directions, respectively. The above equation which sets the numerical stability bound for classical FDTD was first derived and published in [13

A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 23, 623–630 (1975). [CrossRef]

]. Note that if the spatial steps are small, this restriction may result in undesirably long simulation time.

The conventional alternating-direction-implicit finite-difference time-domain (ADI-FDTD) method [14

F. Zheng, Z. Chen, and J. Zhang, “Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method,” IEEE Trans. Microw. Theory Tech. 48, 1550–1558 (2000). [CrossRef]

, 15

T. Namiki, “3-D ADI-FDTD method-Unconditionally stable time domain algorithm for solving full vector Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 48, 1743–1748 (2000). [CrossRef]

], on the other hand, has unconditional stability feature where the time step is no longer restricted by the CFL stability criterion. This has circumvented the above limitation in its own right. However, due to the fact that the algorithm is implicit in nature, its implementation is complicated. Not only are there tridiagonal systems that need to be solved, the right-hand-side (RHS) of each equations involve many update coefficients and field variables, which calls for considerable arithmetic operations and memory indexing. Recently, these overheads has been overcome by the introduction of the fundamental alternating-direction-implicit finite-difference time-domain (FADI-FDTD) method [16

E. L. Tan, “Efficient algorithm for the unconditionally stable 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 7–9 (2007). [CrossRef]

, 17

E. L. Tan, “Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods,” IEEE Trans. Antennas Propag. 56, 170–177 (2008). [CrossRef]

]. In FADI-FDTD, its algorithm is included within a family of fundamental implicit schemes, which feature similar fundamental updating structures that are in simplest forms with most efficient matrix-operator-free RHS. This leads to substantial reduction in arithmetic operations and memory indexing. Nevertheless, while the evolution of computational electromagnetic in time domain from explicit FDTD method to FADI-FDTD has been encouraging, there is still a lack of appropriate complex permittivity models for biological media (in this case, hemoglobin) at optical frequency. Even then, the implementation of such appropriate complex permittivity models into the FADI-FDTD is also highly sought after.

In this paper, we seek to address the aforementioned problems and bridge the necessary gaps. First, we propose a model based on complex conjugate pole-residue pairs to describe the electrical behaviour of hemoglobin, where it is fitted against the complex refractive index data. Once an accurate model is found, it is incorporated into FADI-FDTD method for solving electromagnetic problems involving light interactions with biological media. Finally, some numerical examples are provided to validate the proposed method and to demonstrate some electromagnetic wave interactions with human blood at optical frequency.

2. Complex relative permittivity model for hemoglobin

2.1. Complex refractive index data

We first denote the wavelength dependent complex refractive index of hemoglobin as

n (λ) = n_{r} (λ) - j k (λ)

(2)

where

j = \sqrt{- 1}

, n_r is the real part of the refractive index and k is known as the extinction coefficient which is related to the absorption coefficient, μ_a of the medium by

k = \frac{μ_{a} λ}{4 π}

(3)

with λ being the free space wavelength. Note that the convention of time dependence e^jωt is used throughout the paper. k is always positive to yield negative imaginary part in Eq. (2) for losses or absorption in the media.

We make use of the tabulated molar extinction coefficient, e, data in [5

S. A. Prahl, “Tabulated molar extinction coefficient for hemoglobin in water,” http://omlc.ogi.edu/spectra/hemoglobin/summary.html (1998).

]. The absorption coefficient of various hemoglobin concentrations can be obtained from the molar extinction coefficient by means of the following relation:

μ_{a} = ln (10) \times e \times \frac{C}{64500}

(4)

where C is the hemoglobin concentration in g/L. Note that e has the unit of L mol⁻¹ cm⁻¹ or M⁻¹cm⁻¹ and 64500 is the gram molecular weight of hemoglobin in g/mol. Once the absorption coefficient is obtained, the imaginary part of the complex refractive index can be extracted through Eq. (3). On the other hand, the real part of the complex refractive index of hemoglobin at optical frequency has been determined and it is found to be related to the real part of the complex refractive index of water through a model function proposed in [10

] as

n_{r} = n_{w r} (β c + 1)

(5)

where n _wr is the real part of complex refractive index of water and c is the hemoglobin concentration in g/dL (1 dL = 0.1 L). β is the specific refractive increment which is also wavelength dependent and tabulated in [10

Therefore, using Eq. (4) and Eq. (5), the complex refractive index of hemoglobin at optical frequency can be extracted for various concentrations conveniently. Figure 1 plots the (a) real part of complex refractive index, n_r and (b) extinction coefficient, k, respectively, versus wavelength, λ for various hemoglobin concentrations. It can be seen that n_r increases with higher hemoglobin concentration across all wavelengths. The same also holds true for k, which indicates that higher hemoglobin concentration results in higher absorption (more loss). Another point to note is that when the hemoglobin concentration is at 0 g/dL, its n_r corresponds to that of water [18

G. M. Hale and M. R. Querry, “Optical constants of water in the 200 nm to 200 um wavelength region,” Appl. Opt. 12, 555–563 (1973). [CrossRef] [PubMed]

]. At this 0 g/dL concentration the k value is also reduced to zero. This is based on the fact that water has a “transparency window” at the visible light and near infrared region, in which it has negligible absorption or losses.

Fig. 1 (a) Real part of complex refractive index, n_r and (b) extinction coefficient, k versus wavelength, λ for different hemoglobin concentrations. Both n_r and k increase with higher hemoglobin concentration across all wavelengths.

2.2. Complex conjugate pole-residue pair model

In electromagnetic analysis, we are more interested in the complex relative permittivity rather than the refractive index. The complex relative permittivity is related to the refractive index by

\begin{array}{l} ɛ (ω) & = & n {(ω)}^{2} \\ = & n_{r} {(ω)}^{2} - k {(ω)}^{2} - j 2 n_{r} (ω) k (ω) \end{array}

(6)

where the wavelength dependence has been changed into angular frequency, ω = 2πf dependence. The complex relative permittivity values at various frequencies (wavelengths) can then be extracted using Eq. (6) from the values of n_r and k earlier. We now propose a complex relative permittivity model for hemoglobin. It can be described based on complex conjugate pole-residue pairs as follows:

ɛ (ω) = ɛ_{\infty} + \sum_{p} (\frac{r_{p}}{j ω - a_{p}} + \frac{r_{p}^{*}}{j ω - a_{p}^{*}})

(7)

where ε _∞ is the relative permittivity at infinite frequency, a_p and r_p are the p-th pole and residue, respectively.

The model based on complex conjugate pole-residue pairs in Eq. (7) is then fitted against the complex relative permittivity data extracted from Eq. (6) using the methodology detailed in [19

B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Trans. Power Delivery 14, 1052–1061 (1999). [CrossRef]

] to obtain the values of ε _∞, a_p and r_p . The values are fitted within the frequency range of 300 THz to 1200 THz (250 nm–1000 nm wavelength). In practical medical applications we are often interested in only two hemoglobin concentrations. The first one is the hemoglobin concentration in erythrocytes (red blood cells), which is also commonly known as mean corpuscular hemoglobin concentration (MCHC). Typical MCHC value for a healthy adult is around 33 g/dL. The second one is the overall hemoglobin concentration in human blood, which consists of not only erythrocytes, but also leukocytes (white blood cells), thrombocytes (platelets) and plasma. The overall hemoglobin concentration for a healthy adult is around 15 g/dL, and this corresponds to hematocrit (which is the proportion of blood volume that is occupied by red blood cells) of around 45 percent. It is found that using 6 pairs (order 12) of complex conjugate pole and residue pairs can sufficiently fit the experimental data. The values of ε _∞, a_p and r_p for hemoglobin concentrations of 15 g/dL and 33 g/dL are provided in Tables 1 and 2, respectively. Since all hemoglobin concentrations have rather the same trends of complex relative permittivity or complex refractive index, we assume that they all share the same poles a_p . That is, once the poles are identified for one hemoglobin concentration, they can be reused as poles of other hemoglobin concentrations. Only the new values of residues are calculated to avoid tedius searching of new poles.

Table 1 Fitted Values of a_p , r_p and ε _∞ for Hemoglobin Concentration of 15 g/dL

p	a_p	r_p
1	−3.1751 × 10¹⁴ + j2.4717 × 10¹⁵	−2.0191 × 10¹⁰ + j1.6882 × 10¹⁰
2	−7.1750 × 10¹³ + j3.2611 × 10¹⁵	1.1223 × 10¹¹ − j2.0688 × 10¹¹
3	−3.6183 × 10¹³ + j3.4679 × 10¹⁵	3.6477 × 10¹⁰ −j9.6911 × 10¹⁰
4	−1.4538 × 10¹⁴ + j4.5108 × 10¹⁵	7.1508 × 10¹¹ − j3.6455 × 10¹²
5	−1.5535 × 10¹⁴ + j5.3385 × 10¹⁵	2.0872 × 10¹¹ − j2.2396 × 10¹¹
6	−3.2153 × 10¹³ + j9.5959 × 10¹⁵	−4.2183 × 10¹² − j3.9864 × 10¹⁴

ε _∞ = 1.7844

Table 2 Fitted Values of a_p , r_p and ε _∞ for Hemoglobin Concentration of 33 g/dL

p	a_p	r_p
1	−3.1751 × 10¹⁴ + j2.4717 × 10¹⁵	−3.8401 × 10¹⁰ + j5.8207 × 10¹⁰
2	−7.1750 × 10¹³ + j3.2611 × 10¹⁵	2.8608 × 10¹¹ − j5.0495 × 10¹¹
3	−3.6183 × 10¹³ + j3.4679 × 10¹⁵	1.0080 × 10¹¹ − j2.2246 × 10¹¹
4	−1.4538 × 10¹⁴ + j4.5108 × 10¹⁵	1.6344 × 10¹² − j8.0730 × 10¹²
5	−1.5535 × 10¹⁴ + j5.3385 × 10¹⁵	3.6171 × 10¹¹ − j5.2711 × 10¹¹
6	–3.2153 × 10¹³ + j9.5959 × 10¹⁵	−1.4155 × 10¹³ − j4.5947 × 10¹⁴

ε _∞ = 1.8994

Using the values provided in Tables 1 and 2, we plot the real and imaginary part of the fitted model based on complex conjugate pole-residue pairs and compare them with the experimental data obtained in the previous subsection. Figure 2 plots the (a) real and (b) imaginary parts of the complex relative permittivity versus wavelength for both our model and experimental data. The absorption coefficient is then plotted in Fig. 3. It is evident that our model in Eq. (7), along with the values given in Tables 1 and 2, provide a good fit to the experimental data for hemoglobin concentrations of 15 g/dL and 33 g/dL. The proposed model is thus expected to describe the electromagnetic behaviour of hemoglobin accurately in the wavelength regions of 250 nm to 1000 nm.

Fig. 2 (a) Real part, ε_r and (b) imaginary part, ε_i of complex relative permittivity versus wavelength, λ : model and experimental data for hemoglobin concentrations of 15 g/dL and 33 g/dL. A good fit is observed between our proposed model and experimental data across all wavelengths.

Fig. 3 Absorption coefficient, μ_a versus wavelength, λ : model and experimental data for hemoglobin concentrations of 15 g/dL and 33 g/dL.

3. Implementation using fundamental ADI-FDTD (FADI-FDTD) method

To incorporate the hemoglobin model into FADI-FDTD update equations, we first write down the Maxwell’s curl equations as

\begin{array}{l} \nabla \times H (t) & = & ɛ_{0} ɛ_{\infty} \frac{\partial}{\partial t} E (t) + σ E (t) + ɛ_{0} \sum_{p} ℱ^{- 1} {j ω (\frac{r_{p}}{j ω - a_{p}} + \frac{r_{p}^{*}}{j ω - a_{p}^{*}}) E (ω)} \\ = & ɛ_{0} ɛ_{\infty} \frac{\partial}{\partial t} E (t) + σ E (t) + \sum_{p} (J_{p} (t) + {\hat{J}}_{p} (t)) \end{array}

(8a)

\begin{array}{l} - \nabla \times E (t) & = & μ_{0} μ_{\infty} \frac{\partial}{\partial t} H (t) + σ * H (t) + μ_{0} \sum_{p} ℱ^{- 1} {j ω (\frac{q_{p}}{j ω - b_{p}} + \frac{q_{p}^{*}}{j ω - b_{p}^{*}}) H (ω)} \\ = & μ_{0} μ_{\infty} \frac{\partial}{\partial t} H (t) + σ * H (t) + \sum_{p} (M_{p} (t) + {\hat{M}}_{p} (t)) \end{array}

(8b)

where ε ₀ and μ ₀ are the permittivity and permeability in free space and σ and σ^* are the electric and magnetic conductivities, respectively. For completeness, we have also included the implementation of magnetic complex conjugate pole-residue pairs model with permeability at infinite frequency μ _∞, p-th pole b_p and residue q_p . The symbol ℱ⁻¹ denotes inverse Fourier transform operation and two sets of auxiliary variables J _p and Ĵ _p , M _p and M̂ _p have been introduced. These two sets of variables are related to electric field E and magnetic field H by the following equations:

\frac{\partial}{\partial t} J_{p} (t) - a_{p} J_{p} (t) = ɛ_{0} r_{p} \frac{\partial}{\partial t} E (t)

(9a)

\frac{\partial}{\partial t} {\hat{J}}_{p} (t) - a_{p}^{*} {\hat{J}}_{p} (t) = ɛ_{0} r_{p}^{*} \frac{\partial}{\partial t} E (t)

(9b)

\frac{\partial}{\partial t} M_{p} (t) - b_{p} M_{p} (t) = μ_{0} q_{p} \frac{\partial}{\partial t} H (t)

(10a)

\frac{\partial}{\partial t} {\hat{M}}_{p} (t) - b_{p}^{*} {\hat{M}}_{p} (t) = μ_{0} q_{p}^{*} \frac{\partial}{\partial t} H (t) .

(10b)

It is worth noting that since E(t) and H(t) are always of real value, Ĵ(t) = J ^* (t) and M̂(t) = M ^* (t) always hold and Eq. (8) can be simplified to

\nabla \times H (t) = ɛ_{0} ɛ_{\infty} \frac{\partial}{\partial t} E (t) + σ E (t) + \sum_{p} 2 Re {J_{p} (t)}

(11a)

- \nabla \times E (t) = μ_{0} μ_{\infty} \frac{\partial}{\partial t} H (t) + σ * H (t) + \sum_{p} 2 Re {M_{p} (t)} .

(11b)

Therefore, only variables J and M are needed to be solved. We then discretize Eq. (9a) and Eq. (10a) in time using central difference approximation (at time index n = 1/4) to yield

J_{p}^{n + \frac{1}{2}} = \frac{1 + a_{p} \frac{Δ t}{4}}{1 - a_{p} \frac{Δ t}{4}} J_{p}^{n} + \frac{r_{p}}{1 - a_{p} \frac{Δ t}{4}} (E^{n + \frac{1}{2}} - E^{n}) = k_{1 p} J_{p}^{n} + k_{2 p} (E^{n + \frac{1}{2}} - E^{n})

(12a)

M_{p}^{n + \frac{1}{2}} = \frac{1 + b_{p} \frac{Δ t}{4}}{1 - b_{p} \frac{Δ t}{4}} M_{p}^{n} + \frac{q_{p}}{1 - b_{p} \frac{Δ t}{4}} (H^{n + \frac{1}{2}} - H^{n}) = l_{1 p} M_{p}^{n} + l_{2 p} (H^{n + \frac{1}{2}} - H^{n}) .

(12b)

The time step used in the central difference approximation is Δt/2 because FADI-FDTD method has two substeps with each substep having time step of Δt/2. Note that apart from the central difference approximation, the corrected impulse invariance method in Z-transform theory [20

D. Y. Heh and E. L. Tan, “Corrected impulse invariance method in z-transform theory for frequency-dependent FDTD methods,” IEEE Trans. Antennas Propag. 57, 2683–2690 (2009). [CrossRef]

] can also be utilized. Along with Eq. (12) , we first discretize Eq. (11) into ADI-FDTD of conventional form. The E_x and H_z components of first substep is written as

E_{x}^{n + \frac{1}{2}} - c_{2} \frac{\partial}{\partial y} H_{z}^{n + \frac{1}{2}} = c_{1} E_{x}^{n} - c_{2} \frac{\partial}{\partial z} H_{y}^{n} - c_{2} \sum_{p} Re {(1 + k_{1 p}) J_{x p}^{n}}

(13a)

H_{z}^{n + \frac{1}{2}} - d_{2} \frac{\partial}{\partial y} E_{x}^{n + \frac{1}{2}} = d_{1} H_{z}^{n} - d_{2} \frac{\partial}{\partial x} E_{y}^{n} - d_{2} \sum_{p} Re {(1 + l_{1 p}) M_{x p}^{n}}

(13b)

where

c_{1} = \frac{1 - \frac{σ Δ t}{4 ɛ_{0} ɛ_{\infty}} + \frac{Δ t}{4 ɛ_{0} ɛ_{\infty}} \sum_{p} 2 Re (k_{2 p})}{1 + \frac{σ Δ t}{4 ɛ_{0} ɛ_{\infty}} + \frac{Δ t}{4 ɛ_{0} ɛ_{\infty}} \sum_{p} 2 Re (k_{2 p})}, c_{2} = \frac{\frac{Δ t}{2 ɛ_{0} ɛ_{\infty}}}{1 + \frac{σ Δ t}{4 ɛ_{0} ɛ_{\infty}} + \frac{Δ t}{4 ɛ_{0} ɛ_{\infty}} \sum_{p} 2 Re (k_{2 p})},

(14a)

d_{1} = \frac{1 - \frac{σ * Δ t}{4 μ_{0} μ_{\infty}} + \frac{Δ t}{4 μ_{0} μ_{\infty}} \sum_{p} 2 Re (l_{2 p})}{1 + \frac{σ * Δ t}{4 μ_{0} μ_{\infty}} + \frac{Δ t}{4 μ_{0} μ_{\infty}} \sum_{p} 2 Re (l_{2 p})}, d_{2} = \frac{\frac{Δ t}{2 μ_{0} μ_{\infty}}}{1 + \frac{σ * Δ t}{4 μ_{0} μ_{\infty}} + \frac{Δ t}{4 μ_{0} μ_{\infty}} \sum_{p} 2 Re (l_{2 p})} .

(14b)

To solve for E_x , we need to substitute Eq. (13b) into Eq. (13a) to obtain the implicit electic field update equation as

\begin{array}{l} E_{x}^{n + \frac{1}{2}} - c_{2} \frac{\partial}{\partial y} (d_{2} \frac{\partial}{\partial y} E_{x}^{n + \frac{1}{2}}) & = & c_{1} E_{x}^{n} - c_{2} \frac{\partial}{\partial z} H_{y}^{n} - c_{2} \sum_{p} Re {(1 + k_{1 p}) J_{x p}^{n}} + c_{2} \frac{\partial}{\partial y} (d_{1} H_{z}^{n}) \\ - c_{2} \frac{\partial}{\partial y} (d_{2} \frac{\partial}{\partial x} E_{y}^{n}) - c_{2} \frac{\partial}{\partial y} (d_{2} \sum_{p} Re {(1 + l_{1 p}) M_{x p}^{n}}) . \end{array}

(15)

When the spatial differential operators are approximated using the second-order central differencing operated on a standard Yee’s cell [11

K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]

], the electric field components can be obtained by solving a set of tridiagonal matrices.

The full update equations of conventional ADI-FDTD can be summarized in compact matrix form (in two substeps) as

(I_{18} - D_{2} A + F_{l}) u^{n + \frac{1}{2}} = (D_{1} + D_{2} B + F_{r}) u^{n} + D_{2} s^{n + \frac{1}{2}}

(16a)

(I_{18} - D_{2} B + F_{l}) u^{n + 1} = (D_{1} + D_{2} A + F_{r}) u^{n + \frac{1}{2}} + D_{2} s^{n + \frac{1}{2}}

(16b)

where

u = {[E_{x}, E_{y}, E_{z}, H_{x}, H_{y}, H_{z}, J_{x}^{'}, J_{y}^{'}, J_{z}^{'}, J_{x}^{″}, J_{y}^{″}, J_{z}^{″}, M_{x}^{'}, M_{y}^{'}, M_{z}^{'}, M_{x}^{″}, M_{y}^{″}, M_{z}^{″}]}^{T}

(17)

s = {[- s_{e x}, - s_{e y}, - s_{e z}, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]}^{T}

(18)

A = [\begin{matrix} O_{3} & K_{1} & O_{3} & O_{3} & O_{3} & O_{3} \\ K_{2} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \end{matrix}], B = [\begin{matrix} O_{3} & - K_{2} & O_{3} & O_{3} & O_{3} & O_{3} \\ - K_{1} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \end{matrix}]

(19)

K_{1} = [\begin{array}{r} 0 & 0 & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} & 0 & 0 \\ 0 & \frac{\partial}{\partial x} & 0 \end{array}], K_{2} = [\begin{array}{r} 0 & \frac{\partial}{\partial z} & 0 \\ 0 & 0 & \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} & 0 & 0 \end{array}]

(20)

D_{1} = diag (c_{1}, c_{1}, c_{1}, d_{1}, d_{1}, d_{1}, k_{1}^{'}, k_{1}^{'}, k_{1}^{'}, k_{1}^{'}, k_{1}^{'}, k_{1}^{'}, l_{1}^{'}, l_{1}^{'}, l_{1}^{'}, l_{1}^{'}, l_{1}^{'}, l_{1}^{'})

(21)

D_{2} = diag (c_{2}, c_{2}, c_{2}, d_{2}, d_{2}, d_{2}, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

(22)

F_{l} = [\begin{matrix} O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ - k_{2}^{'} I_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ - k_{2}^{″} I_{3} & O_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & - l_{2}^{'} I_{3} & O_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & - l_{2}^{″} I_{3} & O_{3} & O_{3} & O_{3} & O_{3} \end{matrix}]

(23)

F_{r} = [\begin{matrix} O_{3} & O_{3} & - c_{2} {(1 + k_{1})}^{'} I_{3} & c_{2} k_{1}^{″} I_{3} & O_{3} & O_{3} \\ O_{3} & O_{3} & O_{3} & O_{3} & - d_{2} {(1 + l_{1})}^{'} I_{3} & d_{2} l_{1}^{″} I_{3} \\ - k_{2}^{'} I_{3} & O_{3} & O_{3} & - k_{1}^{″} I_{3} & O_{3} & O_{3} \\ - k_{2}^{″} I_{3} & O_{3} & k_{1}^{″} I_{3} & O_{3} & O_{3} & O_{3} \\ O_{3} & - l_{2}^{'} I_{3} & O_{3} & O_{3} & O_{3} & - l_{1}^{″} I_{3} \\ O_{3} & - l_{2}^{″} I_{3} & O_{3} & O_{3} & l_{1}^{″} I_{3} & O_{3} \end{matrix}] .

(24)

Note that variables χ ^′ = Re(χ) and χ ^″ = Im(χ). For simplicity, we have assumed single complex conjugate pole-residue pair for now. I _r and O _r represent identity and null matrices with dimension r × r. We have also included electric current sources s_ex , s_ey and s_ez in the formulation.

A quick glance at Eq. (16) shows that its right-hand-side (RHS) still involves matrix operators A and B which incur considerable arithmetic operations and memory indexing, especially when one is solving the implicit update equations. Efforts have been made to further reduce such overheads by means of a more efficient ADI-FDTD method in lossless media proposed in [16

E. L. Tan, “Efficient algorithm for the unconditionally stable 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 7–9 (2007). [CrossRef]

, 17

E. L. Tan, “Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods,” IEEE Trans. Antennas Propag. 56, 170–177 (2008). [CrossRef]

]. Such efficient ADI algorithm is included within a family of fundamental implicit schemes, which feature similar fundamental updating structures that are in simplest forms with most efficient matrix-operator-free RHS. Following [17

E. L. Tan, “Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods,” IEEE Trans. Antennas Propag. 56, 170–177 (2008). [CrossRef]

], it can be shown that Eq. (16) is reducible to the fundamental form (FADI-FDTD), which calls for the following update procedures:

{\tilde{v}}^{n} = [(I_{12} + D_{1}) + F_{l} + F_{r}] u^{n} - {\tilde{v}}^{n - \frac{1}{2}}

(25a)

(I_{12} - D_{2} A + F_{l}) u^{n + \frac{1}{2}} = {\tilde{v}}^{n} + D_{2} s^{n + \frac{1}{2}}

(25b)

{\tilde{v}}^{n + \frac{1}{2}} = [(I_{12} + D_{1}) + F_{l} + F_{r}] u^{n + \frac{1}{2}} - {\tilde{v}}^{n}

(25c)

(I_{12} - D_{2} B + F_{l}) u^{n + 1} = {\tilde{v}}^{n + \frac{1}{2}}

(25d)

where Eq. (25a) and Eq. (25b) form the first substep while Eq. (25c) and Eq. (25d) form the second. It is worth noting that the FADI-FDTD has matrix-operator-free RHS, i.e., there is no matrix operator A and B on the RHS of Eq. (25) . The source excitation is also only needed at the first substep [21

E. L. Tan, “Concise current source implementation for efficient 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 748–750 (2007). [CrossRef]

] which further simplifies the implementation compared to the conventional ADI-FDTD method.

We now provide the full update equations of FADI-FDTD from Eq. (25) . Letting

\tilde{v} = {[{\tilde{e}}_{x}, {\tilde{e}}_{y}, {\tilde{e}}_{z}, {\tilde{h}}_{x}, {\tilde{h}}_{y}, {\tilde{h}}_{z}, {\tilde{j}}_{x}^{'}, {\tilde{j}}_{y}^{'}, {\tilde{j}}_{z}^{'}, {\tilde{j}}_{x}^{″}, {\tilde{j}}_{y}^{″}, {\tilde{j}}_{z}^{″}, {\tilde{m}}_{x}^{'}, {\tilde{m}}_{y}^{'}, {\tilde{m}}_{z}^{'}, {\tilde{m}}_{x}^{″}, {\tilde{m}}_{y}^{″}, {\tilde{m}}_{z}^{″}]}^{T},

(26)

we have for the first substep,

{\tilde{e}}_{ξ}^{n} = (1 + c_{2}) E_{ξ}^{n} - {\tilde{e}}_{x}^{n - \frac{1}{2}} - c_{2} \sum_{p} Re {(1 + k_{1 p}) J_{ξ p}^{n}}, ξ = x, y, z

(27a)

{\tilde{h}}_{ξ}^{n} = (1 + d_{2}) H_{ξ}^{n} - {\tilde{h}}_{x}^{n - \frac{1}{2}} - d_{2} \sum_{p} Re {(1 + l_{1 p}) M_{ξ p}^{n}}, ξ = x, y, z

(27b)

{\tilde{j}}_{ξ p}^{n} = (1 + k_{1 p}) J_{ξ p}^{n} - {\tilde{j}}_{ξ p}^{n - \frac{1}{2}} - 2 k_{2 p} E_{ξ}^{n}, \forall p, ξ = x, y, z

(27c)

{\tilde{m}}_{ξ p}^{n} = (1 + l_{1 p}) M_{ξ p}^{n} - {\tilde{m}}_{ξ p}^{n - \frac{1}{2}} - 2 l_{2 p} H_{ξ}^{n}, \forall p, ξ = x, y, z

(27d)

E_{x}^{n + \frac{1}{2}} - c_{2} \frac{\partial}{\partial y} H_{z}^{n + \frac{1}{2}} = {\tilde{e}}_{x}^{n} - c_{2} s_{e x}

(27e)

E_{y}^{n + \frac{1}{2}} - c_{2} \frac{\partial}{\partial z} H_{x}^{n + \frac{1}{2}} = {\tilde{e}}_{y}^{n} - c_{2} s_{e y}

(27f)

E_{z}^{n + \frac{1}{2}} - c_{2} \frac{\partial}{\partial x} H_{y}^{n + \frac{1}{2}} = {\tilde{e}}_{z}^{n} - c_{2} s_{e z}

(27g)

H_{x}^{n + \frac{1}{2}} - d_{2} \frac{\partial}{\partial z} E_{y}^{n + \frac{1}{2}} = {\tilde{h}}_{x}^{n}

(27h)

H_{y}^{n + \frac{1}{2}} - d_{2} \frac{\partial}{\partial x} E_{z}^{n + \frac{1}{2}} = {\tilde{h}}_{y}^{n}

(27i)

H_{z}^{n + \frac{1}{2}} - d_{2} \frac{\partial}{\partial y} E_{x}^{n + \frac{1}{2}} = {\tilde{h}}_{z}^{n}

(27j)

J_{ξ p}^{n + \frac{1}{2}} = k_{2 p} E_{ξ}^{n + \frac{1}{2}} + {\tilde{j}}_{ξ p}^{n}, \forall p, ξ = x, y, z

(27k)

M_{ξ p}^{n + \frac{1}{2}} = l_{2 p} H_{ξ}^{n + \frac{1}{2}} + {\tilde{m}}_{ξ p}^{n}, \forall p, ξ = x, y, z .

(27l)

The second substep can be similarly written down by expanding Eq. (25c) and Eq. (25d). To solve for E_x implicitly, Eq. (27j) is substituted into Eq. (27e) to yield

E_{x}^{n + \frac{1}{2}} - c_{2} \frac{\partial}{\partial y} (d_{2} \frac{\partial}{\partial y} E_{x}^{n + \frac{1}{2}}) = {\tilde{e}}_{x}^{n} - c_{2} s_{e x} + c_{2} \frac{\partial}{\partial y} {\tilde{h}}_{x}^{n} .

(28)

The advantages of the FADI-FDTD over the conventional ADI-FDTD can be summarized as follows:

Matrix-operator-free RHS in FADI-FDTD, c.f. Eq. (25) . On the other hand, the RHS of conventional ADI-FDTD still involves matrix operators A and B, c.f. Eq. (16) .
Source excitation is only needed in first substep of FADI-FDTD, c.f. Eq. (25b), compared to both substeps of conventional ADI-FDTD, c.f. Eq. (16a) and Eq. (16b).
Implicit update equation of FADI-FDTD, c.f. Eq. (28) is much simpler compared to that of conventional ADI-FDTD, c.f. Eq. (15).
Implicit update equation of FADI-FDTD, c.f. Eq. (28) has much lesser spatial differential operators than that of conventional ADI-FDTD, c.f. Eq. (15). Thus, applying finite difference approximation to these spatial differential operators will incur less overheads in FADI-FDTD compared to conventional ADI-FDTD.
No dispersive terms J and M is required in implicit update equation of FADI-FDTD, c.f. Eq. (28). On the other hand, they are required in implicit update equation of conventional ADI-FDTD, c.f. Eq. (15).
The overall number of RHS terms in FADI-FDTD are lesser than that of conventional ADI-FDTD. This reduces the amount of arithmetic operations, memory indexing, and results in higher efficiency.

4. Numerical results

4.1. Transmission and reflection coefficients

To demonstrate the accuracy of our hemoglobin model in the FADI-FDTD method, we first compute the transmission and reflection coefficients of a half space human blood numerically. A modulated Gaussian pulse is launched from free space and upon hitting the human blood interface, undergoes transmission and reflection. The modulated Gaussian pulse has a bandwidth of 1500 THz around centre frequency of 750 THz, thus covering the entire frequency spectrum of the human blood model (300 THz to 1200 THz). The model of 15 g/dL hemoglobin concentration (Table 1) is used. Cell size is set at 3.125 nm which corresponds to cell per wavelength (CPW) of around 50 for the shortest wavelength in the human blood (wavelength in human blood, λ ≈ λ ₀/1.5, where λ ₀ is the wavelength in free space. The shortest wavelenght considered in free space is 250 nm). The time step Δt is always specified as Courant-Friedrich-Lewy number (CFLN), which is relative to the Courant-Friedrich-Lewy time step limit in the explicit Yee’s FDTD method, i.e. CFLN = Δt/Δt_CFL where Δt_CFL is the maximum allowed time step of explicit FDTD method in Eq. (1).

The transmission and reflection coefficents (in frequency domain) can be obtained numerically by taking the ratio of the discrete Fourier transform of the transmitted and reflected wave, respectively to the discrete Fourier transform of the incident wave. The numerical experiment is run for various CFLNs to validate the accuracy of the hemoglobin model under different time steps. Figure 4 graphs the magnitude of the computed (a) transmission and (b) reflection coefficients, respectively versus frequency for various CFLNs. They are compared against the reference transmission and reflection coefficients of our analytical hemoglobin model. Also shown in the figures are the transmission and reflection coefficients computed using the explicit FDTD method. Since the explicit FDTD method is not unconditionally stable, it is run at its maximum time step limit of CFLN=1. Any time step larger than that will render the algorithm unstable. The FADI-FDTD method is unconditionally stable, and therefore we can see that CFLN beyond 1 still yield considerably accurate results. Nevertheless, the accuracy of the FADI-FDTD method deteriorates as the CFLN increases, especially towards higher frequency region. This is because as the frequency increases, higher CFLN (and thus larger time step) fails to resolve these higher frequency components sufficiently. Therefore, larger chosen time step increases the simulation speed, but at the expense of sacrificing a certain degree of accuracy. Based on these results, we can see that the chosen time step in FADI-FDTD is always a compromise between efficiency and accuracy.

Fig. 4 Magnitude of (a) transmission and (b) reflection coefficients versus frequency for various CFLNs.

4.2. Specific absorption rate (SAR)

The specifc absorption rate (SAR) is an important measure for biological media exposed to electromagnetic radiation, which describes the electromagnetic energy absorption rate. It is defined as

SAR = \frac{s_{l}}{ρ}

(29)

where s_l is the power loss density in W/m³ and ρ is the density of the biological media in kg/m³, which gives the SAR the overall unit of W/kg. In most of the previous SAR analysis literatures [22

C. M. Furse and O. P. Gahdhi, “A memory efficient method of calculating specific absorption rate in CW FDTD simulations,” IEEE Trans. Biomed. Eng. 43, 558–560 (1996). [CrossRef] [PubMed]

–25

S. G. Garcez, C. F. Galan, L. H. Bonani, and V. Baranauskas, “Estimating the electromagnetic field effects in biological tissues through the finite-difference time-domain method,” SBMO/IEEE MTT-S Int. Microw. and Optoelectronics Conf. Proc. , 43, 58–62 (2007).

], the power loss density is usually computed using

s_{l} = \frac{1}{2} σ {| \tilde{E} |}^{2}

(30)

where Ẽ is the phasor form (without time dependence) of the electric field components. In time domain methods, such as the FDTD methods, the significant benefit is such that a narrow pulse can be injected to provide a wide range of frequency points using only a single run of simulation. In this case, Ẽ can be the Fourier transform of the recorded electric field components, which gives us the information across a range of frequency resulted from illumination by any sinusoidal (monochromatic) source. This however, can only be performed if the medium of interest is modeled as the simplest case of conductive media, described by only a single conductivity value across the frequency range of interest. This is usually not the case in most biological media, such as the human blood described in this paper, where the conductivity changes with frequency. To use Eq. (30), the simulation has to be carried out point by point across the frequency range of interest using sinusoidal excitation, along with the corresponding conductivity values at different frequency. The advantage of FDTD method is thus not being exploited fully.

Since we have successfully modeled human blood across a range of optical frequency based on complex conjugate pole-residue pairs along with its implementation using the FADI-FDTD method, the power loss density formula in Eq. (30) can be modified to better suit our scenario. For a dispersive medium which has a complex permittivity of ε = ε ^′ – jε ^″, where ε ^′ and ε ^″ are the real and negative imaginary (dissipative) parts, respectively, the power loss density can be found using [26

J. D. Jackson, Classical Electrodynamics , 3rd ed. (John Wiley & Sons, 1998).

, 27

D. M. Pozar, Microwave Engineering , 3rd ed. (Wiley, 2005).

]

s_{l} = \frac{1}{2} ω ɛ^{″} {| \tilde{E} |}^{2} .

(31)

Note that for conductive media with single conductivity value, ε ^″ = σ/ω and Eq. (31) conforms to Eq. (30), which is considered as a special case for simplest conductive media. Hence, Eq. (31) is the generalization of power loss density for dispersive media. Using Eq. (31), we are now able to obtain the information across a range of frequency with only a single run of simulation, which is highly efficient.

In the following numerical example, we consider a capillary (blood vessel) in our body, and has a typical diameter of 8 μm. The straight long capillary is oriented along z-direction and is surrounded by tissue, which is assumed to be made up of mostly water and thus, having similar dielectric properties of water. Since the water has negligible absorption in the visible light and near infrared region, the refractive index is taken as a constant value of 1.33. Both models of 15 g/dL and 33 g/dL hemoglobin concentrations are applied to the capillary. The cell size is chosen as 11.1 nm which corresponds to CPW of 15 at shortest wavelength in the human blood and the time step is chosen at CFLN=4. A z-directed line current source of the same modulated Gaussian pulse as the previous subsection illuminates upon the capillary, exciting the frequency contents from 300 THz to 1200 THz. The perfectly matched layer [28

W. C. Tay, D. Y. Heh, and E. L. Tan, “GPU-accelerated fundamental ADI-FDTD with complex frequency shifted convolutional perfectly matched layer,” PIER M 14, 177–192 (2010). [CrossRef]

] terminates the computational boundaries.

At steady state, the SAR values at different locations of the capillary are computed using Eq. (29) and Eq. (31) and the density, ρ of the human blood is 1060 kg/m³ [25

]. The SAR is computed at frequencies of 440 THz, 550 THz and 720 THz. These frequencies correspond to free space wavelengths of 680 nm, 545 mm and 417 nm. Refering to Fig. 3, wavelength at 417 nm is where the highest absorption peak occurs, followed by 545 nm. The lowest absorption occurs at 680 nm. Figs. 5, 6 and 7 illustrate the SAR values (in W/kg) at different locations of the capillary at 440 THz, 550 THz and 720 THz, respectively. All the SAR values are plotted in logarithmic scale (log₁₀ SAR). All figures (a) refer to hemoglobin concentration of 15 g/dL while figures (b) refer to hemoglobin concentration of 33 g/dL. The SAR values at each location have been scaled to correspond to incident power density

S = \frac{1}{2} Re (\tilde{E} \times \tilde{H} *)

of 1 mW/m² in free space. Note that the source illuminates from the left to right of the figures. We first notice that frequency at 720 THz (417 nm), c.f. Fig. 7 gives rise to the highest SAR values, followed by 550 THz (545 nm), c.f. Fig. 6 and 440 THz (680 nm), c.f. Fig. 5. This is in consistent with the absorption coefficient profile in Fig. 3. Another interesting observation is that at only 720 THz (417 nm), the SAR in the capillary shows a descending trend from the left to right along the wave travelling direction, while other frequencies show rather uniform SAR across the whole capillary region. This indicates that at the absorption peak of 720 THz (417 nm), the bulk of the power is absorbed near the left side, which results in lesser power penetration across the capillary. On the other hand, the absorbed power is well distributed across the whole region of the capillary at other frequencies due to weaker absorption. The standing wave effect may be attributed to the closed volume of the capillary. At the same time, it is observed that capillary with hemoglobin concentration of 33 g/dL generally exhibits higher SAR values across all frequencies. This concentration represents that of only the red blood cells, and could possibly model blood clot inside a capillary, which is formed as a result of entrapment of red blood cells [29

K. Boryczko, W. Dzewinel, and D. A. Yuen, “Modeling fibrin aggregation in blood flow with discrete particles,” Comput. Methods Prog. Biomed. 75, 181–194 (2004). [CrossRef]

]. Thus, the SAR calculation reveals that blood clot has a higher absorption compared to normal blood, which could be useful for medical applications and diagnostics.

Fig. 5 SAR (in logarithmic scale) at different locations of capillary at 440 THz (680 nm) for hemoglobin concentrations of (a) 15 g/dL and (b) 33 g/dL. The wave illuminates from left to right.

Fig. 6 SAR (in logarithmic scale) at different locations of capillary at 550 THz (545 nm) for hemoglobin concentrations of (a) 15 g/dL and (b) 33 g/dL. The wave illuminates from left to right.

Fig. 7 SAR (in logarithmic scale) at different locations of capillary at 720 THz (417 nm) for hemoglobin concentrations of (a) 15 g/dL and (b) 33 g/dL. The wave illuminates from left to right.

We have also performed similar numerical example using the explicit FDTD method. Tables 3 and 4 show the maximum and average SAR of the capillary computed using the FADI-FDTD method and explicit FDTD method at hemoglobin concentrations of 15 g/dL and 33 g/dL. It can be seen that the calculated SAR values are quite close to each other. This shows that using the FADI-FDTD method at CFLN=4 does not suffer too much from accuracy deficiency. It is also found that the simulation time incurred by the FADI-FDTD method at CFLN=4 in Matlab codes is around 3 times faster than that of the explicit FDTD method. These further show that while maintaining accuracy, the FADI-FDTD method is promising in achieving higher efficiency for modeling biological media.

Table 3 Maximum and Average SAR Values Computed using FADI-FDTD and Explicit FDTD at Hemoglobin Concentration of 15 g/dL

Frequency/Wavelength	FADI-FDTD		explicit FDTD
Frequency/Wavelength	SAR max (W/kg)	SAR avg (W/kg)	SAR max (W/kg)	SAR avg (W/kg)
440 THz (680 nm)	2.4870 × 10⁻⁴	1.5662 × 10⁻⁴	2.5134 × 10⁻⁴	1.5659 × 10⁻⁴
550 THz (545 nm)	0.0438	0.0243	0.0445	0.0247
720 THz (417 nm)	0.2783	0.0969	0.2758	0.0964

Table 4 Maximum and Average SAR Values Computed Using FADI-FDTD and Explicit FDTD at Hemoglobin Concentration of 33 g/dL

Frequency/Wavelength	FADI-FDTD		explicit FDTD
Frequency/Wavelength	SAR max (W/kg)	SAR avg (W/kg)	SAR max (W/kg)	SAR avg (W/kg)
440 THz (680 nm)	0.0012	4.5712 × 10⁻⁴	0.0012	4.5704 × 10⁻⁴
550 THz (545 nm)	0.1574	0.0491	0.1654	0.0510
720 THz (417 nm)	0.6021	0.0914	0.5973	0.0906

5. Conclusion

This paper has presented the modeling of hemoglobin at optical frequency (250 nm – 1000 nm) using the unconditionally stable fundamental alternating-direction-implicit finite-difference time-domain (FADI-FDTD) method. An accurate model based on complex conjugate pole-residue pairs has been proposed to model the complex permittivity of hemoglobin at optical frequency. Two hemoglobin concentrations at 15 g/dL and 33 g/dL have been considered. The model is then incorporated into the FADI-FDTD method for solving electromagnetic problems involving interaction of light with hemoglobin. The computation of transmission and reflection coefficients of a half space hemoglobin medium using the FADI-FDTD has validated the accuracy of our model and method. The specific absorption rate (SAR) distribution of human capillary at optical frequency has also been shown. While maintaining accuracy, the unconditionally stable FADI-FDTD method has exhibited high efficiency in modeling hemoglobin.

References and links

1.	W. G. Zijlstra, A. Buursma, and W. P. Meeuwsen-van der Roest, “Absorption spectra of human fetal and adult oxyhemoglobin, deoxyhemoglobin, carboxyhemoglobin, and methemoglobin,” Clin. Chem. 37, 1633–1668 (1991).
2.	W. G. Zijlstra, A. Buursma, H. E. Falke, and J. F. Catsburg, “Spectrophotometry of hemoglobin: absorption spectra of rat oxyhemoglobin, deoxyhemoglobin, carboxyhemoglobin, and methemoglobin,” Comput. Biochem. Physiol. 107B, 161–166 (1994). [CrossRef]
3.	W. G. Zijlstra, A. Buursma, and O. W. van Assendelft, Visible and Near Infrared Absorption Spectra of Human and Animal Haemoglobin: Determination and Application , (VSP, Zeist, 2000).
4.	M. Cope, “The application of near infrared spectroscopy to non invasive monitoring of cerebral oxygenation in the newborn infant,” PhD Dissertation, (1991).
5.	S. A. Prahl, “Tabulated molar extinction coefficient for hemoglobin in water,” http://omlc.ogi.edu/spectra/hemoglobin/summary.html (1998).
6.	J. G. Kim, M. Xia, and H. Liu, “Extinction coefficients of hemoglobin for near-infrared spectroscopy of tissue,” IEEE Eng. Med. Biol. Mag. 24, 118–121 (2005). [CrossRef] [PubMed]
7.	J. G. Kim and H. Liu, “Variation of haemoglobin extinction coefficients can cause errors in the determination of haemoglobin concentration measured by near-infrared spectroscopy,” Phys. Med. Biol. 52, 6295–6332 (2007). [CrossRef] [PubMed]
8.	M. Meinke and M. Friebel, “Complex refractive index of hemoglobin in the wavelength range from 250 to 1100 nm,” Proc. SPIE 5862, 1–7 (2005).
9.	M. Meinke and M. Friebel, “Determination of the complex refractive index of highly concentrated hemoglobin solutions using transmittance and reflectance measurements,” J. Biomed. Opt. 10, 064019 (2005). [CrossRef]
10.	M. Meinke and M. Friebel, “Model function to calculate the refractive index of native hemoglobin in the wavelength range of 250 nm to 1100 nm dependent on concentration,” Appl. Opt. 45, 2838–2842 (2006). [CrossRef] [PubMed]
11.	K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]
12.	A. Taflove and S. C. Hagness, Computational Electrodynamics: the Finite-Difference Time-Domain Method , (Artech House, 2005).
13.	A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 23, 623–630 (1975). [CrossRef]
14.	F. Zheng, Z. Chen, and J. Zhang, “Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method,” IEEE Trans. Microw. Theory Tech. 48, 1550–1558 (2000). [CrossRef]
15.	T. Namiki, “3-D ADI-FDTD method-Unconditionally stable time domain algorithm for solving full vector Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 48, 1743–1748 (2000). [CrossRef]
16.	E. L. Tan, “Efficient algorithm for the unconditionally stable 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 7–9 (2007). [CrossRef]
17.	E. L. Tan, “Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods,” IEEE Trans. Antennas Propag. 56, 170–177 (2008). [CrossRef]
18.	G. M. Hale and M. R. Querry, “Optical constants of water in the 200 nm to 200 um wavelength region,” Appl. Opt. 12, 555–563 (1973). [CrossRef] [PubMed]
19.	B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Trans. Power Delivery 14, 1052–1061 (1999). [CrossRef]
20.	D. Y. Heh and E. L. Tan, “Corrected impulse invariance method in z-transform theory for frequency-dependent FDTD methods,” IEEE Trans. Antennas Propag. 57, 2683–2690 (2009). [CrossRef]
21.	E. L. Tan, “Concise current source implementation for efficient 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 748–750 (2007). [CrossRef]
22.	C. M. Furse and O. P. Gahdhi, “A memory efficient method of calculating specific absorption rate in CW FDTD simulations,” IEEE Trans. Biomed. Eng. 43, 558–560 (1996). [CrossRef] [PubMed]
23.	S. Paker and L. Sevgi, “FDTD evaluation of the SAR distribution in a human head near a mobile cellular phone,” Elektrik 6, 227–242 (1998).
24.	I. Laakso, “FDTD method in assessment of human exposure to base station radiation,” Masters Dissertation, (2007).
25.	S. G. Garcez, C. F. Galan, L. H. Bonani, and V. Baranauskas, “Estimating the electromagnetic field effects in biological tissues through the finite-difference time-domain method,” SBMO/IEEE MTT-S Int. Microw. and Optoelectronics Conf. Proc. , 43, 58–62 (2007).
26.	J. D. Jackson, Classical Electrodynamics , 3rd ed. (John Wiley & Sons, 1998).
27.	D. M. Pozar, Microwave Engineering , 3rd ed. (Wiley, 2005).
28.	W. C. Tay, D. Y. Heh, and E. L. Tan, “GPU-accelerated fundamental ADI-FDTD with complex frequency shifted convolutional perfectly matched layer,” PIER M 14, 177–192 (2010). [CrossRef]
29.	K. Boryczko, W. Dzewinel, and D. A. Yuen, “Modeling fibrin aggregation in blood flow with discrete particles,” Comput. Methods Prog. Biomed. 75, 181–194 (2004). [CrossRef]

OCIS Codes
(170.3660) Medical optics and biotechnology : Light propagation in tissues
(050.1755) Diffraction and gratings : Computational electromagnetic methods

ToC Category:
Optics of Tissue and Turbid Media

History
Original Manuscript: January 4, 2011
Revised Manuscript: March 18, 2011
Manuscript Accepted: March 22, 2011
Published: April 12, 2011

Citation
Ding Yu Heh and Eng Leong Tan, "Modeling hemoglobin at optical frequency using the unconditionally stable fundamental ADI-FDTD method," Biomed. Opt. Express 2, 1169-1183 (2011)
http://www.opticsinfobase.org/boe/abstract.cfm?URI=boe-2-5-1169

Sort: Author | Year | Journal | Reset

References

W. G. Zijlstra, A. Buursma, and W. P. Meeuwsen-van der Roest, “Absorption spectra of human fetal and adult oxyhemoglobin, deoxyhemoglobin, carboxyhemoglobin, and methemoglobin,” Clin. Chem. 37, 1633–1668 (1991).
W. G. Zijlstra, A. Buursma, H. E. Falke, and J. F. Catsburg, “Spectrophotometry of hemoglobin: absorption spectra of rat oxyhemoglobin, deoxyhemoglobin, carboxyhemoglobin, and methemoglobin,” Comput. Biochem. Physiol. 107B, 161–166 (1994). [CrossRef]
W. G. Zijlstra, A. Buursma, and O. W. van Assendelft, Visible and Near Infrared Absorption Spectra of Human and Animal Haemoglobin: Determination and Application , (VSP, Zeist, 2000).
M. Cope, “The application of near infrared spectroscopy to non invasive monitoring of cerebral oxygenation in the newborn infant,” PhD Dissertation, (1991).
S. A. Prahl, “Tabulated molar extinction coefficient for hemoglobin in water,” http://omlc.ogi.edu/spectra/hemoglobin/summary.html (1998).
J. G. Kim, M. Xia, and H. Liu, “Extinction coefficients of hemoglobin for near-infrared spectroscopy of tissue,” IEEE Eng. Med. Biol. Mag. 24, 118–121 (2005). [CrossRef] [PubMed]
J. G. Kim and H. Liu, “Variation of haemoglobin extinction coefficients can cause errors in the determination of haemoglobin concentration measured by near-infrared spectroscopy,” Phys. Med. Biol. 52, 6295–6332 (2007). [CrossRef] [PubMed]
M. Meinke and M. Friebel, “Complex refractive index of hemoglobin in the wavelength range from 250 to 1100 nm,” Proc. SPIE 5862, 1–7 (2005).
M. Meinke and M. Friebel, “Determination of the complex refractive index of highly concentrated hemoglobin solutions using transmittance and reflectance measurements,” J. Biomed. Opt. 10, 064019 (2005). [CrossRef]
M. Meinke and M. Friebel, “Model function to calculate the refractive index of native hemoglobin in the wavelength range of 250 nm to 1100 nm dependent on concentration,” Appl. Opt. 45, 2838–2842 (2006). [CrossRef] [PubMed]
K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]
A. Taflove and S. C. Hagness, Computational Electrodynamics: the Finite-Difference Time-Domain Method , (Artech House, 2005).
A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 23, 623–630 (1975). [CrossRef]
F. Zheng, Z. Chen, and J. Zhang, “Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method,” IEEE Trans. Microw. Theory Tech. 48, 1550–1558 (2000). [CrossRef]
T. Namiki, “3-D ADI-FDTD method-Unconditionally stable time domain algorithm for solving full vector Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 48, 1743–1748 (2000). [CrossRef]
E. L. Tan, “Efficient algorithm for the unconditionally stable 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 7–9 (2007). [CrossRef]
E. L. Tan, “Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods,” IEEE Trans. Antennas Propag. 56, 170–177 (2008). [CrossRef]
G. M. Hale and M. R. Querry, “Optical constants of water in the 200 nm to 200 um wavelength region,” Appl. Opt. 12, 555–563 (1973). [CrossRef] [PubMed]
B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Trans. Power Delivery 14, 1052–1061 (1999). [CrossRef]
D. Y. Heh and E. L. Tan, “Corrected impulse invariance method in z-transform theory for frequency-dependent FDTD methods,” IEEE Trans. Antennas Propag. 57, 2683–2690 (2009). [CrossRef]
E. L. Tan, “Concise current source implementation for efficient 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 748–750 (2007). [CrossRef]
C. M. Furse and O. P. Gahdhi, “A memory efficient method of calculating specific absorption rate in CW FDTD simulations,” IEEE Trans. Biomed. Eng. 43, 558–560 (1996). [CrossRef] [PubMed]
S. Paker and L. Sevgi, “FDTD evaluation of the SAR distribution in a human head near a mobile cellular phone,” Elektrik 6, 227–242 (1998).
I. Laakso, “FDTD method in assessment of human exposure to base station radiation,” Masters Dissertation, (2007).
S. G. Garcez, C. F. Galan, L. H. Bonani, and V. Baranauskas, “Estimating the electromagnetic field effects in biological tissues through the finite-difference time-domain method,” SBMO/IEEE MTT-S Int. Microw. and Optoelectronics Conf. Proc. , 43, 58–62 (2007).
J. D. Jackson, Classical Electrodynamics , 3rd ed. (John Wiley & Sons, 1998).
D. M. Pozar, Microwave Engineering , 3rd ed. (Wiley, 2005).
W. C. Tay, D. Y. Heh, and E. L. Tan, “GPU-accelerated fundamental ADI-FDTD with complex frequency shifted convolutional perfectly matched layer,” PIER M 14, 177–192 (2010). [CrossRef]
K. Boryczko, W. Dzewinel, and D. A. Yuen, “Modeling fibrin aggregation in blood flow with discrete particles,” Comput. Methods Prog. Biomed. 75, 181–194 (2004). [CrossRef]

Baranauskas, V.
Bonani, L. H.
Boryczko, K.
Brodwin, M. E.
Buursma, A.
Catsburg, J. F.
Chen, Z.
Dzewinel, W.
Falke, H. E.
Friebel, M.
Furse, C. M.
Gahdhi, O. P.
Galan, C. F.
Garcez, S. G.
Gustavsen, B.
Hagness, S. C.
Hale, G. M.
Heh, D. Y.
Jackson, J. D.
Kim, J. G.
Liu, H.
Meeuwsen-van der Roest, W. P.
Meinke, M.
Namiki, T.
Paker, S.
Pozar, D. M.
Querry, M. R.
Semlyen, A.
Sevgi, L.
Taflove, A.
Tan, E. L.
Tay, W. C.
van Assendelft, O. W.
Xia, M.
Yee, K. S.
Yuen, D. A.
Zhang, J.
Zheng, F.
Zijlstra, W. G.

Appl. Opt.

Clin. Chem.

W. G. Zijlstra, A. Buursma, and W. P. Meeuwsen-van der Roest, “Absorption spectra of human fetal and adult oxyhemoglobin, deoxyhemoglobin, carboxyhemoglobin, and methemoglobin,” Clin. Chem. 37, 1633–1668 (1991).

Comput. Biochem. Physiol.

W. G. Zijlstra, A. Buursma, H. E. Falke, and J. F. Catsburg, “Spectrophotometry of hemoglobin: absorption spectra of rat oxyhemoglobin, deoxyhemoglobin, carboxyhemoglobin, and methemoglobin,” Comput. Biochem. Physiol. 107B, 161–166 (1994). [CrossRef]

Comput. Methods Prog. Biomed.

K. Boryczko, W. Dzewinel, and D. A. Yuen, “Modeling fibrin aggregation in blood flow with discrete particles,” Comput. Methods Prog. Biomed. 75, 181–194 (2004). [CrossRef]

Elektrik

S. Paker and L. Sevgi, “FDTD evaluation of the SAR distribution in a human head near a mobile cellular phone,” Elektrik 6, 227–242 (1998).

IEEE Eng. Med. Biol. Mag.

J. G. Kim, M. Xia, and H. Liu, “Extinction coefficients of hemoglobin for near-infrared spectroscopy of tissue,” IEEE Eng. Med. Biol. Mag. 24, 118–121 (2005). [CrossRef] [PubMed]

IEEE Microw. Wireless Comp. Lett.

E. L. Tan, “Efficient algorithm for the unconditionally stable 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 7–9 (2007). [CrossRef]
E. L. Tan, “Concise current source implementation for efficient 3-D ADI-FDTD method,” IEEE Microw. Wireless Comp. Lett. 17, 748–750 (2007). [CrossRef]

IEEE Trans. Antennas Propag.

E. L. Tan, “Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods,” IEEE Trans. Antennas Propag. 56, 170–177 (2008). [CrossRef]
D. Y. Heh and E. L. Tan, “Corrected impulse invariance method in z-transform theory for frequency-dependent FDTD methods,” IEEE Trans. Antennas Propag. 57, 2683–2690 (2009). [CrossRef]
K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]

IEEE Trans. Biomed. Eng.

C. M. Furse and O. P. Gahdhi, “A memory efficient method of calculating specific absorption rate in CW FDTD simulations,” IEEE Trans. Biomed. Eng. 43, 558–560 (1996). [CrossRef] [PubMed]

IEEE Trans. Microw. Theory Tech.

A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 23, 623–630 (1975). [CrossRef]
F. Zheng, Z. Chen, and J. Zhang, “Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method,” IEEE Trans. Microw. Theory Tech. 48, 1550–1558 (2000). [CrossRef]
T. Namiki, “3-D ADI-FDTD method-Unconditionally stable time domain algorithm for solving full vector Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 48, 1743–1748 (2000). [CrossRef]

IEEE Trans. Power Delivery

B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Trans. Power Delivery 14, 1052–1061 (1999). [CrossRef]

J. Biomed. Opt.

M. Meinke and M. Friebel, “Determination of the complex refractive index of highly concentrated hemoglobin solutions using transmittance and reflectance measurements,” J. Biomed. Opt. 10, 064019 (2005). [CrossRef]

Phys. Med. Biol.

J. G. Kim and H. Liu, “Variation of haemoglobin extinction coefficients can cause errors in the determination of haemoglobin concentration measured by near-infrared spectroscopy,” Phys. Med. Biol. 52, 6295–6332 (2007). [CrossRef] [PubMed]

PIER M

W. C. Tay, D. Y. Heh, and E. L. Tan, “GPU-accelerated fundamental ADI-FDTD with complex frequency shifted convolutional perfectly matched layer,” PIER M 14, 177–192 (2010). [CrossRef]

Proc. SPIE

M. Meinke and M. Friebel, “Complex refractive index of hemoglobin in the wavelength range from 250 to 1100 nm,” Proc. SPIE 5862, 1–7 (2005).

SBMO/IEEE MTT-S Int. Microw. and Optoelectronics Conf. Proc.

S. G. Garcez, C. F. Galan, L. H. Bonani, and V. Baranauskas, “Estimating the electromagnetic field effects in biological tissues through the finite-difference time-domain method,” SBMO/IEEE MTT-S Int. Microw. and Optoelectronics Conf. Proc. , 43, 58–62 (2007).

Other

J. D. Jackson, Classical Electrodynamics , 3rd ed. (John Wiley & Sons, 1998).
D. M. Pozar, Microwave Engineering , 3rd ed. (Wiley, 2005).
I. Laakso, “FDTD method in assessment of human exposure to base station radiation,” Masters Dissertation, (2007).
A. Taflove and S. C. Hagness, Computational Electrodynamics: the Finite-Difference Time-Domain Method , (Artech House, 2005).
W. G. Zijlstra, A. Buursma, and O. W. van Assendelft, Visible and Near Infrared Absorption Spectra of Human and Animal Haemoglobin: Determination and Application , (VSP, Zeist, 2000).
M. Cope, “The application of near infrared spectroscopy to non invasive monitoring of cerebral oxygenation in the newborn infant,” PhD Dissertation, (1991).
S. A. Prahl, “Tabulated molar extinction coefficient for hemoglobin in water,” http://omlc.ogi.edu/spectra/hemoglobin/summary.html (1998).

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper