Bioimpedance Calculation Considering Spherical Anomalies

Abstract

This paper presents a solution for sphere-shaped heterogeneities applied to the bioimpedance test. Such an approach has applicability for electrical impedance tomography, considering the impedance response for electrodes at the skin level, a healthy tissue composing a half-space, and a spherical anomaly representing a malignant or even benign tumor. The presented methodology consists of an inverse solution; however, it has high computational performance and can be later applied to an optimization procedure for generating medical images. The proposed formulation does not require a domain discretization and is validated in the paper using Finite Element Analysis.

Index Terms
bioimpedance; partial differential equations; point source; electrical impedance tomography

I. Introduction

According to [ ¹ ], the Electrical impedance tomography (EIT) is a non-invasive method to monitor conductivity changes in regions of the human body. Such a method reconstructs biological tissues’ electrical properties (conductivity and permittivity) by electric potential measurements at the body surface generated from an electric current injection. The electric current is alternated and has a small magnitude, being injected near a region of interest. There are pairs of active and passive electrodes, where active electrodes consist of current injection and passive consist of voltage readings. The EIT is a low-cost imaging method, which is also free of ionizing radiation [ ² ]. Despite the low-resolution of the generated images, a tomographer of this type is portable and has low-cost.

There is a current injection at a pair of active electrodes and a voltage reading at passive electrodes to determine electric impedance used to estimate the material’s electrical properties. For [ ² ], this process is repeated for several excitation and response sets until it is collected as much information as possible about conductivity distribution.

The solution for EIT is to determine the size and position of a possible anomaly trough voltage measurements. This solution consists of using measurement data to estimate spatial conductivity and permittivity distribution. This paper approach does the inverse; it calculates the voltage response of an anomaly when its size and position are determined, i.e. a known spatial conductivity and permittivity distribution. The spherical anomaly is a model of a tumor that can be placed anywhere in the threedimensional space. The EIT uses the discrepancy of living tissues’ electric properties (e.g. healthy tissues, a benign tumor, a malign tumor, bones, vessels) for a non-invasive imaging method.

This problem’s usual approach is not to solve analytically, which is limited to simple geometries but using dedicated Finite Element Method (FEM) based software, as in [ ³ ]. Such an approach is not limited to simple geometries and calculates the spatial conductivity distribution based on the passive electrodes’ response to electric current injection. Our approach has the advantage of no need for meshing and gives an analytical solution for a simplified geometry. The authors believe the current paper approach has promising applicability for mammography due to its characteristics. A further step, out of this paper’s scope, is developing an optimization algorithm to reconstruct the conductivity distribution and generate medical images.

II. Electric potential solution

According to [ ³ ], at low frequencies and small field strengths, the electromagnetic properties of living tissue allow its solution trough the generalized Laplace’s equation. Time-varying potentials describe the propagation in the matter for time-varying sources. On the other hand, the retardation effect can be negligible if the wavelength in the matter is much higher than the analyzed physical dimensions; such an approach is called quasi-static.

Based on the assumption above, this paper model considers a uniformly conducting medium of conductivity σ and of an infinite extent within a harmonic current source ${\overrightarrow{J}}_{\text{e}}$ lies. The angular frequency is ω and ε = ε_rε ₀ is the medium permittivity.

(1) $$\overrightarrow{J}=\sigma \overrightarrow{E}+{\overrightarrow{J}}_e$$

The Gauss law, where ρ_v is the space volumetric charge density, is also:

(2) $$\nabla \cdot \overrightarrow{D}={\rho }_{\upsilon }$$

By applying the time-harmonic divergence equation:

(3) $$\nabla \cdot \overrightarrow{J}=\nabla \cdot \left(\sigma \overrightarrow{E}+{\overrightarrow{J}}_{\text{e}}\right)=-j\omega {\rho }_{\upsilon }$$

Combining (2) and (3) one achieves (4) , which is a generic form of the current density with conduction and displacement currents terms.

(4) $$\overrightarrow{J}=\sigma \overrightarrow{E}+j\omega \overrightarrow{D}+{\overrightarrow{J}}_{\text{e}}$$

The quasi-static electric field (neglecting the $-\partial \overrightarrow{A}/\partial t$ term) is derived from scalar potential V , which is

(5) $$\overrightarrow{E}=-\nabla \tilde{V}$$

