On elliptic partial differential equations in bioimpedance

This paper deals with mathematical models in electrical bioimpedance fields that are described by means of elliptic partial differential equations (PDEs). To find solutions that have practical significance and value, it is necessary to gain a deep understanding of the underlying physical phenomena with the parameter details of PDE models as well as the data acquisition systems. This paper describes practical details of boundary conditions and effective coefficients of elliptic PDEs occurring in bioimpedance.


Introduction
Elliptic partial differential equations (PDEs) are used to describe a large variety of physical phenomena. This paper focuses on PDE models in electrical bioimpedance area, which can be derived from suitable arrangements of Maxwell's equations at frequency below 1 MHz. Bioelectrical impedance methods are low cost, noninvasive and portable techniques for real-time clinical status monitoring (e.g. cardiopulmonary monitoring for patients with sleep apnea/obesity hypoventilation syndrome) [29,57], estimating body composition (e.g. body fat and muscle mass) [25,36,37], 1 3 and others. Hence, it has been one of important research topics during the past four decades.
Bioelectrical impedance imaging is based on Ohm's law = A , where A = (a ij ) is 3 × 3 matrix representing an electrical admittivity tensor, is the electric field, and is the current density. To generate and into a biological object occupying a three dimensional domain ⊂ ℝ 3 (with its diameter of less than 2 m), we inject an electrical current by attaching several surface electrodes on the boundary . Since ∇ × ≈ 0 at low frequency (below 1 MHz), the resulting electric field is = −∇u where u ∈ H 1 ( ) is the electrical potential governed by the following elliptic PDE: where = (x 1 , x 2 , x 3 ) denotes position vector, = (n 1 , n 2 , n 3 ) is the unit outer normal vector to the boundary , and g ∈ H −1∕2 ( ) represents Neumann boundary data corresponding to the injected current. Note that g must satisfy ∫ gds = 0 for the existence and u is usually normalized by ∫ uds = 0 for the uniqueness, where ds denotes the surface element.
The coefficient A should be understood as a homogenized admittivity tensor that depends on scale, position, cell structure including molecular compositions of cells, shape and direction of cells, cellular membranes, intra-and extra-cellular fluids, concentrations and mobilities of ions [60]. For the existence and uniqueness of a weak solution according to Lax-Milgram theorem [11], the coefficient A in must satisfy where min and max are positive constants. For stability, the ratio max ∕ min should not be too large. The admittivity tensor A can be derived by asymptotic homogenization techniques with knowledge of pointwise admittivity distribution [3,18,19,44]. However, the pointwise admittivity may not exist in real world, because the admittivity is passive property [9]. What we can observe is its effective or equivalent property [14,43]. The coefficient A in muscle region usually exhibits a strong anisotropy depending on muscle fiber orientation and scale [10]. The coefficient A of tissue is affected by functional and pathological conditions such as ischemia, hemorrhage, edema, inflammation, cancers, and neural activities [60]. Hence, the problem of estimating A of tissue has been an important research topic in biomedical field [25,29,57,58].
In practice, we cannot know the Neumann data g in (1) accurately. The g is determined by the configuration of the array of surface electrodes (that are attached on ) and the current injection pattern. If we use a pair of surface electrodes to inject a dc current of I mA into the subject , the induced electrical field −∇u inside roughly satisfies the following boundary conditions: where E + and E − denote the skin-electrode contact area as shown in Fig. 1. Here, the skin-electrode contact impedance [61] is ignored for mathematical simplicity. Note that the corresponding Neumann data g = ∇u ⋅ is living in the Sobolev space H −1∕2 ( ) and g is unbounded due to its singularity on the edge of E ± . This boundary condition will be discussed in detail.
We need to choose the domain such that the boundary value problem (1) is well-posed in the sense of Hadamard [26]. Consequently, it is necessary to exclude regions with almost extreme (zero or infinite) admittivity from . For example, the air region in stomach has almost zero admittivity and electrodes has almost infinite admittivity. We should note that the Neumann data g is very limited depending on the placement and number of electrodes. In order to apply rapidly oscillating boundary data g, we need to attach many electrodes, resulting in the electrode-skin contact area occupying a large portion of . In this case, unfortunately, most of electrical current flows significantly through the highly conducting electrodes. Hence, in the case when the electrode-skin contact area occupying a large portion of , it would be desirable to include electrode regions to for accurate analysis. There was often a misunderstanding of this issue with respect to and g in the both mathematics and engineering community. This paper describes this issue in detail.
Understanding the practical limitations and fundamental physical phenomena mentioned above is important for solving various problems in bioimpedance including electrical impedance tomography (EIT) and electric impedance myography(EIM). EIT and EIM are to estimate the electrical tissue property A in (1) of biological subject at various frequency below 1 MHz. With some ideal assumptions, the mathematical formulation of EIT system can be described as the well known Calderón problem which aims to identify an equivalent isotropic conductivity A from the Neuman-to-Dirichlet map  (1). The Neumann boundary data is determined by the injection current applied using the pair of electrodes E + and E − value problem (1). There have been obtained novel theoretical results guaranteeing a unique identification of isotropic A from the NtD data [4,8,33,35,45,46,62]. The analysis of PDE models often plays an important role in the achievement of major advances in these areas. However, it should not be studied in isolation, and practical limitations associated with the measurement sensitivity and specificity, noise, interface between and instrument and so on must be properly understood and analyzed [57].
This paper is organized as following. Section 2 is devoted to explain elliptic PDE as a governing equation in bioimpedance and effective coefficient A of the PDE. In Sect. 3, we introduce inverse problems of finding coefficient of the PDE; electrical impedance myography and electrical impedance tomography. Finally, we conclude the manuscript with discussion in Sect. 4.

