A High-Order, Symplectic, Finite-Difference Time-Domain Scheme for Bioelectromagnetic Applications within the Mother/Fetus Model

An explicit high-order, symplectic, finite-difference time-domain (SFDTD) scheme is applied to a bioelectromagnetic simulation using a simple model of a pregnant woman and her fetus. Compared to the traditional FDTD scheme, this scheme maintains the inherent nature of the Hamilton system and ensures energy conservation numerically and a high precision. The SFDTD scheme is used to predict the specific absorption rate (SAR) for a simple model of a pregnant female woman (month 9) using radio frequency (RF) fields from 1.5 T and 3 T MRI systems (operating at approximately 64 and 128 MHz, respectively). The results suggest that by using a plasma protective layer under the 1.5 T MRI system, the SAR values for the pregnant woman and her fetus are significantly reduced. Additionally, for a 90 degree plasma protective layer, the SAR values are approximately equal to the 120 degree layer and the 180 degree layer, and it is reduced relative to the 60 degree layer. This proves that using a 90 degree plasma protective layer is the most effective and economical angle to use.


Introduction
Because the protection of a pregnant woman and her fetus is very important, the algorithm that is used in mother/fetus modeling and simulation must have a high precision and must be numerically stable.
The Finite Different Method (FDM), Finite Element Method (FEM) and Finite Volume Method (FVM) can be used for calculating the SAR values of the mother/ fetus model. The FEM is more precise caused by the possibility to use high-order approximations. And since the methods for the high-order calculations of the FVM are developed by various groups recently [1,2], the statement that highorder calculations are not available for the FVM does not hold true anymore at present time. Nevertheless, the FDM is better for processing the simplified geometries, the FEM and FVM are not necessary for the calculation of the mother/ fetus model. Furthermore, the FDM can give simpler linear equations systems which could be solved faster. Also the finite-difference time-domain (FDTD) scheme is one of the FDM based on the time-domain. It can cover a wide frequency range with a single simulation and treat nonlinear material properties in a natural way. So the FDTD scheme is considered for our mother/fetus modeling and simulation.
To this end, our proposed high-order, symplectic, finite-difference timedomain (SFDTD) scheme is compared to the traditional FDTD scheme and the high-order FDTD scheme.
The traditional FDTD scheme has two shortcomings. First, it cannot accurately model curved, complex surfaces and material discontinuities by using the staircasing approach with structured grids. Second, this scheme accumulates significant errors from numerical instability, dispersion and anisotropies.
To reduce the numerical dispersion in the traditional FDTD scheme, a variety of high-order spatial discretization strategies have been proposed [3,4]. Based on orthonormal Harr wavelet expansions, a multi-resolution time-domain (MRTD) method [5,6] was proposed. Another approach is the high-order FDTD scheme, which uses the Fang(2,4) [7] and Ty (2,4) [8] formats, retains the simplicity of the original Yee algorithm and can conserve computational resources by using coarse grids compared to the traditional FDTD scheme. In addition, the discrete singular convolution [9,10] scheme was proposed, which uses delta cores, such as the Shannon core, Poisson core and the Lagrange core. However, these high-order approaches must use lower Courant-Friedrichs-Levy (CFL) numbers to comply with the stability criterion that subsequently destroys the symplectic structure of the electromagnetic system. Because Maxwell's equations can be written as an infinite-dimension Hamiltonian system, a stable and accurate solution can be obtained by using the symplectic integrator [11,12], which conserves energy in the Hamiltonian system. Although the SFDTD scheme has been used to solve the guided-wave [12], electromagnetic radiation, penetration and scattering problems [13], little research has been performed on bioelectromagnetic simulations. Additionally, the traditional FDTD scheme has been used to calculate SARs using the pregnant woman/fetus model [14,15], but no one has used the SFDTD algorithm to this end. In this article, we apply the high-order SFDTD scheme to a bioelectromagnetic simulation using a simple model of a pregnant woman and her fetus.
In addition, many researchers [14,16] model SARs, which are limited below the safety guidelines to ensure the safety of the patient, in the human body. A standard [17] has been developed to limit the maximum energy deposition within human subjects undergoing an MRI scan. However, few researchers are concerned with protective measures to reduce the risks posed by RF radiation. For lower SAR distributions at 64 MHz, we choose a 2 cm-deep plasma protective layer to reduce the mother's/fetus's SARs from a 1.5 T MRI system. By simulating the mother/ fetus model with the SFDTD scheme used, we find an optimal angle that not only best protects the patient, but it also significantly reduces the raw material costs.