The applying of (5) in (4) yields the following equation

(6) $$-\nabla \cdot \left(\left(\sigma +j\omega {\epsilon }_r{\epsilon }_0\right)\nabla V-{\overrightarrow{J}}_{\text{e}}\right)=0$$

Assuming a complex conductivity as $\overline{\sigma }=\sigma +j\omega {\epsilon }_r{\epsilon }_0$ . For a souce located at $\overrightarrow{r^{\prime }}$ and a calculation point at $\overrightarrow{r}$ , according to [ ⁴ ], one has $\nabla \cdot {\overrightarrow{J}}_e=-\tilde{I}\delta \left(\overrightarrow{r},\overrightarrow{r^{\prime }}\right)$ . The Poisson equation is then obtained for the scalar potential.

(7) $${\nabla }^2\tilde{V}=-\frac{\tilde{I}\delta \left(\overrightarrow{r},\overrightarrow{r^{\prime }}\right)}{\overline{\sigma }}$$

The solution of (6) for a homogeneous space is

(8) $$\tilde{V}\left(\overrightarrow{r},\overrightarrow{r^{\prime }}\right)=\frac{\tilde{I}}{4\pi \overline{\sigma }R},$$

where R′ is the radial distance from the source point.

A. Spherical anomaly

In the literature, an often-used model for inhomogeneity is a sphere. Such a model, as illustrated in Fig. 1 , where the current electrode, injecting a current I into the medium, is located at z = z ₀ on the z-axis of an (x, y, z) coordinate system. The conductivity of the whole space is σ . Therefore the potential V at a point in space can be expressed with spherical symmetry. The sphere external material has a conductivity σ ₁ and permittivity ε ₁, while the internal material σ ₂ and ε ₂. This problem has symmetry around the ϕ variable of a ( r , θ , ϕ ) coordinate system.

Fig. 1
Spherical anomaly model and current injection.

The V ( r , θ ) function is obtained by speration of variables and has the following general solution [ ⁵ ]

(9) $$V\left(r,\theta \right)={\displaystyle \sum _{n=0}^{\infty }\left(A_n{r}^n+B_n{r}^{-n-1}\right)}{P}_n\left(\cos \theta \right),$$

where A_n and B_n are functions that satisfies the boundary conditions for the layers interface and P_n is the n-th degree Legendre polynomial.

According to Fig. 1 , the point source solution has a radial distance from the calculation point

(10) $$R=\sqrt{r^2+z_0^2-2r{z}_0\cos \theta }\cos \theta =\frac{z}{r}.$$

The electric potential for the outer region consists of the superposition of the point source solution for a homogeneous space with the addition of the general solution. Since rⁿ diverges when the distance tends to infinity, the A_n terms must be zero to the solution be consistent with the physical phenomena. The resulting superposition is then

(11) $$V_1=\frac{I}{4\pi {\overline{\sigma }}_{1}R}+{\displaystyle \sum }_{n=0}^{\infty }{B}_n\frac{1}{r^{n+1}}{P}_n\left(\cos \theta \right)$$

By expanding R ⁻¹ in terms of spherical harmonics, which is valid for r ₀ < r < z ₀, one achieves

(12) $$V_1=\frac{I}{4\pi {\overline{\sigma }}_1}{\displaystyle \sum _{n=0}^{\infty }\frac{r^n}{z_0^{n+1}}{P}_n\left(\cos \theta \right)+{\displaystyle \sum _{n=0}^{\infty }{B}_n\frac{1}{r^{n+1}}{P}_n\left(\cos \theta \right)}}.$$

This form is necessary for applying the boundary conditions, where the characteristic functions are defined. This equation is only used as a comparison basis for the sphere interface. Equation (11) is necessary to calculate the potential without restrictions.

For the interior of the sphere, one needs to calculate the potential at r = 0. That solution diverges unless B_n = 0, so for the physical solution, all terms of this kind are null.

(13) $$V_2={\displaystyle \sum _{n=0}^{\infty }{A}_n{r}^n{P}_n\left(\cos \theta \right)}$$

At the sphere shell, the electric potential and current density normal component must be equal, which are the boundary conditions for the problem.

(14) $$\widetilde{V_1}=\widetilde{V_2}\left(r=r_0\right)$$

