Mechanistic multilayer model for non-invasive bioimpedance of intact skin

The potential distribution, current density and measured impedance in cylinder coordinates [3] in the skin can be written as

1r∂∂rr∂Φ∂r+∂2Φ∂z2=0,$$\begin{array}{} \displaystyle \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial\Phi}{\partial r}\right)+\frac{\partial^{2}\Phi}{\partial_{z^{2}}}=0, \end{array}$$(1)

Jr=−σi,eff∂Φ∂r,$$\begin{array}{} \displaystyle J_{r}=-\displaystyle \sigma_{i,eff}\frac{\partial\Phi}{\partial r}, \end{array}$$(2)

Jz=−σi,eff∂Φ∂z,$$\begin{array}{} \displaystyle J_{z}=-\displaystyle \sigma_{i,eff}\frac{\partial\Phi}{\partial z}, \end{array}$$(3)

Z=V02π∫0R1Jz(0≤r≤R1,Hn)rdr−1,$$\begin{array}{} \displaystyle Z=V_{0}\left[2\displaystyle \pi\int\limits_{0}^{R_{1}}J_{z}(0\leq r\leq R_{1},\mathcal{H}_{n})rdr\right]^{-1}, \end{array}$$(4)

subject to the following external boundary conditions for j = 1, …, m electrodes:

Φ(R2j−2≤r≤R2j−1,Hn)=Vj,$$\begin{array}{} \displaystyle \Phi(R_{2j-2}\leq r\leq R_{2j-1},\, \mathcal{H}_{n})=V_{j}, \end{array}$$(5)

∂Φ∂z(R2j−1≤r≤R2j,Hn)=0,(j≠m)$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial z}(R_{2j-1}\leq r\leq R_{2j},\mathcal{H}_{n})=0,\, (j\neq m) \end{array}$$(6)

∂Φ∂z(R2m−1≤r,Hn)=0,$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial z}(R_{2m-1}\leq r,\, \mathcal{H}_{n})=0, \end{array}$$(7)

limr→∞Φ(r,z)=0,$$\begin{array}{} \displaystyle \lim_{r\rightarrow\infty}\Phi(r,z)=0, \end{array}$$(8)

∂Φ∂z(r,0)=0,$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial z}(r,0)=0, \end{array}$$(9)

∂Φ∂r(0,z)=0;$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial r}(0,z)=0; \end{array}$$(10)

and constitutive relations for the individual skin layers, i = 1, …, n:

σi,eff=σ(i)+iωε0εr(i),$$\begin{array}{} \displaystyle \sigma_{i,eff}=\sigma^{(i)}+i\omega\varepsilon_{0}\varepsilon_{r}^{(i)}, \end{array}$$

The radii, R_j, for the electrodes are defined for j = 1, …, m:

R0=0,$$\begin{array}{} \displaystyle R_{0}=0, \end{array}$$(11)

R1=w1,$$\begin{array}{} \displaystyle R_{1}=w_{1}, \end{array}$$(12)

Rj=Rj−1+wj,(j=2,…,m).$$\begin{array}{} \displaystyle R_{j}=R_{j-1}+w_{j},~~ (j=2,\ldots,m) . \end{array}$$(13)

For convenience, we also define the cumulative height, 𝓗_i, of skin layer i = 1, …, n:

Hi=∑k=1ihk.$$\begin{array}{} \displaystyle \mathcal{H}_{i}=\sum_{k=1}^{i}h_{k}. \end{array}$$(14)

In the above equations, r and z are the cylindrical coordinates (see Fig. 2), Φ is the complex-valued potential, J_r and J_z are the current densities in the r and z direction respectively, m and n are the number of electrodes in the EIS probe and layers of skin respectively; σ_i,eff, εr(i)$\begin{array}{} \displaystyle \varepsilon_{r}^{(i)} \end{array}$ and

σ⁽ⁱ⁾ are the effective, complex-valued electrical conductivity, the relative permittivity and conductivity of layer i, respectively; ε₀ is the relative permittivity of vacuum, i is the imaginary unit, ω = 2πν is the angular frequency, ν is the frequency, Z is the impedance, V_j is the prescribed potential at a given electrode, j; w_j is the width of or between an electrode as illustrated in Fig. 2 and h_k is the height of skin layer k.

Note that we do not solve for the potential distribution in the electrodes, because the potential drop there is negligible [3].

