Cole-Cole constants for tissues of muscle and fat.
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.
- computational electromagnetics
- numerical methods
- transmission line matrix
- biological system modeling
- Cole-Cole medium
- fractional derivative
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 , ultra-wide band temperature dependent dielectric spectroscopy  and EM dosimetry , 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  through a time decaying polarization current resulting from the time variation of the polarization vector , the permittivity for this model is given in Eq. (1) and the corresponding argand diagram is depicted in Figure 1.
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  each with its own relaxation time and its amplitude . 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 terms corresponding to the 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  proposed an Argand diagram in which in which the imaginary part of the complex permittivity is plotted as a function of its real part , following the empirical relation:
where is the optical relative permittivity, is the amplitude of the pole, is the static relative permittivity, is the angular frequency is the relaxation time, the parameter that indicates the broadening of the dispersion for this pole, which in the Cole-Cole model must satisfy , in the case of the model is simplified to a Debye dispersion problem and (Figure 2).
The difficulty in the numerical implementation of such a model arises from the 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  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 , 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.  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  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 -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 pole and the electric field is given by:
where is the power law function of frequency which in time domain results in a fractional derivative of order :
One could consider an analytic solution for this equation, but due to the fractional derivative this can only be done for particular values of , 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 , the polarization current vector can be approached using a Lagrange interpolation in the time interval with cardinal functions:
the interpolated expression of is then obtained:
and by applying (5) in the time interval we obtain the cardinal functions:
which after simplification can be rewritten as:
The generalization of the derivation operator to arbitrary non-integer orders, has been subject to intensive research from mathematicians [13, 14] and in electromagnetism . One of the definitions in  for the fractional derivative, symbolized by the operator is given for the power series with fractional exponents by considering a function defined as a power series for and a positive integer, so each term could have a noninteger exponent , then according to Reimann  a fractional derivative of order for this term is given by:
where in our case is a noninteger number .
where the update parameters , and are:
To get the update equation of the electric field components we start from the Maxwell-Ampere equation, and by including the conductivity term:
where is the static ionic conductivity, Eq. (16) formulated at time step and by approaching by the average of its values at time steps and
Finally, the updated electric field is given by:
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 . 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 and at the center of the node are formulated as follows .
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 . 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.
The relationship between the electromagnetic field and the voltage and current waves at the time step and at the center of the node are formulated as follows :
In this algorithm as in [11, 19] the dispersive properties of the medium are accounted for by adding voltage sources 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 , becomes:
where 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:
where the subscript indicates the incident pulses on the indicated stubs.
In the TLM formalism the update equation issued from the ADE in the Eq. (12) becomes:
and the electric field update equation is also formulated in the TLM formalism as:
and the update equation for voltage sources at the center of the node:
that can be simplified to
with update constants:
A final step to establish the TLM formulation is to express the scattering process on the capacitive stubs:
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 , 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 axis using ports 1,5,7 and 12, to account for the propagation along the axis ports 3,6,10 and 11 must be connected to adjacent cells as in Figure 5.
The incident wave is a derivative Gaussian pulse given by where and polarized in the direction, and propagates along the 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 satisfies the stability condition , the time step is set as . The simulation domain is truncated by an unsplit PML layer  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 . 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.
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 :
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 in a by 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 , and for the healthy thyroid tissue , 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.
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 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.
Martellosio A, Pasian M, Bozzi M, Perregrini A, Mazzanti F, Svelto PE, et al. L: 0.5-50 ghz dielectric characterisation of breast cancer tissues. Electronics Letters. 2015; 51(13):974-975
Gavazzi S, Limone P, De Rosa G, Molinari F, and Vecchi G: Comparison of microwave dielectric properties of human normal, benign and malignant thyroid tissues obtained from surgeries: A preliminary study, Biomedical Physics Engineering Express, vol. 4, no. 4, p. 047 003, 2018
Ruvio G, Eaton-Evans J, Shahzad A, OHalloran M. Numerical evaluation of microwave thermal ablation to treat small adrenocortical masses. International Journal of RF and Microwave Computer-Aided Engineering. 2018; 28(3):e21236
Ley S, Schilling S, Fiser O, Vrba J, Sachs J, Helbig M. Ultra-wideband temperature dependent dielectric spectroscopy of porcine tissue and blood in the microwave frequency range. Sensors. 2019; 19(7):1707
Chakarothai J, Wake K, Watanabe S. Convergence of a single-frequency FDTD solution in numerical dosimetry. IEEE Transactions on Microwave Theory and Techniques. 2016; 64(3):707-714
Debye P. Part i. dielectric constant. energy absorption in dielectrics with polar molecules. Transactions of the Faraday Society. 1934; 30:679-684
Cole KS, Cole RH. Dispersion and absorption in dielectrics i. alternating current characteristics. The Journal of Chemical Physics. 1941; 9(4):341-351
Rekanos IT, Yioultsis TV. Approximation of grunwald-letnikov fractional derivative for FDTD modeling of Cole-Cole media. IEEE Transactions on magnetics. 2014; 50(2):181-184
Guo B, Li J, Zmuda H. A new FDTD formulation for wave propagation in biological media with Cole-Cole model. IEEE Microwave and Wireless Components Letters. 2006; 16(12):633-635
Rekanos IT, Papadopoulos TG. An auxiliary differential equation method for FDTD modeling of wave propagation in Cole-Cole dispersive media. IEEE Transactions on Antennas and Propagation. 2010; 58(11):3666-3674
Barba I. Cabeceira A C L, Panizo M, and Represa J: Modelling dispersive dielectrics in TLM method, International Journal of Numerical Modelling: Electronic Networks. Devices and Fields. 2001; 14(1):15-30
Samko SG, Kilbas AA, Marichev OI, et al. Fractional integrals and derivatives. Gordon and Breach Science Publishers, Yverdon Yverdon-les-Bains, Switzerland, 1993. In: vol. 1993
Techniques in the fractional calculus,” in The Fractional Calculus, ser. Mathematics in Science and Engineering, K. B. Oldham and J. Spanier, Eds., vol. 111, Elsevier, 1974, pp. 133–160
N. Engheta: n the role of fractional calculus in electromagnetic theory,” IEEE Antennas and Propagation Magazine, vol. 39, no. 4, pp. 35–46, Aug. 1997
Engheta N. On fractional calculus and fractional multipoles in electromagnetism. IEEE Transactions on Antennas and Propagation. Apr. 1996; 44(4):554-566
Wharmby AW, Bagley RL. The application of the fractional calculus model for dispersion and absorption in dielectrics in terahertz waves. International Journal of Engineering Science. 2015; 93:1-12
Jin H, Vahldieck R. Direct derivations of TLM symmetrical condensed node and hybrid symmetrical condensed node from maxwell’s equations using centered differencing and averaging. IEEE Transactions on Microwave Theory and Techniques. Dec. 1994; 42(12):2554-2561
Christopoulos C, The Transmission-Line Modeling (TLM) Method in Electromagnetics. Morgan Claypool, 2006
Cabeceira ACL, Barba I, Grande A, Represa J. A 2D-TLM model for electromagnetic wave propagation in chiral media. Microwave and Optical Technology Letters. 2005; 46(2):180-182
Yaich MI, Khalladi M. The far-zone scattering calculation of frequency-dependent materials objects using the tlm method. IEEE Transactions on Antennas and Propagation. Nov. 2002; 50(11):1605-1608
Yaich MI, Kanjaa M, Adraoui S, Mounirh K, Khalladi M. An unsplit formulation of the 3D-PML absorbing boundary conditions for TLM-method in time domain. In: 2018 6th International Conference on Multimedia Computing and Systems (ICMCS). May 2018. pp. 1-5
Taflove A, Hagness SC. Computational electrodynamics: the finite-difference time-domain method. Artech house. 2005
Gabriel S, Lau RWand Gabriel C1996b The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues Phys. Med. Biol. 41 227193