(15) $${\overline{\sigma }}_1\frac{\partial \widetilde{V_1}}{\partial r}={\overline{\sigma }}_2\frac{\partial \widetilde{V_2}}{\partial r}\left(r=r_0\right)$$

The applying of the boundary conditions for this kind of problem, adapted from [ ⁶ ], results in the following functions.

(16) $$A_n=\frac{\tilde{I}}{4\pi }\frac{1}{z_0^{n+1}}\frac{{\overline{\sigma }}_1\left(1+2n\right)}{{\overline{\sigma }}_1\left(n+1\right)+n{\overline{\sigma }}_2}$$

(17) $$B_n=\frac{\tilde{I}}{4\pi {\overline{\sigma }}_1}\frac{r_0^{2n+1}}{z_0^{n+1}}\frac{n\left({\overline{\sigma }}_1-{\overline{\sigma }}_2\right)}{{\overline{\sigma }}_1\left(n+1\right)+n{\overline{\sigma }}_2}$$

The solution for the potential functions arises from the substitution of (16) and (17) at (11) and (13) . Figure 2 shows equipotential lines for a section of the domain. The current source is placed at a distance of p from the centre of a spherical anomaly with contrasting electrical properties. In the figure, one can see the equipotential lines’ deformation due to the anomaly; such a deformation is a function of the electrical properties, the sphere radius and location.

Fig. 2
The equipotential lines for a current source which is outside a spherical anomaly.

III. Impedance calculation

The impedance is the relation of the harmonic potential and the injected current. This paper considers a four-electrode array to simulate the real measurement. Two electrodes describe the current flow for the electrode array, and the other two, the measured potential difference. The specialized literature uses bidimensional arrays to compute higher resolution images [ ⁷ ]. This strategy is also viable with this paper equation set; however, a linear array is considered at the simplified study.

The model of Fig. 1 is also symmetrical for any rotation of the z-axis. Consider a rotation of an angle α for the z-axis, the Fig. 3 illustrates the equipotential lines due to a current source for such a scenario. For the same case, if one considers the current source at an interface with the skin-air plan, the same solution is also viable, calculating a mapped depth z ₀ in terms of the sphere depth and horizontal distance from the centre of the sphere.

Fig. 3
Rotated z-axis potential

According to [ ⁸ ], for a half-space such as this model, the full space solution’s equipotential lines only cross half of the solution. The potential for such a problem is then twice the entire space solution. The authors used this assumption to remove this additional heterogeneity.

As an illustrative case, when the sphere is at the centre of a measurement line, where the distance between adjacent electrodes is a , the current source’s position is then 1.5 a from the centre of the measurement line (cf. Fig. 4 ). The z ₀ value for this new configuration is not equal to the sphere’s depth, but proportional to it and half of the array length.

Fig. 4
Linear electrode array

(18) $$z_0=\frac{1}{2}\sqrt{9a^2+4p^2}$$

The equally spaced linear array, where a is the electrode spacing for the current excitation at external electrodes and measured at internal electrodes, is presented in Fig. 4 . The figure considers a sphere of radius r ₀ at depth p from the measurement line for a α degree counterclockwise rotated z-axis for positive source and − α for the negative source.

According to Fig. 4 , the angle for the positive current source is α , being the angle between the rotated z-axis and the first measurement point β and the angle between this axis and the second measurement point is 2( α − β ). For a counterclockwise reference, the negative current source’s displacement angle is − α , −2( α − β ) between this new axis and the first measurement point and − β for the second measurement point. The rotation angle α for the array is calculated for the depth and the array length. The angle β between the rotating axis and the first potential calculation point also depends on the referred parameters.

(19) $$\alpha ={\tan }^{-1}\left(\frac{3a}{2p}\right)\beta ={\tan }^{-1}\left(\frac{a}{p}\right)$$

The solution for Fig. 4 configuration is then a superposition of the solution for a current source rotated counterclockwise, and a negative source rotated clockwise. There are two potential values to be calculated, and a practical measurement consists of this voltage difference.

If one intended to perform all calculations in spherical coordinates, the corrected radial distance for the calculation points, i.e., both potential electrodes is

(20) $$r^{\prime }=\frac{1}{2}\sqrt{a^2+4p^2}$$

One can also use a hybrid coordinate system (spherical and cartesian) to determine the electric potential calculation point. When used the cartesian system, the calculation points are determined trough a rotation linear transform on coordinates. The simple relation between measured voltage and injected current is the impedance, a complex number, which is