Let’s start by exploiting that the current density distributions underneath the m electrodes are near-to constant because of the high resistance to current of the stratum corneum in the considered frequency range [3]. This allows us to rewrite the boundary conditions, Eqs. 5-7, between the probe and the uppermost skin layer n, stratum corneum, as (we drop the subindex `eff’ for notational convenience in the analysis)

−σn∂Φ(r,Hn)∂z=∑j=1mIjAj[U(R2j−1−r)−U(R2j−2−r)],$$\begin{array}{} \displaystyle -\sigma_{n}\frac{\partial\Phi(r,\mathcal{H}_{n})}{\partial z}=\sum_{j=1}^{m}\frac{I_{j}}{A_{j}}\bigg[U(R_{2j-1}-r)\\ \displaystyle\qquad\qquad\quad\,-U(R_{2j-2}-r)\bigg], \end{array}$$(15)

where U() denotes the Heaviside function and I_j is the total current at the j^th electrode; the area of the j^th electrode is

Aj=π(R2j−12−R2j−22).$$\begin{array}{} \displaystyle A_{j}=\pi(R_{2_{j-1}}^{2}-R_{2_{j}-2}^{2}). \end{array}$$(16)

Strictly speaking, we are no longer solving exactly the same set of equations with the rewritten boundary condition, Eq. 15, whence we will refer to our generalized solution as approximate in nature.

From the conservation of current for the entire system, we know that

∑j=1mIj=0.$$\begin{array}{} \displaystyle \sum_{j=1}^{m}I_{j}=0. \end{array}$$(17)

We can now proceed with a Hankel transform of zeroth order, H₀ {}, for the potential,

Ψ(ξ,z)=H0{Φ(r,z)}=∫0∞rΦ(r,z)J0(ξr)dr,$$\begin{array}{} \displaystyle \Psi(\xi,z)=H_{0}\big\{\Phi(r,z)\big\}=\int_{0}^{\infty}{r\Phi}(r,z)J_{0}(\xi r)dr, \end{array}$$(18)

which is ideally suited for the axisymmetric Laplace equation, Eq. 1, on our unbounded domain:

H0{ΔΦ(r,z)}=−ξ2Ψ(ξ,z)+∂2∂z2Ψ(ξ,z)=0,$$\begin{array}{} \displaystyle H_{0} \big\{\Delta\Phi(r,z)\big\}=-\xi^{2}\Psi(\xi,z)+\frac{\partial^{2}}{\partial z^{2}}\Psi(\xi,z)=0, \end{array}$$(19)

with the solution for each skin layer, i,

Ψi(ξ,z)=Ai(ξ)e−ξz+Bi(ξ)eξz,$$\begin{array}{} \displaystyle \Psi_{i}(\xi,\, z)=A_{i}(\xi)e^{-\xi z}+B_{i}(\xi)e^{\xi z}, \end{array}$$(20)

and its z-derivative,

Ψi′(ξ,z)=−Ai(ξ)ξe−ξz+Bi(ξ)ξeξz;$$\begin{array}{} \displaystyle \Psi'_{i}(\xi,z)=-A_{i}(\xi)\xi e^{-\xi z}+B_{i}(\xi)\xi e^{\xi z}; \end{array}$$(21)

here, A_i (ξ) and B_i (ξ) are two functions that we will soon determine, Δ is the Laplace operator, and J₀ is the 0^th - order Bessel function of the first kind. Between layers, we know that the flux and potential are continuous; i.e.

σi∂Ψi(ξ,Hi)∂z=σi+1∂Ψi+1(ξ,Hi)∂z,$$\begin{array}{} \displaystyle \sigma_{i}\frac{\partial\Psi_{i}(\xi,\mathcal{H}_{i})}{\partial z}=\sigma_{i+1}\frac{\partial\Psi_{i+1}(\xi,\mathcal{H}_{i})}{\partial z}, \end{array}$$(22)

Ψi(ξ,Hi)=Ψi+1(ξ,Hi),$$\begin{array}{} \displaystyle \Psi_{i}(\xi,\mathcal{H}_{i})=\Psi_{i+1}(\xi,\mathcal{H}_{i}), \end{array}$$(23)

for 1 ≤ i ≤ n − 1. In addition, we need to satisfy the two boundary conditions at z = 0,

∂Ψ1(ξ,0)∂z=0,$$\begin{array}{} \displaystyle \frac{\partial\Psi_{1}(\xi,0)}{\partial z}=0, \end{array}$$(24)

and at z =𝓗_n,

−σn∂Ψn(ξ,Hn)∂z=∑k=1mIkκk(ξ),$$\begin{array}{} \displaystyle -\sigma_{n}\frac{\partial\Psi_{n}(\xi,\mathcal H_{n})}{\partial z} =\sum_{k=1}^{m}I_{k}\kappa_{k}(\xi) , \end{array}$$(25)

with

κj(ξ)=R2j−1J1(R2j−1ξ)−R2j−2J1(R2j−2ξ)π(R2j−12−R2j−22)ξ,$$\begin{array}{} \displaystyle \kappa_{j}(\xi)=\frac{R_{2{j-1}}J_{1}(R_{2j-1}\xi)-R_{2j-2}J_{1}(R_{2j-2}\xi)}{\pi(R_{2j-1}^{2}-R_{2j-2}^{2})\xi}, \end{array}$$(26)

here we have used that the Hankel transform of zeroth order of a Heaviside function is defined as

H0{U(a−r)}=aξJ1(aξ),$$\begin{array}{} \displaystyle H_{0}\displaystyle \big\{U(a-r)\big\}=\frac{a}{\xi}J_{1}(a\xi), \end{array}$$(27)

with J₁ and a (> 0) denoting the 1^st-order Bessel function of the first kind and a constant respectively [10].

Note that we have also switched to the subscript, k, in the summation over the electrodes in Eq. 25 to avoid confusion with the matrix notation we will introduce later. Let us solve the potential of the first layer and then work our way up to the uppermost n-th layer. Starting with the boundary condition, Eq. 24 at z = 0, we have

Ψ1′(ξ,0)=−A1(ξ)ξ+B1ξ(ξ)=0,$$\begin{array}{} \displaystyle \Psi_{1}'(\xi,0)=-A_{1}(\xi)\xi+B_{1}\xi(\xi)=0, \end{array}$$(28)

which results in

A1(ξ)=B1(ξ)e2ξH0Λ1,$$\begin{array}{} \displaystyle A_{1}(\xi)=B_{1}(\xi)e^{2\xi \mathcal{H}_{0}}\Lambda_{1}, \end{array}$$(29)

with

Λ1=1$$\begin{array}{} \displaystyle \Lambda_{1}=1 \end{array}$$(30)

H0=0;$$\begin{array}{} \displaystyle \mathcal{H}_{0}=0; \end{array}$$(31)

the role of e^2ξ𝓗₀ Λ₁ (=1) will become apparent as we proceed with our analysis. Proceeding to the other boundary condition for the first layer, Eq. 23, we have for z = 𝓗₁ that

Ψ1(ξ,H1)=Ψ2(ξ,H1),$$\begin{array}{} \displaystyle \Psi_{1}(\xi,\mathcal{H}_{1})=\Psi_{2}(\xi,\mathcal{H}_{1}), \end{array}$$(32)

i.e.,

B1(ξ)(Λ1e−ξH1+eξH1)=A2(ξ)e−ξH1+B2(ξ)eξH1,$$\begin{array}{} \displaystyle B_{1}(\xi)(\Lambda_{1}e^{-\xi \mathcal{H}_{1}}+e^{\xi \mathcal H_{1}})=A_{2}(\xi)e^{-\xi \mathcal{H}_{1}}+B_{2}(\xi)e^{\xi \mathcal{H}_{1}}, \end{array}$$(33)

from which we find (N.B.: 𝓗₁ = h₁)

B1(ξ)=A2(ξ)e−2ξH1+B2(ξ)Λ1e−2ξh1+1.$$\begin{array}{} \displaystyle B_{1}(\displaystyle \xi)=\frac{A_{2}(\xi)e^{-2\xi \mathcal{H}_{1}}+B_{2}(\xi)}{\Lambda_{1}e^{-2\xi h_{1}}+1}. \end{array}$$(34)

Moving on to the second layer, we find after satisfying the continuity of flux at z = 𝓗₁,

σ1∂Ψ1(ξ,H1)∂z=σ2∂Ψ2(ξ,H1)∂z,$$\begin{array}{} \displaystyle \sigma_{1}\frac{\partial\Psi_{1}(\xi,\mathcal H_{1})}{\partial z}=\sigma_{2}\frac{\partial\Psi_{2}(\xi,\mathcal{H}_{1})}{\partial z}, \end{array}$$(35)

that

A2(ξ)=B2(ξ)e2ξH1Λ2(ξ),$$\begin{array}{} \displaystyle A_{2}(\xi)=B_{2}(\xi)e^{2\xi\mathcal H_{1}}\Lambda_{2}(\xi), \end{array}$$(36)

with

Λ2(ξ)=σ1(1−Λ1e−2ξh1)−σ2(1+Λ1e−2ξh1)σ1(Λ1e−2ξh1−1)−σ2(1+Λ1e−2ξh1).$$\begin{array}{} \displaystyle \Lambda_{2}(\xi)=\frac{\sigma_{1}(1-\Lambda_{1}e^{-2\xi h_{1}})-\sigma_{2}(1+\Lambda_{1}e^{-2\xi h_{1}})}{\sigma_{1}(\Lambda_{1}e^{-2\xi h_{1}}-1)-\sigma_{2}(1+\Lambda_{1}e^{-2\xi h_{1}})}. \end{array}$$(37)

Between the second and third layer, z = 𝓗₂, we require continuity of the potential,

Ψ2(ξ,H2)=Ψ3(ξ,H2),$$\begin{array}{} \displaystyle \Psi_{2}(\xi,\mathcal{H}_{2})=\Psi_{3}(\xi,\, \mathcal{H}_{2}), \end{array}$$(38)

which reveals

B2(ξ)=A3(ξ)e−2ξH2+B3(ξ)Λ2(ξ)e−2ξh2+1.$$\begin{array}{} \displaystyle B_{2}(\displaystyle \xi)=\frac{A_{3}(\xi)e^{-2\xi \mathcal{H}_{2}}+B_{3}(\xi)}{\Lambda_{2}(\xi)e^{-2\xi h_{2}}+1}. \end{array}$$(39)

Let us explicitly continue with the third layer and then generalize the solution for layer i, where 1 ≤ i ≤ n − 1. In the third layer, we find

A3(ξ)=B3(ξ)e2ξH2Λ3(ξ),$$\begin{array}{} \displaystyle A_{3}(\xi)=B_{3}(\xi)e^{2\xi \mathcal{H}_{2}}\Lambda_{3}(\xi), \end{array}$$(40)

and

B3(ξ)=A4(ξ)e−2ξH3+B4(ξ)Λ3(ξ)e−2ξh3+1,$$\begin{array}{} \displaystyle B_{3}(\displaystyle \xi)=\frac{A_{4}(\xi)e^{-2\xi \mathcal{H}_{3}}+B_{4}(\xi)}{\Lambda_{3}(\xi)e^{-2\xi h_{3}}+1}, \end{array}$$(41)

with

Λ3(ξ)=σ2(1−Λ2e−2ξh2)−σ3(1+Λ2e−2ξh2)σ2(Λ2e−2ξh2−1)−σ3(1+Λ2e−2ξh2),$$\begin{array}{} \displaystyle \Lambda_{3}(\xi)=\frac{\sigma_{2}(1-\Lambda_{2}e^{-2\xi h_{2}})-\sigma_{3}(1+\Lambda_{2}e^{-2\xi h_{2}})}{\sigma_{2}(\Lambda_{2}e^{-2\xi h_{2}}-1)-\sigma_{3}(1+\Lambda_{2}e^{-2\xi h_{2}})}, \end{array}$$(42)

after satisfying the boundary conditions with the second and fourth skin layer.

We can now generalize the solution as

Ai(ξ)=Bi(ξ)e2ξHi−1Λi(ξ),$$\begin{array}{} \displaystyle A_{i}(\xi)=B_{i}(\xi)e^{2\xi\mathcal H_{i-1}}\Lambda_{i}(\xi), \end{array}$$(43)

and

Bi(ξ)=Ai+1(ξ)e−2ξHi+Bi+1(ξ)Λi(ξ)e−2ξhi+1,$$\begin{array}{} \displaystyle B_{i}( \xi)=\frac{A_{i+1}(\xi)e^{-2\xi \mathcal{H}_{i}}+B_{i+1}(\xi)}{\Lambda_{i}(\xi)e^{-2\xi h_{i}}+1}, \end{array}$$(44)

for 1 ≤ i ≤ n − 1, with

Λ1=1,$$\begin{array}{} \displaystyle \Lambda_{1}=1, \end{array}$$(45)

for i = 1, and the recursive formula for Λ_i,

Λi(ξ)=σi−1(1−Λi−1e−2ξhi−1)−σi(1+Λi−1e−2ξhi−1)σi−1(Λi−1e−2ξhi−1−1)−σi(1+Λi−1e−2ξhi−1),$$\begin{array}{} \displaystyle \Lambda_{i}(\xi)=\frac{\sigma_{i-1}(1-\Lambda_{i-1}e^{-2\xi h_{i-1}})-\sigma_{i}(1+\Lambda_{i-1}e^{-2\xi h_{i-1}})}{\sigma_{i-1}(\Lambda_{i-1}e^{-2\xi h_{i-1}}-1)-\sigma_{i}(1+\Lambda_{i-1}e^{-2\xi h_{i-1}})}, \end{array}$$(46)

for 2 ≤ i ≤ n − 1.

In the n-th layer, we find after satisfying the continuity of potential with the n − 1 layer and the boundary condition, Eq. 25, between the electrodes and the uppermost skin layer,

An(ξ)=Bn(ξ)e2ξHn−1Λn(ξ),$$\begin{array}{} \displaystyle A_{n}(\xi)=B_{n}(\xi)\mathrm{e}^{2\xi\mathcal H_{n-1}}\Lambda_{n}(\xi), \end{array}$$(47)

by setting

σn−1∂Ψn−1(ξ,Hn−1)∂z=σn∂Ψn(ξ,Hn−1)∂z;$$\begin{array}{} \displaystyle \sigma_{n-1}\frac{\partial\Psi_{n-1}(\xi,\mathcal{H}_{n-1})}{\partial z}=\sigma_{n}\frac{\partial\Psi_{n}(\xi,\mathcal{H}_{n-1})}{\partial z}; \end{array}$$(48)

and

Bn(ξ)=e−ξHn∑k=1mIkκk(ξ)σnξΛn(ξ)e−ξhn−1,$$\begin{array}{} \displaystyle B_{n}(\xi)=\frac{e^{-\xi\mathcal H_{n}}\sum\limits_{k=1}^{m}I_{k}\kappa_{k}(\xi)}{\sigma_{n}\xi\left[ \Lambda_{n}(\xi)e^{-\xi h_{n}}-1\right]}, \end{array}$$(49)

from the uppermost boundary condition at z = 𝓗_n for layer n:

−σn∂Ψn(ξ,Hn)∂z=∑k=1mIkkk(ξ).$$\begin{array}{} -\displaystyle \sigma_{n}\frac{\partial\Psi_{n}(\xi,\mathcal{H}_{n})}{\partial z}=\sum_{k=1}^{m}I_{k}k_{k}(\xi). \end{array}$$(50)

At this stage, we can write the local, pointwise transformed potential in layer i as

Ψi(ξ,z)=Bi(ξ)Λi(ξ)eξ(2Hi−1−z)+eξz.$$\begin{array}{} \displaystyle \Psi_{i}(\xi,z)=B_{i}(\xi)\left[\Lambda_{i}(\xi)e^{\xi(2H_{i-1}-z)}+e^{\xi z}\right]. \end{array}$$(51)

In this expression, the coordinate, z, is bounded for each layer as

Hi−1≤z(in layeri)≤Hi.$$\begin{array}{} \displaystyle \mathcal{H}_{i-1}\leq z (\text{in layer}~ i) \leq \mathcal{H}_{i}. \end{array}$$(52)

Now, taking the inverse Hankel transform,

Υi(r,z)=H0−1{Ψi(ξ,z)}=∫0∞ξΨi(ξ,z)J0(ξr)dξ,$$\begin{array}{} \displaystyle \Upsilon_{i}(r,z)=H_{0}^{-1}\{\Psi_{i}(\xi,z)\}=\int_{0}^{\infty}\xi\Psi_{i}(\xi,z)J_{0}(\xi r)d\xi, \end{array}$$(53)

we find the local potential distribution throughout the viable skin as

Υi(r,z)=∫0∞ξBi(ξ)Λi(ξ)eξ(2Hi−1−z)+eξZJo(ξr)dξ$$\begin{array}{} \displaystyle \Upsilon_{i}(r,z)=\int_{0}^{\infty}\xi B_{i}(\xi)\left[\Lambda_{i}(\xi)e^{\xi(2\mathcal{H}_{i-1}-z)}+e^{\xi_{Z}}\right]J_{o}(\xi r)d\xi \end{array}$$(54)

and

Φi(r,z)=Υi(r,z)+c.$$\begin{array}{} \displaystyle \Phi_{i}(r,z)=\Upsilon_{i}(r,z)+c. \end{array}$$(55)

We have m + 1 unknows: I₁ to I_m and c; the latter constant provides a reference point for the potential since we have in the approximate solution only solved for Neuman boundary conditions. In order to determine the unknowns, we return to our pointwise Dirichlet boundary conditions for the potentials underneath the electrodes, Eq. 5, and rewrite Eqs. 54 and 55 in terms of the average potential underneath the j^th electrode in the uppermost skin layer, n:

Φn,javg=2πAj∫R2j−2R2j−1rΥn(r,Hη)dr+c=Vj,$$\begin{array}{} \displaystyle \Phi_{n,j}^{avg}=\frac{2\pi}{A_{j}}\int_{R_{2j-2}}^{R_{2j-1}}r\Upsilon_{n}(r,\mathcal{H}_{\eta})dr+c=V_{j}, \end{array}$$(56)

or, after integration once, as

Φn,javg=2πσn∫0∞∑k=1mIkκk(ξ)κj(ξ)Ω(ξ)dξ+c=Vj,$$\begin{array}{} \displaystyle \Phi_{n,j}^{avg}=\frac{2\pi}{\sigma_{n}}\int_{0}^{\infty}\left[\sum_{k=1}^{m}I_{k}\kappa_{k}(\xi)\right]\kappa_{j}(\xi)\Omega(\xi)d\xi+c=V_{j}, \end{array}$$(57)

where

Ω(ξ)=Λn(ξ)e−2ξhn+1Λn(ξ)e−2ξhn−1.$$\begin{array}{} \displaystyle \Omega(\xi)=\frac{\Lambda_{n}(\xi)e^{-2\xi h_{n}}+1}{\Lambda_{n}(\xi)e^{-2\xi h_{n}}-1}. \end{array}$$(58)

During the integration, we used that

∫J0(ar)rdr=raJ1(ar)+constant;$$\begin{array}{} \displaystyle \int J_{0}(ar)rdr=\frac{r}{a}J_{1}(ar)+ \text{constant}; \end{array}$$(59)

J1(0)=0.$$\begin{array}{} \displaystyle J_{1}(0)=0. \end{array}$$(60)

The overall conservation of current, Eq. 17, allows us to reduce the number of unknowns by one, as we can write

Im=−∑j=1m−1Ij.$$\begin{array}{} \displaystyle I_{m}=-\displaystyle \sum_{j=1}^{m-1}I_{j}. \end{array}$$(61)

Introducing the expression for the current of the outermost electrode, I_m, into Eq. 57, yields the following system of m linear equations:

Ax=b,$$\begin{array}{} \displaystyle \mathbf{Ax}=\mathbf{b}, \end{array}$$(62)

with

A=A1,1A1,2⋯A1,m−11A2,1A2,2⋯A2,m−11⋮⋮⋱⋮⋮Am−1,1Am−1,2⋯Am−1,m−11Am,1Am,2⋯Am,m−11,$$\begin{array}{} \displaystyle \mathbf{A}=\left[\begin{array}{lllll} A_{1,1} & A_{1,2} & \cdots & A_{1,m-1} & 1\\ A_{2,1} & A_{2,2} & \cdots & A_{2,m-1} & 1\\ ~~~~\vdots & ~~~~\vdots & \ddots & ~~~~~\vdots & \vdots\\ A_{m-1,1} & A_{m-1,2} & \cdots & A_{m-1,m-1} & 1\\ A_{m,1} & A_{m,2} & \cdots & A_{m,m-1} & 1 \end{array}\right], \end{array}$$(63)

x=IlI2⋮Im−lc,$$\begin{array}{} \displaystyle \mathbf{x}=\left[\begin{array}{l} ~~I_{\mathrm{l}}\\ ~~I_{2}\\ ~~~\vdots\\ I_{m-\mathrm{l}}\\ ~~~c \end{array}\right], \end{array}$$(64)

b=V1V2⋮Vm−lVm,$$\begin{array}{} \displaystyle \mathbf{b}=\left[\begin{array}{l} ~~V_{1}\\ ~~V_{2}\\ ~~~~\vdots\\ V_{m-\mathrm{l}}\\ ~~V_{m} \end{array}\right], \end{array}$$(65)

and

Ai,j=2πσn∫0∞κj(ξ)−κm(ξ)κi(ξ)Ω(ξ)dξ.$$\begin{array}{} \displaystyle A_{i,j}=\displaystyle \frac{2\pi}{\sigma_{n}}\int_{0}^{\infty}\left[\kappa_{j}(\xi)-\kappa_{m}(\xi)\right]\kappa_{i}(\xi)\Omega(\xi)d\xi. \end{array}$$(66)

Once the linear system of equations together with the improper integrals given by Eq. 66 have been solved, the impedance can be found from

ZH=V0I1,$$\begin{array}{} \displaystyle Z_{H}=\displaystyle \frac{V_{0}}{I_{1}}, \end{array}$$(67)

assuming that V₀ is set as reference potential for the measurements and that the first electrode is used as the sense. We’ll give one example of this in the next subsection. If the electrode system is different, one needs to couple the right current with the given potential to determine the impedance.

This solution corresponds to EIS measurements in our earlier experimental study [9] and the model thereof with skin encompassing three layers (n=3): stratum corneum, viable skin and adipose tissue, as illustrated in Fig. 3.

Fig. 3

Schematic of a 4-electrode probe and skin modeled as a three-layer entity.

Here, the first electrode is the sense, the second serves as the guard, the third is the second inject with a potential set as αV₀ and the fourth is the primary inject with a potential of V₀; α is referred to as a depth setting on the probe; i.e. m=4. The boundary conditions for the complexvalued Laplace equation reduce to

Φ(0≤r≤R1,H3)=0,$$\begin{array}{} \displaystyle \Phi(0\leq r\leq R_{1},\ \mathcal{H}_{3})=0, \end{array}$$(68)

Φ(R2≤r≤R3,H3)=0,$$\begin{array}{} \displaystyle \Phi(R_{2}\leq r\leq R_{3},\, \mathcal{H}_{3})=0, \end{array}$$(69)

Φ(R4≤r≤R5,H3)=αV0,$$\begin{array}{} \displaystyle \Phi(R_{4}\leq r\leq R_{5},\, \mathcal{H}_{3})=\alpha V_{0}, \end{array}$$(70)

Φ(R6≤r≤R7,H3)=V0,$$\begin{array}{} \displaystyle \Phi(R_{6}\leq r\leq R_{7},\,\mathcal{H}_{3})=V_{0}, \end{array}$$(71)

∂Φ∂z(R1≤r≤R2,H3)=0,$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial z}(R_{1}\leq r\leq R_{2},\, \mathcal{H}_{3})=0, \end{array}$$(72)

∂Φ∂z(R3≤r≤R4,H3)=0,$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial z}(R_{3}\leq r\leq R_{4},\,\mathcal{H}_{3})=0, \end{array}$$(73)

∂Φ∂z(R5≤r≤R6,H3)=0,$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial z}(R_{5}\leq r\leq R_{6},\,\mathcal{H}_{3})=0, \end{array}$$(74)

∂Φ∂z(R7≤r,H3)=0,$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial z}(R_{7}\leq r,\,\mathcal{H}_{3})=0, \end{array}$$(75)

and

limr→∞Φ(r,z)=0,$$\begin{array}{} \displaystyle \lim_{r\rightarrow\infty}\Phi(r,z)=0, \end{array}$$(76)

∂Φ∂z(r,0)=0,$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial z}(r,0)=0, \end{array}$$(77)

∂Φ∂r(0,z)=0.$$\begin{array}{} \displaystyle \frac{\partial\Phi}{\partial r}(0,z)=0. \end{array}$$(78)

Our generalized solution becomes

Λ1=1,$$\begin{array}{} \displaystyle \Lambda_{1}=1, \end{array}$$(79)

Λ2(ξ)=σ11−Λ1e−2ξh1−σ21+Λ1e−2ξh1σ1Λ1e−2ξh1−1−σ21+Λ1e−2ξh1$$\begin{array}{} \displaystyle \Lambda_{2}(\xi)=\frac{\sigma_{1}\left(1-\Lambda_{1}e^{-2\xi h_{1}}\right)-\sigma_{2}\left(1+\Lambda_{1}e^{-2\xi h_{1}}\right)}{\displaystyle\sigma_{1}\left(\Lambda_{1}e^{-2\xi h_{1}}-1\right)-\sigma_{2}\left(1+\Lambda_{1}e^{-2\xi h_{1}}\right)} \end{array}$$(80)

Λ3(ξ)=σ21−Λ2e−2ξh2−σ31+Λ2e−2ξh2σ2Λ2e−2ξh2−1−σ31+Λ2e−2ξh2,$$\begin{array}{} \displaystyle \Lambda_{3}(\xi)=\frac{\sigma_{2}\left(1-\Lambda_{2}e^{-2\xi h_{2}}\right)-\sigma_{3}\left(1+\Lambda_{2}e^{-2\xi h_{2}}\right)}{\displaystyle\sigma_{2}\left(\Lambda_{2}e^{-2\xi h_{2}}-1\right)-\sigma_{3}\left(1+\Lambda_{2}e^{-2\xi h_{2}}\right)}, \end{array}$$(81)

where

Ω(ξ)=Λ3(ξ)e−2ξh3+1Λ3(ξ)e−2ξh3−1.$$\begin{array}{} \displaystyle \Omega(\xi)=\frac{\Lambda_{3}(\xi)e^{-2\xi h_{3}}+1}{\Lambda_{3}(\xi)e^{-2\xi h_{3}}-1}. \end{array}$$(82)

The individual average currents for the four electrodes are as follows (we have dropped the comma notation for the subscripts of the elements in A for notational convenience; b = [0, 0, αV₀, V₀]^T):

I1=V0D(−A12A23+A33A12−A32A13−A23αA42+A13αA42−A13αA22−A12A43α+A12A23α+A43αA22+A32A23−A33A22+A13A22),$$\begin{array}{} I_{1}=\displaystyle \frac{V_{0}}{D}(-A_{12}A_{23}+A_{33}A_{12}-A_{32}A_{13}-A_{23}\alpha A_{42}\\ \quad+A_{13}\alpha A_{42}-A_{13}\alpha A_{22}-A_{12}A_{43}\alpha+A_{12}A_{23}\alpha\\ \quad~~+A_{43}\alpha A_{22}+A_{32}A_{23}-A_{33}A_{22}+A_{13}A_{22})\ , \end{array}$$(83)

I2=−V0D(A11A33−A11A43α−A11A23+A11A23α−A13A21α−A31A13+A41A13α+A43A21α−A23αA41+A31A23+A13A21−A33A21),$$\begin{array}{} I_{2}=-\displaystyle \frac{V_{0}}{D}(A_{11}A_{33}-A_{11}A_{43}\alpha-A_{11}A_{23}+A_{11}A_{23}\alpha\\ \quad~-A_{13}A_{21}\alpha-A_{31}A_{13}+A_{41}A_{13}\alpha+A_{43}A_{21}\alpha\\ \qquad-A_{23}\alpha A_{41}+A_{31}A_{23}+A_{13}A_{21}-A_{33}A_{21}) , \end{array}$$(84)

I3=V0D(αA41A12−αA41A22−A11αA42+A32A11−A31A12+αA42A21−αA12A21−A32A21+A12A21+A31A22−A11A22+A11αA22),$$\begin{array}{} I_{3}=\displaystyle \frac{V_{0}}{D}(\alpha A_{41}A_{12}-\alpha A_{41}A_{22}-A_{11}\alpha A_{42}+A_{32}A_{11}\\ \qquad-A_{31}A_{12}+\alpha A_{42}A_{21}-\alpha A_{12}A_{21}-A_{32}A_{21}\\ \qquad+A_{12}A_{21}+A_{31}A_{22}-A_{11}A_{22}+A_{11}\alpha A_{22}), \end{array}$$(85)

and

I4=−I1−I2−I3,$$\begin{array}{} I_{4}=-I_{1}-I_{2}-I_{3}, \end{array}$$(86)

with the constants

c=V0D(−A13A21αA42+A23A11αA42−A22A11A43α−A33A12A21−A31A13A22−A11A32A23+A11A33A22+A31A12A23+A32A13A21+A41A13αA22+A43αA12A21−A23αA41A12),$$\begin{array}{} c=\displaystyle\frac{V_{0}}{D}(-A_{13}A_{21}\alpha A_{42}+A_{23}A_{11}\alpha A_{42}-A_{22}A_{11}A_{43}\alpha\\ \qquad\quad-A_{33}A_{12}A_{21}-A_{31}A_{13}A_{22}-A_{11}A_{32}A_{23}\\ \qquad\quad+A_{11}A_{33}A_{22}+A_{31}A_{12}A_{23}+A_{32}A_{13}A_{21}\\ \quad~+A_{41}A_{13}\alpha A_{22}+A_{43}\alpha A_{12}A_{21}-A_{23}\alpha A_{41}A_{12}), \end{array}$$(87)

and

D=A33A41A12−A33A41A22−A13A21A42+A31A13A42−A31A13A22−A23A41A12−A31A23A42−A33A12A21+A32A41A23−A32A41A13+A12A21A43+A31A12A23−A31A12A43+A31A22A43+A22A41A13+A33A42A21+A11A33A22+A11A32A43−A11A33A42+A32A13A21−A32A43A21−A11A22A43+A11A23A42−A11A32A23).$$\begin{array}{} D=A_{33}A_{41}A_{12}-A_{33}A_{41}A_{22}-A_{13}A_{21}A_{42} \\ ~~+A_{31}A_{13}A_{42}-A_{31}A_{13}A_{22}-A_{23}A_{41}A_{12} \\ ~~-A_{31}A_{23}A_{42}-A_{33}A_{12}A_{21}+A_{32}A_{41}A_{23} \\ ~~-A_{32}A_{41}A_{13}+A_{12}A_{21}A_{43}+A_{31}A_{12}A_{23} \\ ~~-A_{31}A_{12}A_{43}+A_{31}A_{22}A_{43}+A_{22}A_{41}A_{13} \\ ~~+A_{33}A_{42}A_{21}+A_{11}A_{33}A_{22}+A_{11}A_{32}A_{43} \\ ~~-A_{11}A_{33}A_{42}+A_{32}A_{13}A_{21}-A_{32}A_{43}A_{21} \\ ~-A_{11}A_{22}A_{43}+A_{11}A_{23}A_{42}-A_{11}A_{32}A_{23}). \end{array}$$(88)

Note that the entries, A_ij, in the above expressions entail solving the improper integral, Eq. 66, twelve times.

The parameters are given in Table 1.

For the skin, we can write the constitutive relations for the relative permittivity and the conductivity (S m^-1) as follows [9] based on the dimensionless logarithmic frequency, N = log₁₀(ν/1 Hz).

In the stratum corneum, i=3, we have

log10⁡(σ(3))=−0.011516×N5+0.25509×N4−2.2537×N3+9.9955×N2−21.404×N+11.803,$$\begin{array}{} \quad~\log_{10}(\sigma^{(3)})=-0.011516\times N^{5}+0.25509\times N^{4}\\ -2.2537\times N^{3}+9.9955\times N^{2}-21.404\times N+11.803, \end{array}$$(89)

log10⁡(εr(3))=−0.011319×N5+0.23953×N4−2.0255×N3+8.5278×N2−17.961×N+17.57,$$\begin{array}{} \quad\,\log_{10}(\varepsilon_{r}^{(3)})=-0.011319\times N^{5}+0.23953\times N^{4}\\ -2.0255\times N^{3}+8.5278\times N^{2}-17.961\times N+17.57, \end{array}$$(90)

In the viable skin, i=2, we have

log10⁡(σ(2))=0.015178×N5−0.34167×N4+3.0088×N3−12.952×N2+27−471×N−23.688,$$\begin{array}{} \quad~~\log_{10}(\sigma^{(2)})=0.015178\times N^{5}-0.34167\times N^{4}\\ +3.0088\times N^{3}-12.952\times N^{2}+27_{-}471\times N-23.688, \end{array}$$(91)

log10⁡(εr(2))=−0.044465×N5+0.92429×N4−7.7649×N3+32.847×N2−70.466×N+67.61,$$\begin{array}{} \quad~\log_{10}(\varepsilon_{r}^{(2)})=-0.044465\times N^{5}+0.92429\times N^{4}\\ -7.7649\times N^{3}+32.847\times N^{2}-70.466\times N+67.61, \end{array}$$(92)

And finally in the adipose tissue, i=1, we have

σ(1)=1/54.9236,$$\begin{array}{} \sigma^{(1)}=1/54.9236, \end{array}$$(93)

log10⁡(εr(1))=−0.00052423×N5−0.00040926×N4+0.11797×N3−0.77214×N2+0.083142×N+7.584∠$$\begin{array}{} \quad\,\log_{10}(\varepsilon_{r}^{(1)})=-0.00052423\times N^{5}-0.00040926\times N^{4}\\ +0.11797\times N^{3}-0.77214\times N^{2}+0.083142\times N+7.584{\small\angle} \end{array}$$(94)

We solve the model, Eqs. 1-13, in the finite-element solver Comsol Multiphysics 5.2a [13] and the approximate semianalytical model, Eqs. 79-88, and its pointwise solution in the skin layers, Eq. 55, in the general-purpose numerical software Matlab R2015 [14].