Stochastic Bioimpedance-Based Channel Model of The Human Body for Galvanic Coupling

Abstract Human body communication (HBC) uses the human body as the channel to transfer data. Extensive work has been done to characterize the human body channel for different HBC techniques and scenarios. However, statistical channel bioimpedance characterisation of human body channels, particularly under dynamic conditions, remains relatively understudied. This paper develops a stochastic fading bioimpedance model for the human body channel using Monte Carlo simulations. Differential body segments were modelled as 2-port networks using ABCD parameters which are functions of bioimpedance based body parameters modelled as random variables. The channel was then modelled as the cascade of these random 2-port networks for different combinations of probability distribution functions (PDFs) assumed for the bioimpedance-based body parameters. The resultant distribution of the cascaded body segments varied for the different assumed bioimpedance based body parameter distributions and differential body segment sizes. However, considering the distribution names that demonstrated a best fit (in the top 3 PDF rankings) with highest frequency under the varying conditions, this paper recommends the distribution names: Generalized Pareto for phase distributions and Log-normal for magnitude distributions for each element in the overall cascaded random variable ABCD matrix.


Introduction
Human body communication (HBC) uses the human body as the channel to transfer data [1]. Appropriate channel models are extremely important for analysis and optimization of BAN (Body Area Network) performance. The state of the human body (e.g., the glucose/other analyte content, position of the electrodes on the body and action being performed by body parts) may be reflected by a correlation with the channel response. This data is invaluable for monitoring patients' conditions in medical use cases as well as monitoring athletes' conditions to optimize their performance [2,3,4]. The channel response may be used to develop application-specific physical layer technologies, such as Channel Equalization and Automatic Gain Control (AGC) techniques, which aim to improve the communication performance of transceivers used in BANs. Hence, channel models can be leveraged to improve BANs through better error performance, better receiver sensitivity, higher data rate, better channel resilience and power efficiency [5,6].
Extensive work has been done to characterize HBC channel models for different HBC communication techniques: eHBC (or E-field HBC which involves modulating electric fields through the methods: galvanic coupling and capacitive coupling) and mHBC (or magnetic HBC which involves modulating magneto-quasistatic fields through the method of magnetic induction) [1,5,7,8]. For example, [8] investigated the human body channel for 3 major electrode configurations: on-body, off-body and in-body. There has been notable work as well to develop channel circuit models for the different HBC communication techniques: eHBC (or E-field HBC, which involves modulating electric fields through the methods: galvanic coupling and capacitive coupling) and mHBC (or magnetic HBC which involves modulating magneto-quasistatic fields through the method of magnetic induction) [1,5,7,8]. However, relatively limited emphasis has been placed on the effect statistical properties of the body parameters have on the channel model. Some have suggested PDFs (Probability Density Func-tions) for characterizing fading for limited channel variations (sitting, standing and walking) [6,9] with no single PDF or channel model accounting for dynamic channel state; while others proposed fading margins to account for dynamic channels (e.g., [8,10]). In many instances the investigations focused upon the channel gain, for limited channel variations and electrode positions, with no reference to statistical characterization of fading as a result of these dynamic channel scenarios [11,10,12,13,14]. Finally, [13,14] suggested coherence bandwidths based on these limited channel variations. Consequently, in this paper, stochastic channel modelling for the human body is further investigated using Monte Carlo simulations, with focus on fading for the galvanic coupling HBC technique. It is noted however, that the results equally apply to other HBC techniques as well, due to the common underlying factors which would affect the channel response. The HBC band (i.e., 18.375MHz-23.625MHz) for galvanic communication as HBC is standardised in that region [15]. Consequently, for standardised communication, this model is developed focusing on the HBC band.