(21) $$\overline{Z}=\frac{\text{Δ}\tilde{V}}{\tilde{I}}=\frac{{\tilde{V}}_A-{\tilde{V}}_B}{\tilde{I}}=\Re \left(\overline{Z}\right)+j\Im \left(\overline{Z}\right)$$

For the Fig. 4 potential at voltage probes, denominating Ṽ _{1 p} the Ṽ ₁ value for a positive current source and Ṽ _{1 m} its value considering a negative current source of the same magnitude is for electrode A

Fig. 5
Equipotential lines deformation due to anomaly for the positive and negative current sources.

(22) $${\tilde{V}}_A=2\left({\tilde{V}}_{1p}\left(r^{\prime },\beta \right))+{\tilde{V}}_{1m}(r^{\prime },2\left(\alpha -\beta \right)\right),$$

and the value for electrode B is then

(23) $${\tilde{V}}_B=2\left({\tilde{V}}_{1p}(r^{\prime },2\left(\alpha -\beta \right))+{\tilde{V}}_{1m}\left(r^{\prime },\beta \right)\right),$$

the measured impedance is therefore Z̅ = ( Ṽ_A − Ṽ_B )/ Ĩ . Note that our model only solves linear problems because it is based on the superposition principle.

A. Finite Element Analysis

This paper validates the proposed model by a Finite Element Analysis (FEA). Such analysis consists of transforming a system of differential equations describing an electromagnetic phenomenon in a system of algebraic equations whose solution is close to resolving a differential equations system at nodes in meshes of finite elements. Such methodology requires a closed domain surrounded by a surface for a three-dimensional case. According to [ ⁹ ], for problems where the fields extend to infinite, we need to delimit the study domain by applying a closed surface, separated far enough from the device, from which we can consider the field small enough. For field computation at a conductive medium, the pellicular effect occurs, and we can truncate the tissue correctly for a determined depth without loss of accuracy.

The system of equations has no solution without imposing boundary conditions, also called constraints. At the solution limits, part of the boundary has a fixed potential, and the other part an imposed derivative for electric potential. The fixed-voltage condition can model the sources, which can be a point on the three-dimensional space. The zero derivative condition models the air-skin interface, and forces the electric field to be parallel to such an interface also representing symmetry conditions.

The current paper uses the proprietary software Altair Flux 3D 20.1 [ ¹⁰ ] steady-state AC currents application for modelling the potential generated by the probes at a person’s skin and stimulate the impedance response. There is no current injection condition for this software application, therefore the signal injection where modelled by punctual fixed-voltages at the current probe locations. The electric current is calculated by the imposed voltage difference and its resulting complex power. The impedance is the reason of this voltage difference at the central electrodes and the calculated current. Due to its symmetry, we simulate only half of the domain, which was also truncated based on the probes’ pellicular depth and spacing.

Figure 6 presents the mesh generated by the Altair Flux. Only half of the domain is necessary due to the problem symmetry, aiming to reduce computer efforts. The model represents the electrodes at the skin level, located at the top of the device. There is a refinement for the mesh close to the electrodes, which is necessary due to higher fields magnitude at the region, improving the model’s accuracy for the impedance calculation. Figure 7 shows the simulated mesh cross-section near the spherical anomaly and electrodes. Every electrode spacing requires a reformulation of this mesh, i.e. to compute a complete impedance profile, we need multiple FEM simulations. One can observe the refinement near the measurement line, which comprehends the injected signal points and voltage difference reading points.

Fig. 6
a) The half device, b) the full device. Only half of the domain is simulated due to its symmetry by the imposition of a boundary condition. The post-processing calculations consider the full domain.

Fig. 7
The simulated FEM mesh near the spherical anomaly and electrodes.

This paper intends to reduce computational effort compared to the FEA since there is no need for the domain discretization and multiple simulations. The proposed analytical formula calculates every point in space isolated and requires no meshing for the entire domain. The FEA has the advantage of simulating complex geometries; on the other hand, our paper model uses a simplified approach specialized in the studied geometry. The authors believe this novel approach has promising applicability for EIT imaging due to its computational performance.

B. Example application

