A Novel Efficient FEM Thin Shell Model for Bio-Impedance Analysis

In this paper, a novel method for accelerating eddy currents calculation on a cell model using the finite element method (FEM) is presented. Due to the tiny thickness of cell membrane, a full-mesh cell model requires a large number of mesh elements and hence intensive computation resources and long time. In this paper, an acceleration method is proposed to reduce the number of mesh elements and therefore reduce the computing time. It is based on the principle of replacing the thin cell membrane with an equivalent thicker structure. The method can reduce the number of mesh elements to 23% and the computational time to 17%, with an error of less than 1%. The method was verified using 2D and 3D finite element methods and can potentially be extended to other thin shell structures. The simulation results were validated by measurement and analytical results.


Introduction
β-dispersion analysis has been widely investigated for medical [1,2] and industrial [3,4] applications. For example, it has been proposed for the detection of tumours [5,6] and cerebral stroke [7], and it has been applied to the quality inspection of food such as meat, fruit and vegetables [3]. The mechanism of β-dispersion was first described by Schwan [8,9]. β-dispersion takes place at radio frequency, mainly from hundreds of hertz to megahertz. This dispersion is caused by Maxwell-Wagner effect, an interfacial polarisation of cell membrane blocking the ion-flow of intra and extracellular dielectrics.
Analytical solutions are well developed to calculate the β dispersion [10][11][12]. However, the analytical solution is only designed to analyze the models with a regular shape (i.e., spherical cell model) [13,14]. In reality, most cell shapes are anomalous. The finite element method is a methodology that discretises an entire continuous domain into discrete domains. An irregular cell shape can be simulated in a high accuracy using the finite element method. Therefore, FEM is a feasible and effective way to simulate the dielectric dispersions for models with irregular shapes [15,16].
The difficulty of FEM is that the number of meshing elements is significantly large due to the tiny thickness of the cell membrane; it takes a long time to compute with numerous meshing elements. In this paper, a novel FEM acceleration method was proposed to simulate β dispersion. The acceleration method is to reduce the number of meshing elements and computing time by replacing the full-mesh cell membrane model with an equivalent reduced-mesh model. In this paper, the FEM The two-dimensional simulation model is a single cell model, as shown in Figure 1. The spherical cell is put in a suspension with an applied alternating electric field. The voltage of the upper electrode was set to be 10V, and the voltage of the lower electrode was set to be 0V [14]. The intracellular and extracellular fluid were conductive while the cell membrane blocked current flow between intracellular and extracellular fluid at low frequency. Dielectric relaxation describes a phenomenon that the permittivity decreases when the frequency of applied electric field increases at a certain range. The relaxation occurs when the electric dipoles are not able to reverse at the same rate with the frequency of the applied field. -dispersion was observed from 10 kHz to 1 MHz. The dispersion was described by Schwan and could be explained by Maxwell-Wagner effect [8]. The effect occurs at the boundary of layered or inhomogeneous dielectrics. Applying an electric field would build up charges at the boundaries created by different dielectrics. The double-layer dielectrics system can be represented by one capacitance in parallel with one another. Applying a static electric field would charge both capacitances and creates charges at the boundary of the dielectrics. This phenomenon is also called interfacial polarisation.
The equivalent complex relative permittivity of the cell model in Figure 1 can be calculated according to Asami [12]: The two-dimensional simulation model is a single cell model, as shown in Figure 1. The spherical cell is put in a suspension with an applied alternating electric field. The voltage of the upper electrode was set to be 10 V, and the voltage of the lower electrode was set to be 0 V [14]. The intracellular and extracellular fluid were conductive while the cell membrane blocked current flow between intracellular and extracellular fluid at low frequency. Dielectric relaxation describes a phenomenon that the permittivity decreases when the frequency of applied electric field increases at a certain range. The relaxation occurs when the electric dipoles are not able to reverse at the same rate with the frequency of the applied field. β-dispersion was observed from 10 kHz to 1 MHz. The dispersion was described by Schwan and could be explained by Maxwell-Wagner effect [8]. The effect occurs at the boundary of layered or inhomogeneous dielectrics. Applying an electric field would build up charges at the boundaries created by different dielectrics. The double-layer dielectrics system can be represented by one capacitance in parallel with one another. Applying a static electric field would charge both capacitances and creates charges at the boundary of the dielectrics. This phenomenon is also called interfacial polarisation.
The equivalent complex relative permittivity of the cell model in Figure 1 can be calculated according to Asami [12]: where n , d m is the cell membrane thickness, R is the radius of the cell, and n stands for the factor of dimension. When n = 2, the equation is applied to the two-dimension model. When n = 3, the equation is applied to three-dimension model. ε * p is the complex relative permittivity of the cell model, ε * m is the complex relative permittivity of the cell membrane and ε * c is the complex relative where P is the volume fraction, ε * is the complex relative permittivity of cell suspension, ε * a is the complex relative permittivity of the external medium.
The conductivity of the suspension can be derived: where σ r and ε r are the conductivity and permittivity of the suspension, respectively. σ * is the complex conductivity of the suspension, ε 0 is the permittivity of vacuum, ω = 2πf where f is the frequency of the applied field and j is the imaginary unit. According to Schwan [9], when d m R 1, k m k a = k c , the dielectric relaxation of the cell suspension can be represented by the Debye relaxation model: where ε h and ε l are the high-frequency limit and low-frequency limit of the relative permittivity, τ is the relaxation time and k l is the low-frequency limit of the conductivity.

