Fractional calculus based FDTD modeling of layered biological media exposure to wideband electromagnetic pulses

: Electromagnetic ﬁelds are involved in several therapeutic and diagnostic applications such as hyperthermia and electroporation. For these applications, pulsed electric ﬁelds (PEFs) and transient phenomena are playing a key role for understanding the biological response due to the exposure to non-ionizing wideband pulses. To this end, the PEF propagation in the six-layered planar structure modeling the human head has been studied. The electromagnetic ﬁeld and the speciﬁc absorption rate (SAR) have been calculated through an accurate ﬁnite-difference time-domain (FDTD) dispersive modeling based on the fractional derivative operator. The temperature rise inside the tissues due to the electromagnetic ﬁeld exposure has been evaluated using both the non-thermoregulated and thermoregulated Gagge’s two-node models. Moreover, additional parametric studies have been carried out with the aim to investigate the thermal response by changing the amplitude and duration of the electric pulses.


Introduction
The modeling of electromagnetic field propagation through dispersive materials is a subject of increasing research activities since these dielectric media are found in a growing number of applications [1][2][3][4][5][6][7].In particular, biological tissues are an important class of dispersive media especially taking into account that their interaction with the electromagnetic fields has a great influence on the behavior of living systems.In fact, the increasing number of power and telecommunication systems leads to human exposure to non-ionizing radiation, causing an increased public concern about the potential health hazards [8][9][10].On the other hand, biological tissues are involved in a variety of therapeutic and diagnostic applications of electromagnetic fields such as hyperthermia, electroporation and the treatment of specific diseases [11][12][13][14][15].
Pulsed electric fields and transient phenomena are playing a key role in understanding the biological response to electromagnetic field exposure.As a result, the development of theoretical models and computational techniques describing the interaction between pulsed electric field (PEF) and dispersive materials is an invaluable research activity to increase the knowledge of the electromagnetic field distribution inside biological tissues, as well as to impose the development of specific therapeutic approaches.In fact, the efficiency of PEF therapies and treatments strongly depends on the electric field strength, the duration of the applied field, the shape of the electric pulse, the polarity and the number of intervals between pulses [12,13,[16][17][18].
Generally, accurate modeling of the nanosecond PEF propagation inside biological tissues needs to account for the effects of dielectric relaxation.As revealed by a number of experiments, the dielectric response of a broad variety of biological media cannot be described by a simple exponential expression with a single relaxation time [19].Therefore, exponential laws based on Debye-type dispersion or the combination of such dispersions cannot fully model the relaxation processes in such materials.As a consequence, a number of empirical dispersion functions (Cole-Cole, Cole-Davidson, Havriliak-Negami, Raicu) have been proposed to fit such types of dielectric spectra.
The finite-difference time-domain (FDTD) method has been widely used in electromagnetic modeling due to its straightforward implementation and its ability to manage a broad range of exposure conditions [20,21].Since Cole-Cole, Cole-Davidson, Havriliak-Negami and Raicu dispersion functions include fractional powers of angular frequency jω, mathematical models adequately representing this kind of response have to be incorporated into FDTD simulators.Recently, the authors have proposed a novel FDTD formulation based on the Riemann-Liouville theory of fractional differentiation where the fractional operator is directly approximated using finite differences [22][23][24][25][26].The reflectance, transmittance, absorbance and penetration depth were computed and compared with the corresponding analytical solution in a wideband frequency range.
Considering that common protocols for therapeutic applications of PEF involve the delivering of multiple pulses and taking into account that a high potential drop occurs across the skin, where the electric field is the largest, deleterious thermal damage due to Joule heating can occur.To this aim, the proposed FDTD scheme has been extended to incorporate further physics for evaluating the temperature distribution inside a layered heterogeneous biological system.In particular, a multiphysics model based on the fractional calculus-based FDTD algorithm, bioheat equation and thermoregulated Gagge's two-node model [27] has been developed.As a test case, a six-layered planar structure modeling the human head and consisting of air, skin, fat, bone, dura, cerebro-spinal fluid (CSF) and brain was taken into account [28].The computation of the electromagnetic field, power density and specific absorption rate (SAR), as well as the temperature rise inside the tissues when exposed to an incident plane wave was performed.Moreover, the temperature profiles with and without thermoregulation phenomenon have been calculated and compared.In particular, they have been evaluated by changing electric parameters such as the amplitude and duration of the electric pulses.

Complex Permittivity
The dielectric properties of biological tissues result from the interaction of electromagnetic energy with the tissue constituents at the cellular and molecular level.As a result, the dielectric permittivity and conductivity values vary from tissue to tissue and depend on frequency.In order to model realistic electromagnetic wave propagation over a broad frequency range, an exhaustive and reliable analytic model of dielectrics is needed.To this aim, the relative permittivity, ε r , of the biological tissues exhibiting an N-th order relaxation process has been modeled by using the general Havriliak-Negami (HN) relationship [29,30]: where τ i and ∆ε ri are the relaxation time and the amplitude change of the i-th relaxation process, respectively, N is the number of dielectric relaxation processes, r ∞ is the relative permittivity at ωτ → ∞, σ is the static ionic conductivity, 0 < α i , β i < 1 are heuristically-derived fitting parameters, ω = 2π f , and 0 is the free space permittivity.
To test the developed numerical code, a planar layered system is considered.It consists of skin, fat, bone, dura, CSF and brain, and a schematic diagram is shown in Figure 1.This structure well describes the human head locally and does not require extensive computational efforts.In fact, considering the occurrence of tissues having small thickness, the existing anatomical computer models can provide reliable results if huge computational resources are available.Despite its simplicity, it was proven that this simple planar model produces SAR values very close to those obtained by using more complex ones [31].Moreover, it has been demonstrated that the developed fractional derivative-based FDTD scheme allows an accurate space-time evaluation of electromagnetic field profiles in a broad frequency range [22][23][24][25][26].The reflection and transmission of a electromagnetic wave at the tissue interface depends on the frequency, permittivity and conductivity of the involved biological materials.Tissues are very inhomogeneous materials, and the fluid content of the material generates large differences in their electric properties.Dielectric data of several bulk tissues in a wide frequency range have been reported in the literature [32,33].Therefore, to estimate the dielectric properties in a desired frequency range, these experimental results have been fitted using (1) and minimizing the error function, E , by means of a dedicated numerical procedure based on the enhanced weighted quantum particle swarm optimization (EWQPSO) [34]:

CSF
In Equation (2), δ is the maximum tolerable error, exp r is the measured permittivity and r represents the general HN dielectric response.
By using (2) in the frequency bandwidth ranging from f min = 100 MHz to f max = 10 GHz and considering N = 2, the HN parameters recovered for each medium are listed in Table 1.

Electromagnetic Analysis
In order to evaluate the electromagnetic field propagation inside the layered structure involving HN dielectric materials, the FDTD scheme proposed in [22][23][24][25][26] has been used.In particular, it implements a more general series representation of the Riemann-Liouville fractional derivative operator, and it takes into account multiple relaxation times, as well as ohmic losses occurring in common biological media displaying different dispersion mechanisms.Dedicated uniaxial perfectly matched layer boundary conditions have been also derived and implemented in combination with the basic time-marching scheme.Moreover, the conventional total field/scattered field approach and the Courant stability criterion are considered.Thus, applying a second order accurate finite-difference scheme and the procedure detailed in [22][23][24][25][26] for the finite-difference discretization, the FDTD update equations for the electric field, E, the magnetic field H and the auxiliary displacement current density, J i , can be calculated within the HN media and UPML regions.In particular, by considering the standard Cartesian Yee cell grid for the FDTD scheme and applying the central difference discretization, a second order O (∆z) 2 truncation error is achieved.

Bioheat Equation with Thermoregulation
It is well known that the interaction of an electromagnetic field with living biological tissues generates an increase of the local temperature [31].For this reason, it is important to know the amplitude of the electromagnetic wave propagating inside the considered biological structure to evaluate the temperature rise.The electromagnetic power transferred inside the biological tissue is related to the SAR, and it is defined as follows: where σ m (S/m) and ρ m (kg/m 3 ) are the m-th tissue conductivity and density, respectively, and E is the electric field at the generic z coordinate.The temperature distribution inside each tissue exposed to radio frequency (RF) power can be evaluated solving the Pennes bioheat equation.It is based on the classical Fourier law of heat conduction, and it incorporates the effects of metabolism and blood perfusion on the energy balance within tissue in terms of distributed heat sinks or sources.The generalized 1D Pennes equation for the m-th tissue can be written as [35]: where T is the tissue temperature, T art is the time-dependent temperature of the arterial blood, ρ b = 937 kg/m 3 is the blood density, C b = 3889 J/kgK is the blood specific heat, W b,m (1/s) is the blood perfusion rate, q m (W/m 3 ) is the metabolic heat generation rate, k m (W/mK) and C m (J/kgK) are the thermal conductivity and the specific heat of the m-th tissue, respectively.The Pennes model does not take into account the thermoregulation mechanisms of the human body such as shivering, regulatory sweating and vasomotion affecting the thermal response of the human body.In order to overcome this drawback and to calculate a more accurate temperature field, a thermoregulated bioheat model combining Pennes' equation and Gagge's two-node model has been coupled to the fractional calculus-based FDTD code.In this approach, the layered structure, called the outer layer, is sandwiched between the inner layer, called the core, and the surrounding air (see Figure 1).
A part of the heat produced in the core by metabolic and perfusion rates (when it is treated as the source) is lost via conduction and mass transfer to the outer layer.A portion of this energy may be expended as external work done by the muscles and as dispersion at the skin/air interface.The rest of the core energy is stored, and it causes temperature to rise.All these phenomena are controlled by hypothalamus via a complex thermoregulation mechanism.In particular, the warm and cold receptors located in the core and outer layers send signals to the hypothalamus which, depending on the integration of such stimuli, suitably activates the vasomotion, shivering and regulatory sweating (see Figure 2).By using this concept, Gagge's two-node mathematical model is based on the energy balance equation involving the outer and inner layers.In particular, the energy balance equation for the core is [27]: where α by is the fraction of the body mass concentrated in the skin, ρ by = 985 kg/m 3 is the body density, C by = 3490 J/kgK is the specific heat of the body [36], T cr is the core temperature, r m is the remaining metabolic coefficient, K eff = 5.28 W/m 2 K is the effective conductance between the core and outer layer, ṁb (kg/m 2 s) is the rate of blood flow due to the vasomotion of the blood vessels caused by the cold/warm thermal conditions and l by = V by /A D (m) is the characteristic length of the body where V by (m 3 ) is the volume of the human body and A D (m 2 ) is the nude body surface area given by the DuBois formula [37]: where m and l are the body mass and height, respectively.Moreover, taking into account that the external mechanical efficiency is insignificant in many human thermal response applications, the remaining metabolic coefficient can be defined as [38]: where T A ( • C) is the air temperature and P A (kPa) is the water vapor pressure in the air given by [39]: Finally, the parameter α by can be defined as [40]: Discretizing ( 5) gives: which is used as thermal boundary condition for the core in Pennes' equation.The thermal boundary condition at the outer layer surface z * = z sk , z sk being the z coordinate of the skin/air interface, is:

Air
convective heat loss z=z * , where k sk = 0.42 W/mK is the thermal conductivity of the skin, h = 10.5 W/m 2 K is the convective heat transfer coefficient at the air temperature, sk = 0.98 is the skin emissivity, σ s = 5.67 × 10 −8 W/m 2 K 4 is the Stefan-Boltzmann constant and w sk is the total skin wettedness.Moreover, considering that only a part of the total skin surface, d SAR , is irradiated by the electromagnetic field, the skin temperature T sk is computed using the equation: where T (z, t) z=z * is the temperature evaluated when the electromagnetic field source is applied (E = 0) and T noSAR (z, t) z=z * is the temperature for the part of the body not irradiated by the electromagnetic field (E = 0).In Equations ( 10) and (11), some physiological parameters such as q m , ṁb , w sk , α by depend on the thermoregulation mechanisms.The metabolic heat production mainly depends on the physical activity of the human body (q m,act ), as well as it can be increased by shivering against cold (q m,shiv ) as follows [38]: where: In (14), CSIG sk and CSIG cr are the cold signals of the human body for the outer layer and core compartments, respectively, given by [38]: and: where T sk,n = 33.7 • C and T cr,n = 36.8• C are the neutral outer layer and core temperature, respectively.Vasomotion of the blood vessels changes both parameters ṁb and α by .In particular, ṁb can be calculated using the equation [40]: where WSIG cr is the warm signal for the core compartments, and it is given by: The regulatory sweating increases the skin wettedness.In particular, the total skin wettedness is due to both diffusion through the skin (w dif ) and regulatory sweating (w rsw ) [27]: where: and: In Equation ( 21), h fg = 2260 × 10 3 J/kg is the heat of water vaporization, ṁrsw ( kg/m 2 s) is the rate of the sweat production per unit of skin area and q max evap is the maximum evaporative potential.The parameter q max evap can be estimated by using the equation [40]: where P sk is the water vapor pressure in the saturated air at the skin temperature given by: and R et = 0.017 m 2 kPa/W is the total evaporative resistance between the body and the air.
The parameter ṁrsw can be estimated by the following equation [40]: ṁrsw = 4.7 × 10 −5 WSIG by exp{WSIG sk /10.7} (24) where: is the warm signal for the outer layer compartments.In (24), WSIG by is the warm signal of the body given by: WSIG by = max 0, T by − T by,n .
Moreover, the body temperature, T by ( • C), and the neutral body temperature, T by,n ( • C), are given by: T by,n = α by T sk,n + (1 − α by )T cr,n .

Numerical Results
In order to validate the fractional calculus-based FDTD model, the electromagnetic field propagation in the six-layered dielectric slab shown in Figure 1 and over a frequency band spanning from f min = 100 MHz to f max = 10 GHz has been simulated.In the simulation, the geometrical and dielectric parameters listed in Table 1 have been considered.The system is irradiated by a plane wave propagating along the z-direction with the electric field linearly polarized along the x-axis.In particular, the excitation signal is a sinusoidally time-modulated Gaussian pulse having the same parameters reported in [22].The chosen time and spatial steps are ∆t = 0.3 ps and ∆z = 0.1 mm, respectively.Moreover, the time and spatial computational domain are t tot = 10 ns and z tot = 6 cm, respectively.The stability of the proposed FDTD numerical scheme has been tested by following von Neumann's spectral approach and assuming the harmonic time dependence of the electromagnetic field quantities.At any time step, the instantaneous distributions of the E, H and J i fields can be Fourier transformed with respect to the space variable z to provide the spectrum of monochromatic plane-wave modes propagating along the computational lattice.In this way, it has been checked that the resulting characteristic polynomial equation has zeros inside the stability circle [20,24].
The numerical results given by the FDTD scheme are compared with an analytical frequency-domain technique.In particular, the reflectance |R| 2 and transmittance |T | 2 have been calculated.The transmission coefficient is evaluated as the ratio of the Fourier transform of the electric field computed in two different points, z and z , across the biological structure, where E is the amplitude of the electric field.The frequency-domain behavior of reflectance and transmittance calculated using the developed FDTD and a rigorous fully-analytical technique based on the transfer matrix approach for layered lossy media [41] is illustrated in Figure 3.The excellent agreement with the analytical technique validates the proposed methodology.
The bioheat Equation ( 4) and the boundary conditions ( 5) and (11) were discretized using the standard finite difference method in both the spatial and temporal domain.Considering that the temporal evolution of the heat is slower than the electromagnetic field one, the time step used to determine the SAR is not suitable for the heat transfer simulation.To this aim, a new time step has been assigned to solve the bioheat equation.More specifically, to fulfill the constrains due to the numerical stability and accuracy requirements, it has been chosen according to the von Neumann's condition: where C (J/kgK) is the specific heat, ρ (kg/m 3 ) is the medium density and k (W/mK) is the thermal conductivity.In particular, the spatial and temporal step size considered in the bioheat simulations are 0.1 mm and 3.5 ms, respectively, as well as the blood temperature is assumed to be 37 • C.Moreover, the thermal properties of the involved biological tissues are listed in in Table 2 [31].The temperature rise, ∆T, due to the PEF interaction with the biological tissues is calculated as: where T noSAR is the steady-state temperature distribution calculated solving the bioheat equation with SAR = 0 and assuming that the core and skin temperatures are 36.8• C and 33.7 • C, respectively, and T SAR is the temperature calculated with SAR = 0 and using the model including the thermoregulation mechanisms.Moreover, in order to quantify the thermoregulation effect, the temperature rise ∆T 1 has been calculated as: where T noreg is the temperature evaluated neglecting the thermoregulation mechanisms.The SAR distribution inside the stratified system has been evaluated using (3), where the electric field is computed solving the Maxwell equations.The proposed model has been validated by comparison with the literature results [42].In particular, in Table 3 is reported the maximum temperature rise inside the whole head calculated by considering an electromagnetic source with the power equal to 0.6 W at 900 MHz and 0.27 W at 1500 MHz.The good agreement between the two models highlights the numerical algorithm's ability to model the biological media exposure to electromagnetic pulses well.Table 3. Thermal increase comparison inside the brain layer by considering the thermal model published in [42] and the proposed thermoregulated model.

Model
Thermal model in [42] 0.160 K 0.132 K Proposed thermoregulated model 0.156 K 0.117 K In the following investigations, a PEF burst characterized by a hold time t hold = 1 ns, rise time t r = t hold /5, repetition time t prt = 10 ns and maximum amplitude E max has been considered as the source.Figure 4a shows the time evolution of the PEF source by fixing E max = 1600 V/m, and Figure 4b illustrates the corresponding spatial distribution of the SAR.The different dielectric losses characterizing the biological tissues result in a non-homogeneous SAR profile.In particular, an higher value of SAR can be noticed inside the biological tissues exhibiting higher losses.Figure 5a illustrates the increase in temperature ∆T versus the time in the skin, bone and brain layers.In the simulations, a PEF exposure of 30 min followed by 30 min without exposure have been used.Moreover, the thermoregulated model and the external temperature equal to 23 • C are considered.By an inspection of the figure, it can be inferred that, with respect to the skin and bone layers, the brain is more protected by the temperature increase.In fact, the temperature grows more than 1.3 • C in the outer layer of the head and around 0.7 • C inside the brain.The temperature rise in each biological medium is directly connected to the corresponding thermal properties such as the thermal conductivity and blood perfusion rate.In particular, higher thermal conductivity enhances the thermal diffusion resulting in a lowering of the mean temperature.Meanwhile, higher blood perfusion rate increases the heat exchange with the blood ensuring a more stable temperature around the steady-state value.On the basis of these remarks, it is possible to justify the highest increase in temperature occurring inside the bone layer since it exhibits the lowest values of both the thermal conductivity and blood perfusion rate.Instead, the higher values of such parameters result in a more limited temperature rise in the brain layer.On the other hand, the temperature improvement due to the PEF exposure is mitigated by the thermoregulation process, which reduces the metabolic heating and improves, via vasodilatation, the blood flow, as well as it enables the production of sweat.Moreover, when the sources is turned off, the reduced metabolic activity and the starting of the sweating process allow a temporary temperature reduction below the steady-state value.

t (ns)
Figure 5b highlights the temperature changes due to the model including the thermoregulation mechanisms.It can be noticed that when the PEF signal is turned on, ∆T 1 ≥ 0 in each biological layer.This means that the temperature rise calculated using the thermoregulated model is lower than that calculated using the non-thermoregulated one.This occurrence can be explained by considering that the thermoregulation phenomenon modifies the blood perfusion and the metabolic heating, as well as it activates the sweating process.When the PEF signal is turned off, the sweating process and the reduced metabolic heating generate a faster temperature decrease compared to the one calculated using the non-thermoregulated model.As a result, a temporary negative value of ∆T 1 occurs.In response to such a negative value, the thermoregulation phenomenon should seek to stabilize the temperature as much as possible closer to the steady-state value.However, the noticeable discordance is a compelling argument for necessarily using the more complex model employing the thermoregulation mechanisms.In order to gain insights into the interaction of PEF with the multilayered human head, additional parametric studies have been carried out.Two different test cases have been considered.In the first one, the PEF amplitude is changed, and the hold time is equal to t hold = 1 ns.The second test case pertains to the parametric analysis performed by varying the hold time and fixing E max = 1600 V/m. Figure 6 displays the spatial distribution of the SAR concerning (a) the first and (b) the second test case.As it would be expected, the electromagnetic power enhancement due to the PEF amplitude rise generates higher SAR in each biological layers (see Figure 6a).The SAR also increases for longer pulse durations (see Figure 6b).In the latter case, the SAR improvement depends on the spectral band associated with the pulsed signal as well as on the penetration depth of the PEF inside the head.For shorter PEF signals, the SAR value is higher in the outer layer because the electromagnetic wave is strongly attenuated by the skin.Instead, by increasing the pulse duration, the electromagnetic wave propagates deeper, generating higher values of SAR inside the CSF layer, which shields the brain from the electromagnetic radiation.
Figure 7 shows the temperature versus the time for three different PEF amplitudes (a) at the air-skin boundary and (b) at the outer surface of the brain layer.In the simulations, the SAR profiles illustrated in Figure 6a have been used, as well as the source is active for the first 30 min and turned off for the following 30 min.The temperature improvement is directly correlated with the amplitude of the PEF source.In fact, for E max = 800 V/m, a negligible temperature rise can be observed.On the other hand, for E max = 3200 V/m, a strong thermal effect occurs inside the whole multilayered structure.Such behavior is due to the electromagnetic power associated with the pulse propagating inside the biological system.The power lost inside the biological tissues increases as the PEF amplitude increases, producing a higher heat load.In particular, the temperature increases more than 5 • C on the skin surface and around 3.5 • C inside the brain, reaching values that would result in damage to the biological tissues.Figure 8 shows the temperature versus the time for three different PEF durations (a) at the air-skin boundary and (b) at the outer surface of the brain layer.The SAR profiles illustrated in Figure 6b have been used, as well.It is clear that the temperature increases by increasing the pulse duration.In particular, with respect to the steady state value, the temperature improvement at the air-skin boundary changes from 0.7 • C, for t hold = 0.5 ns, to 2.8 • C, for t hold = 2 ns.However, at the outer surface of the brain, a more limited temperature change occurs ranging from 0.2 • C, for t hold = 0.5 ns to 2 • C, for t hold = 2 ns.The temperature rise exhibits a non-linear dependance respect to the PEF duration.In fact, it is higher for a longer pulse duration since the electromagnetic field propagates deeper.As a result, a more efficient heating occurs.For this reason, the temperature at the air-skin surface linearly increases for longer PEF durations, and a more evident changes inside the brain can be observed.In fact, for shorter pulses, the high frequency spectral components exhibit a strong attenuation before reaching the brain, thus producing a lower temperature rise.On the other hand, by extending the pulse duration, the spectral energy is mainly associated with lower frequency spectral components, which can reach the brain layer.

Conclusions
In this paper, a novel mathematical model based on an electromagnetic FDTD fractional scheme coupled with the bioheat equations and thermoregulated Gagge's two-node model has been illustrated.In the presented model, the fractional derivatives appearing in the permittivity function characterizing the dispersive materials are directly incorporated into the FDTD scheme by using the Riemann-Liouville definition.The developed numerical tool has been applied to the analysis of pulse wave propagation in a multilayered human head structure.After a validation carried out by comparison with a fully-analytical approach, the numerical code has been applied to study various test cases.The obtained results show that the presented multiphysics model provides the space-time temperature distribution inside the biological layers forming the human head when it is exposed to PEF radiation.Moreover, in order to investigate the effects of the amplitude and pulse duration on the temperature rise inside the head, several parametric simulations have been carried out.The obtained results demonstrate that the proposed numerical model could be a useful tool to study the coupled thermal and electromagnetic problems occurring in biological tissues.

Figure 1 .
Figure 1.Sketch of the layered structure considered in the simulations.

Figure 2 .
Figure 2. Sketch illustrating the thermoreceptor signals and the thermoregulation mechanisms involved in Gagge's two-node model.

FrequencyFigure 3 .
Figure 3. Reflectance and transmittance spectrum of the six-layered dielectric slab shown in Figure 1.

Figure 4 .
Figure 4. (a) PEF source having E max = 1600 V/m and t hold = 1 ns; (b) the corresponding SAR profile inside the multilayered head structure.

Figure 5 .
Figure 5. (a) Time-dependent temperature rise in the skin, bone, and brain layers; (b) comparison between non-thermoregulated and thermoregulated models.

Figure 6 .Figure 7 .
Figure 6.SAR profile obtained by considering a PEF source with (a) three different pulse amplitudes; (b) three different hold times.

Figure 8 .
Figure 8. Temperature versus time for three different PEF durations; (a) at the air-skin boundary and (b) at the outer surface of the brain layer.

Table 1 .
Dielectric and geometrical parameters of the simulated layered structure.

Table 2 .
Bioheat parameters of the layered media.