Figure 5 illustrates the equipotential lines due to spherical anomaly; in such a figure, medium 1 is healthy tissue, and medium 2 is a malignant tumor model. Reference [ ¹¹ ] measured the electrical properties of breast tissues for healthy tissue, cancerous, and benign tumor properties—the different kinds of tissue present distinct electrical properties, which are also frequency-dependent. The example parameters and physical dimensions are as follows:

• Electrode spacing ( a ): 1 cm

• Sphere depth ( p ): 2 cm

• σ₁ = 0.245 S/m ε _{r 1} = 9.0170 (500 kHz)

• σ₂ = 1.697 S/m ε _{r 2} = 66.696 (500 kHz)

The simulated frequency is 500 kHz, which wavelength at the healthy tissue is 9.03 m. The simulation studies a few centimetres resolution; for such, the authors neglected the propagation effect for this frequency.

Figure 5 shows positive equipotential lines, represented as continuous lines, and the negative lines are dotted. One can observe the deformation of such lines due to the spherical anomaly. This figure is a sectional view of the electric potential; however, the potential formula contemplates a three-dimensional space. The impedance is the voltage difference between passive electrodes and the current injected at the active electrodes. Higher electrode spacings produce a higher penetration of the electromagnetic fields on the tissue. On the other hand, a smaller electrode spacing results in shallower field penetration.

A Finite Element Model for the same simplified geometry was performed using Altair Flux software. The simulation uses several values for electrodes spacing a, and the impedance is measured through voltage electrodes. Table I presents the FEA results compared to the calculated trough the proposed solution. As one can see in the table, the results are in high accordance, which indicates the paper’s model accuracy. Our model excels in calculating the potential for specific points when FEA needs a three-dimensional domain discretization for such.

Table I
Comparison with FEA

a (cm)	Z̅ proposed		Z̅ FEA		deviation
	ℜ( Z̅ )	ℑ( Z̅ )	ℜ( Z̅ )	ℑ( Z̅ )
1	63.424	-0.065	64.744	-0.067	-2.04%
2	31.289	-0.032	32.327	-0.033	-3.21%
3	20.977	-0.022	21.594	-0.022	-2.86%
4	15.864	-0.016	16.180	-0.017	-1.95%

• ^*
  The impedance is expressed in Ω

C. Parametric study

This section performs a parametric study for the linear array measurement. Our model is not limited to this particular configuration, but we use it as a comparison basis. All the parametric studies in this section utilize several electrode spacings. A small electrode spacing produces an impedance reading representing better shallower tissue conductivities and an ample electrode spacing deeper ones. As the figure illustrates, the real part of the impedance is much more significant than the imaginary part for the analyzed frequency. Although less significant, the imaginary part can aid the anomaly localization through a proper parameter estimation routine.

This section intends to demonstrate that several factors imply a different result for the measured impedance. All the figures present the real and imaginary part of the impedance. Several measurements at different locations are the basis for constructing medical images to determine a tumor and its location and size. For a comparison purpose, all figures also have the impedance response for the homogeneous medium composed of only healthy tissue.

Figure 8 presents a simulated measurement, where the spacing of the electrodes is varied, for a parametric variation of the sphere radius and a fixed depth. The sphere and outer media are the same as Fig. 5 , for the same frequency and electrical properties. All studies of this section use the same electric properties for the media, a simulation of healthy and cancerous tissues. Note the impedance response is more contrasting for a more prominent anomaly as expected, and the impedance becomes smaller for a greater radius.

Fig. 8
Measured impedance for several tumor (sphere) radius. a) Real part, b) Imaginary part.

Figure 9 also shows a simulated measurement, but for a parametric variation of the sphere depth and a fixed radius. Both Figures 8 and 9 use a measurement line aligned with the sphere for z = 0, representing the skin level. The figure shows that impedance response is inversely proportional to the anomaly’s depth. If the tumor size is small and deep, the impedance response tends to the healthy tissue response. We have to consider the size and depth ratio for a higher impedance response. Note the impedance is more contrasting for a shallower anomaly, and its value is smaller for a better response as in the first case.

Fig. 9
Measured impedance for several tumor (sphere) depth. a) Real part, b) Imaginary part.

The simulated measurement of Fig. 10 shows a parametric variation of the measurement line offset from the sphere center. The measurement line is at the skin level but has no alignment with the sphere. This result is crucial because the sphere location in practice is unknown, as its depth and radius. Note the impedance response has more significant values for small offset, and there is a curve shape deformation due to offset. This study result is significant because it indicates a characteristic for locating the anomaly based on the curve shape.

