A TLM Formulation Based on Fractional Derivatives for Dispersive Cole-Cole Media

Abstract

An auxiliary differential equation (ADE) transmission line method (TLM) is proposed for broadband modeling of electromagnetic (EM) wave propagation in biological tissues with the Cole-Cole dispersion Model. The fractional derivative problem is surmounted by assuming a linear behavior of the polarization current when the time discretization is short enough. The polarization current density is approached using Lagrange extrapolation polynomial and the fractional derivation is obtained according to Riemann definition of a fractional α-order derivative. Reflection coefficients at an air/muscle and air/fat tissues interfaces simulated in a 1-D domain are found to be in good agreement with those obtained from the analytic model over a broad frequency range, demonstrating the validity of the proposed approach.

1. Introduction

In the last two decades there has been a growing interest in the interaction between biological tissues and electromagnetic field at microwave frequencies. New promising applications of this technology in biomedical engineering like microwave imaging [1, 2], minimal invasive cancer therapies as thermal ablations [3], ultra-wide band temperature dependent dielectric spectroscopy [4] and EM dosimetry [5], rely heavily on an accurate mathematical model of the response of these tissues to an external electromagnetic field. Numerical resolution of the propagation problem within these tissues requires a robust mathematical model and the previous incorporation of the dielectric data that at all the working frequencies.

A first model of the time response in a time varying electric field of biological tissues was formulated by Debye [6] through a time decaying polarization current $j\left(t\right)=\frac{\partial P}{\partial t}$ resulting from the time variation of the polarization vector $P\left(t\right)$, the permittivity for this model is given in Eq. (1) and the corresponding argand diagram is depicted in Figure 1.

Figure 1.

Argand diagram for a unique pole Debye model obtained by representing Eq. (1) in the complexe plane.

$${\epsilon }_r^{\ast }\left(\omega \right)-{\epsilon }_{r\infty }={\epsilon }_r\left(\omega \right)-{\epsilon }_{r\infty }-j{\epsilon }_r^{\prime }\left(\omega \right)=\frac{\mathrm{\Delta }{\epsilon }_p}{1+j{\mathit{\omega \tau }}_p}$$E1

This model even if it fits the experimental results in liquids it loses its accuracy when applied over a large band of frequency or in the presence of more than one type of polar molecule. This non-Debye relaxation is attributed to the existence of different relaxation processes [7] each with its own relaxation time $\tau$ and its amplitude $\mathrm{\Delta }\epsilon$. A more accurate approximation to the behavior of this category of dielectrics is given by the sum of the individual relaxation processes leading to a multipole model, the permittivity is given as the sum of the $p$ terms corresponding to the $p$ poles each with its own amplitude. In the case of the biological tissue there is a multiple contribution of each relaxation process resulting in a broadening of the relaxation zone. Cole-Cole [8] proposed an Argand diagram in which in which ${\epsilon }_r$ the imaginary part of the complex permittivity is plotted as a function of its real part ${\epsilon }_r$, following the empirical relation:

$${\epsilon }_r^{\ast }\left(\omega \right)-{\epsilon }_{r\infty }={\epsilon }_r\left(\omega \right)-{\epsilon }_{r\infty }-j{\epsilon }_r^{\prime }\left(\omega \right)=\sum _{p=1}^p\frac{\mathrm{\Delta }{\epsilon }_p}{1+{\left(j{\mathit{\omega \tau }}_p\right)}^{{\alpha }_p}}$$E2

where ${\epsilon }_{r\infty }$ is the optical relative permittivity, $\mathrm{\Delta }{\epsilon }_p={\epsilon }_{\mathit{rs}}-{\epsilon }_{r\infty }$ is the amplitude of the $p-\mathit{th}$ pole, ${\epsilon }_s$ is the static relative permittivity, $\omega$ is the angular frequency $,{\tau }_p$ is the relaxation time, ${\alpha }_p$ the parameter that indicates the broadening of the dispersion for this pole, which in the Cole-Cole model must satisfy $0<{\alpha }_p<1$, in the case of ${\alpha }_p=1$ the model is simplified to a Debye dispersion problem and $j=\sqrt{-}1$ (Figure 2).