Elliptic PDE in bioimpedance
The admittivity distribution A in an imaging domain (a part of human body) is the passive electrical property which can be estimated by Ohm's law, using the relationship between the Neumann data g and the Dirichlet data u| with u being the solution of the elliptic PDE (1). This boundary value problem is derived from Maxwell's equations [57]. In order to measure A, we inject a sinusoidal current I sin( t) between a pair of electrodes ( E + and E − ) to produce the current density ℜ( e i t ) and the electric field ℜ( e i t ) inside . Here, i = √ −1 , ℜ(f ) is the real part of f, is angular frequency, t is time, is the time-harmonic electric field, and is the timeharmonic current density. Throughout this paper, we assume that the boundary is smooth, electrodes E ± are rectangular shape as shown in Fig. 1, the diameter of is less than 1 m, and the applied frequency ∕2 is below 1 MHz.
For the given injection current, the corresponding Neumann data g = ⋅ | satisfies From the time-harmonic Maxwell's equation, we have where is the time-harmonic magnetic field and 0 = 4 × 10 −7 H∕m is the magnetic permeability of free space. The human body is assumed to be the constant 0 . Assuming that 0 ≈ 0 at frequency below 1 MHz, there exists a complex potential u such that −∇u ≈ . Since ∇ ⋅ (∇ × ) = 0 , it follows from (5) that u satisfies the elliptic PDE (1) with g being ⋅ | .
Note that A can be expressed as A = + i , where and are the conductivity and permittivity, respectively. The A = A( , ) depends not only on the position On elliptic partial differential equations in bioimpedance = (x 1 , x 2 , x 3 ) but also on the angular frequency . The effective admittivity tensor A (macroscopic scale) in satisfies the usual ellipticity condition: where c is a positive constant being away from zero. We should note that A in a very microscopic scale may not satisfy the ellipticity condition (6) with c ⪊ 0 . It is because cell membranes are almost electrically insulating ( A ≈ 0 ) at = 0.
For the ease of explanation, we focus on A at dc current ( = 0 ). Let us understand the Neumann data g and = A∇u on the electrode area E + ∪ E − . Since electrodes are highly conductive materials such as copper and carbon, we have × | E + ∪E − ≈ 0 , under the assumption that the contact impedance underneath the electrode is ignored. Hence, the complex potential u approximately satisfies Setting ∫ u ds = 0 , we obtain a unique solution u ∈ H 1 ( ).
On the edge of E ± , there exists a singularity of g. The following observation explains the edge singularity on the electrodes.
The proof is based on the square root edge singularity [23,51] on the electrodes. For the ease of explanation, assume that A is the identity matrix and = ℝ 3 − = { ∶ x 3 < 0} (the lower half space). We also assume that both E + and E − are rectangular domains lying on . Let w be the extension of u such that . Moreover, w| E + and w| E − are constants. Hence, ∇w has the square root edge singularity on E + ∪ E − , resulting in Now, it remains to prove ‖∇u‖ L 2 (E + ∪E − ) = ∞ . It can be proven the well-known Rellich type identity [50]. For any smooth vector field F = (F 1 , F 2 , F 3 ) , we have Integrating the above identity over , we obtain is used to express a term bounded by ‖∇u‖ 2 L 2 ( ) . The identity (10) follows from Since the smooth vector field F is arbitrary, it follows from (11) This proves Observation 2.1 for the case where A is the identity matrix and = { ∶ x 3 < 0} (the lower half space). The general case can be similarly proved by a proper modification. Next, we explain how the potential u changes as additional electrodes ( E 1 , E 2 , ⋯ , E n ) are attached between E + and E − . See Fig. 2. Let all electrodes be the same size and shape, and let the spacing between electrodes be a constant s n . We denote The following observation shows that the boundary value problem (7) changes as the additional electrodes n increase.  (7)) is almost the same as the middle image because s n (the spacing between electrodes described in Observation 2.2) is not so small. However, the right image is very different from the left image because s n ≈ 0 1 3 On elliptic partial differential equations in bioimpedance Observation 2.2 Let u n be the potential corresponding to u in (7) with adding electrodes ( E 1 , E 2 , ⋯ , E n ) between E + and E − , as shown in Fig.2. The potential u and u n satisfy the same first two conditions in (7) (because of the same A, , and the same injection current), but the last two boundary conditions are different. The last boundary conditions in (7) should be replaced by: Moreover, as s n → 0 (the spacing between electrodes converges to zero), the potential difference u n | E + − u n | E − converges to zero and where is any positive constant and dist ( , ) is the distance between and .
In electrical impedance tomography, several researchers have tried to use as many electrodes as possible in order to get a rough version of Neumann-to-Dirichlet map [8,35]. However, if the spacing between electrodes is very small due to many electrodes, it is almost impossible to get any information of A in a remote internal region from the boundary, according to (14) in Observation 2.2.
For the well-posedness of the boundary value problem (1) from Lax-Milgram theorem, we usually exclude very high and low conducting region such as electrode and air to satisfy the condition (6). Indeed, in the case when a small number of electrodes are attached on , solution u of (1) can be approximated by that of (1) with excluding electrode region. However if electrodes almost surround the object then the approximation fails, because electric currents mostly flow boundary as shown in Fig. 2.
For the determination of A, according to the Calderón problem [8], it is necessary to have full boundary data which requires closely packed electrodes on the boundary. Most of practical studies (e.g. EIT, EIM) of finding A in (1) from a partial information of Neumann-to-Dirichlet map ∶ g ↦ u| use less than or equal to 16 small electrodes.

