An ADE-TLM Algorithm for Modeling Wave Propagation in Biological Tissues with Debye Dispersion

In this paper, we present a Transmission Line Matrix (TLM) algorithm for the simulation of electromagnetic wave interaction with a Debye dispersive medium. This new formulation is based on the use of the polarization currents in the medium. The auxiliary differential equation (ADE) method is considered to deal with dispersion after the classical discretization. The accuracy and efficiency of this approach were tested on 1D Debye medium by calculating the reflection coefficient on an air-dielectric interface. The potential of the developed algorithm to model the existence of tumors in a human breast is also demonstrated. The obtained results compared with the analytic model show a good agreement. The number of operations needed for each iteration has been reduced, hence the computational time in comparison with time convolution techniques, while maintaining a comparable numerical accuracy. INDEX TERMS Computational electromagnetics, numerical methods, Transmission Line Matrix, biological system modeling, Debye medium, microwave imaging.


I. INTRODUCTION
I N In the last two decades, the interest in modeling and simulating electromagnetic wave propagation in dispersive media has considerably grown, fostered by the fact that many materials exhibit frequency dependent electromagnetic properties like polymers, dielectric fluids and more recently of particular interest are biological tissues. Modelling such media by means of time domain numerical methods implies incorporating the dispesive effect to the formulation.
The numerical solutions proposed for modeling such dispersive media could be grouped in three main categories: recursive convolution methods based on a recursive computation of a convolution integral between the frequencydependent susceptibility and the electric field [1] [2], the Z-transform method [3] [4], and the auxiliary differential equation method (ADE) [5]. The latter is of particular interest owing to its easier arithmetic implementation [6], first introduced by Kishawa to model Debye media, Lorentz media, and media obeying to the Cole-Cole circular relaxation arc model [7] [8]. It was then extended to cover non linear dispersive media by Goorjian and Taflove [9],and more recently extended to cover 2-D Kerr and Raman nonlinear dispersive media response to an optical pulse [10]. Association of the ADE method with previous cited techniques was also prooven to be efficient for dealing with Drude critical point model [11].
In the TLM method, a first modelisation and simulation of a non linear dispersive medium in a 2D TLM mesh was described in [12] where the behaviour of the dielectric was modeled as an anharmonic oscillator and by connecting a one-port nonlinear network to the nodes of the TLM mesh allowing the simulation of the propagation of a soliton in the despersive medium. In [13] the dispersion is illustrated by incorporating a nonlinear capacitor represented by an open circuit stub to the TLM node. The authors in [14] extended the TLM method to take into account the presence of anisotropic and dispersive first order Debye medium by adding voltage sources on each node of the TLM mesh. This scheme required the calculation of a convolution product between the electric suceptibility and the electric field on each iteration in order to obtain the value of the electric flux density field D. Later on, Barba et al. [15] extended the same method to a Cole-Davidson and a Cole-Cole dispersive medium.
Extension of these simulation methods to the biological tissues have gained relevance due to their application in the emergent field of microwave tomography and bioelectromagnetic dosimetry. An accurate model would give a reliable prediction of these tissues physical response when exposed to electromagnetic waves [16] [17], leading to a more efficient design of wearable or implantable devices and the design of a more efficient imaging systems [18] [19]. An efficient model for these tissues in the microwave range is given by the Cole-Cole model but with a major flaw: its computational complexity, mainly due to the fractional power of the poles, which in time domain results in a partial derivative of the polarisation current [20] [21] or the polarisation vector [22]. An accurate approximation to the Cole-Cole model but with less computational complexity is given by the debye model [23]. This model was first proposed by Debye [24] for materials with permanent dipole moments, it is mostly used to model electromagnetic wave interactions with water-based substances. It is worth mentioning that the ADE technique has not been applied before to model the Debye despersive media in the TLM method. In this work, an ADE-TLM approach to numerically model a double pole Debye medium of the human skin, fat and tumorous breast tissue is implemented. Unlike [14], in this work voltage sources are added on each TLM node [25] following the differential equation of the polarisation currents without requiring the calculation of a convolution product between the electric susceptibility and the electric field on each iteration. The numerical results are compared to those obtained by the analytic model.