Materials and Methods
Electromagnetic-Field-Solver General Formulations for the SFDTD Scheme A helicity Hamiltonian [18,19] for Maxwell's equations in a homogeneous, lossless, and sourceless medium can be introduced as Where H and E are the magnetic and electric fields. e and m are the permittivity and permeability of the medium. + is Curl, which is a vector operator that describes the infinitesimal rotation of a 3-dimensional vector field.
Based on the variational principle, (2) can be rewritten as Where 0 f g 3|3 is the 3|3 null matrix and R is the 3|3 matrix representing the three-dimensional curl operator. Now, Maxwell's equations with the FDTD scheme applied can be written as an infinite-dimension Hamiltonian system. A stable and accurate solution can be obtained by using the symplectic scheme.
If x~(x 1 , ::: x n , j 1 , ::: j n ),y~(y 1 , ::: y n , j 1 , ::: j n ) [ R 2n , the symplectic inner then v(wx,wy)~v(x,y). So, w is called a symplectic transformation or a canonical transformation. In other words, the symplectic scheme is just the transformation, in which the symplectic inner products are changeless. For the temporal direction, the electromagnetic field solution of formula (3) from t~0 to t~Dt is expressed by using the exponential operator as: Where exp Dt UzV ð Þ ð Þis the time evolution matrix of Maxwell's equations Using the product of elementary symplectic mappings, the exact solution (6) can be approximately constructed as [20]: Where c l and d l are the constant coefficients of the symplectic integrator. m m §1 ð Þ is the stages of the approximation and r is the order of the approximation. The coefficients can be found by using the Baker-Campbell-Hausdorff (BCH) formula [12,20,21,22].
For the operator U and V containing the curl operator R, Maxwell's equations must be discreted in the spatial direction by the higher-order difference for the numerical solutions.
For the spatial direction, the explicit, fourth-order-accurate difference expressions in conjunction with the staggered Yee lattice are used to discretize the first-order spatial derivatives as follows: Where the coefficients 9 8 and { 1 24 are derivated from the Tayor expanded formula.
So, the relationship between variable E and H is derived. As each tissues in the mother/fetus are the dispersive medium, the formula (1) is changed as: Where D is the electric displacement vector, e 0 is the permittivity of vacuum, v~2pf , f is the frequency of the incident wave and e Ã r v ð Þ is the complex permittivity of each tissues in the mother/fetus. : D are introduced into the formula (9) as: Where m 0 is the permeability of vacuum For the complex permittivity of each tissues, e Ã r v ð Þ is expressed as Where e r is the averaged relative permittivity of each tissues, s is the conductivity of each tissues.
Substituting the formula (11) into the formula (10): So, as the integration algorithm is substituted, the formula (13) can be changed with partial summation as: E n can be calculated from the formula (14) as With I n is determined by I n~s Dt e 0 X n i~0 E i , the formula (15) becomes: So, the relationship between variable D and E is derived. From the two-dimensional Maxwell's equations, the SFDTD scheme, which is fourth-order-accurate in space and fifth-stage-accurate in time, can be obtained by using the discretization approaches above. When m~5, a detailed expression of all the components of the scaled electric field from the s-th stage to (s{1)-th stage can be written as [12]: H nzs=5 gbz(i,j)~s : Dt=e 0 ð22bÞ Where e~e r e 0 and e r denotes the averaged relative permittivity. Dt is the timesteps and D is the length of space-steps. CFL is Courant-Friedrichs-Lewy condition, which is a necessary condition for stability while solving certain partial differential equations numerically by the method of finite differences. Time-steps and space-steps must satisfy the CFL condition. The coefficients s and m obey the symmetry relations c s~cm{sz1 1ƒsƒm ð Þ , d s~dm{s 1ƒsƒm{1 ð Þand d m~0 . D z can be calculated form the formula (17) and I z can be calculated form the formula (19). Then D z and I z are substituted into the formula (18) to get E z . Finally, E z is substituted into the formula (20) and the formula (21) to get H x and H y to complete the whole solution of Maxwell's equations.

Mother/Fetus Model (Month 9)
A significant amount of research has shown that only a small amount of the body's tissues (such as the brain) contains magnetic material. In addition, most of the body's tissues are made of non-magnetic material. The magnetic permeability m of these tissues is almost 1 and does not need to be considered. For B~mH, where B is the magnetic flux density and m<1, the relationship between B and H is linear. So only the electric properties of D is needed to be considered. Thus, we only address the dielectric properties of the tissues listed in Table 1. To study protective measures for a pregnant woman and her fetus, we developed a simple computational model of a mother and her fetus using the SFDTD scheme because of its high precision and numerical stability. The model used in this article consists of three different types of tissues shown in Table 1: the ''uterus'', ''placenta'' and the ''fetus''.
In Fig. 1, from the numerical model, the radius of the uterus, placenta and fetus are set to r uterus~1 1:28cm, r placenta~1 1:08cm, and r fetus~1 0:08cm to coincide with their physical dimensions as closely as possible when the fetus is nine months old.
The SFDTD scheme is applied to Maxwell's equations to calculate the SAR distributions within our model. We establish a computational domain of 450|450 cells for this numerical model. Each mother/fetus model consists of a lattice of cubic cells with a length of space-steps D~0:0006 m, which is sufficient to provide acceptable accuracy for the calculation of the SARs for the fetus. D~0:0006 m also is the maximum length of the grid, which is used to split the mother/fetus model into a lot of small calculation area. For the D is smaller, the 2r Where E x and E ŷ are the peak values of the electric field components and s and r denote the conductivity and the mass density of each tissues, respectively. Because this mother/fetus model is just two-dimensional now, the z field of E is not to be considered. In the future, the three-dimensional SFDTD scheme will be expanded to the three-dimensional electromagnetic modeling and simulation of the mother/fetus. A plane wave is incident perpendicular to the protective layer, with frequencies of approximately 64 and 128 MHz to simulate the RF fields of 1.5 T and 3 T MRI systems, respectively.
To ensure the safety of the mother and the fetus, we can add a protective layer where appropriate to reduce the body's SAR distributions [27,28]. In our realistic model, we choose a plasma material that has a thickness of 2 cm, plasma frequency of 3 GHz and collision frequency of 10 GHz as the protective layer. As shown in Fig. 1, the angle of the plasma protective layer is 60 degree. In addition, the angle may be 0 degree (no protection), 90 degree, 120 degree or 180 degree.
The parameter values of the numerical model can be modified according to different requirements. By using 800 time-steps iterations, with an incident 64 MHz plane wave and a protective layer of 60 degree, the results obtained for the E z electric field are shown in Fig. 2

Results and Discussions
Demonstrating the SFDTD scheme's benefits A one-dimensional hard source that uses a Gaussian pulse can be given as.
A long-term simulation is performed with CFL~0:48 and the length of one space-step D~1 cm, for which is the 10th part of a wavelength. One time-step is D 6|10 10 s. By using the perfect electric conductor (PEC) boundary condition, an one-dimensional resonant cavity is constructed. In Fig. 3, the propagation of the one-dimensional Gaussian pulse is simulated with the analytical solution, the traditional FDTD scheme and the SFDTD scheme respectively at the 436 th and  3146 th time-steps iteration. Note that Fig. 3 is the one-dimensional solution to prove the SFDTD scheme's advantage, because only the one-dimensional solution has the analytic solutions to be compared. And the analytic solutions are the curves of the ''PLUSE'' in Fig. 3. With the number of the time-steps iterations is larger form 436 to 3146, the propagation distance of the pulse is longer.
In Fig. 3, all the algorithms give acceptable results for the pulse calculation with the same waveforms after 436 time-steps iterations. However, at higher iterations, the waveform from the traditional FDTD scheme is greatly distorted after 3146 time-steps iterations with the curve of the traditional FDTD method has been shocked before about 50 relative cells. Where ''shocked'' means that the curve is not stable and the traditional FDTD method cannot preserve the constant energy of the Hamiltonian system. And relative cells are the relative number of grid. 120 relative cells in Fig. 3 are the size of the computational domain in this simulation.
In addition, the energy of the electromagnetic system computed by the traditional FDTD scheme gradually becomes attenuated. In contrast, under the long-term simulation at the 3146 th iteration, the conservation of energy by the SFDTD scheme is verified, as the curve is smooth and in good agreement with the curve of the Gaussian pulse. The peak value of the curve by the SFDTD scheme is always greater than 1 V/m at the 436 th iteration or the 3146 th iteration. Nevertheless, the peak value of the curve has been less than 1 V/m at the 3146 th iteration by the traditional FDTD scheme.
Above all, this demonstrates that the SFDTD scheme can remain stable and accurate by using the symplectic integrator, which preserves the constant energy of the Hamiltonian system. Therefore, the SFDTD scheme can meet the requirements for modeling the mother/fetus model and can calculate the local SARs below with a high precision and in a numerically stable manner.
The Local SARs using the SFDTD scheme at approximately 64 MHz and 128 MHz with a 60 degree protective layer The SAR distributions found on the surface of the mother/fetus model are shown in Fig. 4. Different SAR patterns are observed for 64 and 128 MHz. This is mainly caused by the MRI operating frequency as well as the conductivities used in the mother/fetus models. Our simulation results suggest that higher local SARs are found at 128 MHz rather than at 64 MHz. At 128 MHz, the peak SAR is greater than 5:10|10 {3 (W/kg), a fourfold increase over the SAR of 1:2|10 {3 (W/kg) at 64 MHz. In these cases, we choose the 1.5 T MRI system operating at approximately 64 MHz to perform a scan on a pregnant woman because this is safer as a result of the lower SAR distributions.
The Local SARs using the SFDTD scheme, under the 64 MHz controlled mode, with different plasma protective layer angles In Fig. 5, in comparison to the case of no protective layer, the peak SAR is significantly reduced using a 60 degree plasma protective layer. The maximum value of the SARx is reduced from 1:96|10 {3 (W/kg) to 1:75|10 {3 (W/kg), and the maximum value of the SARy is reduced from 1:76|10 {3 (W/kg) to 1:55|10 {3 (W/kg) with a 60 degree plasma protection with a reduction of 10 to 12%. It illustrates that the plasma protective layer plays an important role in reducing the electromagnetic radiation received by the mother and fetus.  Additionally, compared to using a 90 degree plasma protection layer, the peak SAR is nearly equal to the 120 degree layer and the 180 degree layer, and it is lower at 60 degree. By increasing the angle of the plasma protective layer from 60 degree to 90 degree, the maximum value of SARx is reduced from 1:75|10 {3 (W/kg) to 1:69|10 {3 (W/kg), and the maximum value of SARy is reduced from 1:55|10 {3 (W/kg) to 1:42|10 {3 (W/kg) with a reduction of 3.4 to 8.4 per cent. This proves that using a 90 degree plasma protection layer is most effective and economical for the radiation protection of a pregnant woman and her fetus.

Conclusion
This article derives the SFDTD differential equations and demonstrates the superiority of the SFDTD scheme for researching the mother/fetus model. The SFDTD scheme can make the results more stable and accurate with high precision and meet the requirements for calculating the local SARs of the mother/fetus model. Also, with the SFDTD scheme used, a better frequency of the MRI system and an optimal angle of the plasma protective layer, are found to reduce the values of the local SARs and protect the pregnant women and fetus better.
In addition, not only the SFDTD scheme is better than the traditional FDTD scheme and the high-order FDTD scheme, but also the geometries of the mother/ fetus model are relatively simple, therefore the SFDTD scheme is better than the FEM and the FVM for researching the SARs and the optimal angle of the protective layer. Also with the mother/fetus model becoming more and more complex, the symplectic integrator is considered to be added into the FEM or the FVM based on unstructured grids for its advantages in the future research. And the massively parallel computing, which is based on the new algorithm, such as the SFDTD scheme and the FEM above, is considered too. Of course, a more realistic mother/fetus model for electromagnetic simulation will be constructed in the future research. Also in the future, the algorithm, which is used in this article or will be used in the future, will be improved continuously. And some clinical experiments will need to be done to prove our results.
All the forward simulations and calculations are executed on a personal computer with a 2.93 GHz Intel Core2 Duo processor with 2 GBs of RAM using the C/C++ languages. After the calculated data are obtained by the C/C++ codes, the figures above are generated by the Matlab software using these data.

Author Contributions
Conceived and designed the experiments: YJG HWY. Performed the experiments: YJG HWY. Analyzed the data: YJG. Contributed reagents/materials/analysis tools: YJG HWY. Wrote the paper: YJG.