D FEM Simulation
Normally, the finite element problems could be solved by two kinds of methods which are Galerkin's method and Ritz's method [19]. In this paper, Galerkin's method is used to build up simulations.
According to the Gaussian equation, D is the electric flux density, → E is the electrical field intensity, ε = ε r ε 0 , ε r is the relative permittivity of the material, ρ is the quantity of electric charge, which is zero in the cell model. So the equation could be rewritten as: where σ is the conductivity of the material and ε = ε r ε 0 . ε r is the relative permittivity of the material. ε 0 is the permittivity of vacuum. U is the electric potential which is the objective value to be solved [20]. To obtain the equivalent conductivity of the total cell suspension, we calculate the total current flowing through the cell suspension from the top electrode to the bottom electrode [15]. According to Kirchhoff's current law, the current flow through any plane, which is in parallel with the electrodes is the same with the total current flow through the cell suspension. So the current I flow through a chosen plane in parallel with electrodes can be calculated as: Biosensors 2020, 10, 69 where j (i) is the current density of the element on the chosen plane, σ (i) is the conductivity of the element on the chosen plane, is the background electrical field on the chosen plane. Then the complex conductivity of the cell suspension can be calculated as: where σ * and ρ * are the complex conductivity and resistivity, respectively. L is the distance between the top electrode and the bottom electrode. S is the area of the electrodes.

D FEM Simulation
An induction model was built for the three-dimension FEM simulation. As shown in Figure 2, the cylinder stands for the cell suspension, and the spherical model stands for the cell. Imaginary transmitter and receiver coils are put on the top of the suspension to provide an alternating magnetic field as an exciting signal and detect induced secondary magnetic field. The single-cell model in a three-dimension situation is shown in Figure 2b. The applied field is a magnetic field excited by the transmitter in Figure 2a. The radius of the sphere cell model is R and the cell membrane thickness of the sphere cell model is d m . All the electric parameters of the sphere cell model are set to be the same with the parameters in Figure 1. The Mxwell-Wagner effect can be calculated by calculating the eddy current distribution induced by the applied magnetic field. where j ( ) is the current density of the element on the chosen plane, σ ( ) is the conductivity of the element on the chosen plane, E ⃗ ( ) is the background electrical field on the chosen plane. Then the complex conductivity of the cell suspension can be calculated as: where * and * are the complex conductivity and resistivity, respectively. L is the distance between the top electrode and the bottom electrode. S is the area of the electrodes.