II. FORMULATION AND EQUATIONS
In the Debye medium, absorption is caused by dipole relaxation [26] and for a medium fitted with multiple relaxation poles, each pole has an amplitude ∆ p = sp − ∞ and a relaxation time τ p and the complex permittivity is written: where ∞ is the infinite frequency relative permittivity, 0 is the dielectric constant of the free space, sp is the static permittivity and τ p represents the relaxation time for the pth pole, and σ 0 is the electric conductivity at DC. In this work the Debye dispersion will be fitted with a 2 Debye pole model ( p=1,2). We derive a set of auxiliary differential equations for the polarization currents by transforming the frequency representation of the permittivity into a time domain differential operator. The polarization current vector for the pth pole is related to the electric field through : Taking the inverse Fourier transform we obtain the differential equation describing the temporal relation between J and E: Applying central difference approximation at the time instance t = (n + 1 2 )∆t with the discrete time step ∆t the update equation for the polarization current vector is obtained by [5] : with update coefficients : Once the updating relation for the polarization current vector is obtained we can consider the evolution of the time dependent Maxwell equation: which after dicretization yields the relation between E n+1 and E n and the polarization currents attributed to the different debye poles: The update equation for the electric field can then be obtained as follows: where

III. THE TLM FORMALISM
The TLM method is based on the analogy between the Maxwell's equations and the signal propagation through a network of multiport circuits. Hence instead of electric and magnetic field components the electromagnetic field is represented by voltage and current waves. Each circuit is represented by a symmetrical condensed node (SCN), consisting of interconnected transmission lines. As shown in Fig.1, this structure models a unit cell of the propagation medium and each face of the cell corresponds to two ports orthogonal to each other and labeled from 1 to 12 [29],with the supplementary ports from 13 to 15 referring to the capacitive stubs. To take into account the dispersive properties of the medium voltage sources s V are included on ports 16 to 18 [15] [31].The SCN nodes are connected to each other to form the simulation domain given in Figure 2 for a 1-D propagation.

FIGURE 1: The Symmetrical condensed node (SCN)
In the case of a uniform TLM meshing of a 3-D domain, the relationship between the electromagnetic field and the voltage and current waves at the time step n∆t at the center of the node are formulated as follows: In this scheme, the update equation obtained from the ADE in the Eq.(9) becomes: compared to the voltage update equation in the TLM scheme: where s V n are the sources added to the 16,17,18 symmetrical Condensed Node ports to take into account the dispersive property of the debye medium and Y ocx,y,z is the normalized admittance of the capacitive stub.
From the analogy between eqations (12) and (13) the expression of the normalized admittance of the stub added to the SCN node is obtained straightforward : and the update equation for the voltage sources : where The final update equation for the voltage in the center of the SCN node after developing the curl operator by means of incident pulses: where the exponent i indicates the incident pulses, the subscripts (1,...,12) the port number in the standard SCN, (13,14,15) denote the capacitive ports,while at the ports (16,17,18) voltage sources are injected .  Figure 3 illustrates the successive steps leading to the update of the electric field.The reduced number of arrows in this chart shows that the method requires less arithmetic operations than other algorithms such as PLRC [27]. The dashed lines refer to the operations particular to this ADE-TLM model. Starting from J n and E n and s V n the update of the electric sources s V n is obtained using Eq. (15). The latter is used alongside with the value of E n to retrieve the value of E n+1 through Eq. (18) and finally the update of the polarisation current vector is obtained using Eq.(4). It is worth mentioning that the electric field components are obtained from the voltage components in each iteration VOL. 9, NO. 3, NOVEMBER 2020