Channel Model
Consider the human body channel where currents flow between the electrodes placed on it as shown in Fig. 1. Galvanic coupling involves a differential, modulated signal applied between the transmitter (Tx) electrodes direct in contact with the body. Galvanic currents, are induced in the body and collected by the receiver (Rx) electrodes which are also in contact with the body [5]. The signal path is effected primarily through the skin [1]. Hence, the channel model is dependent upon skin properties.
The literature demonstrates that fading is as a result of channel variations (e.g., [8,11,10,12,13,14]). Fading is the change in channel gain, and can be related to bioimpedance as this affects the channel response. [1,16,17,18] all show signals sent in the HBC band propagate through the skin. In fact, the electric field, of the propagated signal, in the HBC band, is distributed through different layers of cells; where the skin and fat layers showing the greatest average electric field relative to the other tissue layers [17]. In [1] the authors propose a validated circuit-based static model for the skin for the galvanic coupling HBC technique. Circuit parametersskin resistance (R(ω)) and skin admittance (Y (ω) = G(ω)+jB(ω)) which comprises of a shunt circuit between the skin conductance (G(ω)) and skin susceptance (B(ω); where ω is the angular frequency of the signal propagating through the skin-represent the transfer function of the channel.
Bioimpedance is affected by the body's biological condition. As investigated in [19], circuit-based parameters depend on electrical properties of the body channel (i.e., skin conductivity and permittivity, in particular), which were observed to be roughly frequency flat for the HBC band (18.375MHz-23.625MHz). However, skin circuit parameters also depend on various non-deterministic, time-varying factors such as chemical composition of the body and channel geometry (e.g., skin thickness, skin deformation due to movement). Thus, non-deterministic, time-varying factors affect the body channel [20,21,22,23,24]. For example, [25] showed skin impedance changes when subject to externally-applied forces. In [26] the authors demonstrated skin impedance can vary when skin is damaged through needle punctures. Furthermore, [19] showed skin properties, which affect skin impedance and admittance [1], change when skin is subject to different conditions such as wet or dry scenarios. Bioimpedance is influenced by biological conditions. For example, glucose concentration affects body analyte composition [27]). [28] showed bioimpedance is related to biological state of liver ischemia. [29,30] showed bioimpdeance body parameters such as relative permittivity differ in the various body parts and subsections such as the dermis and epidermis skin layers.
Consequently, the human body channel is taken as a non-deterministic, time-varying channel. In this work, the human body is modelled as a lossy transmission line subdivided into sections, ∆z i (i.e., the i th body segment), as seen in Fig. 2. Thus, each differential/discretized section, of length ∆z i , can be modelled as an ABCD 2port network given by: This follows from [1], where the i th body differential segment parameters in (1) are: • Y i = G i + jB i , the admittance of the channel composed of the tissue conductance (G i ) and susceptance (B i ), which account for coupling between the conductive pathways of the tissue; and • R i , the shunt resistance of the tissue which accounts for signal propagation between cells.
Since the skin and body are inhomogeneous [19], this model simplifies the incorporation of this property by dividing the skin into homogeneous, differential sections; each with non-deterministic time-varying parameters. This segmented transmission line model form allows modelling of random variations due to causes of fading (e.g., movement and biological conditions in particular). Hence, the on-body human body channel can be represented as a random process matrix, Γ i , with probability density functions (PDFs) for the parameters that each of its elements depend on -R i , G i and B i (denoted by f R i , f G i and f B i respectively). Each ∆z i is initially, assumed to be independent and identically distributed (IID). Hence, the PDFs for each R i , G i and B i become R, G and B respectively with PDFs: f R , f G and f B respectively. To the best of the authors' knowledge, based on the extent of the literature surveyed, these PDFs have not been characterized.
Thus, Γ i can be considered as a continuous-valued, discrete-parameter random matrix, containing corresponding random processes: where N is the number of ∆z i segments that constitute the channel. Hence, the channel length, l = N i=1 ∆z i where ∆z i is the length of the i th segment. Since each ∆z i segment is assumed to be IID, then this process can be considered stationary as well [31]. Furthermore, this implies that one realization can be used to characterize at least the mean and auto-correlation, since it satisfies the index invariance (in this case, space invariance for short time intervals) condition, which points to Γ i being ergodic in the mean and auto-correlation [31].
Given the use of ABCD network modelling, the entire human body channel can be considered as an aggregate 2-port network equivalent to the cascade of N differential 2-port networks, to give the overall HBC network representation under matrix multiplication as follows: This research aims to statistically characterise Γ through numerical simulations based on circuit-based electrical body parameters.

Methodology
Even under the simplest assumptions for the channel segments, statistical models for the overall ABCD matrix for the channel may not be analytically tractable. Recall work has yet been done investigating the effect of random variations in biological state on the channel response. The authors' approach involves a Monte Carlo simulationbased investigation of various candidate models based upon assumptions for the PDFs for the channel segments which make sense practically. Distributions were assumed for f R , f G and f B with parameters determined through the literature surveyed. Segments were assumed to be independent and identically distributed. It is noted, however, that correlation can be factored in.
Through (1) and (2), the cascaded ABCD parameters of the human body channel may tend towards a distribution based on the assumed PDFs. The following PDFs for these parameters are considered: • Normal Distribution: This was chosen as a candidate for segments based upon scenarios where there may be a central tendency with a spread that can be approximated by the Gaussian distribution. Thus: where each PDF contain finite means (µ R , µ G and µ B ) and variances (σ 2 R , σ 2 G and σ 2 B ). The mean values, for R, G and B, were selected using the formulae for these parameters in [1]; with the values for the skin's permittivity (ε(ω)) and conductivity (σ(ω)) taken from the study done by [30] for the HBC band. The variances were set to different multiples (×1, ×2 and ×3 times) of the means to account for parameter changes seen in the literature surveyed (e.g., [20,21,22,23,24] which showed skin resistance changes for different scenarios).
• Uniform Distribution: This was chosen as a candidate for segments based upon scenarios where the range of values for R, G and B may be equally likely.
where the minimum (R min , G min and B min ) and maximum (R max , G max and B max ) values for each PDF are also finite. The minimum values were selected as 0 for each parameter as these are possible based on the literature surveyed. The maximums were set to different multiples of the means determined from the Normal distribution case as done for the variances.
• Gamma Distribution: Since f R , f G and f B may take on different shapes, the Gamma distribution was considered. This distribution is typically used to model wait times between Poisson distributed events; i.e. it models sums of exponentially distributed random variables. Thus, it can be linked to this use case. Furthermore, this use case takes advantage of this distribution's versatility in shape through manipulation of the shape parameters: a and b.
Thus Gamma(a B , b B ), where the shape parameters are finite.
The means and variances derived for the Normal distribution were used to derive the shape parameters of this distribution where: mean = ab and v ar i ance = ab 2 .
Note, these PDFs were chosen as they are commonly used for phenomena that have the traits stated above. However, the approach adopted can be generalized to accommodate other distributions as well (e.g., Rayleigh, Log-normal). Now consider each differential section of skin to be of dimensions: t × ∆z L × ∆z L , where t is skin thickness and ∆z L is the length of a differential segment (∆z i ). R, G and B are determined based upon [1]. Since each R, G and B are random variables, and the size ∆z L does not affect the values of R, G and B [1], the same PDFs for these parameters can be assumed for all values of ∆z L . Table 1 shows the PDFs assumed for R, G and B whose parameters were estimated based on the literature surveyed. Fig. 3 shows the algorithm followed to plot the PDFs of the elements of the resultant cascade, Γ, using the different combinations of the PDFs for R, G and B in Table1 for different ∆z L lengths (0.01mm to 10mm). These PDFs can be used to represent the dynamic channel.

Results and Analysis
This section presents the findings from evaluation of the stochastic model through simulation and empirical analysis.

Simulation of Stochastic Model
A repository of the results obtained from the simulation process can be accessed online [32]. Fig. 4 shows an example PDF plot approximation for the magnitude of element A of the resultant cascade, Γ, for test case A with ∆z L = 0.1mm. No single distribution was found to be the best fit for all variations of the inputted parameters. Table 2 and Fig. 5 show the frequency of distribution names (considering the top 3 best-fit distributions) for Γ's elements considering both magnitude and phase. Distribution names can be assumed for the magnitude and phase of each element of Γ based on the distribution names that occur with the highest frequency under the different PDF test cases and ∆z L sizes. Thus, based on this criteria, the following distribution names can be assumed for the magnitude and phase of each element of the cascaded random variable ABCD matrix, Γ: • Magnitude of Elements: Log-normal.
Inverse Gaussian demonstrated a marginally better fit than with the Log-normal distribution only in the case of the element D. However, this is small enough to validate the Log-normal best fit for element D.
NOTE: Only distribution names are suggested since in cases where they rank in the top 3 best fit PDFs, their shape parameters vary under the different conditions investigated.
Log-normal distributions are used in cases where a large number of individual effects, not strictly independent, act together on a signal [33]. Multiple factors, outlined in  The mean of the channel losses for the distributions (for the magnitude of element A where ∆z L = 0.1mm in test cases A-C) fall within measured static channel loss values in [1] and [12] under similar conditions simulated. Thus, to some extent, the PDFs derived for the magnitude of element A where ∆z L = 0.1mm in test cases A-C are validated. This gives credence to the channel length: ∆z L = 0.1mm being the best candidate to be used in the stochastic channel model of the human body for galvanic coupling.
Some additional findings from the simulated results are: • Increasing ∆z L decreases the mean of the distribu-tions for the magnitude of the elements of the cascaded random variable ABCD matrix, Γ. This is intuitive since a higher ∆z L length corresponds to a smaller cascade length and hence a smaller summation term for the cascaded element.
• Increasing ∆z L induces slight changes in the distribution parameters for the phases of the elements of Γ.
• Increasing channel length changes the distributions by increasing their mean. Since the distribution of the parameters for R, G and B are the same for all ∆z i segments for different ∆z L sizes, increasing the channel length will be equivalent to increasing the number of ∆z i samples. In other words, the results obtained for a ∆z L of 1mm will be the equivalent to increasing the channel for a ∆z L case of 10mm by a factor of 10.
• The magnitude and phase of each element of Γ are dependent. This was determined through finding the Pearson correlation coefficient between the magnitude and phase. This is a non zero value for all PDF test cases-of R, G and B-and ∆z L sizes except for 3 cases: for cases D,E and F for element C where ∆z L = 0.05m. Thus, magnitude and phase are correlated for each element in the cascaded random variable ABCD matrix except for 3 cases: for cases D,E and F for element C where ∆z L = 0.05m. A counter example was found for each of these uncorrelated cases that showed no independence; i.e. the joint PDF of magnitude and phase did not equate to the product of the marginal PDFs of magnitude and phase. Hence, the magnitude and phase are dependent for all elements in Γ under all conditions. Consequently, there is no further need to investigate independence between magnitude and phase.
• For all test cases, except test cases D-F (which considers uniform distributions for R, G and B), the magnitude and phase demonstrate high negative correlation.
For cases D-F there is relatively low correlation between the magnitude and phase. The magnitude and phase became more negatively correlated as the ∆z L size increases for cases D-F. Thus, the magnitude and phase of each element become less correlated as ∆z L increased only for cases where there are uniform distributions for R, G and B. Since the distribution of the parameters for R, G and B are the same for the ∆z i segments for different ∆z L sizes, increasing the channel length will be equivalent to increasing the number of ∆z i samples. In other words, the results obtained for a ∆z L of 1mm will be the equivalent to increasing the channel for a ∆z L case of 10mm by a factor of 10. Thus, as the channel length increases, the magnitude and phase become more correlated for each of the elements in Γ if uniform distributions are assumed for R, G and B.
• Distributions typically used for fitting tail ends of other distributions-extreme value and gereralized Pareto -fit better for the phase elements when compared to the other distribution fits. The Rayleigh distribution fit relatively well distributions associated with the phase elements for cases where ∆z L = 0.01mm. Rayleigh distributions are used to model complex Gaussian random events making it a suitable candidate for modelling the phase elements in Γ [33]. These 3 distributions have similar shapes which points to an empirically derived phase element distribution resembling that typically used to measure tail-ended distributions.

Conclusion
Stochastic models of the channel were produced for different channel model configurations via a simulation-based approach (∆z L size and input PDF combination cases that describe Γ i outlined in Table 1). These stochastic models changed as the channel model configurations changed suggesting that there was no single universal approximate distribution fading model. Further, the investigations focussed upon scenarios without channel segment correlation. Considering the distribution names that demonstrated a best fit (in the top 3 PDF rankings) with highest frequency under the different channel model configurations, this paper recommends the following distributions for consideration when modelling channel fading for dynamic HBC scenarios: Generalized Pareto for phase distributions and Log-normal for magnitude distributions for each element in the cascaded random variable ABCD matrix. Additionally, there is a need to investigate these findings empirically. Furthermore, the method for determining the stochastic fading model for the HBC technique, galvanic coupling, in this paper can be adopted in future work for other communication techniques (i.e., capacitive coupling or magnetic induction).

Ethical approval
The conducted research is not related to either human or animal use.