D FEM Simulation
An induction model was built for the three-dimension FEM simulation. As shown in Figure 2, the cylinder stands for the cell suspension, and the spherical model stands for the cell. Imaginary transmitter and receiver coils are put on the top of the suspension to provide an alternating magnetic field as an exciting signal and detect induced secondary magnetic field. The single-cell model in a three-dimension situation is shown in Figure 2b. The applied field is a magnetic field excited by the transmitter in Figure 2a. The radius of the sphere cell model is R and the cell membrane thickness of the sphere cell model is . All the electric parameters of the sphere cell model are set to be the same with the parameters in Figure 1. The Mxwell-Wagner effect can be calculated by calculating the eddy current distribution induced by the applied magnetic field. According to Oszkar Biro [20], the basic formulas for the three-dimension edge elements simulation of eddy current problems, Galerkin's equation is shown as following: where L is the elemental interpolation of i node corresponding to the nth element of it; N stands for the vector interpolation of i edge relating to the nth edge element of it; A ( ) stands for the excited edge vector potential of the nth element; A represents the original edge vector potential of the nth element; V ( ) represents the electrical potential excited by the nth element on the receiving coil; v stands for the reluctivity of the model; v is the reluctivity of the air; σ is the conductivity of the model [22−24]. The eddy current density is: According to Oszkar Biro [21], the basic formulas for the three-dimension edge elements simulation of eddy current problems, Galerkin's equation is shown as following: where L i is the elemental interpolation of i th node corresponding to the nth element of it; N i stands for the vector interpolation of i th edge relating to the nth edge element of it; A (n) stands for the excited edge vector potential of the nth element; A s represents the original edge vector potential of the nth element; V (n) represents the electrical potential excited by the nth element on the receiving coil; v stands for the reluctivity of the model; v 0 is the reluctivity of the air; σ is the conductivity of the model [22][23][24]. The eddy current density is: where E (i) denotes the vector sum of the electrical field on all the edges of each tetrahedral element, σ (i) is the complex conductivity parameter of the of each tetrahedral element. The model is shown in Figure 2. Assuming there is an equivalent model with uniform dielectric and the background suspension has exactly the same shape with the simulation model in Figure 2.
Since the normal component of the electrical field relative to each cross-section surface of suspension (shown in Figure 2 where, E Since the equivalent model is uniform, the electrical background field E (i) s is always vertical to every cylindrical cross-section surface of the suspension model shown in Figure 2. Then the equivalent complex conductivity of the original suspension (as shown in Figure 2) can be derived from (8) and (9) that s is the eddy current density through each cylindrical cross-section of equivalent suspension model, J (i) is the eddy current density through each cylindrical cross-section of the original suspension model with arranged cells and → n is the current flow direction.

Acceleration Method
The complicated calculation progress and the long computing time are caused by a large number of meshing elements around the thin cell membrane. The number of meshing elements could be reduced by enlarging the cell membrane thickness. However, changing the cell membrane thickness affects the conductivity spectrum of the cell suspensions. The aim of the proposed acceleration method is to keep the accuracy of the simulation result while enlarging the cell membrane thickness (reducing the number of elements). This acceleration method replaces the full-mesh model with an equivalent thicker cell membrane (reduce-mesh model) with an equivalent complex conductivity, as shown in Figure 3. The M 1 in Figure 3 is the full-mesh model cell membrane with normal thickness l 1 . The M 2 in Figure 3 is the enlarged part of the reduce-mesh model cell membrane with thickness l 2 , R is the radius of the cell model. The electrical parameter of the cell membrane M 2 is exactly the same with the intracellular fluid. The cell membrane M 1 and M 2 combine as one thicker equivalent cell membrane (reduce-mesh model cell membrane), and the thickness of the reduce-mesh model cell membrane is l 1 + l 2 . The number of meshing elements is reduced due to a thicker equivalent cell membrane.
In this case, the factor v = 1 − d m  (2), as the parameter P and ε * a of the equivalent model are fixed, the dielectric behaviour of the suspension model ε * only depends on the dielectric property of the cell model ε * p . According to Equation (1), ε * p is determined by ε * m and the factor v. In order to retain the dielectric behaviour of the suspension model, the equivalent dielectric parameters of the equivalent cell membrane ε * m can be calculated by Equations (14)- (18). With the equivalent complex conductivity, the behaviour of Maxwell-Wagner effect, which leads to β-dispersion remains the same.
Biosensors 2020, 10, x FOR PEER REVIEW 6 of 15 With the equivalent complex conductivity, the behaviour of Maxwell-Wagner effect, which leads to β -dispersion remains the same.
where σ * and σ * are the complex conductivity of M cell membrane and M cell membrane in Figure 3, respectively. σ and σ are the conductivity of cell membrane M and cell membrane M in Figure 3, respectively. ε and ε are the relative permittivity of the cell membrane M and cell membrane M , respectively. Considering the two cell membranes M and M are connected as two electrolytes in series. The total impedance Z is: where Z and Z are the impedance of the cell membrane M and M , respectively. σ * and σ * are the complex conductivity of cell membrane M and M respectively. l and l are the thickness of the cell membrane M and the cell membrane M in Figure 3, respectively. S and S are the area of cell membrane M and the cell membrane M , respectively. As l and l are negligible comparing with radius R, the area S = S . Then the equivalent complex conductivity can be derived as: σ * = (l + l )σ * σ * l σ * + l σ * (18) * is the complex conductivity of the equivalent cell membrane (reduce-mesh model cell membrane).
The replacement of the thin cell membrane with an equivalent thicker structure appears to have reflected an equivalent Maxwell-Wagner effect from the computation results. Theoretically, less contrasting EM properties between the boundary of two materials cause less Maxwell-Wagner effect (i.e., less charge building-up). At the same time, thicker materials between the two charged surfaces would mean the equivalent capacitance is smaller. These two effects would mean the electric field E = Q/C would be kept to a constant; this is possibly the reason why we can use the replacement to carry out the calculation of the dispersion caused by the MW effect. Equation (17) where σ * 1 and σ * 2 are the complex conductivity of M 1 cell membrane and M 2 cell membrane in Figure 3, respectively. σ 1 and σ 2 are the conductivity of cell membrane M 1 and cell membrane M 2 in Figure 3, respectively. ε 1 and ε 2 are the relative permittivity of the cell membrane M 1 and cell membrane M 2, respectively. Considering the two cell membranes M 1 and M 2 are connected as two electrolytes in series. The total impedance Z is: where Z 1 and Z 2 are the impedance of the cell membrane M 1 and M 2 , respectively. σ * 1 and σ * 2 are the complex conductivity of cell membrane M 1 and M 2 respectively. l 1 and l 2 are the thickness of the cell membrane M 1 and the cell membrane M 2 in Figure 3, respectively. S 1 and S 2 are the area of cell membrane M 1 and the cell membrane M 2 , respectively. As l 1 and l 2 are negligible comparing with radius R, the area S 1 = S 2 .
Then the equivalent complex conductivity can be derived as: σ * e is the complex conductivity of the equivalent cell membrane (reduce-mesh model cell membrane). The replacement of the thin cell membrane with an equivalent thicker structure appears to have reflected an equivalent Maxwell-Wagner effect from the computation results. Theoretically, less contrasting EM properties between the boundary of two materials cause less Maxwell-Wagner effect (i.e., less charge building-up). At the same time, thicker materials between the two charged surfaces would mean the equivalent capacitance is smaller. These two effects would mean the electric field E = Q/C would be kept to a constant; this is possibly the reason why we can use the replacement to carry out the calculation of the dispersion caused by the MW effect. Equations (17) and (18) considered the complex conductivity, and therefore, the physics would be masked to a certain degree, but the principle can be explained as above.

Meshing Details
The suspension model shown in Figure 1 was discretised into three regions, cell membrane and intra and extracellular fluid, for meshing. The commercial software Comsol 5.0 was used for meshing. The convergent criterion was accelerated BICGS with an optimal precondition. The initial precondition was based on optimised initial guess: the final solution from the previous frequency [25][26][27][28][29][30]. As shown in Table 1, the minimum and maximum element sizes are set to be the same to maintain the meshing accuracy. The maximum element growth rate controls the maximum changing rate of the dimension of the element among the adjacent subdomains. The narrow factor controls the mesh density on the narrow region. Increasing the cell membrane thickness reduces the narrow regions. Therefore, the meshing density is larger in the cell model with smaller cell membrane thickness. Table 2 shows the number of meshing element in each specific region, including total suspension, cell membrane, intracellular fluid and extracellular fluid. In addition, a specific region surrounding the cell membrane was created to obtain the number of meshing elements around the cell membrane. It is obvious from Table 2 that the most meshing elements concentrate at the region around the cell membrane. The narrow region is reduced when the cell membrane thickness is enlarged, thus, the number of meshing elements is decreased when increasing cell membrane thickness.

D Single Spherical Cell Model
The full-mesh model cell membrane thickness is 5 nm. There are two simulations with equivalent thicker cell membrane (reduce-mesh model) thickness of 10 nm and 20 nm, respectively. The ratio between the full-mesh model cell membrane (cell membrane M_1 in Figure 3) thickness and the enlarged part of reduce-model cell membrane (cell membrane M_2 in Figure 3) thickness is 1:1 and 1:3, respectively. The simulation results are the progress of β-dispersion which reflects dispersions on relative permittivity and conductivity, as shown in Figures 4 and 5. The main characteristics of β-dispersion are the dispersion frequency range and the magnitude of the relative permittivity and conductivity [27,30,31]. As shown in Figure 4, the dispersion frequency is ranging from 100 kHz to 10 MHz. The magnitude of relative permittivity at low frequency and high frequency are 960 and 80, respectively. The magnitude of conductivity at low frequency and high frequencies were 0.925 S/m and 1.04 S/m, respectively. This work was aimed at accelerating the computing progress with high accuracy,  Table 3.
Biosensors 2020, 10, x FOR PEER REVIEW 8 of 15 from 100 kHz to 10 MHz. The magnitude of relative permittivity at low frequency and high frequency are 960 and 80, respectively. The magnitude of conductivity at low frequency and high frequencies were 0.925 S/m and 1.04 S/m, respectively. This work was aimed at accelerating the computing progress with high accuracy, 15 sample points are selected in Figure 4 to describe the main feature of an entire β -dispersion and verify the accuracy of the acceleration method. The results of acceleration models are validated by the full-mesh model to verify the accuracy of the acceleration method. The errors between the acceleration model and full-mesh model at each point in Figure 4 and Figure 5 have been calculated. The calculated error and computing time for acceleration model and full-mesh model are shown in Table 3.   Table 3, the number of elements of the full-mesh model can be reduced to 24% when the equivalent cell membrane (reduce-mesh model) is four times as the thickness of the full-mesh model cell membrane ( M cell membrane to M cell membrane thickness ratio is 1:3). The computing time was reduced from 73 minutes to 13 minutes, with only a 0.2% error on the simulation result. The computing time of the reduce-mesh model is 18% of the full-mesh model. The spectroscopy of conductivity and relative permittivity is shown in Figure 4 and Figure 5, the results of the reduce-model and full-mesh model show the same magnitude over the same frequency range with an error less than 0.2%. The β dispersion of reduce-model and full-mesh model both starts at 100 kHz and both ends at 10 MHz. This shows that the reduce-mesh model can significantly accelerate the computing progress with only a tiny error on the result in the 2D-FEM spherical model simulation. Table 3. Results of the 2D spherical cell model. from 100 kHz to 10 MHz. The magnitude of relative permittivity at low frequency and high frequency are 960 and 80, respectively. The magnitude of conductivity at low frequency and high frequencies were 0.925 S/m and 1.04 S/m, respectively. This work was aimed at accelerating the computing progress with high accuracy, 15 sample points are selected in Figure 4 to describe the main feature of an entire β -dispersion and verify the accuracy of the acceleration method. The results of acceleration models are validated by the full-mesh model to verify the accuracy of the acceleration method. The errors between the acceleration model and full-mesh model at each point in Figure 4 and Figure 5 have been calculated. The calculated error and computing time for acceleration model and full-mesh model are shown in Table 3.   Table 3, the number of elements of the full-mesh model can be reduced to 24% when the equivalent cell membrane (reduce-mesh model) is four times as the thickness of the full-mesh model cell membrane ( M cell membrane to M cell membrane thickness ratio is 1:3). The computing time was reduced from 73 minutes to 13 minutes, with only a 0.2% error on the simulation result. The computing time of the reduce-mesh model is 18% of the full-mesh model. The spectroscopy of conductivity and relative permittivity is shown in Figure 4 and Figure 5, the results of the reduce-model and full-mesh model show the same magnitude over the same frequency range with an error less than 0.2%. The β dispersion of reduce-model and full-mesh model both starts at 100 kHz and both ends at 10 MHz. This shows that the reduce-mesh model can significantly accelerate the computing progress with only a tiny error on the result in the 2D-FEM spherical model simulation. Table 3. Results of the 2D spherical cell model.  As shown Table 3, the number of elements of the full-mesh model can be reduced to 24% when the equivalent cell membrane (reduce-mesh model) is four times as the thickness of the full-mesh model cell membrane (M 1 cell membrane to M 2 cell membrane thickness ratio is 1:3). The computing time was reduced from 73 min to 13 min, with only a 0.2% error on the simulation result. The computing time of the reduce-mesh model is 18% of the full-mesh model. The spectroscopy of conductivity and relative permittivity is shown in Figures 4 and 5, the results of the reduce-model and full-mesh model show the same magnitude over the same frequency range with an error less than 0.2%. The β dispersion of reduce-model and full-mesh model both starts at 100 kHz and both ends at 10 MHz. This shows that the reduce-mesh model can significantly accelerate the computing progress with only a tiny error on the result in the 2D-FEM spherical model simulation.

D Cell Deformation Model
This simulation is to verify that the acceleration method not only works on spherical cell model but also works on the deformation model (oval model).
The deformation model is an oval model in order to simulation the cell deformation [28]. The electrical parameters of the oval model are the same as that of the custom spherical model. The model is shown in Figure 6, where a and b are the length of the semi-major and semi-minor axis, respectively. The length of the semi-major axis is set to be a = 12 µm and the length of the semi-minor axis is set to be b = 2 µm, k a and k c are the conductivity of the extracellular and intracellular fluid, respectively, ε c and ε a are the relative permittivity of the intracellular and extracellular fluid respectively, ε m is the relative permittivity of cell membrane, d m is the cell membrane thickness. The full-mesh model cell membrane thickness is still 5 nm and the electrical properties are calculated using equation 18. The simulation result of β dispersion is shown in Figure 7, Figure 8 and the parameters of the computing time and error are shown in Table 4.

D Cell Deformation Model
This simulation is to verify that the acceleration method not only works on spherical cell model but also works on the deformation model (oval model).
The deformation model is an oval model in order to simulation the cell deformation [28]. The electrical parameters of the oval model are the same as that of the custom spherical model. The model is shown in Figure 6, where a and b are the length of the semi-major and semi-minor axis, respectively. The length of the semi-major axis is set to be a = 12 and the length of the semi-minor axis is set to be b = 2 , and are the conductivity of the extracellular and intracellular fluid, respectively, and are the relative permittivity of the intracellular and extracellular fluid respectively, is the relative permittivity of cell membrane, is the cell membrane thickness. The full-mesh model cell membrane thickness is still 5 and the electrical properties are calculated using equation 18. The simulation result of dispersion is shown in Figure 7, Figure 8 and the parameters of the computing time and error are shown in Table 4.

D Cell Deformation Model
This simulation is to verify that the acceleration method not only works on spherical cell model but also works on the deformation model (oval model).
The deformation model is an oval model in order to simulation the cell deformation [28]. The electrical parameters of the oval model are the same as that of the custom spherical model. The model is shown in Figure 6, where a and b are the length of the semi-major and semi-minor axis, respectively. The length of the semi-major axis is set to be a = 12 and the length of the semi-minor axis is set to be b = 2 , and are the conductivity of the extracellular and intracellular fluid, respectively, and are the relative permittivity of the intracellular and extracellular fluid respectively, is the relative permittivity of cell membrane, is the cell membrane thickness. The full-mesh model cell membrane thickness is still 5 and the electrical properties are calculated using equation 18. The simulation result of dispersion is shown in Figure 7, Figure 8 and the parameters of the computing time and error are shown in Table 4.  The number of elements of full-mesh model can be reduced to 24.1% as shown in Table 4 and the computing time is reduced from 2h 24mins to 25mins. The simulation error between reduce-mesh model and full-mesh model is only 0.11%. As shown in Figure 7 and Figure 8, the results of the full-mesh model and reduce-mesh model exhibits dispersion with tiny error in the same frequency range. The results show that the acceleration method is feasible to simulate not only the spherical model but also deformation model (oval model). The computing time is significantly reduced with acceptable error.   The number of elements of full-mesh model can be reduced to 24.1% as shown in Table 4 and the computing time is reduced from 2h 24mins to 25mins. The simulation error between reduce-mesh model and full-mesh model is only 0.11%. As shown in Figures 7 and 8, the results of the full-mesh model and reduce-mesh model exhibits β dispersion with tiny error in the same frequency range. The results show that the acceleration method is feasible to simulate not only the spherical model but also deformation model (oval model). The computing time is significantly reduced with acceptable error.

D Spherical Cell Model
The simulation model used for 3D FEM simulation is shown in Figure 2. 3-D single sphere cell model is used. The radius of the cell model is set to be 5 µm, and the cell membrane thickness is 5 nm. The top view of eddy current distributions at lower frequency 1 kHz and higher frequency 10 MHz are shown in Figure 9. The single-cell model blocks eddy currents at a lower frequency and becomes conductive at a higher frequency which meets the expectations of maxwell-wagner effect and other published works [9,11,15,29]. The full-mesh model cell membrane thickness is 5 nm. There are two reduce-models with equivalent thicker cell membrane (reduce-mesh model) thickness of 10 nm and 20 nm, respectively. The ratio between the full-mesh model cell membrane thickness and the reduce-mesh cell membrane thickness is 1:1 and 1:3, respectively. The relative permittivity and conductivity result are shown in Figures 10 and 11. β dispersion can be observed from Figures 10 and 11. The β dispersion of the full-mesh model and reduce-mesh model has an error of 2% on the magnitude. The frequency range of the β dispersion of full-mesh model and reduce-mesh model are the same, ranging from 100 kHz to 10 MHz. The results, including the number of elements, the computing time and error of full-mesh model and reduce-mesh model, are shown in Table 5. The number of elements of the full-mesh model can be reduced to 28%, as shown in Table 5, and the computing is reduced from 4 h 37 min to 71 min. The simulation error between reduce-mesh model and full-mesh model is 2%. The error is larger than the two-dimensional simulation result, but it is acceptable for β-dispersion simulations. The results show that the acceleration method is feasible to simulate a three-dimensional spherical model. The computing time is significantly reduced with acceptable error.

MHz.
The full-mesh model cell membrane thickness is 5 nm. There are two reduce-models with equivalent thicker cell membrane (reduce-mesh model) thickness of 10nm and 20nm, respectively. The ratio between the full-mesh model cell membrane thickness and the reduce-mesh cell membrane thickness is 1:1 and 1:3, respectively. The relative permittivity and conductivity result are shown in Figure 10 and Figure 11. β dispersion can be observed from Figure 10 and Figure 11. The β dispersion of the full-mesh model and reduce-mesh model has an error of 2% on the magnitude. The frequency range of the β dispersion of full-mesh model and reduce-mesh model are the same, ranging from 100 kHz to 10 MHz. The results, including the number of elements, the computing time and error of full-mesh model and reduce-mesh model, are shown in Table 5. The number of elements of the full-mesh model can be reduced to 28%, as shown in Table 5, and the computing is reduced from 4 h 37 mins to 71 mins. The simulation error between reduce-mesh model and full-mesh model is 2%. The error is larger than the two-dimensional simulation result, but it is acceptable for β -dispersion simulations. The results show that the acceleration method is feasible to simulate a three-dimensional spherical model. The computing time is significantly reduced with acceptable error.

Experimental Validation by Contact-Electrode Measurement
The four-electrode measurement system is used to measure the impedance spectroscopy of biological samples. The electrode is contacted directly to the samples to measure the impedance over the samples using an impedance analyser. The measurement system work in the following way: the electrodes generate an input signal which applies an alternating electric field and current in the samples. By measuring the current and voltage across the samples, the conductivity and permittivity of the samples can be calculated. The impedance analyser used in this work is Solatron 1260, which is efficient and accurate over the frequency ranging from hundreds of hertz to 32 MHz. The effective area of the measurement electrodes is fixed to S = 2 cm * 3 cm = 6 cm , and the distance between the electric field measuring electrodes is fixed to L = 6.6 cm. The conductivity can then be calculated by σ = 1/ρ = L/RS. R is the real part of the measured impedance, ρ is the resistivity. The measurement sample is fresh potato (obtained from a local market). This measurement is to validate the FEM simulation result. All FEM simulation parameters are set to be the realistic value obtained from the measurement.
In Figure 12, the FEM and measurement results show similar β -dispersion with the same dispersion frequency range. There is some error on the magnitude of the conductivity between the Figure 11. The conductivity of 3D single spherical cell model.

Experimental Validation by Contact-Electrode Measurement
The four-electrode measurement system is used to measure the impedance spectroscopy of biological samples. The electrode is contacted directly to the samples to measure the impedance over the samples using an impedance analyser. The measurement system work in the following way: the electrodes generate an input signal which applies an alternating electric field and current in the samples. By measuring the current and voltage across the samples, the conductivity and permittivity of the samples can be calculated. The impedance analyser used in this work is Solatron 1260, which is efficient and accurate over the frequency ranging from hundreds of hertz to 32 MHz. The effective area of the measurement electrodes is fixed to S = 2 cm * 3 cm = 6 cm 2 , and the distance between the electric field measuring electrodes is fixed to L = 6.6 cm. The conductivity can then be calculated by σ = 1/ρ = L/RS. R is the real part of the measured impedance, ρ is the resistivity. The measurement sample is fresh potato (obtained from a local market). This measurement is to validate the FEM simulation result. All FEM simulation parameters are set to be the realistic value obtained from the measurement.
In Figure 12, the FEM and measurement results show similar β -dispersion with the same dispersion frequency range. There is some error on the magnitude of the conductivity between the measurement result and the simulation result. However, the error is acceptable for the simulation approach of the measurement result. The conductivity of measurement and simulation result is low at a lower frequency, and the dispersion starts at around 50 kHz. The dispersion ends at around 2 MHz and the conductivity of measurement and simulation result is increased to 0.15 S/m. The spectroscopy curve is flat over the higher frequency range. The measurement result validated the full-mesh and reduce-mesh FEM simulation result and agreed with other published works [9,13,15,21,25,32].
Biosensors 2020, 10, x FOR PEER REVIEW 13 of 15 Figure 12. Impedance spectroscopy of FEM and measurement result.

Verification by Analytical Result
The analytical result can be calculated according to Maxwell-Wagner Equations (1) and (2). The analytical result is compared with the 3D finite element simulation result. So the factor v in equation (1) is defined as v = (1 − ) . According to Schwan [9], the impedance spectroscopy of the cell model can be described by Debye relaxation based on equation (4). In Figure 13 and Figure  14, both analytical and FEM results show the same magnitude and frequency range of beta dispersion with an error less than 0.1%. The volume fraction of the cell is set to be the same as P=3.5%. The magnitude and frequency range of -dispersion on the full-mesh model and reduce-mesh model agree with the analytical result. The result validates that the three-dimension FEM simulation on cell suspension is accurate and the improved acceleration method is also validated. The proposed acceleration method can be applied to further FEM simulations on irregular shape cell models and thin shell models.

Verification by Analytical Result
The analytical result can be calculated according to Maxwell-Wagner Equations (1) and (2). The analytical result is compared with the 3D finite element simulation result. So the factor v in equation (1) is defined as v = 1 − d m . According to Schwan [9], the impedance spectroscopy of the cell model can be described by Debye relaxation based on Equation (4). In Figures 13 and 14, both analytical and FEM results show the same magnitude and frequency range of beta dispersion with an error less than 0.1%. The volume fraction of the cell is set to be the same as P = 3.5%. The magnitude and frequency range of β-dispersion on the full-mesh model and reduce-mesh model agree with the analytical result. The result validates that the three-dimension FEM simulation on cell suspension is accurate and the improved acceleration method is also validated. The proposed acceleration method can be applied to further FEM simulations on irregular shape cell models and thin shell models.
Biosensors 2020, 10, x FOR PEER REVIEW 13 of 15 Figure 12. Impedance spectroscopy of FEM and measurement result.

Verification by Analytical Result
The analytical result can be calculated according to Maxwell-Wagner Equations (1) and (2). The analytical result is compared with the 3D finite element simulation result. So the factor v in equation (1) is defined as v = (1 − ) . According to Schwan [9], the impedance spectroscopy of the cell model can be described by Debye relaxation based on equation (4). In Figure 13 and Figure  14, both analytical and FEM results show the same magnitude and frequency range of beta dispersion with an error less than 0.1%. The volume fraction of the cell is set to be the same as P=3.5%. The magnitude and frequency range of -dispersion on the full-mesh model and reduce-mesh model agree with the analytical result. The result validates that the three-dimension FEM simulation on cell suspension is accurate and the improved acceleration method is also validated. The proposed acceleration method can be applied to further FEM simulations on irregular shape cell models and thin shell models.

Conclusions
This paper proposed a method to accelerate FEM calculation with bio-cell models. The idea is to replace the thin cell membrane (full-mesh model) with an equivalent thicker cell membrane (reduce-mesh model). Then the number of meshing elements of the full-mesh model is reduced, and thus the computing time is reduced.
According to the simulation results, the reduce-mesh model can be used in both 2D and 3D FEM simulations. The amount of computing time is significantly reduced with an error no more than 0.5% in 2D simulation. The error on 3D simulation result is no more than 2%. All simulation results are validated by measurements and analytical results.
The proposed acceleration method is validated to be fast and accurate on cell models and the acceleration method has potential for all other thin shell FEM models.