IV. SIMULATION AND DISCUSSION
To test the accuracy of the proposed TLM-ADE algorithm, several cases have been studied, the choices made are driven by the aforementioned interest in the capability of the ADE-TLM Debye approximation to describe the behavior of biological tissues in the perspective of implementing it in a microwave imaging technique. Under this scope, the parameters of the simulated two pole Debye model of the human skin ,breast tissue in both healthy and tumorous cases taken as in [28] are listed in Tab.1. The computational domain is discretised as a one-dimensional grid of 800 cells aligned on the z-axis with a dimension of ∆l = 75µm each , the Debye dielectric medium is located at 300 ≤ z ≤ 700, the remaining of the domain is occupied by air. The domain is truncated by an unsplit PML layer [32] of 10 cells at the beginning and the end. The TLM grid is illuminated by a modulated Gaussian wave E x (t) = sin(2πf ∆t).exp( −(n∆t−tmax) 2 ) at a point 20 cells from the limit of the grid. The wave propagates in the z direction and the electric field is polarized in the x direction. The simulation is performed for 10000 time steps corresponding to a total time of 1.25ns, the transient reflected field values were retrieved at the output point 10∆l TLM nodes before the air-dielectric interface. Fourier transform is then used to transpose the results to the frequency domain. Final results are compared to the analytic prediction given by: where r is the complex permittivity of the Debye medium. it is worth to mention that in order for the ADE-TLM to be computationally stable, the following condition must be satisfied ( c.∆t ∆l ≤ 1) where ∆l is the spacial discretization and ∆(t) is the time-step and c is the speed of light, satisfying this requirement the values of the time and space steps are ∆l = 75µm and ∆(t) = 0.125ps. First, a simulation of an air water and air muscle incidence is carried out. Air parameters are taken as = 1, µ = 1 and σ 0 = 0S.m −1 . the water is modeled as a single pole Debye medium with ∞ = 1.8, s = 81 and τ w = 9.4ps and the two pole Debye parameters of the muscle are ∞ = 32, ∆ 1 = 2159, ∆ 2 = 24.99,τ 1 = 46.25ns,τ 2 = 0.0907ns with σ 0 = 0.106S.m −1 . figure 4 shows the incident impulse on the dielectric slab at iteration 173,while figure 5 shows both reflected and transmitted pulses through the slab at iteration 677. Figure  6 depicts the time domain simulation where we can see the time domain incident and reflected waveforms.     figure 8 show the simulation results and the theoretical reflection coefficients for an air-water interface and air-muscle interface respectively. Very good agreement between the two data sets for both cases over the whole frequency range is obtained, showing the accuracy of the proposed algorithm.     Figure 11. Also in this case we can see that a good agreement is achieved, with a maximum error of 0.4 dB at frequencies below 20 GHz and 0.1dB at the upper limit of the frequency range of the simulation. IT is worth mentioning that the above listed Debye fitting parameters are dependent on the hydration of the biological tissue , an approach to this dependence was presented by Robinson et al. [33] for static properties. In [34] authors used a perturbation method to demonstrate that for an average dehydration of 2% (measured in body weight loss) there is a 1% decrease in permittivity.

V. REFLECTION ON THREE LAYERED MEDIUM WITH DIFFERENT DEBYE RELAXATION
To accurately assess the microwave imaging techniques the simulation domain must be of a heterogeneous nature. In this work the heterogeneity is modeled by a layered setup, each layer has the Debye parameters of the corresponding tissue in human breast morphology, the computational domain is bounded by Unsplit PML layers. The Debye medium constants are taken as in Table 1 [35]. The computational domain is 800 cells long and as in the previous simulation the propagation is along the z-axis and the electric field is polarized in the x direction. The Debye medium is defined first as a skin layer of 10 cells thickness followed by healthy breast tissue . A first simulation with a malignant tumor located at a variable depth in the healthy tissue is carried out. Figure 13 depicts the results of a parametric simulation study, where the parameter is the depth of the tumor ranging from 50 to 250 cells with a 50 cell step. The obtained results show that there is a difference in the reflection coefficient even for a minimal variation of the tumor depth of 50 cells corresponding to a real variation of 3.7 mm. Furthermore, we have calculated the reflection coefficient with the width of the tumor as a parameter at a fixed depth. The results are presented in figure 14. We can clearly see that there is a variation of the reflection magnitude with the increase of the tumor size.In both cases the reflection coefficient is sensitive to the position of the tumor.

VI. CONCLUSIONS AND FUTURE WORK
The authors have proposed an ADE-TLM method to simulate a multi-pole Debye medium. This model was compared to the theoretical results for the skin, muscle, brain, and breast tissues. Furthermore, the authors applied the model to simulate the response of a breast with tumor. A parametric study was conducted to show the ability of the model to discriminate the tumor depth and width. The proposed model has the advantage of being more simple because it doesn't have to deal with the convolution product as in the recursive convolution methods. Future work will incorporate the effect of the heterogeneity of the breast in order to extract a more precise information about the location and characteristics of the tumor from the back scatter.