Fig. 10
Measured impedance for several measurement line offsets. a) Real part, b) Imaginary part.

One can observe by results analysis that the impedance has smaller magnitudes for a contrasting response in comparison with the solution without anomaly because of the higher conductivity and permittivity of the cancerous tissue. The response is a function of the depth and size of the tumor and the electrode configuration. Higher electrode spacings are better for more profound readings and a smaller spacing for shallower readings. The results also presented the deformation of measurement lines due to the electrodes offset, which is an excellent indicator of the anomaly’s location.

IV. Conclusions

The paper presents a model for spherical anomalies in a half-space of uniform electrical properties. It admits a known size and location for a sphere and calculates the potential function for points within the interest region. There is a frequency limitation for the model once the authors neglect the propagation effects. A point current source’s potential function is generalized for any position on the air-skin interface, and the superposition principle is then applied for a medium of linear electric properties reaching the impedance solution for a four-electrode configuration. As illustrated in the parametric studies, there is a contrasting response for a shallower and bigger tumor, which results in a lower impedance reading. The electrode spacing influences the depth of tissue fields, which are higher for a larger spacing. The offset between the measurement line and the anomaly implies a deformed impedance profile with smaller amplitudes. This paper’s model resulted in good accordance with the FEA, which is an indicator of its accuracy. Our model excels the FEA because it has an analytical solution and needs no domain discretization. The authors believe the model is a promising application for Electrical Impedance Tomography as an alternative approach. Further development is an appropriate interpretation of the measurements (bidimensional array impedance) to generate medical images for determining a possible tumor location.

Acknowledgments

Authors would like to thank Altair Engineering inc. for kindly supporting the Laboratory of Applied Electromagnetics with their software packages, particularly the Altair Flux used in this paper.

References

• ^[1]
  B. Brandstätter, H. Scharfetter, and C. Magele, “Multi frequency electrical impedance tomography,” COMPEL-The international journal for computation and mathematics in electrical and electronic engineering, 2001.
• ^[2]
  F. Dalvi-Garcia, M. N. d. Souza, and A. V. Pino, “Algoritmo de reconstrução de imagens para um sistema de tomografia por impedância elétrica (tie) baseado em configuração multiterminais,” Revista Brasileira de Engenharia Biomédica, vol. 29, pp. 133–143, Junho 2013.
• ^[3]
  J. L. Mueller, D. Isaacson, and J. C. Newell, “A reconstruction algorithm for electrical impedance tomography data collected on rectangular electrode arrays,” IEEE Transactions on Biomedical Engineering, vol. 46, no. 11, pp. 1379–1386, Nov 1999.
• ^[4]
  J. D. Bronzino, The Biomedical Engineering Handbook, 2nd ed. CRC Press, Dec. 1999, vol. 1.
• ^[5]
  D. Griffiths, Introduction to Electrodynamics. Cambridge University Press, 2017.
• ^[6]
  J. Wait, Geo-electromagnetism. Academic Press, 1982.
• ^[7]
  X. Zhang, W. Wang, G. Sze, D. Barber, and C. Chatwin, “An image reconstruction algorithm for 3-d electrical impedance mammography,” IEEE Transactions on Medical Imaging, vol. 33, no. 12, pp. 2223–2241, Dec 2014.
• ^[8]
  J. R. Cardoso, Engenharia Eletromagnética, 1st ed. Elsevier Editora Ltda., 2010.
• ^[9]
  J. R. Cardoso, Electromagnetics through the Finite Element Method: A Simplified Approach Using Maxwell’s Equations. Boca Raton: CRC Press, 2017.
• ^[10]
  Altair Engineering, Altair Flux™ 2020.1 User Guide, 2020.
• ^[11]
  F. M. Cheng Y, “Dielectric properties for non-invasive detection of normal, benign, and malignant breast tissues using microwave theories,” Thorac Cancer, vol. 9, no. 4, pp. 459–465, April 2018.

Publication Dates

• Publication in this collection
  29 Sept 2021
• Date of issue
  Sept 2021

History

• Received
  21 Jan 2021
• Reviewed
  21 Jan 2021
• Accepted
  29 Apr 2021