The effective coefficient A of heterogeneous media and homogenization
The effective admittivity tensor A in a cubic voxel Q centered at a position can be determined by Ohm's law: where and are the time-harmonic electric field and current density, respectively, in (5). Since the time-harmonic fields change with the angular frequency , the effective admittivity tensor A of biological tissue changes with . Hence, A = A( , ) depends on position , the applied frequency , and the size of the voxel. The effective tensor A is decomposed into its real and imaginary parts: The effective coefficient A can be derived from the pointwise admittivity, denoted by A pt = pt + i pt , on the basis of the two-scale homogenization theory [2]. The frequency dependency behavior of A is owing to Maxwell-Wagner polarization effect [20], that is related to the geometry of biological tissue structure with cells, extracellular matrix, intra-and extra-cellular fluids [21,24]. The frequency-dependent A is mainly influenced by the cell membranes with their thickness being in the order of several nm [34]. Let us understand intuitively how the effective conductivity A is affected by the tissue structure. At low frequency, cell membranes are insulators and current can flow around them. As cell swelling occurs, the extra-cellular space within the region decreases and this results in a reduce effective conductivity. If cells are densely packed, the effective conductivity becomes smaller for the same reason. If cells are tightly packed along the horizontal direction and the gaps between cells in the vertical direction are bigger, then the effective conductivity will be different depending on the direction of the applied electric field. This means that we should consider the anisotropy as well as the inhomogeneity in the effective conductivity.
In 1924, Fricke [18,19] described a mathematical form of the homogenized coefficient in the simplest case of dilute single suspension. Ammari et al [3] analyzed the role of the membrane in terms of the frequency-dependent behavior of a solution u in the framework of the elliptic PDE in two dimensional case where cell membranes are immersed in a domain . For a rigorous analysis, they assume the followings: (1) is divided periodically in each direction in identical very small squares and a cell lives in each small squares. (2) The medium outside the cell membranes is a homogeneous isotropic medium with A 0 pt = 0 pt + i 0 pt and each thin membrane has the uniform thickness of d ≈ 0 with the isotropic admittivity A m pt = m pt + i m pt with m pt 0 pt ≈ 0 . Fig. 3 illustrates the above assumptions more clearly, where the normalized square [0, 1] 2 contains the reference cell membrane given by { ∶ dist ( , ) < d} for a cell contour . Let us denote by c the area of inside , which inside is the domain inside the cell contour , as shown in Fig. 3. We assume that d∕c ≈ 0 . Given an applied sinusoidal current at the angular frequency , the resulting harmonic potential u in term of its pointwise coefficient A pt approximately satisfies Although the concept of effective admittivity has been studied deeply using homogenization concepts, its practical and intuitive definition is not clear due to the ideal assumptions of periodic and dilute suspension. Seo et al [59] provided a practical way to measure the effective admittivity A for a given cube from measurable boundary current-voltage data. Given a small cubic subject Q ⊂ , the effective admittivity A can be computed by its pointwise admittivity A pt , as in homogenization approaches. Since A( ) as a function of is regarded as an ensemble average of pointwise admittivity A pt over the cube Q, the A must satisfy where u j is the solution of (17) Here, we apply a sinusoidal current of I mA at angular frequency . For each pair (j, k) ∈ {1, 2, 3} 2 , we denote the voltage difference between two sides of the cube Q: Then, the effective A in the cube Q can be evaluated as in Fig 4.