Figure 2.

Argand diagram for a simple pole Cole-Cole model obtained by representing Eq. (2) in the complexe plane.

The difficulty in the numerical implementation of such a model arises from the ${\alpha }_p$ parameter which is a noninteger. The consequence is a fractional order differentiation in the time domain which is more challenging compared to the Debye time domain solution. To address this problem, authors in [9] used the Letnikov fractional derivative by introducing a Stirling asymptotic formula and a recursive relation for the polarization vector. The computational demands of the FDTD scheme were considerably reduced, but nevertheless it required the storage of a large number of previous values of the polarization vector. In [10], authors used the Z transform to formulate the frequency dependence between the electric flux density and the electric field, the fractional derivative was approximated by using a polynomial method. Rekanos et al. [11] used the Pade approximation, where the fractional power term is approximated using a Pade technique. All the latter approximations were made in the context of the FDTD method.

The TLM method, even if it is flexible and wide band and being a time domain method, cannot deal with the dispersive aspect of the Cole-Cole medium directly. In [12] an approach based on a convolution product between the susceptibility and the electric field and the temporal behavior is deduced by a DFT and a nonrecursive summation leading to a considerable computational cost and a nonnegligible error at high frequencies. The causality principle is used to justify the minor dependence of the recent susceptibility values on previous ones, but the problem of the fractional Differentegration wasn’t addressed.

In this work an ADE-TLM algorithm to model the Cole-Cole dispersion. The polarization current density is approached using an extrapolation with Lagrange polynomial method and the fractional derivation is obtained using the Riemann definition of the $\alpha$-order derivative. The auxiliary differential equation is used to establish the update equation of the polarization currents which are included later in the general structure of the SCN-TLM node.

2. Formulation of the method

In the Cole-Cole model the relationship between the polarization current related to the $p_{\mathit{th}}$ pole and the electric field is given by:

$${\mathbf{J}}_p\left(\omega \right)={\epsilon }_0\mathrm{\Delta }{\epsilon }_p\frac{\mathit{j\omega }}{1+{\left(j{\mathit{\omega \tau }}_p\right)}^{{\alpha }_p}}\mathbf{E}\left(\omega \right)$$E3

where ${\left(j{\mathit{\omega \tau }}_p\right)}^{{\alpha }_p}$ is the power law function of frequency which in time domain results in a fractional derivative of order $\alpha$:

$${\mathbf{J}}_p\left(t\right)+{\tau }^{{\alpha }_p}{D}^{{\alpha }_p}{\mathbf{J}}_p\left(t\right)={\epsilon }_0\mathrm{\Delta }{\epsilon }_p\frac{\partial \mathbf{E}\left(t\right)}{\partial t}$$E4

One could consider an analytic solution for this equation, but due to the fractional derivative this can only be done for particular values of ${\alpha }_p$, hence a numerical solution is necessary with a previous discretization of the temporal values of the vectors.

In order to obtain the discretized expression of ${\mathbf{J}}_p\left(t\right)$, the polarization current vector can be approached using a Lagrange interpolation in the time interval $t_i<t<t_j$ with cardinal functions:

$$L_i\left(t\right)=\prod _{j=1;i\ne j}^n\frac{\left(t-t_j\right)}{\left(t_i-t_j\right)}$$E5

the interpolated expression of ${\mathbf{J}}_p\left(t\right)$ is then obtained:

$${\mathbf{J}}_p\left(t\right)=\sum _{i=1}^n{\mathbf{J}}^i{L}_i\left(t\right)$$E6

and by applying (5) in the time interval $n\mathrm{\Delta }t<t<\left(n+1\right)\mathrm{\Delta }t$ we obtain the cardinal functions:

$$L_n\left(t\right)=\frac{t-\left(n+1\right)\mathrm{\Delta }t}{n\mathrm{\Delta }t-\left(n+1\right)\mathrm{\Delta }t}{L}_{n+1}\left(t\right)=\frac{t-n\mathrm{\Delta }t}{\left(n+1\right)\mathrm{\Delta }t-n\mathrm{\Delta }t}$$E7

By substituting (7) in (6) the polynomial extrapolation of the polarization current density between time steps $n\mathrm{\Delta }t$ and $\left(n+1\right)\mathrm{\Delta }t$ can be found in a straightforward manner as follow:

$${\mathbf{J}}_p\left(t\right)=\frac{t-\left(n+1\right)\mathrm{\Delta }t}{n\mathrm{\Delta }t-\left(n+1\right)\mathrm{\Delta }t}{\mathbf{J}}_p^n+\frac{t-n\mathrm{\Delta }t}{\left(n+1\right)\mathrm{\Delta }t-n\mathrm{\Delta }t}{\mathbf{J}}_p^{n+1}$$E8

which after simplification can be rewritten as:

$${\mathbf{J}}_p\left(t\right)=\frac{{\mathbf{J}}_p^{n+1}-{\mathbf{J}}_p^n}{\mathrm{\Delta }t}t+\left(n+1\right){\mathbf{J}}_p^n-n{\mathbf{J}}_p^{n+1}$$E9

Furthermore, by substituting (8) into (3) we derive the auxiliary differential equation with a fractional order derivative given by:

$$\frac{{\mathbf{J}}_p^{n+1}+{\mathbf{J}}_p^n}{2}+{\tau }^{{\alpha }_p}\frac{{\mathbf{J}}_p^{n+1}-{\mathbf{J}}_p^n}{\mathrm{\Delta }t}{D}^{{\alpha }_p}t+D^{{\alpha }_p}\left(\left(n+1\right){\mathbf{J}}_p^n-n{\mathbf{J}}_p^{n+1}\right)={\epsilon }_0\mathrm{\Delta }{\epsilon }_p\frac{\partial \mathbf{E}}{\partial t}$$E10

The generalization of the derivation operator to arbitrary non-integer orders, has been subject to intensive research from mathematicians [13, 14] and in electromagnetism [15]. One of the definitions in [16] for the fractional derivative, symbolized by the operator ${\mathbf{D}}_x^{\alpha }$ is given for the power series with fractional exponents by considering a function $h\left(x\right)$ defined as a power series $h\left(x\right)={\sum }_{i=1}^n{A}_i{\left(x-a\right)}^{\nu +i/n}$ for $\nu >-1$ and $n$ a positive integer, so each term ${\left(x-a\right)}^p$ could have a noninteger exponent $p=\nu +i/n$, then according to Reimann [15] a fractional derivative of order $\alpha$ for this term is given by:

$$D_x^{\alpha }{\left(x-a\right)}^p=\frac{d^{\alpha }{\left(x-a\right)}^p}{d{\left(x-a\right)}^{\alpha }}=\frac{\mathrm{\Gamma }\left(p+1\right)}{\mathrm{\Gamma }\left(p-\alpha +1\right)}{\left(x-a\right)}^{p-\alpha }\mathit{\text{for}}x>a$$E11

where in our case $\alpha$ is a noninteger number [16].

Hence by applying the fractional derivation operator given in (11) to (10) and for $p=1$ and $a=0$ we obtain the update equation for the polarization current:

$${\mathbf{J}}_p^{n+1}=-\frac{{\varphi }_p}{{\psi }_p}{\mathbf{J}}_p^n+\frac{{\zeta }_p}{{\psi }_p}\left({\mathbf{E}}^{n+1}-{\mathbf{E}}^n\right)$$E12

where the update parameters ${\psi }_p$, ${\varphi }_p$ and ${\zeta }_p$ are:

$${\psi }_p=\frac{1}{2}+\frac{{\tau }_p^{{\alpha }_p}}{\mathrm{\Delta }t}\frac{\mathrm{\Gamma }\left(2\right)}{\mathrm{\Gamma }\left(2-{\alpha }_p\right)}{t}^{1-{\alpha }_p}-n\frac{{\tau }_p^{{\alpha }_p}}{\mathrm{\Gamma }\left(1-{\alpha }_p\right)}{t}^{-{\alpha }_p}$$E13

$${\varphi }_p=\frac{1}{2}-\frac{{\tau }_p^{{\alpha }_p}}{\mathrm{\Delta }t}\frac{\mathrm{\Gamma }\left(2\right)}{\mathrm{\Gamma }\left(2-{\alpha }_p\right)}{t}^{1-{\alpha }_p}+\left(n+1\right)\frac{{\tau }_p^{{\alpha }_p}}{\mathrm{\Gamma }\left(1-{\alpha }_p\right)}{t}^{-{\alpha }_p}$$E14

$${\zeta }_p=\frac{{\epsilon }_0\mathrm{\Delta }{\epsilon }_p}{\mathrm{\Delta }t}$$E15

To get the update equation of the electric field components we start from the Maxwell-Ampere equation, and by including the conductivity term:

$$\nabla \times \mathbf{H}={\epsilon }_0{\epsilon }_{\infty }\frac{\partial \mathbf{E}}{\partial t}+{\sigma }_0\mathbf{E}+\sum _{p=1}^p{\mathbf{J}}_p$$E16

where ${\sigma }_0$ is the static ionic conductivity, Eq. (16) formulated at time step $n+\frac{1}{2}$ and by approaching ${\mathbf{J}}^{n+\frac{1}{2}}$ by the average of its values at time steps $n$ and $n+1$

$${\mathbf{J}}_p^{n+\frac{1}{2}}=\frac{{\mathbf{J}}_p^{n+1}+{\mathbf{J}}_p^n}{2}$$E17

Finally, the updated electric field is given by:

$$E_{x,y,z}^{n+1}=E_{x,y,z}^n\left(1-\frac{2{\sigma }_0\mathrm{\Delta }t}{D}\right)-\sum _p\frac{\left(1-\frac{{\varphi }_p}{{\psi }_p}\right)\mathrm{\Delta }t}{D}{J}_{\mathit{px},y,z}^n+\frac{2\mathrm{\Delta }t}{D}{\left(\nabla \times \mathbf{H}\right)}_{x,y,z}^{n+\frac{1}{2}}$$E18

$$D=2{\epsilon }_0{\epsilon }_{\infty }+\sum _p\frac{{\zeta }_p\mathrm{\Delta }t}{{\psi }_p}+{\sigma }_0\mathrm{\Delta }t$$E19

3. The TLM Formalisme

The TLM method is a numerical technique based on the discretization of the computational domain according to Huygens principle as an alternative to the Maxwell equations used in the FDTD method [17]. In this method the simulation domain is discretized in cells where a series of uniform transmission lines, parallel and series stubs and additional sources are used to take account of the real characteristics of the propagation through the medium as given by Maxwell equations. Therefore, instead of electric and magnetic field components the electromagnetic field is represented by voltage and current waves, propagating through the unit cell circuit referred to as symmetrical condensed node (SCN). The relationship between the electromagnetic field and the voltage and current waves at the time step $n\mathrm{\Delta }t$ and at the center of the node are formulated as follows [18].

As can be seen in Figure 3(b), the SCN consists of interconnected transmission lines. This structure models a unit cell of the propagation medium as in Figure 3(a). Each face of the cell corresponds to two ports orthogonal to each other and labeled from 1 to 12, these stubs model the propagation through free space and have therefore its characteristic impedance $Z_0=\sqrt{{\epsilon }_0{\mu }_0}$. Three additional open circuit (13,14,15) stubs must be added to the node to account for the dispersive behavior of the medium. The SCN nodes are connected to each other to form the simulation domain given in Figure 4 for a 1-D propagation.

