Abstract
The interaction of electromagnetic fields and biological tissues has become a topic of increasing interest for new research activities in bioelectrics, a new interdisciplinary field combining knowledge of electromagnetic theory, modeling, and simulations, physics, material science, cell biology, and medicine. In particular, the feasibility of pulsed electromagnetic fields in RF and mm-wave frequency range has been investigated with the objective to discover new noninvasive techniques in healthcare. The aim of this contribution is to illustrate a novel Finite-Difference Time-Domain (FDTD) scheme for simulating electromagnetic pulse propagation in arbitrary dispersive biological media. The proposed method is based on the fractional calculus theory and a general series expansion of the permittivity function. The spatial dispersion effects are taken into account, too. The resulting formulation is explicit, it has a second-order accuracy, and the need for additional storage variables is minimal. The comparison between simulation results and those evaluated by using an analytical method based on the Fourier transformation demonstrates the accuracy and effectiveness of the developed FDTD model. Five numerical examples showing the plane wave propagation in a variety of dispersive media are examined.
1. Introduction
Recently, multifrequency and ultrawideband microwave technologies have receiving increasing attention due to their promising applications in biological and medical domains [1–3]. Several studies concern microwave dielectric spectroscopy and imaging since their represent potential solutions for early–stage cancer detection and treatment as well as for cancer investigations and diagnostic [4, 5]. Many research activities are focused on the use of electromagnetic waves for discriminating malignant tissues from healthy ones. Moreover, recent applications employing electromagnetic technology include medical implant communication service, wireless medical telemetry service, and medical body area network [6–8]. Electromagnetic wave technologies have also been introduced as a rapid method of delivering high temperatures to destroy the cancer cells during thermal ablation [9]. Finally, continuous and pulsed electric field are successfully used in a variety of therapeutic and diagnostic applications such as hyperthermia, electroporation, and treatment of specific diseases [10].
All these technologies involve the interaction of electromagnetic fields with living tissues. As a result, detailed theoretical modeling and computational techniques are essential to gain insight into the several phenomena occurring within the biological materials subject to an imposed electromagnetic field. For example, the right knowledge of the electromagnetic field distribution inside biological tissues is necessary to demonstrate the effects of electromagnetic radiation absorption. Furthermore, numerical simulations under various conditions can be used to identify the fundamental parameters involved in noninvasive diagnosis and medical sensors as well as to provide guidance for computational dosimetry and for the development of specific therapeutic approaches.
In recent decades, fractional calculus has been a fruitful field of research in science and engineering. It has found use in studies of viscoelastic materials and anomalous diffusion processes, as well as in many fields of science and engineering including signal processing, mathematical biology, electrochemistry, electrical networks, electromagnetic theory, and probability [11–15]. As well known, the concept of fractional exponents is an outgrowth of exponents with integer value. In the same way, nonintegral order of integration is a generalization of the mathematical operations of differentiation and integration to arbitrary, general, noninteger order. Although it better models the higher complexity by nature, it is still fairly easy to physically represent its meaning. However, just as fractional exponents may find their way into innumerable equations and applications, it will become apparent that integrations of fractional order can find practical use in many modern problems. Taking in consideration that the study of fractional calculus opens the mind to entirely new branches of thought, we applied such concept in a specific electromagnetic problem in order to demonstrate that the fractional derivatives operator can be an interesting and useful tools in electromagnetic theory.
Biological tissues are heterogeneous mixtures of various materials and quantities such as water, ions, membranes, and macromolecules with a broad variety of profiles. Therefore, their interaction with electromagnetic radiation takes place in different relaxation processes including (i) reorientation of water and protein-bound water molecules, (ii) interfacial polarization arising by the presence of two or more regions with different electrical properties, (iii) ionic diffusion, (iv) tumbling motion of the protein molecules, and (v) relaxations due to the nonspherical shape. These phenomena cause a frequency dispersion pattern of permittivity and conductivity which is characterized by a series of steps as the frequency increases [16, 17].
The Debye model has been widely used to describe the dielectric response arising by the dipolar relaxation [18, 19]. However, this simple exponential law cannot describe a wide class of relaxation processes occurring in biological materials since the intrinsic structural disorder and the ion-ion interaction induce anomalous dynamics based on nonexponential temporal decay. As a consequence, a number of empirical dispersion functions, such as Cole–Cole, Cole–Davidson, Havriliak–Negami, and Raicu spectral models, have been developed to accurately fit the fractional power-law decays of the dielectric response [17].
Finite-Difference Time-Domain (FDTD) method has been widely used in electromagnetic simulation due to its straightforward implementation and ability to model a broad range of initial and boundary values problems. Taking into account that the dielectric response includes fractional powers of , the design of the FDTD algorithm requires special treatments aimed at embedding the approximation of fractional derivatives into the simulator [20–23].
Besides the requirements linked to the frequency behavior of dispersive materials, a complete numerical code should be able to correctly incorporate, inside its kernel, the spatial dependence of the material parameters. This is a fundamental skill to meet especially when biological tissues have to be employed in the simulation. In fact, from the macroscopic point of view, a biological tissue can be treated as a random mixture of different medium compounds having a spatial-dependent concentration.
Considering these premises and the requirement to model a more general permittivity function, in this paper we illustrate a more accurate study of the time-domain electromagnetic field propagation in arbitrary dispersive media. The study is fully motivated to seek for an extended model flexibility enabling a better parametrization of the arbitrary dispersive media properties as well as a better fitting, over broad frequency ranges, of the experimental dielectric response. In particular, using a general fractional polynomial series approximation and the fractional calculus theory both the spatial and frequency dispersion characteristics of the dielectric response have been incorporated into the developed FDTD scheme. Dedicated uniaxial perfectly matched layer boundary conditions have been also derived and implemented in combination with the basic time-marching scheme. In particular, the iterative -field and -field updating while applying boundary conditions results in a marching-in-time procedures that simulate the continuous actual electromagnetic waves in a finite spatial region by sampled-data numerical analogs propagating in a computer data space [24]. Moreover, ohmic losses and total field/scattered field approach are taken into account. The accuracy of the proposed FDTD scheme has been tested by considering several test cases involving different spatial and time-dependent permittivity profiles.
2. Formulation of the Proposed FDTD Method
FDTD method is a powerful and efficient numerical technique used to solve both simple and complicated electromagnetic problems. This method has attracted much attention due to its simplicity, accuracy, robustness, low computational footprint, and its capability in treating in a straightforward and effective way many types of dispersive media. FDTD technique is based on the direct numerical integration of Maxwell’s two curl equations [25]. In dielectric materials, such equations are as follows:where is the electric flux density, is the electric field, is the magnetic field, is the magnetic permeability of free space, is the source current density, and is the static conductivity due to the ionic or electronic charge transport and modeling the ohmic losses. To calculate the spatial and temporal distributions of the electromagnetic fields in a medium, the constitutive relation between the fields and has to be combined with the equations (1). In frequency domain, the relative electric permittivity of a linear and nonmagnetic media can be expressed as [25, 26]where is the permittivity of free space, is the relative permittivity at high frequency limit, and is the electric susceptibility.
FDTD methods have been developed for incorporating Debye, Cole–Cole, Cole–Davidson, and Havriliak–Negami electric susceptibility models [27–30]. The common approach is to approximate the relative complex permittivity by means of rational or polynomial functions resulting in auxiliary differential equation FDTD models. When implemented to predict the electromagnetic field propagation inside quite simple dispersive media, such numerical models provide reliable results. On the other hand, their effectiveness falls when more complex and multirelaxed dispersive media are employed since the approximation does not accurately fit the frequency behavior of the electric permittivity. In order to overcome this drawback, we have successfully implemented in the FDTD formulation a more general electric susceptibility function having the following fractional polynomial dispersion response [31]:where is the relaxation strength with being the relative permittivity at low frequency limit, , is the relaxation time, and is a real coefficients related to the material parameters. The condition is recommended to meet the requirement on the numerical stability of the developed time-marching scheme. Moreover, implementing an optimization algorithm based on swarm intelligence the parameters and can be recovered to suitably fit a generic frequency domain permittivity function [32]. In particular, if is the maximum expansion order in (3) and is a given small positive threshold to be used for controlling the accuracy of the approximation, a dedicated numerical procedure based on the enhanced weighted quantum particle swarm optimization (EWQPSO) was used to solve the minimization problemwhere is the dispersion response to approximate. This technique has proved to feature superior effectiveness in terms of convergence rate and accuracy in comparison with alternative evolutionary stochastic search methods such as genetic algorithm [31, 32].
In linear dispersive materials, the frequency domain equation relating the electric and polarization fields isor by substituting (3),where .
Taking the inverse Fourier transform of (6),where Moreover, considering the gamma function, the following equation could be achieved (see Appendix A):where . Differentiating both sides of (7) with respect to time and using Leibniz’s rule for differentiation under the integral sign the following equation can be obtained:whereis the polarization current. By using the constitutive relationAmpere’s law in time domain can be written asApplying a second-order accurate finite-difference scheme, one readily obtains at the time instant where the vectors terms appearing on the right-hand side of (14) are evaluated using the semi-implicit approximation:Evaluating (10) at the following equation can be obtained (see Appendix B):Finally, by carrying out a second-order accurate finite-difference approximation of Faraday’s law in the time domain, the following update equation for the magnetic field is obtained:In (14) and (17) the curl of magnetic and electric fields are approximated by the standard Yee algorithm. In our computations, a uniform space-time grid is employed. Moreover, a flowchart summarizing the main steps of the proposed FDTD algorithm is shown in Figure 1, where and represent the time and space indexes, respectively.
The arbitrary biological tissue can be generally visualized as a heterogeneous material. As a result, its dielectric properties are governed by a spatial average which can be expressed through an appropriate spatial dependence of the material parameters. To this aim, both Maxwell–Garnet and Bruggeman models have been taken into account [33]. In particular, in (2)-(3) the electric permittivity at both low and high frequency limits is expressed asfor Maxwell–Garnet model, andfor Bruggeman model. In (18)-(19), is the permittivity of inclusions having ellipsoidal shape, is the permittivity of environment material, is the depolarization factor, and is the space-dependent filling factor function.
3. Numerical Results
To test the developed numerical code, the slab model illustrated in Figure 2 is implemented. It consists of a half space filled up with dispersive and nonhomogeneous biological material. The system is irradiated by a plane wave propagating along the positive -direction, with electric field linearly polarized along the -axis and located at a given point . In particular, the time-domain signal source iswhere GHz and GHz are selected in order to obtain a spectral bandwidth ranging from GHz to GHz.
The time and spatial extension of the computational domain are ns and cm, respectively, the thickness of the biological media is cm, and the time and spatial steps are ps and mm, respectively. Moreover, the computational domain is closed from both sides by dedicated uniaxial perfectly matched layers [21].
The relative electric permittivity considered in the simulations iswhere ps, and are evaluated according (18)-(19), and , , , and . Therefore, a linear and exponential filling factor functions are considered as follows:where is the initial value at the points and .
In order to test the effectiveness of the presented fractional power series approach, a benchmarking against the multiterm Debye approximation and Padé approximation has been carried out. In particular, the relative permittivity function (21) with and has been approximated as follows:To quantify the accuracy of the different approaches, the root-mean-square deviation on the complex electric permittivity has been evaluated aswhere the index identifies the specific type of approximation. As it can be noticed in Table 1, higher order of both Debye and Padé approximations results in a smaller fitting errors but at the same time the computational model becomes more complex. On the other hand, fractional power expansions lead to all-pole approximation of permittivity function. In this way, benefits in terms of memory usage and FDTD numerical implementation are achieved.
In order to validate the proposed approach, the numerical results given by the FDTD scheme are compared with those derived using a fully analytical frequency domain technique [34]. To this aim, the total field/scattered field (TF/SF) formulation, detailed in [24], has been implemented. TF/SF boundary separates the computational domain into a total field region, containing both the incident and scattered fields, and a scattered field region, containing only the scattered field. The TF/SF boundary can, in principle, be used to introduce any type of incident field into the FDTD grid, but, in practice, it is used almost exclusively to introduce plane wave excitation.
The simulations have been performed by considering a biological media described by (18), (21), and (23), where and and (needles). In particular, the reflectance, , and transmittance, , have been calculated, where the reflection, , and transmission, , coefficients have been evaluated as follows:where is the electric field reflected at the air-biological tissue interface, is the electric field excited at the location cm, with being inside the biological tissue, and is the electric field excited at the source location .
Figure 3(a) shows the space-time distribution of the electric field resulting from the developed FDTD-based procedure. It worthwhile to note the transmission process as well as the reflection phenomenon occurring at the air-biological tissue interface. Moreover, the pulse spreading due to the propagation inside the dispersive biological media as well as the spatial variation of the group velocity due to the spatial dependence of the biological tissue permittivity is evident. Figure 3(b) illustrates the reflectance and transmittance spectra as evaluated by using the proposed FDTD procedure and the analytical technique. The excellent agreement between them fully validates the developed FDTD formulation.
(a)
(b)
To further investigate the accuracy and flexibility of the proposed FDTD scheme, the pulse propagation in complex media involving several spatial-dependent permittivity profiles has been studied. The numerical simulations have been performed by changing the filling factor function, the depolarization factor, and the mixing model. In particular, besides the function profiles (22)-(23) and the Maxwell–Garnet and Bruggeman model, the depolarization factor sets (sphere), , (needle), and , (disc) are taken into account. In addition to the transmittance and reflectance, the penetration depth spectrum is calculated, too. This parameter is used to denote the depth at which the power density has decreased to 37 of its initial value at the surface.
The first test case regards the pulse propagation inside a dispersive medium described by the Maxwell–Garnet model, the linear filling factor function, and the depolarization factor . The calculated numerical results are illustrated in Figure 4. By an inspection of Figure 4(a) it can be observed that the reflectance increases by increasing the initial value . In fact, considering that both and are higher than and , rising generates higher permittivity values resulting in an improvement of the permittivity variation at the air-biological tissue interface. On the other hand, lower values of enable a better penetration of the electromagnetic field inside the biological media and as shown in Figure 4(b), the penetration depth increases. Moreover, the different relative variations of the transmittance and penetration depth at low and high frequency can be explained by considering that .
(a)
(b)
The second test case regards the pulse propagation inside a dispersive medium described by both the Maxwell–Garnet and Bruggeman models. The linear filling factor function, with , and the depolarization factor are considered to perform the simulations. The obtained reflectance and transmittance as well as the penetration depth are illustrated in Figures 5(a) and 5(b), respectively. Differences in the propagation characteristics can be noticed especially at low frequency. Such result highlights a critical point regarding the choice of the suitable model describing the inhomogeneous permittivity function. In fact, for denser composites Bruggeman’s formula is better suited than the Maxwell–Garnett one.
(a)
(b)
The third test case pertains the pulse propagation inside a dispersive medium described by both the linear and exponential filling factor functions having . The Maxwell–Garnet model and the depolarization factor are considered, too. Figures 6(a) and 6(b) depict the transmittance/reflectance and penetration depth spectra, respectively. Since is the same in linear and exponential functions, negligible reflectance variations occur. More significant transmittance and penetration depth changes can be observed since their trend mainly depends on the filling factor function.
(a)
(b)
The forth test case regards the pulse propagation inside a dispersive medium modeled by the depolarization factor sets (sphere), (needle), and (disc). Linear filling factor function having and the Maxwell–Garnet model are considered, too. The calculated reflectance/transmittance and penetration depth spectra are reported in Figures 7(a) and 7(b), respectively. By an inspection of these figures it is clear that the depolarization factors strongly affect the propagation properties of the dispersive media since the geometry and the surface area change the material response. In fact, the spherical shape creates the smallest dipole moment and when the spherical symmetry is strongly broken, as it happen for extremely squeezed ellipsoids (discs, needles), a quite strong deviation from the polarizability of the spherical shape occurs.
(a)
(b)
The final test case is the more general one since it pertains the electromagnetic analysis of a dispersive media characterized by a piecewise linear filling factor function:where cm and cm. In particular, the coefficient sets and are simulated. Moreover, the Maxwell–Garnet model, the permittivity values , , and , and the depolarization factor are also considered. Figure 8 shows the space-time distribution of the electric field as calculated by means of the developed FDTD technique. Figure 8(a) pertains the coefficient set which generates a growing piecewise function composed by two linear functions having different slopes. In this case, the resulting permittivity decreases from its maximum value at as the coordinate increases. Moreover, a quite high permittivity gap at the interface occurs. As a result, the incident plane wave is strongly reflected by the air-biological tissue interface and the speed of the transmitted wave increases as its penetration inside the biological media increases. Figure 8(b) pertains the coefficient set which generates a decreasing piecewise function. In this case, the permittivity of the biological medium at is well matched to the air permittivity and the incident plane wave is poorly reflected at the interface . Moreover, because of the permittivity increasing as a function of the coordinate, the transmitted wave slowly penetrates inside the dispersive media.
(a)
(b)
4. Conclusion
An accurate FDTD model for simulating the pulse propagation in arbitrary dispersive media is proposed. Using the fractional derivative operator, the fractional part resulting from a general series expansion of the permittivity function of dispersive materials is directly incorporated into the FDTD scheme, avoiding the use of auxiliary differential equations. Moreover, dedicated uniaxial perfectly matched layer boundary conditions are derived and total field/scattered field formulation is adopted in combination with the basic time-marching scheme. As a result, the newly proposed algorithm gets excellent accuracy with minimum auxiliary variables and thus minimum storage spaces. Its application in modeling the electromagnetic pulse propagation in several test cases is presented and simulation results demonstrate the reliability and efficiency of the proposed method.
Such a technique can be considered as a general scheme for modeling various types of dispersions including Debye, Cole–Cole, Cole–Davidson, Havriliak–Negami, Drude–Lorentz, and universal dielectric response. The presented FDTD scheme can be considered as a general-purpose tool useful to address complex electromagnetic problems in the field of bioengineering. It could be successfully applied in the field of remotely powered implantable devices for optimizing the trade-off between the received power and tissue absorption. It could be also useful to evaluate in detail the temperature elevation in the biological tissues as well as to perform an accurate investigation of the spatial and temporal distribution of the transmembrane voltage and electroporation phenomenon in cells exposed to high voltage electric pulses. Finally, it could be extended to perform an extensive numerical analysis of pulsed electromagnetic fields propagation in colloids, emulsions, and random dielectric materials.
Appendix
A. Fourier Transform of the Power-Law Function
For a real number , the gamma function is defined as The Fourier transform of the power-law function iswhere for , and using the variable change , it is straightforward to obtainThe inverse Fourier transform of (A.3) givesand comparing (A.4) with (8), (9) can be easily achieved.
B. Time-Marching Scheme
Evaluating (10) at yieldswhereButwhereRearranging (14) and using (15) and (B.1)–(B.4) givewhereRearranging (B.1) gives
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work has been partially supported within the “Aiuti a Sostegno dei Cluster Tecnologici Regionali” plan: LPIJ9P2 PERFOrm TECH-Puglia Emerging FOOD Technology.