Estimation of A in electrical bioimpedance
In bioimpedance, the simplest way of measuring an overall version of the effective admittivity A of the biological subject (occupying a three dimensional domain ) is the four-electrode method [24], as shown in Fig. 5; the pair of external electrodes (denoted by E 1 + and E 1 − ) is used to inject a sinusoidal current of I mA at angular frequency and the induced voltage difference is measured the pair of inner electrodes (denoted by E 2 + and E 2 − ). et u 1 denote the induced potential due to the injection current using E 1 + and E 1 − . Then u 1 is governed by On elliptic partial differential equations in bioimpedance In the above model, skin-electrode contact impedance [61] is ignored for simplicity. The principle of reciprocity invented by Helmholtz [28] is a very useful tool in bioimpedance. The reciprocity principle says that the voltage difference of u 1 between the pair of inner electrodes is the same as the voltage difference of u 2 between the pair of outer electrodes, where u 2 is a solution of (23) with interchanging subindex 1 and 2.
Theorem 1 Let u 1 be a solution of (23) and u 2 be a solution of (23) with interchanging subindex 1 and 2. Then, we have the following reciprocity principle: Proof For simplicity, we assume that all electrodes are the same size. The reciprocity principle can be shown by the following integration by parts:

Electric impedance myography (EIM)
Electrical impedance myography (EIM) is a technique for the assessment of muscle admittivity A in order to evaluating neuromuscular diseases both for their diagnosis and for their ongoing assessment of the progress or therapeutic intervention [49,52,54]. Since muscle cells are composed of muscle fibers, they have anisotropic properties in which admittivity A varies depending on the direction of the muscle fibers [16,17]. So muscle admittivity A has two aspects; along and across the muscle fiber. When the electric field is formed along the muscle fiber direction, the admittivity of the muscle that affect the formation of the electromagnetic field is called longitudinal admittivity and denoted by A longi ∶= longi + i 0 longi with vacuum permittivity 0 = 8.85 × 10 −12 F ⋅ m −1 , longitudinal conductivity longi and relative permittivity longi . Similarly, we call transverse admittvity A trans ∶= trans + i 0 trans , conductivity trans , and relative permittivity trans , in the case when electric field is formed across the muscle fiber. Thereby, the admittivity of muscle is anisotropy and it can be written as matrix, for example when muscle fiber lies along x 1 -axis, then a 11 = A longi , a 22 = A trans , a 33 = A trans , and remaining a ij = 0 (Fig. 6). According to the [22] (in vivo study of rat muscle), longitudinal conductivity is 10 times higher than that of transverse, i.e. longi ≈ 10 trans . Note that the coefficient A is complex anisotropy. So, the real and imaginary parts of u [solution of (23)] gives system of equations. If anisotropic ratios of real and imaginary parts of A are same, then the system of equations can be separated.
Observation 3.1 Let u be a solution of (7) and for simplicity muscle fiber assume to be on x 1 -axis. Define the anisotropy ratio of real and imaginary parts of A as 2 ∶= trans ∕ longi and 2 ∶= trans ∕ longi , respectively. In general from the governing equation ∇ ⋅ (( + i )∇(v + iw)) = 0 , we have where = ℜ(A) , = ℑ(A) , and u = v + iw . Then, we have a system of equations as below where = trans ∕ trans . If 2 = 2 , the above system of equations can be separated as

3
On elliptic partial differential equations in bioimpedance In order to estimate A, EIM uses set of voltage differences with dividing current amplitude: which is called impedance with same notations in the Theorem 1. Due to the anisotropic property of A, it is required to have at least two impedances, for example impedance measured on longitudinal and transverse directions. The measured impedance is a some sense of average of both longitudinal and transverse admittivity.