Figure 3.

Structure of the TLM symmetrical condensed node, (a) unit volume of thedielectric material with dimensions Δx, Δy,and Δz, (b) equivalent SCN.

Figure 4.

Structure of the 1-D SCN grid used in this work to model the simulation domain.

The relationship between the electromagnetic field and the voltage and current waves at the time step $n\mathrm{\Delta }t$ and at the center of the node are formulated as follows [18]:

$$E_{x,y,z}^n=\frac{V_{x,y,z}^n}{\mathrm{\Delta }l}{H}_{x,y,z}^n=\frac{I_{x,y,z}^n}{\mathrm{\Delta }l}$$E20

In this algorithm as in [11, 19] the dispersive properties of the medium are accounted for by adding voltage sources ${}_{s}V$ to each node, these sources constitute the counterparts of the terms in the FDTD formulation that have no equivalents in the TLM formalism. Therefore the voltage at the center of the node, as deduced from the conservation laws [20], becomes:

$$V_{x,y,z}^{n+1}=V_{x,y,z}^n+\frac{1}{4+Y_{\mathit{ocx},y,z}}\left({}_s{V}_{x,y,z}^{n+1}{+}_s{V}_{x,y,z}^n\right)+\frac{4}{4+Y_{\mathit{ocx},y,z}}\frac{\mathrm{\Delta }l\mathrm{\Delta }t}{{\epsilon }_0}{\left(\nabla \times \mathbf{H}\right)}_{x,y,z}^{n+\frac{1}{2}}$$E21

where $Y_{\mathit{ocx},y,z}$ indicates the normalized admittance of the open circuit stub added to the node.

When expressed in terms of incident voltages on the corresponding stubs and accounting for the voltage sources injected in the stubs (16,17,18) to complete the model of the Cole-Cole medium:

$$V_{x,y,z}^{n+1}=\frac{2}{4+Y_{\mathit{ox},y,z}}\left[V_{\mathrm{1,3,5}}^i+V_{\mathrm{2,4,6}}^i+V_{\mathrm{9,8,7}}^i+V_{\mathrm{12,11,10}}^i+Y_{\mathit{oc},x,y,z}{V}_{\mathrm{13,14,15}}^i+{\frac{1}{2}}_s{V}_{\mathrm{16,17,18}}^{n+1}\right]$$E22

where the subscript $i$ indicates the incident pulses on the indicated stubs.

In the TLM formalism the update equation issued from the ADE in the Eq. (12) becomes:

$$J_{\mathit{px},y,z}^{n+1}=-\frac{{\varphi }_p}{{\psi }_p}{J}_{\mathit{px},y,z}^n+\frac{{\zeta }_p}{{\psi }_p\mathrm{\Delta }l}\left(V_{x,y,z}^{n+1}-V_{x,y,z}^n\right)$$E23

and the electric field update equation is also formulated in the TLM formalism as:

$$V_{x,y,z}^{n+1}=V_{x,y,z}^n\left(1-\frac{2{\sigma }_0\mathrm{\Delta }t}{D}\right)-\sum _p\frac{\left(1-\frac{{\varphi }_p}{{\psi }_p}\right)\mathrm{\Delta }t\mathrm{\Delta }l}{D}{J}_{\mathit{px},y,z}^n+\frac{2\mathrm{\Delta }t\mathrm{\Delta }l}{D}{\left(\nabla \times \mathbf{H}\right)}_{x,y,z}^{n+\frac{1}{2}}$$E24

The analogy between Eqs. (22) and (24) the expression of the normalized admittance of the stub added to the SCN node is obtained straightforward:

$$Y_{\mathit{oc}}=4\left(\frac{D}{2{\epsilon }_0}-1\right)$$E25

and the update equation for voltage sources at the center of the node:

$${}_s{V}_{x,y,z}^{n+1}={-}_s{V}_{x,y,z}^n-\frac{4{\sigma }_0\mathrm{\Delta }t}{{\epsilon }_0}{V}_{x,y,z}^n-\sum _p\frac{2\left(1-\frac{{\varphi }_p}{{\psi }_p}\right)\mathrm{\Delta }t\mathrm{\Delta }l}{{\epsilon }_0}{J}_{\mathit{px},y,z}$$E26

that can be simplified to

$${}_s{V}_{x,y,z}^{n+1}={-}_s{V}_{x,y,z}^n+C_1{V}_{x,y,z}^n+\sum _p{C}_{2p}{J}_{\mathit{px},y,z}$$E27

with update constants:

$$C_1=-\frac{4{\sigma }_0\mathrm{\Delta }t}{{\epsilon }_0}$$E28

$$C_{2p}=-\frac{2\left(1-\frac{{\varphi }_p}{{\psi }_p}\right)\mathrm{\Delta }t\mathrm{\Delta }l}{{\epsilon }_0}$$E29

A final step to establish the TLM formulation is to express the scattering process on the capacitive stubs:

$${\left[V_{\mathrm{13,14,15}}^r\right]}^n={\left[V_{x,y,z}\right]}^n-{\left[V_{\mathrm{13,14,15}}^i\right]}^n$$E30

4. Simulation and results

In order to verify the new proposed algorithm, and for the sake of simplicity a one dimensional problem is simulated, where a Cole-Cole medium (fat or muscle) occupies the region $z\ge 0$, while the rest of space is air, 2 dimensional and 3 dimensional propagation and reflection models can be naturally deducted by ensuring the connexion between TLM cells in the $y$ axis using ports 1,5,7 and 12, to account for the propagation along the $x$ axis ports 3,6,10 and 11 must be connected to adjacent cells as in Figure 5.

Figure 5.

2-D ADE-TLM simulation of an incident and reflected pulse on the air/muscle, the tissue is modeled by a 4 pole Cole-Cole model, the interface is located at the TLM cell $z=100$, at time steps 160 (a) and 218 (b).

