Analysis of a Mechanistic Model for Non-invasive Bioimpedance of Intact Skin

We start with the mathematical model for EIS measurements derived by Birgersson et al. [12], which treats the skin as a three-layer entity: stratum corneum (SC) as the outermost layer characterized by corneocytes; viable skin (VS) consisting of both living epidermis and dermis; and adipose tissue (AT) to account for the subcutis layer containing a high percentage of fat, as illustrated in Fig. 1a-b While the non-invasive EIS probe in the studies by Birgersson et al. [12] comprised four electrodes – one current detector (sense), one guard and two injects – we will here consider a simpler probe with similar dimensions but only one current detector (at I) and one voltage injector (at II). The dimensions and material properties of the electrodes and skin layers can be found in Table 1.

Fig.1

The mathematical model can be written as

Here, the complex-valued phasors, Φ(r) and J(r) , for the potential and current density are functions of space only (although we note that they can be pulled back into the time domain), _n is the unit normal vector pointing outwards to any given boundary, V₀ is the applied voltage, and the location of the boundary conditions is given by roman numerals, see Fig. 1 for their location; the effective conductivity, σeffk,$\sigma _{eff}^{k},$is given by

where σ ^k is the electric conductivity, ε₀ and εrk$\varepsilon _{r}^{k}$are the relative permittivities of vacuum and the material respectively, i is the imaginary unit, and ω = 2πν is the angular frequency. The superscript, k, denotes either viable skin, stratum corneum soaked with a saline solution with 0.9% NaCl for 1 min, adipose tissue, or the electrode.

The material properties for the three skin layers on the volar forearm are defined as follows [17]:

Here, cjk$\mathfrak{c}_{j}^{k}$and djk$\mathfrak{d}_{j}^{k}$are the material coefficients, given in Table 2; v˜=v/(1Hz$\widetilde{v}={v}/{\left( \text{1}\,\text{Hz} \right.}\;$is a dimensionless frequency and N is the uppercase greek letter of ν , introduced for notational convenience. These constitutive relations for the electrical material properties are valid between 1 kHz and 1 MHz; see Birgersson et al. [17] for the validation of these relations for two studies comprising 26 and 120 volunteers respectively.

The impedance is defined as

with the total current measured at the sense, boundary I in Fig. 1b, defined as

We will refer to this model as the complete model.

We define the following dimensionless variables and dimensionless numbers:

Here, [ ] is the scale (absolute value), R = w₁ + w₂ + w₃ is the overall radius of the outer electrode (inject) and Asen=πw12${{A}^{sen}}=\pi w_{1}^{2}$is the area of the inner electrode (sense); we have taken the scale for the impedance as [Z]=V0/([J]Asen).$\left[ Z \right]={{{V}_{0}}}/{\left( \left[ J \right]{{A}^{sen}} \right)}\;.$The dimensionless number, Π^k , is the ratio of the displacement and ohmic current; Θ^k relates the height of a given skin layer to the overall length in the radial direction; Λ₁ arises in the non-dimensionalization of the total current measured at the sense; and Ξk${{\Xi }^{k}}$gives the ratio of the current flowing through skin layer k to the current flowing through stratum corneum in the z -direction – and will provide valuable quantification of where the currents flow later in the analysis. We note that at this moment, the scale for the current, [J], is unknown.

These dimensionless variables and numbers are substituted into the governing equations for each layer of the skin, which we shall consider separately. Between the different layers, we require continuity of the current density as well as potential – we do not explicitly state those boundary conditions in the coming analysis for the sake of brevity.

In the stratum corneum, we obtain the following dimensionless transport equations:

where

Here, Π^SC ~ 1, as can be inferred from Fig. 2, whence we can identify the current density scale, [J], from Eqs. 16 and 17 as

Fig.2

Substituting the current density scale back into the governing equations, Eqs. 15, 16, 17, gives

where

Furthermore, the stratum corneum is slender with

R ~10^-2 m and h^SC ~10^-5 m, such that

and

In other words, the current density in the radial direction is around a thousand times smaller than the current in the z-direction for the conditions and probe dimensions considered here, whence Eq. 18 reduces to

subject to the boundary conditions adjacent to the electrodes

and a continuous potential at the interface with the viable skin (N.B.: the potential drop in the electrodes is negligible at leading order as we shall show later). The current between the inject and the sense in the radial direction thus has to occur in the viable skin and/or adipose tissue, because the voltage is not enough to drive the current through the SC in the radial direction.

Returning to the current density scale for the frequency range of 1 kHz to 1 MHz that is of interest here, we find that

which provides the scale for the current density in the z-direction for the base-case parameters and conditions; in the radial direction, the scale for current density is [J]ΘSC∼10−5−10−2Am-2$\left[ J \right]{{\Theta }^{SC}}\sim {{10}^{-5}}-{{10}^{-2}}\,\,\text{A}\,{{\text{m}}^{\text{-2}}}$We verify these scales with the numerical solution in the stratum corneum for 1 kHz to 1 MHz, which gives |Jr|∼10−6−10−2Am-2,$\left| {{J}_{r}} \right|\sim {{10}^{-6}}-{{10}^{-2}}\,\text{A}\,{{\text{m}}^{\text{-2}}},$evaluated at r=w1+w2/2,0<z<hSC,$r={{w}_{1}}+{{{w}_{2}}}/{2}\;,\,\,\,0<z<{{h}^{SC}},$which is between the two electrodes; and |Jz|∼10−2−10Am-2,$\left| {{J}_{z}} \right|\sim {{10}^{-2}}-10\,\,\text{A}\,{{\text{m}}^{\text{-2}}},$evaluated at 0 < r < w₁, z = h^SC , which is underneath the sense.

Finally, we confirm that scaling the potential, Φ˜=Φ/V0,$\widetilde{\Phi }={\Phi }/{{{V}_{0}},}\;$in the stratum corneum with V₀ in the entire frequency range is correct by comparing with the numerical simulations of the absolute voltage drop underneath the sense (r=w1/2,0≤z≤hSC)$\left( r={{{w}_{1}}}/{2,0\le z\le {{h}^{SC}}}\; \right)$and the inject (r=w1+w2+w3/2,0≤z≤hSC),$\left( r={{{w}_{1}}+{{w}_{2}}+{{w}_{3}}}/{2,0\le z\le {{h}^{SC}}}\; \right),$which are around 0.033 – 0.024 V and 0.017 – 0.013 V respectively for 1 kHz to 1 MHz. These voltage drops are indeed O(V₀ ). The higher

potential drop underneath the sense compared to the inject originates from the fact that the area of the inject is roughly two times larger than the sense for the base case conditions.

For the viable skin,

where

and

In order for the current to be carried through the viable skin between the inject and the sense, we require that ΞVS∼${{\Xi }^{VS}}\sim $1 in Eqs. 28 and 29 for both the left-hand and righthand sides to be O(1), which is not the case since

as illustrated in Figure 3.

Fig.3

This in turn suggests that we need to revisit the potential drop in the viable skin, which we assumed to be ~ V₀ when we scaled Φ = Φ / V₀ . If we rescale the potential,

with ΔΦ^VS denoting an as-of-yet unknown scale for the potential drop, Eqs. 27, 28, 29 take the form

with

Again, we require that J˜r∼J˜z∼(N.B.ΘVS∼1)${{\widetilde{J}}_{r}}\sim {{\widetilde{J}}_{z}}\sim \left( \text{N}\text{.B}\text{.}\,\,\,\,{{\Theta }^{VS}}\sim 1 \right)$for the current to be carried through the viable skin, i.e.

which gives the scale for the potential drop,

For the base case,

which agrees well with the numerical solutions of the complete model (not shown here).

For the adipose tissue,

where

and

Now, because

the current in the adipose tissue is mainly ohmic, such that we can write

For the current to be carried through the adipose tissue, we require that ΞAT∼1;${{\Xi }^{AT}}\sim 1;$returning to Fig. 3, however, we see that ΞAT∼20−2×10−2${{\Xi }^{AT}}\sim 20-2\times {{10}^{-2}}$at 1 kHz and 1 MHz respectively.

Rescaling the potential drop like what we did for the viable skin,

and requiring J˜r∼J˜z∼1(N.B.ΘAT∼1)${{\widetilde{J}}_{r}}\sim {{\widetilde{J}}_{z}}\sim 1\left( N.B.\,\,{{\Theta }^{AT}}\sim 1 \right)$i.e.

we find the scale for the potential drop in the adipose tissue,

With the base case parameters,

Clearly, the current from the probe cannot be carried through the adipose tissue at 1 MHz since it would require around forty times the applied voltage, V₀; at the other end of the considered spectrum, 1 kHz, the voltage drop in the adipose tissue is around ten times larger than that of the viable skin; i.e.

In other words, most of the current from the inject would choose the easier path through the viable skin before even going deeper through the viable skin to reach the adipose tissue, after which the current would need to pass through the adipose tissue in the radial direction and then upwards through the viable skin again.

In the electrodes, the dimensionless governing equations can be written as

where

Because Π^EL ~ 10^-15 -10^-12 ≪ 1, the displacement current is negligible compared to the ohmic counterpart in the electrodes; in addition, ΞEL~1011−108${{\Xi }^{EL}}\tilde{\ }{{10}^{11}}-{{10}^{8}}$for 1 kHz to 1 MHz, which implies that the voltage drop in the electrodes is negligible compared to the drop across the skin layers. There is thus no need to consider the potential drop in the electrodes.

Now that we have analyzed the currents in the electrodes and three skin layers in detail to deduce the leading-order phenomena, we can combine the findings to arrive at a reduced model, as illustrated in Fig. 1b-1f

In the first step from Fig. 1b to 1c, the electrodes do not need to be resolved due to the negligible potential drops in these; we can thus move the boundaries for the ground, I, and applied potential, II, down to the stratum corneum; i.e. at the sense (I),

and at the inject (II),

where R₁ = w₁, R₂ = w₁ + w₂ , R₃ = w₁ + w₂ + w₃ (= R) for notational convenience.

In the second and third steps from Fig. 1c to 1e, the leading-order transport in the stratum corneum for the frequencies considered here is in the z-direction, whence the PDE reduces to an ordinary differential equation, which we can integrate analytically and so reduce this layer to the modified boundary conditions, I’and II’ adjacent to the viable skin. We can capture these two steps in our reduced model by returning to Eq. 23 and writing it in dimensional form as

which gives us that J _z is constant in the z-direction in the stratum corneum (the radial component is O(Θ^SC ) ≪ 1 for the conditions we consider here); i.e., it is a function of r only. With the definition of the current density, we can write

which can be integrated once to give

With the boundary condition, Eq. 55, underneath the inject, we can determine the integration constant, C, and write the potential as

At the interface between the stratum corneum and viable skin (II'; z=0 in Eq. 59), we can now introduce the modified boundary condition in lieu of the stratum corneum and the inject as

where, again for convenience,

Returning to the sense (I) and the boundary condition, Eq. 54, we can repeat the integration above and write the modified boundary condition for the reduced model as

In the fourth step from Fig. 1e to 1f, the contribution of the adipose tissue to the impedance is negligible at leading order, whence we do not need to resolve this layer.

After these four steps, we can finally formulate the dimensional governing equations and boundary conditions, because it is not necessary to go through the entire nondimensionalization to employ the model:

We note that the initial boundary conditions – two Dirichlet and one Neuman condition – have been replaced by two Robin and one Neuman condition.

The reduced model – a Laplace equation with a Neuman and two Robin conditions – still requires a numerical method to solve. It is thus not as straight-forward to employ as equivalent circuits nor does it readily yield the functional form of the measured impedance and its dependence on the parameters. We therefore proceed with our analysis with a view to secure semi-analytical or analytical solutions (step v in Fig. 1). As we shall see, these solutions will be approximate in nature; i.e. we will make assumptions during the derivation.

Now, because we are primarily interested in the measured impedance (macroscopic, integral quantity) and not the pointwise potential distribution (microscopic quantity) throughout the electrodes and skin layers, we start by integrating the reduced model, Eqs. 63, 64, and invoking the divergence theorem,

where V is the volume and A is the outer area bounding the volume.

Introducing the boundary conditions, Eqs. 65, 67, into the integral formulation, Eq. 69, gives

where Ф^sen and Ф^inj are the area-averaged potentials at the sense and inject respectively:

Here, I is the total current and the areas of the sense and inject are

Let us define the overall potential drop across the viable skin as

and argue based on linearity that ΔV ^VS should be a function of the current density, the effective conductivity of the viable skin and a representative average length for the current flowing through the viable skin, L; i.e.

Assuming that we will later be able to quantify L, we now have four equations for our primary variables, Φ ^sen , Φ^inj , I, and Z:

After some algebra, we find the current and impedance as

with

Returning to the definition of S in Eq. 61, we can express C as

which acts as a dimensionless correction factor for the current and impedance; i.e. if ℭ≪1$\mathfrak{C}\ll 1$we obtain

In particular, if we take A^inj ~ R² and L - R , we see that

based on the potential drops that were defined in the earlier analysis for the skin layers. The condition ℭ≪1$\mathfrak{C}\ll 1$thus corresponds to the scenario when the impedance is dominated by the stratum corneum, ΔΦ^VS ≪ ΔΦ^SC , which is typically the case for frequencies around 1 kHz. Noteworthy is that the closed-form solution for the impedance in Eq. 86 is identical to our analytical solution for a four-electrode probe around 1 kHz [13] when only considering two electrodes.

Let us now try to determine L and its functional form. For this purpose, we will carry out a Hankel transform of zeroth order, H₀ { Φ ( r , z ) } , similar to Cheng et al. [21],

because it allows us to transform our axisymmetric Laplace equation in cylindrical coordinates (Eq. 63),

to an ordinary differential equation,

with the solution

where A(ξ ) and B(ξ ) are two functions that we will soon determine from the transformed boundary conditions; J₀ is the 0^th-order Bessel function of the first kind and Δ is the Laplace operator. At this stage, we note that our transformed solution automatically satisfies the circular symmetry and requires that

which is reasonable since we expect the measured impedance to mainly be determined by the current flowing between the two electrodes; in fact, we have in the original boundary condition, Eq. 67, artificially limited the computational domain by assuming insulation for r = R^skin and choosing R^skin sufficiently large not to affect the impedance significantly.

If we further assume that the current across each electrode is constant, we can simplify Eqs. 65 and 66 to

or, after the Hankel transform,

where we have used that the Hankel transform of zeroth order of a Heaviside function, U(a - r) , is defined as

with J₁ and a denoting the 1^st-order Bessel function of the first kind and a constant respectively. The remaining boundary condition, Eq. 67, becomes after transformation

Substituting the z-derivative of the solution, Eq. 89,

into Eq. 94 gives

Proceeding to the other boundary condition, Eq. 91, we find after substituting the solution thus far,

whence

and

With A(ξ) and B(ξ) in Eq. 89,

or

since

Taking the inverse Hankel transform,

we find the local potential distribution throughout the viable skin as

Note that the local potential, Υ(r, z) , at this stage is not the same as Φ(r, z) in our initial model with a Robin boundarycondition, Eqs. 65 and 66. We need to correct the potential distribution arising from the Hankel transform with the simplified boundary condition, Eq. 91, by shifting it with a constant, c; i.e.

The potential drop across the viable skin can now be found from

where

because

and

It is noteworthy that the integral on the right-hand-side (RHS) of Eq. 110 only depends on the lengths, R₁, R₂, R₃ and h^VS and not on the frequency. One thus only needs to solve the integral once for a given electrode design and thickness of the viable skin. Most importantly, we see that the functional form agrees with Eq. 76 and that

Substituting our expression for 1/ L into Eqs. 81, 82, 83 finally reveals the functional form for the current and impedance based on the Hankel transform (denoted with the subscript H) as

It is clear that this solution requires one numerical step for the integral, ℑ, in the correction factor, ℭ_H, whence we refer to it as a semi-analytical solution (a simple Matlab code can be found in Appendix A). With the current, I_H, we can also calculate the average potentials at the sense and inject:

Let us determine the constant c; this can be done by introducing a reference potential, e.g. Φ^sen from 117 into Eq. 106 together with Eq. 108, such that

whence

If we want to identify the contributions from the individual layers, we can write the impedance as

or in its dimensionless form as

with the dimensionless numbers

The first two terms on the RHS in Eqs. 121 and 122 quantify the impedance contribution of the current passing through the stratum corneum underneath the sense and inject and the third term as it passes through the viable skin. In particular, the dimensionless form of the impedance reveals that the full model for the measured EIS can be reduced down to four dimensionless numbers Λ ₂ , Λ ₃ , Π^SC and Π^VS .

We solve both the complete mathematical model and the reduced counterpart in the finite-element solver COMSOL Multiphysics 5.2a [22]. The axisymmetric two-dimensional geometries of the complete, see Fig. 1b, and reduced model, see Fig. 1f, were resolved with around 10⁴ and 500 mesh elements after a mesh-independent study to ensure mesh-independent solutions. The direct solver MUMPS was selected as the linear solver with a relative convergence tolerance of 10^-6; and a typical run for one frequency required around 1 s (wall-clock time) with around 10⁴ degrees of freedom on a quad-core 3.4 GHz workstation with 16 GB RAM for the complete model; the reduced model with 10³ degrees of freedom was solved almost instantaneously (reported as 0 s in COMSOL). Additionally, the approximate semi-analytical and analytical solutions were calculated and post-processed in Matlab R2015 [23].

At the macroscopic level, we have three expressions for the impedance:

We compare their magnitudes and phases with the numerical solution of the complete model, Eqs. 1, 2, 3, 4, 5, for 1 – 10³ kHz at base-case conditions (Table 1) in Fig. 4. Here, several features are apparent. First and foremost is that the reduced model, Z_red, can predict the magnitude and phase of the EIS in the considered frequency range: the maximum relative error is less than 2% and 0.1% for the magnitude and phase respectively. For the base-case parameters, the key dimensionless number Θ^SC = h^SC / R ~ 10^-3 ≪1 – the key condition for the reduced model to be valid. This condition should be fulfilled for most probe designs if the length scale, R ≪ h^SC ~ 10^-5 m. Second is that the magnitude and phase of the impedance predicted by the approximate solution from the Hankel transform, Z_H, agrees well throughout the frequency range with a maximum relative error of around 2% and 1% for the magnitude and phase respectively. It agrees well despite our approximation of the current being constant underneath each electrode, which raises the question as to why this is.

Fig.4

The answer can be found in the stratum corneum, which acts as a significant resistance to the current and so evens out the current density profile underneath the electrodes – this is in particular the case around 1 kHz where ΞVS≫1.${{\Xi }^{VS}}\gg 1.$Third is the good agreement of Z₁ around 1 kHz with maximum relative errors of 0.8% and 0.3% for the magnitude and phase respectively when ℭ≪1;$\mathfrak{C}\ll 1;$and the large deviation of 22% and 34% around 1 MHz when ℭ~1.$\mathfrak{C}\tilde{\ }1.$We can thus simply use the closed-form expression by Z₁ around 1 kHz, but need to resort to either Z_red or Z_H for frequencies around 10⁴ Hz and higher.

Parameters

To further test the fidelity of the reduced model, Z_red, and approximate solution, Z_H, let us compare it with the full set of equations for a set of parameters in Table 3 that correspond to different electrode designs with smaller and larger electrodes – w₁, w₂ and w₃ – and thicknesses of stratum corneum, h^SC, and viable skin, h^VS, and frequencies, ν . All permutations of these parameters are considered for a total of 7200 parameter combinations. For the two skin layer thicknesses, we have taken the mean values ± two standard deviations [12]; i.e. 14±6 μm for stratum corneum and 1.2 ± 0.4 mm for viable skin.

Starting with the reduced model in Fig. 5a-b, we see that it captures all the 7200 permutations of the parameters well with respect to the full model – the maximum relative error is 3% and 0.1% for the magnitude and phase respectively. For the approximate solution based on the Hankel transform, Z_H, the picture is slightly different: It agrees well in terms of the predicted magnitude with a maximum relative error of 3% for all parameter combinations, but the error for its predicted phase increases for smaller angles and reaches around 24% for negative phase values of around 20^o.

Fig.5

A closer inspection reveals that the specific parameter combinations for these predicted low negative phase values of around 20^o are as follows: 1 MHz frequency, where the width of the sense, w₁, is typically 3 mm, the width of the inject, w₃, is 3 mm and the spacing between the electrodes, w₂, is 1.2 - 2 mm, whilst h^SC is either 8 or 1.4 μm (but not at the higher thickness of 2 μm) and h^VS is at the lowest thickness of 0.8 mm. The key factors here are the large size of the electrodes (the inject is three times wider than in the base case) and the frequency of 1 MHz, for which ∂Φ/∂z|z=0${{\left. {\partial \Phi }/{\partial z}\; \right|}_{z=0}}$varies significantly across the electrodes, whence the assumption that ∂Φ/∂z|z=0${{\left. {\partial \Phi }/{\partial z}\; \right|}_{z=0}}$is constant for the Hankel solution becomes less accurate. For these conditions, we need to solve the reduced model to maintain low relative errors as compared to the full model. From the point of view of experimental measurements, it should be noted that the standard deviation in the measured phase increases as well when we sweep from 1 kHz to 1 MHz with around 20^o±10^o around 1 MHz [12].

We also see that Z₁ captures the measured magnitude reasonably well at high absolute values of the impedance but grows to a maximum relative error of around 70% at lower values of around 10 kΩ, whereas it only captures the phase close to 1 kHz (values in Fig. 5f close the diagonal) and then fails to describe the measured phase at frequencies larger than 1 kHz with a maximum relative error of around 220%. The latter is due to the fact that Z₁ does not account for the current flow through the viable skin.

At the microscopic level throughout the viable skin, we can determine the point-wise potential distribution by solving the reduced model numerically or take the approximate solution arising from the Hankel transform of the reduced model together with the assumption of constant currents at the electrodes. The latter is given by

Let us first look at the potential distributions at the electrodes, as depicted in Fig. 6 for the base-case conditions.

Fig.6

Again, we note the good agreement of the reduced model and approximate solution, Eq. 127, vis-à-vis the complete counterpart with a maximum relative error of around or less than 1%. Throughout the viable skin, the agreement remains good, as can be inferred from the visual comparison of the absolute value of the potential at 1 MHz in Fig. 7.

Fig.7

Finally, we note that we could, although we do not do so here, also calculate and visualize the potential distribution in the stratum corneum by returning to the earlier analysis during which we reduced the stratum corneum to a boundary condition.