Observation 3.2 For the ease of explanation, we consider the half space domain
, same anisotropy ratio 2 = 2 , and four linearly aligned electrodes with spacing a, b, and a. When the electrodes alignment angle is and muscle fiber direction is , the measured impedance Z can be represented by [53] where the geometry factor F geom = F domain ∕F eled , the domain factor F domain = 2 (or F domain = 4 for = ℝ 3 ), and the electrode factor F eled ∶= 2b a(a+b) .
Hence the inverse problem of EIM is to estimate A trans and A longi from a set of measured impedance of several measuring angle [38,39].

Electrical impedance tomography (EIT)
Electrical impedance tomography (EIT) aims to visualize the distribution of (isotropic) coefficient A with meaning that A is human body using finite number of surface electrodes. One of most successful clinical application in EIT is lung ventilation monitoring, i.e. A is lung (Fig. 7). For the reconstruction of the distribution of A, a set of Neumann-to-Dirichlet map is constructed with several choices of current injection electrodes and corresponding voltage measuring electrodes for each pair of current injection electrodes. For the sake of simple explanation, let {E 1 , E 2 , … , E 16 } be a set of surface electrodes attached on the boundary of human body with adjacent pairs of current and voltage electrodes. Then, we have the following Neumann-to-Dirichlet data: 16,14 where z i,k is skin-electrode contact impedance [61]. Hence, the inverse problem of EIT is to reconstruct (equivalent isotropic admittivity) A from Neumann-to-Dirichlet data in (26). There are several EIT algorithm as in [1,5,27,30,32,42].

Discussions and conclusions
In medical imaging, mathematical methodologies have evolved to improve our ability to visualize various physical phenomena and features accurately and reliably. The analysis of PDE models often plays an important role in the achievement of major advances in these areas. Technical advances have been followed by theoretical progress aimed at understanding the solution's structure. However, the subject of partial differential equations should not be studied in isolation, because much intuition comes from a thorough understanding of applications.
Let us briefly discuss the inverse problem of reconstructing the equivalent isotropic conductivity of A. In spite of novel theoretical results guaranteeing a unique identification of A from the Neumann-to-Dirichlet data [8,45,62], numerous experiences have shown that static EIT for recovering isotropic admittivity A has fundamental drawbacks because EIT data depends strongly on the boundary geometry and electrode positions, whereas it is much less sensitive to a local perturbation of isotropic A away from the measuring electrodes [56]. In 2003, Magnetic Resonance Electrical Impedance Tomography (MREIT) was invented to deal with the wellknown ill-posedness of the image reconstruction problem of EIT [55]. Mathematicsoriented research overcomes technical barriers in electrical tissue property imaging. This paper observed that ∇ ln A ⋅ (A∇u × (0, 0, 1)) with u being a solution of (7) is During lung ventilation, pulmonary conductivity A changes due to conversion of air in the alveoli measurable quantity using MRI system so that it can probe changes in the logarithm of the isotropic A along any equipotential curve in each imaging slice. Here, (0, 0, 1) is the direction of the main magnetic field of MRI scanner. This method uses two different Neumann data g 1 , g 2 so that the area of the parallelogram made by these two vector fields ∇u 1 and ∇u 2 is non-zero at every position in the imaging slice. Taking advantage of these mathematical observations, they found a representation formula for the equivalent isotropic A which can offer state-of-the-art conductivity imaging using MRI animal and human experiments [56,60,63]. There have been developed numerous bioimpedance-based health care systems such as body fat assessment [31], stroke volume monitoring, hypoventilation monitoring [7], lung monitoring [15], and so on. All of these are tied to the effective property of A, that is the coefficient of the elliptic PDE. Recently, bioimpedance techniques are used for non-destructive continuous monitoring of in vitro chondrogenesis for the production of high-performance engineered cartilage of clinical use and its quality control [41,47].
Impedance imaging techniques have been used various areas beside the biomedical area such as ground monitoring [13,40,48] and industrial process monitoring [6]. Chipot et al [12] developed a pressure-sensitive conductive fabric sensor, which is based on a design of a composite fabric consisting of an electrically conductive yarn and a sponge-like non-conductive fabric with high pore density. In this model, the conductive yarn is woven in a wavy pattern to possess a pressure-sensitive conductive property (being capable of changing its effective electrical property of the admittivity A due to an applied pressure, in the sense of homogenization theory).
Developing mathematical models with practical significance and value requires the fusing of the knowledge and techniques of traditional engineering fields with pure and applied mathematics [57]. It is necessary to understand practical limitations imposed by measurement methods in mathematical models.