The incident wave is a derivative Gaussian pulse given by ${\overrightarrow{E}}^{\mathit{inc}}=\frac{\left(t-t_0\right)}{{\tau }^2}\mathit{exp}\left(-\frac{4\pi ({\left(t-t_0\right)}^2}{{\tau }^2}\right){\widehat{e}}_x$ where $t_0=200\mathrm{\Delta }t$ and $\tau =220\mathrm{\Delta }t$ polarized in the $x$ direction, and propagates along the $z$ axis. In the TLM implementation, the considered 1-D simulation domain (see Figure 6) consists of a grid of 1000 cells interconnected on the z axis as in Figure 4, 500 of which were used to model the Cole-Cole medium and the remaining cells were used to model the air. The cell dimension $\mathrm{\Delta }l$ satisfies the stability condition $\mathrm{\Delta }t=\mathrm{\Delta }l/2c$, the time step is set as $\mathrm{\Delta }t=0.125\mathit{ps}$. The simulation domain is truncated by an unsplit PML layer [21] of 10 cells at the beginning and the end.

The fourth order Cole-Cole parameters for fat tissue as well as for muscle tissue are listed in Table 1 [9]. The values of the electric field are recorded at a point P located at the center of the SCN node 10 cells before the air-medium interface. Figure 7(a) shows the incident impulse at iteration 399 propagating in the air and Figure 7(b) at iteration 927 depicts the reflected and transmitted impulse on the air/muscle interface.

Tissue	${\epsilon }_{\infty }$	${\sigma }_0$	$\mathrm{\Delta }{\epsilon }_1$	${\tau }_1$ (ps)	${\alpha }_1$	$\mathrm{\Delta }{\epsilon }_2$	${\tau }_2$ (ps)	${\alpha }_2$	$\mathrm{\Delta }{\epsilon }_3$	${\tau }_3$ ($\mathit{\mu s}$)	${\alpha }_3$	$\mathrm{\Delta }{\epsilon }_4$	${\tau }_4$ (ps)	$\mathrm{\Delta }{\epsilon }_4$
Muscle	2.5	0.035	9	7.96	0.8	35	15.92	0.9	3.3E4	159.15	0.95	1.0E7	15.915	0.99
Fat	4	0.2	50	7.23	0.9	7000	353.68	0.9	1.2E6	318.31	0.9	2.5E7	2.274	1

Table 1.

Cole-Cole constants for tissues of muscle and fat.

Figure 6.

2-D ADE-TLM simulation of a reflected and transmitted pulse on the air/muscle , the tissue is modeled by a 4 pole Cole-Cole model, the interface is located at the TLM cell z = 100 , at (a) time step 305, and (b) at time step 446.

In Figures 8 and 9, the simulation results for both cases (fat and muscle) are compared to the theoretical reflection at the air/Cole-Cole medium interface which can be obtained by using the following equation [22]:

$$\mid R\mid =\frac{\mid 1-{\sqrt{\epsilon }}_r\mid }{\mid 1+{\sqrt{\epsilon }}_r\mid }$$E31

Figure 7.

Propagation of an incident derivative Gaussian pulse through Cole-Cole medium model of human muscle (a) at iteration 399 (incident pulse) (b) at iteration 677 (reflected and transmitted pulse).

Figure 8.

Reflection coefficient (dB) at normal incidence on an air/muscle interface.

A good agreement over the whole frequency band is observed. The slight discrepancy between the numerical and theoretical results can be ascribed to the approximation made to the polarization current value when performing the Lagrange extrapolation on each iteration. The propagation and reflection through a 2-D TLM grid modeling the incidence on an air/muscle interface is presented in Figures 5 and 10 at different time steps, the interface is located at $z=100$ in a $10$ by $200$ cells grid, Figure 5(a) shows the Gaussian wave before the arrival on the interface while Figure 10(a) and (b) depicts it after reflection and transmission. A more precise simulation that targeted the Dielectric properties of pathological thyroid tissue types including adenoma, thyroiditis and the properties of the healthy thyroid based on the results of the experimental study presented in a previous published work in [2], and for the healthy thyroid tissue [23], in both works experimental results are used to extract the 2 poles Cole-Cole model parameters. Based upon this parameters we conducted a simulation for the 4 tissue types, the results of this simulation are presented in Figure 11(a) and (b) for the air/adenoma and the air/normal thyroidian interface respectively, Figure 12(a) and (b) depicts the reflection coefficient on an air/thyroyditis affected thyroidian tissue and an air/tumor respectively also in this case a perfect agreement over the working frequency is observed.

Figure 9.

Reflection coefficient (dB) at normal incidence on an air/fat interface.

Figure 10.

2-D ADE-TLM simulation of a reflected and transmitted pulse on the air/muscle, the tissue is modeled by a 4 pole Cole-Cole model, the interface is located at the TLM cell $z=100$, at time steps are 305 and 446.

Figure 11.

Reflection coefficient (dB) at normal incidence on an air/adenoma (a), air/Normal tyroidian tissue interface (b).

Figure 12.

Reflection coefficient (dB) at normal incidence on an air.

5. Conclusion

Numerical methods are an essential part of the modeling process, new propagation media with anomalous relaxation properties impose innovative modeling methods. An effective ADE-TLM formulation way to model electromagnetic wave propagation in biological tissues using Cole-Cole dispersion model. The fractional-order derivative in the Cole-Cole model is tackled by using a polynomial extrapolation of the polarization current. This approximation reduces considerably both the numerical cost of the simulation and its complexity since only two previous values of the $J$ field are needed at each iteration. a Reimann derivation on the obtained polynomial extrapolation of the polarization current is then performed to solve the fractional order derivative. The presented results indicate the accuracy of our approach.

Conflict of interest

The authors declare no conflict of interest.